Pollinator and host sharing lead to hybridization and introgression in Panamanian free‐standing figs, but not in their pollinator wasps

Abstract Obligate pollination mutualisms, in which plant and pollinator lineages depend on each other for reproduction, often exhibit high levels of species specificity. However, cases in which two or more pollinator species share a single host species (host sharing), or two or more host species share a single pollinator species (pollinator sharing), are known to occur in current ecological time. Further, evidence for host switching in evolutionary time is increasingly being recognized in these systems. The degree to which departures from strict specificity differentially affect the potential for hybridization and introgression in the associated host or pollinator is unclear. We addressed this question using genome‐wide sequence data from five sympatric Panamanian free‐standing fig species (Ficus subgenus Pharmacosycea, section Pharmacosycea) and their six associated fig–pollinator wasp species (Tetrapus). Two of the five fig species, F. glabrata and F. maxima, were found to regularly share pollinators. In these species, ongoing hybridization was demonstrated by the detection of several first‐generation (F1) hybrid individuals, and historical introgression was indicated by phylogenetic network analysis. By contrast, although two of the pollinator species regularly share hosts, all six species were genetically distinct and deeply divergent, with no evidence for either hybridization or introgression. This pattern is consistent with results from other obligate pollination mutualisms, suggesting that, in contrast to their host plants, pollinators appear to be reproductively isolated, even when different species of pollinators mate in shared hosts.


| INTRODUC TI ON
Hybridization and introgression have contributed to the evolution of species and clades across the tree of life (Anderson, 1953;Mallet, 2005;Mallet et al., 2016;Stebbins, 1959). Genome-scale data and advances in analytical techniques, in particular, have enabled the detection and documentation of hybridization and introgression in many groups and suggest that these processes are more widespread than previously thought (Taylor & Larson, 2019). These processes can be important contributors to speciation and adaptive radiations, which can be spurred by the introduction of beneficial alleles and multilocus combinations to a recipient lineage (adaptive introgression; e.g., Edelman & Mallet, 2021;Hedrick, 2013). Adaptive introgression has been documented in a diversity of lineages, including cichlid fishes (Malinsky et al., 2018;Meier et al., 2017;Svardal et al., 2020), butterflies (Edelman et al., 2019;Enciso-Romero et al., 2017;Pardo-Diaz et al., 2012), and oaks (Eaton et al., 2015;Leroy et al., 2020;McVay et al., 2017). Interspecific gene flow can also contribute to adaptive responses to rapidly changing environments, such as human-caused environmental modifications (Hamilton & Miller, 2016). Consequently, there is growing appreciation for the importance of hybridization and introgression in generating and maintaining organismal adaptation and diversity.  de Medeiros & Farrell, 2020;Hembry & Althoff, 2016;Pellmyr et al., 2020). Brood pollination mutualisms often exhibit strict host specificity (i.e., only one pollinator species is consistently associated with only one host species), which is thought to promote reproductive isolation for both the host and the pollinator.
Several studies, however, have revealed that two or more pollinator species per host (host sharing), or two or more host species per pollinator (pollinator sharing), are not uncommon in some systems (Cornille et al., 2012;McLeish & Van Noort, 2012;Molbo et al., 2003;Starr et al., 2013;Su et al., 2022;Wang et al., 2016;Yang et al., 2015). These reduced specificity ecological interactions, along with accumulating molecular evidence of historical host switching (e.g., Cruaud et al., 2012;Hembry et al., 2013;Satler et al., 2019), are consistent with opportunities for hybridization in either the host, the pollinator, or both (Arteaga et al., 2020;Berg, 1989;Cornille et al., 2012;Leebens-Mack et al., 1998;Machado et al., 2005;Rentsch & Leebens-Mack, 2012;Wang et al., 2016;. This raises the more general question of how different patterns of host specificity affect the evolutionary dynamics (e.g., hybrid formation and introgression versus reproductive isolation) for each partner species in obligate plant-pollinator mutualisms. These dynamics, in turn, will affect the processes of speciation and diversification in both taxa.

Figs (Ficus, family Moraceae) and their pollinator wasps (fam-
ily Agaonidae) represent an ancient (~80 Ma) and diverse (~900 described species of figs) obligate pollination mutualism (Cook & Rasplus, 2003;Cruaud et al., 2012;Janzen, 1979;Machado et al., 2001;Weiblen, 2002). Pollination in this keystone mutualism results in fruit production that supports diverse frugivores across tropical and subtropical habitats worldwide (Shanahan et al., 2001). When receptive, the enclosed fig inflorescences (syconia) produce volatile chemicals that attract mated pollen-bearing female fig wasps (Cornille et al., 2012;Grison-Pigé et al., 2002;Hossaert-McKey et al., 2010;Van Noort et al., 1989;Wang, Yang, et al., 2021;Ware et al., 1993). These foundress wasps enter the syconia to pollinate flowers and oviposit in a subset of them. The pollinated flowers usually develop as viable seeds, but those flowers that receive wasp eggs usually become galls that support the development of the wasp offspring. After maturing, the pollinator wasp offspring then mate within their natal fig syconia before females collect pollen and disperse-typically several kilometers (Ahmed et al., 2009;Herre, 1989;Nason et al., 1998) Weiblen, 2002), which has contributed to the paradigm of one wasp to one fig. However, with increasingly deep spatial and temporal sampling, and the application of more sophisticated genetic techniques, it is increasingly clear that there is often considerable deviation from this paradigm. For example, several studies (based primarily on mitochondrial DNA) suggest some level of pollinator sharing or host sharing, and also imply host switching (e.g., Cornille et al., 2012;Darwell et al., 2014;Machado et al., 2005;Molbo et al., 2003;Wachi et al., 2016;Wang et al., 2016;Yang et al., 2015;Yu et al., 2019). Further, roughly 30% of fig species has been estimated to host multiple pollinator species at either local or regional scales (Yang et al., 2015). This contemporary host sharing is complemented by detailed studies documenting an evolutionary history of host switching (Satler et al., 2019).  Berg, 1989;Jackson et al., 2008;Machado et al., 2005;Wilde et al., 2020). This appears to be the case even among distantly-related fig species across distinct subgenera (Compton, 1990;Ramírez, 1994;. In contrast to their fig hosts, however, the few studies conducted to date using microsatellites or even deeper genomic tools find little or no evidence of hybridization or successful introgression between co-occurring pollinator wasp species (Molbo et al., 2003(Molbo et al., , 2004Satler et al., 2022;Sutton et al., 2017). Importantly, no study has directly applied genomic tools that can reveal the presence of hybridization and introgression across both the host figs and their associated pollinator species comprising an entire local community. Therefore, it is unclear the degree to which pollinator sharing, host sharing, host switching, or some combination of the three, generate hybridization and introgression in the host figs or pollinator wasps, and what these different ecological interactions mean for the strengthening or weakening of species boundaries.
Here, we use genome-wide sequence data to test for hybridization and introgression in all species comprising a community of . We ask whether hybridization and introgression have been operating within either the host figs, pollinator wasps, or both, and if so, whether currently observed levels of species specificity (or lack thereof) explain these evolutionary processes. For the pollinating wasps, we find no evidence for hybridization or introgression among any of the Tetrapus species.
We do identify, however, at least two pollinator species that are consistently pollinating and reproducing in more than one host fig species. For the host figs, we find that host species that frequently share the same pollinator species exhibit evidence of recent hybridization events (genetically identified F1 hybrids that are morpholog-   (Croat, 1978) and comprise roughly one-quarter of the 22 described species of neotropical free-standing figs (Berg, 2006). The individual trees sampled were located along the shoreline of the Rio Chagres, Lake Gatun, and adjacent seasonally dry tropical forest, typically a few hundred meters to DNA was extracted using a Qiagen DNeasy blood and tissue kit (Qiagen Inc.). Illumina libraries were generated using a KAPA Hyper prep kit with custom indices as described in Glenn et al. (2019).
Samples were sheared using a Covaris sonicator to an average size of 400-500 base pairs. Following library prep, we grouped samples into sets of eight and conducted probe hybridization targeting 2590 ultraconserved element (UCE) loci using the hymenopteran probe v2 set of Branstetter et al. (2017). Size distributions were assessed with a Bioanalyzer, and samples were grouped in equimolar concentrations for sequencing. We sequenced libraries on an Illumina sequencer (HiSeq 3000 and HiSeq 4000) generating 150 bp pairedend reads.
We also generated a phased data set for the UCE loci following the outline of Andermann et al. (2018). Briefly, we aligned our cleaned sequence reads back to aligned loci with BWA-MEM as implemented in bwa v0.7.17 (Li & Durbin, 2010). Data were phased using the phase command in samtools v1.9 (Li et al., 2009) resulting in two alleles per individual per locus. Phased data sets were cleaned as outlined above, and loci with a minimum of 70% of individuals were once again saved for downstream analysis.

| Wasp population structure and hybridization
We used two approaches to test for species boundaries and hybridization in the fig wasps. First, we used principal components analysis (PCA) to determine species groupings. In the absence of recent hybridization and introgression, individuals are expected to form distinct clusters corresponding to species. Hybrid individuals (F1s or subsequent backcrosses), by contrast, are expected to be located equidistant between species clusters while limited introgression is expected to result in intermixed species clusters. We conducted the PCA in R v3.6.3 (R Core Team, 2018) using the dudi.pca command in adegenet v2.1.3 (Jombart, 2008). For our input data set, we subsampled a single biallelic SNP per UCE locus. Missing data were replaced with the global mean allele frequency for that SNP. Because of a high amount of missing data for one individual (FW514), this sample was excluded from the analysis. We visualized results by plotting the first two principal component axes of variation.
Second, we explicitly tested for hybridization and introgression using the population graph and admixture approach as implemented in TreeMix v1.13 (Pickrell & Pritchard, 2012). TreeMix estimates a population graph with an a priori number of migration events between lineages, here species. This method allowed us to test whether a model with migration between wasp species is a better fit to the data than a strictly bifurcating model without migration. We analyzed the data in TreeMix using zero to three interspecific migration events and used the proportion of variance explained by the model to determine the optimal number of migration events. We used the phased data set and randomly subsampled a single biallelic SNP per locus for this analysis. If hybridization and introgression are not operating within this system, we would expect negligible improvement to the model as we add migration events.
This approach allowed us to determine whether individuals sam-

| Wasp mitochondrial DNA
We wanted to compare phylogenetic patterns between the nuclear (UCE) and mitochondrial genomes. If individuals belonged to different clades between the species tree (estimated with nuclear UCE data) and the mitochondrial gene tree, the cytonuclear discordance could be explained by interspecific hybridization and introgression. To generate mtDNA data from our samples, we followed the outline of Satler et al. (2022). Briefly, we used NOVOPlasty v4.3.1 (Dierckxsens et al., 2017) to identify mitochondrial reads and generated haplotypes from off-target reads present in the UCE sequencing files. We used a COI sequence from a Tetrapus species (AY148155) as our seed sequence. After recovering mtDNA haplotypes, we aligned these data with MAFFT v7.471 and trimmed the matrix to match the length of the seed sequence to minimize missing data. We then estimated an ML gene tree with IQ-TREE, used ModelFinder with BIC to select the substitution model of best fit, and generated 1000 bootstrap replicates with the ultrafast bootstrap approximation. Finally, we tested for cytonuclear discordance by comparing the species compositions recovered with the mitochondrial DNA with those recovered with the nuclear (UCE) DNA. (with the exception of F. tonduzii). Genomic DNA was extracted, and COI sequences were generated and aligned following the methods described in Marussich and Machado (2007). These COI mtDNA data provide sufficient information for confirming species identification and for generating host association frequencies.

| Host associations
This independent COI data set was generated from samples collected between February 1997 and May 2005, earlier than the samples collected here for UCE sequencing. Because of potential pollinator turnover, we needed to confirm that the wasp species sampled for the newer UCE data set and the older COI data set were the same. To confirm wasp species identity between the two data sets, we combined the older COI data set of 201 wasps with our newer COI data set recovered from NOVOPlasty, resulting in a total of 257 sequences. We realigned these data with MAFFT, then estimated an ML gene tree with IQ-TREE (as outlined above) to test for continuity of wasp species over these two time periods. Since the two data sets resulted in congruent wasp species inference, we used this information from the combined COI gene tree to determine host-pollinator association frequencies.  (Table S2). This included five trees that were putatively identified as F. glabrata × maxima hybrids and one tree identified as an F. insipida × maxima hybrid. Our initial hybrid identifications were based on intermediate leaf morphology and growth form. Genomic DNA was extracted using a modified CTAB protocol (Doyle & Doyle, 1987). Extractions were sent to Floragenex Inc. for restriction site-associated DNA (RAD) library preparation.

| Fig sampling and sequencing
Single-end RAD libraries were generated with the PstI restriction enzyme following the standard protocol (Baird et al., 2008). Libraries were sequenced on an Illumina HiSeq 3000 using 100 bp single-end sequencing.
DNA sequence reads were processed with ipyrad v0.9.62 (Eaton & Overcast, 2020). No mismatches were allowed in barcodes when demultiplexing samples, with strict filtering used for removing any adapter contamination. Up to five low-quality base calls were allowed in a read. We used a clustering threshold of 85% sequence similarity when assembling reads into loci within species. Within individuals, we allowed up to 5% Ns and 5% heterozygous sites per locus. Alleles were clustered across individuals using an 85% sequence similarity threshold. For clustered loci, we allowed up to 20% SNPs, up to 20% heterozygous sites, and up to eight total indels. Data sets were output varying the amount of missing data depending on the downstream application.

| Fig population structure and hybridization
We conducted a PCA to determine whether fig species cluster in multivariate space and to identify any potential hybrid individuals located between species clusters. To reduce the potential negative effects of missing data, we output loci sampled from at least 90% of individuals and, from these data, selected a single SNP per locus.
We conducted the PCA in R as described above for the wasps and visualized the first two axes of variation.
Based on the PCA of SNP data (see 'Section 3'), we identified six individual fig trees as putative recent hybrids between F. glabrata and F. maxima and identified one individual as a putative hybrid between F. maxima and F. yoponensis. To further evaluate hybridization, we used fastSTRUCTURE (Raj et al., 2014) to estimate population membership between pure species and potential hybrids. fast-STRUCTURE uses a variational Bayesian framework to approximate the structure (Pritchard et al., 2000) model for estimating population membership. We created two data sets for fastSTRUCTURE, one that included individuals of F. glabrata, F. maxima, and putative hybrids between the two species, and one that included individuals of F. maxima, F. yoponensis, and their putative hybrid. We ran each data set under a two-population model (K = 2). If individuals represent hybrids, we would expect them to show population membership in both clusters in their respective analyses. For fastSTRUCTURE, we used unlinked SNPs present in at least 50% of sampled individuals.
Of the hybrid figs identified in our community, we next wanted to know whether they were first-generation hybrids (F1) or first-generation backcrosses (BC1) to either parental species.
Individuals backcrossing to a parental species are of interest because they provide a mechanism for the introgression of genetic material between species. To address this question, we used snapclust (Beugin et al., 2018) implemented in the R package adegenet. Snapclust is a maximum likelihood approach for assigning individuals to clusters, including both pure species and hybrids. Specifically, snapclust can model F1 hybrids as well as first-and second-generation backcrosses, allowing the identification of the specific generation for a sampled hybrid. We partitioned samples into the same two data sets as described for fastSTRUCTURE and used the same genomic data as described for the PCA.

| Fig phylogenetics and introgression
To estimate the phylogenetic relationships among these freestanding figs, we used the coalescent-based approach SVDQuartets (Chifman & Kubatko, 2014) as implemented in PAUP* v4.0a168 (Swofford, 2003). SVDQuartets uses site patterns in the nucleotide data to estimate a phylogeny under the coalescent model. Because hybridization and introgression are not accounted for in this method, we removed hybrid individuals identified by the above analyses from the unlinked SNP data set. Individuals were assigned to species a priori, all quartets were evaluated, and nodal support values were To estimate the phylogenetic network, we first generated concordance factors from our unlinked SNP data sampled from pure species (as in SVDQuartets, hybrids were removed) as outlined in Olave and Meyer (2020). Specifically, we used the R function SNPs2CF (www.github.com/melis aolav e/SNPs2CF), sampled 100 alleles per species quartet (n.quartet = 100), and generated 100 bootstrap replicates. Using these concordance factors as input to SNaQ, we estimated networks allowing a maximum (hmax) number of between zero and three hybrid edges (i.e., hybridization events), doing 10 runs per analysis. For the network analysis with zero hybrid edges, we used the SVDQuartets species tree as the starting tree. For each subsequent analysis, we used the hmax -1 network as our starting network. Pseudolikelihoods were compared across runs with different numbers of hybrid edges to estimate the network model with the best support. We then used the best model to estimate 100 bootstrap replicates to generate support for the presence of the hybrid edge(s) indicating historical gene flow and introgression between species.

| Wasp population structure and hybridization
We generated 180,616,361 raw reads for the 57 fig wasp pollinator samples representing six pollinator species (Table S3). Individuals had on average 3,168,708 (±1,603,597) raw reads. Following data processing, individuals had on average 154,359 (±74,989) contigs with an average length of 379 (±211) base pairs. We generated contig data from 2248 total UCE loci, with each individual being sequenced at 1423 (±151) loci on average (Table S3).
All six wasp species are well differentiated with little intraspecific variation in PCA space (Figure 1). PC1 and PC2 explained 44.19% and 20.11% of the variation, respectively. Because we see tight clusters of individuals within species, and no spread of individuals between species, the PCA supports the pollinators as genetically distinct species with no recent hybridization or introgression.
TreeMix results indicate that a model without hybridization provides the best fit to the data ( Table 1). Although we estimated models with up to three admixture edges, there was essentially no increase in the proportion of variation explained by the model with the addition of admixture edges. Thus, the PCA and TreeMix results both show the pollinator species to be genetically distinct and reproductively isolated, with no evidence of hybridization.

| Wasp phylogenetics
In agreement with the PCA and TreeMix results, a concatenated ML tree of the UCE data detected six well-defined and reciprocally monophyletic species (Figure 1). The tree was characterized by low intraspecific divergence and high interspecific divergence. All interspecific nodes have 100 bootstrap support values.

| Wasp mitochondrial DNA
We were able to generate COI data from 56 of the 57 wasp individuals. After alignment and edge trimming, our COI matrix was composed of 816 base pairs. In agreement with the nuclear UCE phylogeny, our COI gene tree recovered six clades, all supported with bootstrap values of 100 ( Figure S1). Once again, these species were characterized by low intraspecific divergence and high interspecific divergence. Results recovered with mtDNA data mirrored those recovered with UCE data, showing no cytonuclear discordance among species compositions between the two data sets.

| Host associations
Species compositions were identical between the UCE nuclear phylogeny and the mitochondrial gene tree. Pollinator species have also been consistent temporally in this community between the previously generated COI data and the current UCE data set. We therefore combined our current sampling with the 201 samples directly sequenced for COI to quantify associations between host figs and pollinator wasps in this community ( Figure 2). Additionally, because we identified six individual figs as being recent hybrids between F. glabrata and F. maxima (see below), we grouped these individuals as   Figure 2). Additionally, one pollinator species (T. sp. 3) is primarily associated with F. maxima, but individuals were also frequently sampled from F. glabrata × maxima hybrids. We thus have two cases of host sharing: (T. sp. 1) and (T. sp.

2) associated with F. glabrata, and (T. sp. 2) and (T. sp. 3) associated
with maxima and also with F. glabrata × maxima hybrids ( Figure 2). In sum, four of the pollinator species are host species-specific while two are associated with multiple hosts.

| Fig population structure and hybridization
We  For each analysis, the two pure species were recovered as distinct with the putative hybrid samples showing admixed ancestry consistent with recent hybridization (Figure 3). To estimate whether these recent hybrids are first-generation hybrids (F1s) or first-or secondgeneration backcrosses, we analyzed the two data sets in snapclust.
The result for the F. glabrata and F. maxima data set estimated that five individuals are F1 hybrids and one individual is a first-generation backcross (BC1) to F. glabrata ( Figure S2)

| Pollinator sharing leads to hybridization in the figs
Pollinating fig wasps often appear to be highly species-specific in their associations with host figs (Bronstein, 1987;Moe et al., 2011;Ramírez, 1970;Satler et al., 2022). The degree to which this is true constrains the opportunities for hybridization (and subsequent introgression) in both lineages. Strict species specificity by the pollinators necessarily limits potential interspecific pollination in their hosts. For host figs to have opportunities for hybridization, several barriers must be overcome by pollinators. A wasp bearing heterospecific pollen must disperse and recognize a species different from the one in which she developed (Compton, 1990;Nason et al., 1996). Once inside the syconium, the wasp must be able to put viable pollen grains in contact with the stigmatic surfaces of receptive flowers, so that they can potentially produce viable F1 seeds. glabrata and F. maxima. We then used a phylogenetic network approach restricting our data set to only individuals representing pure species. Using this method, the best model placed a hybrid edge indicating introgression between F. glabrata and F. maxima (Figure 3).

However
These results demonstrate that hybridization and introgression be-  Van Noort et al., 1989;Wang, Yang, et al., 2021;Ware et al., 1993),  -Lund et al., 2017;Machado et al., 2005;Renoult et al., 2009;Wilde et al., 2020), there is little evidence these processes affect fig wasp pollinators (see Molbo et al., 2003Molbo et al., , 2004Satler et al., 2022;Sutton et al., 2017). This suggests that the processes governing reproductive isolation op- And because exceptions to the one-to-one fig-to-pollinator association are becoming more evident with increased taxon and genetic sampling (Darwell et al., 2014;Souto-Vilarós et al., 2019;Su et al., 2022;Sutton et al., 2017;Yang et al., 2015;Yu et al., 2019), coupled with a growing appreciation of the role of host switching in these systems (Satler et al., 2019), it is probable that hybridization and introgression in the host figs are more widespread than previously realized. Indeed, an increase in genomic data sets has led to an increase in the detection of hybridization events across the tree of life (Taylor & Larson, 2019). We suggest this will also be the case for figs as more systems are explored with genome-scale data and are explicitly tested for hybridization and introgression .
In our community of Tetrapus wasps, species are genetically distinct and are highly divergent (Figure 1). This result is consis-  Although hybridization has long been considered to be more prevalent in plants than in animals (Stebbins, 1959), there is a growing appreciation of the role hybridization has played in the animal tree of life (Mallet et al., 2016;Taylor & Larson, 2019). This is becoming more apparent as access to genomic data sets and newly

ACK N OWLED G M ENTS
We thank Juan C Penagos Zuluaga, Brian Park, Finn Piatscheck, and Jose Lopez for extracting DNA from the host figs. We thank the Heath Lab for helpful comments and discussion that improved the manuscript. Computational resources were provided by ResearchIT and the College of Liberal Arts and Sciences at Iowa State University.
The Smithsonian Tropical Research Institute provided the research facilities and intellectual community that helped make this study possible.

This work was supported by a grant from the National Science
Foundation [DEB-1556853] to JDN, TAH, and AEH and was supported by Yale University (KCJ).

CO N FLI C T O F I NTE R E S T
There is no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw sequence data are available from the NCBI Sequence Read Archive (SRA) under BioProject ID: PRJNA904825 (BioSample accessions: SAMN31644000-SAMN31644086). NCBI BioSample accession numbers are found in Table S1 for individual wasps and in Table S2 for individual figs. COI mtDNA pollinator barcode data are available at GenBank under accession numbers OP870658-OP870858. All data sets and custom scripts are available on Dryad (https://doi.org/10.5061/dryad.t1g1j wt60).