1. Introduction

Hybridization plays an important role in biological evolution (Mallet, 2007; Rieseberg and Willis, 2007; Abbott et al., 2013; Arnold, 2015), particularly through reorganizing standing genetic variation upon which natural selection may sculpt new species (Grant and Grant, 1992; Seehausen, 2004). However, successful hybrid speciation is likely curtailed as any newly forming "hopeful monster" (sensu Goldschmidt, 1933) of hybridization may have lower fitness than its parental species within the parental distribution range. Thus, it is not clear under what ecological conditions hybridization may be favored. Consistent reports of hybridization in contact zones suggests that many distinct but closely related taxa remain reproductively isolated primarily due to their occupancy of alternate geographic ranges (e.g., Abbott et al., 2013). However, consideration of distinct taxa with extensive contact zones may be more illuminating as we may assume that historic selection has imposed significant barriers to gene flow (e.g., Coyne and Orr, 2004; Hey, 2006) permitting reproductively isolated sympatry. Thus, evaluation of hybridization across species pairs occupying geographically extensive contact ranges may give greater insight as to when barriers to gene flow are likely weakened, and subsequently, where a significant proportion of diversification likely occurs.

Hybridization between closely related species in contact zones (e.g., Mebert, 2008; Garroway et al., 2010; Chunco, 2014; Wachowiak et al., 2016; Kindler et al., 2017) often results from secondary contact after post-glacial range shifts (e.g., Liao et al., 2019; Porto-Hannes et al., 2021). Meanwhile, contact zones are often at range peripheries in marginal habitats where hybridization is hypothesized to occur more frequently (Cruzan and Arnold, 1993; Strelkov et al., 2007; Kawecki, 2008; Hellberg et al., 2016; Chhatre et al., 2018; De Cahsan et al., 2021), because admixed variants may exhibit greater relative fitnesses at locations where parental stock likely endure a more precarious adaptive foothold (Bridle and Vines, 2007). In general, these reports suggest that the tendency to hybridize varies according to ecological conditions and adaptive opportunity (Abbott, 2017). Most data are generated from localized studies with homogeneous conditions and a strict appraisal of ecological influence is not forthcoming. Nevertheless, some evidence suggests that species pairs remain distinct in more benign conditions while introgressing in more marginal habitats (Hellberg et al., 2016; McCauley et al., 2019). Unsurprisingly, there is also evidence that climate warming drives the movement of hybridization zones. For example, in North America the hybrid zones of salamanders (Walls, 2009), and crickets (Britch et al., 2001) have both undergone northwards shifts. A greater understanding of how ecological conditions impact patterns of hybridization may yield greater insight as to how future biodiversity may flux under current climate change models (Chunco, 2014; De Cahsan et al., 2021).

Incidences of hybridization are unevenly distributed across the vascular plants and varies with factors such as life history, breeding system, environmental disturbance, genetic predisposition, and pollination syndrome (e.g., Grant, 1981; Ellstrand et al., 1996; Whitney et al., 2010; Abbott, 2017; Mitchell et al., 2019). It is predicted to occur less frequently among obligately insect-pollinated plants such as pollination mutualisms: e.g., Phyllantheae–Epicephala (Kato et al., 2003), Yucca–yucca moths (Pellmyr, 2003) and fig trees–fig wasps (Cruaud et al., 2012). This is likely due to a tendency for stronger pre-zygotic barriers to develop driven by fine-tuning of co-evolved interacting traits. Host plants are typically pollinated by a single species which reproduces by laying eggs in flowers (Janzen, 1979; Kato et al., 2003; Pellmyr, 2003). The reciprocal adaptations between interacting partners in reproductive phenology (Patel, 1996; Patel and Hossaert-McKey, 2000; Peng et al., 2010), chemical attraction (Ware et al., 1993; Grison-Pigé et al., 2002; Hossaert-McKey et al., 2010; Souto-Vilarós et al., 2018), partner identification and physical compatibility (van Noort and Compton, 1996; Kato et al., 2003), augment species specificity and predicate much stronger pre-zygotic barriers that will discourage hybridization (Jousselin et al., 2003; Pellmyr, 2003; Svensson et al., 2010; Cruaud et al., 2012).

Pollinator species that shift hosts or visit atypical hosts (hereafter referred to as "heterospecific visitation") are one of the main ecological processes that break patterns of species specificity between tightly bound interacting species and may permit host plant hybridization (Machado et al., 2005; Wang et al., 2021; Su et al., 2022). In particular, anthropogenic climate change can cause temporal and spatial mismatches between plants and interacting pollinators (Memmott et al., 2007; Schweiger et al., 2008; Hegland et al., 2009; Gérard et al., 2020), facilitating heterospecific visitation.

Ficus (Moraceae) species are well-known for their intricate, 75 Ma relationship with co-evolved pollinator wasps (Hymenoptera, Chalcidoidea, Agaonidae) (Herre et al., 2008; Cruaud et al., 2012) and regularly follow the one-to-one host-pollinating species rule (Cook and Rasplus, 2003; Cook and Segar, 2010). Given the extreme species-specificity in fig pollination, fig hybridization was historically considered unlikely. However, heterospecific visitation and resulting pollinator sharing have been regularly detected during detailed case studies (e.g., Cornille et al., 2011; Liu et al., 2015a; Su et al., 2022; Yu et al., 2022), or inferred through incongruent fig-pollinator phylogenies (e.g., Machado et al., 2005; Renoult et al., 2009; Cruaud et al., 2012; Wang et al., 2021). Thus challenging the ubiquity of strict host-specificity. Furthermore, heterospecific visitations have been found to mainly occur among closely related fig species due to reduced chemical and physical barriers (Ware and Compton, 1994a; van Noort and Compton, 1996; Grison-Pigé et al., 2002; Renoult et al., 2009; Cornille et al., 2011; Cruaud et al., 2012; Ghana et al., 2015a; Wang et al., 2016).

The frequency of heterospecific visitation is also assumed to vary across sites. Availability of phenologically compatible host figs or pollinators varies temporally and spatially (Khadari et al., 1995; Harrison, 2000; Zhang et al., 2014; Liu et al., 2015a). Pollinator wasps are predicted to visit alternative hosts if their typical host is locally unavailable, especially if the wasp species is incapable of long distances dispersal. Janzen (1979) further hypothesized that heterospecific visitations are likely to occur on islands and in harsh environments where availability of host figs or pollinator wasps may be reduced (e.g., Parrish et al., 2003; Kusumi et al., 2012; Tsai et al., 2015). Marginal habitats are vulnerable to climate change, anthropogenic threats, and stochastic events, and populations are often sparse, fragmented, prone to local extinctions, poor pollinator performance and heterospecific visitations (González-Megías et al., 2005; Kawecki, 2008; Liu et al., 2015a; Chen et al., 2018; Fréjaville et al., 2020).

Heterospecific visitation permitting pollen transfer imparts genomic signatures of hybridization. Cases of natural and experimental hybridization among intra-subgenus fig species have yielded viable hybrid seeds (e.g., Ramírez, 1986, 1994; Parrish et al., 2003; Ghana et al., 2015a, Ghana et al., 2015b; Tsai et al., 2015; Wang et al., 2016), demonstrating fairly weak post-zygotic isolation between closely related figs. Cophylogenetic analyses have suggested hybridization as a major mechanism generating the tremendous extant fig diversity (Machado et al., 2005; Wang et al., 2021). Genetic hybridization is predicted to positively correlate with the frequency of heterospecific visitation and may be more likely to occur among marginal populations in mainland species. Nevertheless, as a tropical genus, most Ficus members are thermophilic and hygrophilous. Low latitude margins with increased temperature and humidity are expected to be more favorable than high latitude margins for their mutualisms and empirical research has confirmed the reduced fitnesses of figs and pollinator wasps in high latitude margins (Chen et al., 2018). Thus, enhanced hybridization signatures may be predicted at high latitude contact margins.

Ficus heterostyla and F. squamosa (section Sycocarpus) are closely related dioecious figs (Cruaud et al., 2012). They occur across geographic ranges from northern India to Vietnam, with contact zones from southern Thailand to southwest China (Fig. 1). Heterospecific colonization of F. squamosa figs by the nominal pollinator of F. heterostyla has been detected in southwest China, and experimental hybrid seeds can be produced bidirectionally (Liu et al., 2015a). A continuous contact belt with gradually reduced climatic suitability from low to high latitude margins provides a critical test for predicting environmental drivers of heterospecific visitation mediated introgressive hybridization. As such, we employed latitudinal sampling using nuclear microsatellites to detect patterns of clinal change in interspecific hybridization and introgression. To our knowledge, our study is the first to evaluate such patterns across a comprehensively sampled wide latitudinal contact range among closely related species. Thus, our data offers a controlled evaluation of whether ecological clines demonstrably alter tendency to hybridize. Focal species do show increased signatures of hybridization at northern climatic margins. Moreover, our data also indicate that range shifts caused by ongoing climate change is likely to increase or greatly alter incidences of hybridization and concomitantly influence future planetary biodiversity outcomes.

2. Materials and methods 2.1. Study species

Ficus heterostyla is a small dioecious tree growing up to 5(–8) m tall typically under the forest canopy from southwest China to southern Vietnam. Figs are located in rooting stolons near or under the soil (Berg and Chantarasuwan, 2007; Liu et al., 2015a). F. squamosa is a short dioecious shrub up to 2(–3) m tall with a creeping stem, typically found growing along riverbanks or near fast flowing streams in tropical forests and occurring from northern India to southern Thailand. Figs are located on branches near to water level or even below (Zhou and Gilbert, 2003; Berg et al., 2011; Pothasin et al., 2016) (Fig. S1). Despite distinct habitat dependencies, these two figs co-occur from southern Thailand to southwest China (Fig. 1) and grow close together in some sympatric areas sometimes only a few metres apart. In Xishuangbanna, southwest China, during months when there is a shortage of typical F. squamosa pollinators, the nominal pollinator of F. heterostyla, Ceratosolen spp., enters receptive figs of sympatric F. squamosa facilitated by complementary pattern of reproductive phenology.

2.2. Sample collection and DNA extraction

Sympatric populations of the two species were collected roughly along longitude 100°E, from 10°N of southern Thailand to 22°N of southwest China. To contrast patterns of interspecific gene flow between sympatric and allopatric populations, four Ficus heterostyla populations (H-tko, H-tpb, H-tch, H-tur) were also sampled from southeast Thailand where no F. squamosa are recorded. Three to 17 individuals for each population were collected at least 10 m apart. Healthy leaves were collected and dried with silica gel separately. Voucher specimens were deposited in the Herbarium of Xishuangbanna Tropical Botanical Garden (HITBC). In total, 176 individuals from 14 F. heterostyla populations and 143 individuals from 13 F. squamosa populations were sampled across southern China, Myanmar, Thailand, and Laos (Table 1; Fig. 1). No morphological intermediates were found.

2.3. Nuclear microsatellite amplification and genotyping

Total genomic DNA was extracted from silica-gel dried material using the Plant Genomic DNA Kit (Tiangen Biotech, Beijing, China). Thirteen microsatellite loci (3-N173, 4-101, 5-N108, 6-N104, 7-N245, 19-N530, 20-N291, 23-N724, 21-N197, 26-N180, 28-N247, 30-N457 and 32-N125) identified by Li et al. (2020) were amplified and genotyped according to described protocols. Post-PCR products were analyzed by capillary electrophoresis on an ABI 3730XL DNA analyzer (Applied Biosystems, Foster City, California, USA) with the GeneScan 500 ROX Size Standard. Microsatellite fragment sizes were determined using GeneMapper v.3.2.

2.4. Genetic diversity, population structure and spatial signatures of hybridization

Classical indices of genetic diversity, including mean number of alleles per locus (N_a), observed heterozygosity (H_O), and expected heterozygosity (H_E), were estimated using GenAIEx 6.5 (Peakall and Smouse, 2012). Genetic differentiation among population pairs of each species were evaluated by analysis of molecular variance (AMOVA) performed in Arlequin v.3.5 (Excoffier and Lischer, 2010) with 10,000 permutations.

To assess the effects of geographic distance (IBD; isolation by distance) and environmental conditions (IBE; isolation by environment) on genetic structure, Mantel and partial Mantel tests were conducted using the R package Vegan (Oksanen et al., 2017) based on Pearson's correlation coefficient method with 9999 permutations. The pairwise geographic distance matrix was generated using the geographic coordinates of populations. While pairwise genetic distances among populations were calculated using (F_ST/(1 − F_ST)). We collected 19 standard bioclimatic variables (for the period 1970–2000 at 30 arc-second resolution) of the 27 population sampling sites as environmental data from WorldClim v.2.1 (http://www.worldclim.org). These standardized environmental variables were converted into a matrix using PCA and the first five principal components were employed to calculate Euclidean distance-based dissimilarity in environmental conditions among populations by the 'dist' function in R.

To estimate inter-population genetic affinity of studied populations, principal coordinate analysis (PCoA) was conducted with GeneAIEx 6.5 based on euclidean distance. PCoA was also performed on the complete set of individual samples and a subset excluding individuals of cluster 1. Scatterplots were visualized in two-dimensional space for the first and second principal coordinates. To further understand patterns of clustering and genetic admixture between the two species, Bayesian clustering of individual genotypes was investigated using STRUCTURE v.2.3.4 (Pritchard et al., 2000). We employed a model with admixture, with a burn-in of 100,000 and run length of 1,000,000 iterations varying K from K = 1 to K = 10. For each value of K, ten runs were performed. The STRUCTURE HARVESTER online program (Earl and Vonholdt, 2012) was used to detect the optimal K value (Evanno et al., 2005). The posterior probability of genetic admixture (q) in STRUCTURE analysis describes the proportion of an individual genotype originating from each cluster. We examined geographic trends of q regressing population admixture proportions obtained from STRUCTURE analysis against latitudes under the assumption that admixture proportions should show enhanced trends towards high latitudes. Furthermore, we estimated the correlation between admixture proportions and each of the 19 bioclimatic variables to evaluate the potential for environmental factors to influence interspecific hybridization.

The program Newhybrids v.1.1 (Anderson and Thompson, 2002) was used to check the assignment of sampled individuals to six genotypic categories (pure Ficus heterostyla, pure F. squamosa, F₁ hybrid, F₂ hybrid, backcrosses to F. heterostyla and backcrosses to F. squamosa, respectively) based on posterior probabilities (PP). Newhybrids analyses were run independently five times without prior species information using Jeffreys-like priors and 500,000 MCMC sweeps after a burnin of 100,000. We used PP ≥ 0.50 to designate genotypic class. If PP < 0.50, the individual was not assigned to any category.

Limited dispersal distances among pollinator wasps is expected to favor heterospecific visitations, rather than search for distant typical hosts when local typical hosts are unavailable. Genetic similarity between hybridizing species should be greater at sympatric locations than for allopatric populations. Thus, we tested for IBD by correlating geographic distance with F_ST/(1 − F_ST) between population pairs of different species using the Vegan package in R. Further, we also examined the potential influence of environmental differences on genetic distance between populations of different species by conducting IBE tests.

2.5. Body size of pollinator wasps

Wasp body size generally correlates with fig size (Renoult et al., 2009). Ficus squamosa is significantly larger than F. heterostyla in both female and male receptive figs (Liu et al., 2015a). We compared body sizes of wasps sampled from the two fig species in order to provide insights regarding the direction of interspecific gene flow. The wasp samples for measuring body size were collected from population H-cyx for F. heterostyla and population S-cym for F. squamosa in Xishuangbanna. Male figs in male phase were sampled randomly and individually placed in mesh bags (200 × 200 mm). After the natural emergence of wasps from figs, ten female wasps from each fig were sampled randomly to measure the lengths of head, front femur and ovipositor, as well as head widths, under a microscope with a micrometer (Liu et al., 2011). Two-sample t test was used to compare paired differences in wasp body size.

3. Results 3.1. Genetic diversity

Across 13 loci, 78 and 73 microsatellite alleles were identified for 176 Ficus heterostyla and 143 F. squamosa individuals, respectively. The genetic diversity parameters for each population (Table 1) indicate an observed species-level H_O of 0.624 and H_E of 0.868 for F. heterostyla, and H_O of 0.495 and H_E of 0.707 for F. squamosa. At the population level, N_a ranged from 1.92 to 10.23 (mean 6.26); H_O ranged from 0.451 to 0.840 (mean 0.643); H_E ranged from 0.423 to 0.802 (mean 0.641) among F. heterostyla populations, while N_a ranged from 1.46 to 5.15 (mean 3.50); H_O ranged from 0.269 to 0.658 (mean 0.489); H_E ranged from 0.173 to 0.591 (0.442) among F. squamosa populations, respectively. Within F. heterostyla, populations collected from southeast Thailand (H-tko, H-tur, H-tpb, H-tch) exhibit lower N_a, H_O, H_E compared to populations in other locations. The results show that F. heterostyla has a higher level of genetic diversity than F. squamosa at both the species and population levels.

3.2. Population structure

STRUCTURE results indicate an optimal K value of three using the deltaK criterion (Fig. S2), rather than an expected value of two corresponding to the two studied species (Fig. 2a; Table S1). Membership proportion of collected samples at K = 2 is also shown in Fig. S3. At K = 3, geographical distributions (see Fig. 1) indicate two genetically divergent and geographically structured clusters among Ficus heterostyla populations: four populations from southeast Thailand formed cluster 1 while the ten remaining populations were assigned to cluster 2. A few F. heterostyla individuals presented genetic admixture between the two intraspecific clusters, suggesting a low level of inter-cluster gene flow. All 13 F. squamosa populations constituted cluster 3. After excluding the four southeast Thailand allopatric F. heterostyla populations, an optimal K value of two was determined; the pattern of genetic admixture between sympatric F. heterostyla and F. squamosa populations remained unchanged by this (Fig. 2b). PCoA results further corroborated the finding of two identified F. heterostyla clusters and supported the lack of a clear boundary defining clusters 2 and 3 (Fig. 3).

A significant correlation between genetic and geographical distances was observed in both Ficus heterostyla (r = 0.321, p = 0.015, Fig. 4a), and F. squamosa (r = 0.794, p < 0.001, Fig. 4b). The partial Mantel test also supported a significant pattern of IBD in F. heterostyla (r = 0.282, p = 0.032) and F. squamosa (r = 0.478, p < 0.001) when controlling the influence of environmental distance. After excluding the inferred hybrids indicated by STRUCTURE analyses (see "Genetic admixture analysis"), the remaining samples of F. heterostyla from cluster 2 (r = 0.495, p = 0.045, Fig. S4a) and F. squamosa (r = 0.759, p < 0.001, Fig. S4b) also revealed with significant patterns of IBD. Mantel correlations between genetic distance and environmental distance were nonsignificant in F. heterostyla (r = 0.160, p = 0.180, Fig. 4c), but was significant in F. squamosa (r = 0.734, p < 0.001, Fig. 4d). However, the significant pattern of IBE among F. squamosa populations was no longer supported (r = 0.194, p = 0.097) by a partial Mantel test controlling the influence of geographic distance. These results indicate that the spatial distribution of genetic variation in our two focal figs is primarily influenced by geographic distance, rather than environmental heterogeneity.

An AMOVA test revealed F_ST of 0.313 (p < 0.001), and 68.7% of genetic variation contained within populations, with 31.3% of genetic variation observed among F. heterostyla populations. 18.1% of genetic variation was contained between clusters 1 and 2 (F_CT = 0.181, p < 0.001), with 16.1% of genetic variation observed among populations within clusters. A slightly higher level of genetic differentiation was revealed among F. squamosa populations (F_ST = 0.341, p < 0.001), and 65.9% genetic variation was found within populations, with 34.1% of genetic variation among populations. When AMOVA analyses were assessed between F. heterostyla and F. squamosa, 14.0% of the genetic variation was accounted for by interspecific differences, lower than the differentiation between the two F. heterostyla clusters, suggesting interspecific gene flow (Table 2).

3.3. Genetic admixture analysis

Ficus heterostyla was posited as a relatively pure species but showing two distinctly intraspecific clusters (Fig. 2, Fig. 3), with an east-west geographic differentiation. The eastern cluster (i.e. cluster 1) is allopatric with F. squamosa, while western cluster (i.e. cluster 2) is sympatric with F. squamosa. In contrast, five F. squamosa populations, including S-cym, S-cyn, S-tph, S-lpn and S-tme, displayed weak to strong signatures of interspecific genetic admixture (Fig. 2). The proportion of admixture (q) is used to identify hybrids, with a threshold value of above or below 0.90 to classify each individual as purebred or hybrid, respectively (Vähä and Primmer, 2006; Starr et al., 2013; Li et al., 2021). In total, 47 of 319 individual samples were inferred as hybrids (Fig. 2a), with two of them morphologically classified as F. heterostyla in populations H-tla and H-mmk, while the remaining 45 inferred hybrids were all morphologically classified as F. squamosa, indicating asymmetrical bidirectional gene flow. STRUCTURE identified the two admixed individuals morphologically classified as F. heterostyla with q value of 0.17 hailing from F. squamosa, while the 45 admixtured individuals among F. squamosa populations have q values ranging from 0.23 to 0.63 hailing from F. heterostyla. When members of cluster 1 were excluded from STRUCTURE analyses, the two F. heterostyla individuals from population H-tta displaying genetic admixture between clusters 1 and 2 were inferred as hybrids (Fig. 2b). Newhybrids analyses showed that both inferred hybrid individuals by STRUCTURE in populations H-tla and H-mmk were assigned as F₂ individuals (PP ≥ 0.60). All other F. heterostyla samples were assigned to parental genotypes. The five inferred hybrids of population S-tph by STRUCTURE were assigned to F₂ (PP > 0.76), while the two inferred hybrids of population S-tme were assigned to backcrosses to F. squamosa (PP > 0.63). However, hybrid assignments to samples of populations S-cym, S-cyn and S-lpn were unstable. For example, only a single sample of population S-cyn was assigned as an F₂ hybrid when all 319 samples were included in the analyses. Other samples of these three populations were assigned to parental F. squamosa with PP values ranging between 0.52 and 0.99 (Fig. S5). However, if populations S-cym and S-lpn were excluded from the analyses, all samples of population S-cyn displayed a high posterior probability (PP > 0.95) of assignment to F₂ hybrids. We suspect that intraspecific differentiation at a level similar to or higher than interspecific differentiation (see the pairwise F_ST between populations in Table S2) may contribute to this instability, but future research is needed. Because of the instability of Newhybrids results, subsequent discussions around genetic admixture are mainly based on STRUCTURE analyses and pays less attention to genotype categories of sampled individuals.

PCoA (Fig. 5) performed on all 319 individuals further showed a clear separation between purebreds and the 47 inferred hybrids displaying varying degrees of kinship among purebreds. Results show a few inferred hybrids overlapped with the Ficus squamosa cluster 3 and F. heterostyla cluster 2 (Fig. 5). PCoA scatterplots of individuals featuring only in clusters 2 and 3 are presented in Fig. S6.

3.4. Spatial pattern of interspecific hybridization

Admixture proportions in Ficus squamosa populations increased significantly towards high latitudes (r = 0.699, p < 0.001, Fig. 6a), while this pattern was not observed for F. heterostyla (r = 0.089, p = 0.176, Fig. 6b). Among 19 bioclimatic variables, 17 correlated significantly with admixture proportions in F. squamosa populations and five variables presented significant correlations with admixture proportions in F. heterostyla (Table S3). Overall, admixture proportions showed significant positive correlations with temperature seasonality and precipitation, but a significant negative correlation with temperature.

Signatures of hybridization were significantly (r = 0.256, p < 0.001) higher among geographically proximate populations (Fig. 7a), likely indicating the importance of atypical host availability. There was no significant correlation between interspecific genetic admixture and environmental difference (r = 0.156, p = 0.204, Fig. 7b), further indicating that environmental heterogeneity may not influence propensity to hybridize.

3.5. Wasp body size

In total, 150 wasps from 15 figs of three male Ficus heterostyla trees and 110 wasps from 11 figs of five male F. squamosa trees were collected to measure wasp body sizes. Wasps from the two fig species show no significant differences in either head length (p = 0.158) or ovipositor length (p = 0.931). However, wasps of F. squamosa are significantly bigger than the wasps of F. heterostyla in both head width (p < 0.001) and front femur length (p = 0.001) (Table 3).

4. Discussion

While hybridization is considered to be a significant contributor to biodiversity patterns (Mallet, 2007; Rieseberg and Willis, 2007; Abbott et al., 2013), including speciation events (Grant and Grant, 1992; Seehausen, 2004), little is known about the ecological conditions in which it is favored. A theoretical factor inhibiting hybridized speciation is the likelihood of a novel hybrid demonstrating superior fitness to its progenitor species in the environmental conditions in which they are adapted to. However, it should be more regularly favored at range or climatic margins where the adaptive fitness of parent species may diminish, and admixed variants may stochastically exhibit higher relative fitnesses (e.g., Kawecki, 2008). Hybridization is often recorded in species' contact zones (e.g., Abbott et al., 2013) and there is some evidence for such differential outcomes at core versus marginal occupancies (e.g., Porites corals; Hellberg et al., 2016). Nevertheless, there is a lack of data systematically evaluating fine-grained patterns of species pairs along extensive geographic contact ranges, that may discern whether observed patterns arise mechanistically, rather than stochastically. Here we report that two closely related figs with a longitudinal contact belt across Southeast Asia demonstrate reproductive isolation across their core contact range before hybridization signatures increase at northern climatic margins. As well as providing insight to ecologically contingent diversity fluxes among co-evolving species, our data also suggest that predicted climate change scenarios will likely alter the outcomes of ecosystem hybridization dynamics, and potential contingent speciation events, over future decades.

Host-specificity in fig pollination appears extreme compared with most other plant–insect interactions (Cruaud et al., 2012). However, increasing discoveries of pollinator wasps visiting multiple host species provoke an interest in fig hybridization and its causal drivers. Cases of fig hybridization have focused on native figs and closely related introduced figs where there is an artificial breakdown of geographic barriers and an absence of pollinator wasps specific to an introduced species (Ramírez and Montero, 1988; Ware and Compton, 1992; Ramírez, 1994; Ghana et al., 2015a, Ghana et al., 2015b). Further, island species have also been studied as their pollinators are more likely to visit atypical hosts due to a reduction in typical host availability (Parrish et al., 2003; Kusumi et al., 2012; Tsai et al., 2015). The three fig species endemic to the isolated Ogasawara Islands, Ficus boninsimae, F. nishimurae and F. iidaiana, were inferred to originate from hybridization events (Kusumi et al., 2012) and the likely local extirpation of parental species due to the harsh environment. Human activity, habitat fragmentation, introduction of exotic species, and increasing anthropogenic climate change, could all increase the occurrence of partner mismatch and heterospecific visitation in fig-wasp mutualisms and other intimate plant–pollinator associations (Ware and Compton, 1992; Kiers et al., 2010; Schweiger et al., 2010; Hagen et al., 2012; Grass et al., 2018; Bernard et al., 2020; Gérard et al., 2020). Thus, tendency towards hybridization may be exacerbated (Muhlfeld et al., 2014).

4.1. Restricted intraspecific gene flow

Both Ficus heterostyla (F_ST = 0.313) and F. squamosa (F_ST = 0.341) showed high population differentiation, more so than that of the Asian dioecious species, F. hirta (F_ST = 0.038, Yu and Nason, 2013) and F. pumila (F_ST = 0.123, Liu et al., 2015b), and slightly higher than that of F. tikoua (F_ST = 0.281, Chen et al., 2011). Alongside significant patterns of isolation by distance (Fig. 4), this indicates restricted gene flow among our focal species. Short distance pollen dispersal among dioecious Ficus mediated by pollinator wasps is commonly considered to likely result from high population densities (Harrison, 2003; Wang et al., 2009) associated with low host growth forms (Chen et al., 2011). Similar to F. tikoua, the figs of F. heterostyla partially submerse in soil or vegetation debris and are unusually hidden from their pollinator wasps meaning volatile attractants are likely to be highly localized, potentially restricting gene flow (Chen et al., 2011). Moreover, F. heterostyla typically grows in the understorey of forests, and tendency to cluster in fragmented localized patches may further restrict gene flow. The lower levels of genetic diversity observed in the four southeast F. heterostyla populations included in cluster 1 may result from the more fragmented habitats of Thailand they occupy. Ecological niche modeling of F. heterostyla indicated that habitat loss and fragmentation have been ongoing since the Last Glacial Maximum and are projected to continue into future decades, especially in southeast Thailand (Fungjanthuek et al., 2022). The figs of F. squamosa are close to the ground or at water level (even below water) and restricted to riverine habitats (Liu et al., 2015a; Pothasin et al., 2016). This may often lead to isolated discontinuity further restricting wasp dispersal and gene flow. As suggested by Harrison (2001, 2003), the slow recovery of dioecious fig-pollinator populations after local extinctions likely explained by their limited dispersal. This is reminiscent of the conditions underlying atypical host colonization between naturally sympatric related figs as driven by local pollinator shortages (Parrish et al., 2003; Tsai et al., 2015). Limited dispersal abilities of nominally host-specific pollinators may select for atypical host colonization and resultant subsequent hybridization when locally receptive typical host figs are unavailable. Notably, temporary shortages of F. squamosa wasps have been recorded among Xishuangbanna F. squamosa populations (Liu et al., 2015a).

The lack of IBE pattern may reflect the passivity of wasp dispersal. Dispersal flights of tiny (2–3 mm) and short-lived (1–2 days) wasps mostly involves floating on air currents (Ware and Compton, 1994b, Ware and Compton, 1994c). After entering the air currents, wasps are unable to exert directional control, potentially meaning that environmental heterogeneity does not influence dispersal.

4.2. Asymmetric interspecific gene flow

A strikingly asymmetric interspecific gene flow pattern was detected between F. heterostyla and F. squamosa, similar to that often observed in other closely related fig pairs (Tsai et al., 2015; Wang et al., 2016). A much higher level of gene flow from F. heterostyla into F. squamosa was observed. Several factors contributing to this pattern are discussed as possible facilitating mechanisms.

4.2.1. Fig and pollinator size

The ostiole exerts strong selective pressure on pollinator wasps by its size and shape and thus forms an important physical barrier preventing the entrance of atypical wasps with poorly matched morphological traits (Janzen, 1979; Verkerke, 1989; van Noort and Compton, 1996). Ostiole size is generally positively correlated with fig diameter (Moe et al., 2011). F. squamosa is significantly broader than F. heterostyla in both female and male receptive figs (see Table 1 in Liu et al., 2015a) while the present results show that wasps of F. squamosa are significantly bigger than wasps of F. heterostyla (Table 3). Thus, smaller F. heterostyla wasps should more easily enter larger F. squamosa figs, while bigger F. squamosa wasps should find access to smaller F. heterostyla figs more restrictive (Renoult et al., 2009; Cook and Segar, 2010; Tsai et al., 2015), leading to asymmetric tendencies of heterospecific visitation and subsequent hybridization.

4.2.2. Fruiting position

Wasp dispersal ability appears strongly correlated with breeding system, population density, fruiting positions, and host growth form (Harrison, 2003; Wang et al., 2009; Chen et al., 2011). F. squamosa is a shrub bearing figs positioned along branches close to the ground while the figs of F. heterostyla grow along the ground with figs often partially covered (Liu et al., 2015a) (Fig. S1). It may be problematic for F. squamosa wasps to locate partially covered F. heterostyla figs among terrestrial substrate.

4.2.3. Complementary fruiting phenology

Pollinator-specificity breakdown may occur if a pollinator of a fig species is locally extinct or temporarily absent from an area as has been recorded among Xishuangbanna F. squamosa populations. Figs of F. squamosa are produced almost exclusively in cold dry months and there is often a gap between wasp emergence and the availability of receptive figs at Xishuangbanna. In contrast, both male and female F. heterostyla trees bear fruit year-round with a well-defined summer peak. Wasps emerge from F. heterostyla figs throughout the period when receptive F. squamosa figs are present (but F. squamosa pollinators are absent). These partially disjunct fruiting patterns may have facilitated the colonization of F. squamosa by the typical pollinator of F. heterostyla (Liu et al., 2015a). This situation is different from that of Ficus auriculata complex in which limited overlap of reproductive phenology prevented hybridization between the sympatric interfertile forms (Wei et al., 2014).

4.2.4. Local host abundance

Species with reduced local abundance may receive heterospecific pollen and are likely to produce a higher proportions of hybrid offspring than species with high local abundance (Burgess et al., 2005; Lepais et al., 2009). A contrasting abundance of host figs (Tsai et al., 2015) or pollinator wasps (Yu et al., 2022) may contribute to asymmetric interspecific gene flow. F. squamosa is found growing abundantly though restricted to stream sides (Pothasin et al., 2014), while F. heterostyla occurs under the forest canopy or along roadsides at low density. The pollinator wasps of F. heterostyla more likely encounter a local reduction in receptive figs and may be forced to enter figs of closely related F. squamosa, while the high local abundance of F. squamosa means its wasps should cycle within their regular host population, potentially resulting in asymmetric hybridization and introgression.

4.3. Enhanced hybridization signatures towards high latitudes

The present study reveals enhanced hybridization signatures between F. heterostyla and F. squamosa among four high latitude northern populations suggesting heterospecific visitation among these populations. The female pollinator wasps emerging from natal figs must find and enter receptive figs. Phenological co-adaptation ensures that partner species can typically locate each other in temporal space (Rafferty et al., 2015). However, the markedly different taxonomic forms typically comprising symbiotic species (e.g., plants and arthropods) may differentially respond to climatic perturbation particularly outside their core ranges where they are most highly adapted. A shortage of matched pollinators or host trees due to mismatched phenologies will increase the likelihood of heterospecific visitation. Reproductive phenology of fig plants is highly sensitive to climatic (Peng et al., 2010; Pothasin et al., 2016; Chen et al., 2018) and environmental change (Bain et al., 2014) and there are large seasonal variations in both rainfall and temperature at the northern margins of the Asian tropics, such as at Xishuangbanna (Peng et al., 2010; Chen et al., 2015, 2018; Liu et al., 2015a). Correspondingly, reproductive phenology of figs has large seasonal fluctuations correlated with temperature and rainfall at northern margins (Spencer et al., 1996; Peng et al., 2010; Pothasin et al., 2016) that reduce towards the equator (Corlett, 1987; Bronstein et al., 1990). Fig production in both F. heterostyla and F. squamosa exhibit marked seasonal variation in Xishuangbanna. A partial disruption in male and receptive phase fig production occurs among F. squamosa Xishuangbanna populations (Liu et al., 2015a), though this disruption is not typical of the plant further south (e.g., in Chiang Mai, northern Thailand; Pothasin et al., 2016). Furthermore, pollinator wasps are more vulnerable to climate change than their hosts. Low temperature and low humidity will greatly reduce the wasp fitness (Dunn et al., 2008; Warren et al., 2010; Chen et al., 2018). Thus, synchronization of mutualistic cycles is greater at lower latitudes, and the occurrence and frequency of heterospecific visitation should correspondingly reduce increasing the likelihood of a clinally heterogeneous pattern of hybridization signatures.

The high latitude tropics are generally the geographic and climatic margins of tropical fig species, often featuring less suitable habitat and resulting in reduced fitnesses (Chen et al., 2018). Compared with populations located nearer to equatorial regions, climatically marginal pollinator populations, similar to island populations (Kusumi et al., 2012), are more vulnerable to extinction or temporary absence, and/or a reduction of host availability. Recent evidence indicates that southwest China represents the northern boundary of suitable habitats for the two targeted Ficus species (Fungjanthuek et al., 2022). The present data corroborate these hypothesized predictions evinced by nuclear microsatellite hybridization peaking in northern marginal areas and decreasing towards low latitudes. In general, our data, examining co-evolved mutualist species predicted to rarely hybridize, support predictions that hybridization is adaptively favored at climatic margins where parental species may exhibit maladaptive performance. Moreover, our examination of natural history data identifies multiple intuitive drivers of this dynamic.

4.4. Maintenance of morphological divergence

Although hybridization and introgression may be responsible for reduced molecular divergence between marginal F. heterostyla and F. squamosa populations, no morphologically identifiable intermediates were observed during field survey despite hybrids being collected among purebred populations. This phenotypic constancy despite introgressive gene flow is also observed between F. pumila and F. thunbergii in Okinoshima Island in Japan (Tsai et al., 2015), where phenotypic differentiation between the species is maintained against the homogenizing effects of interspecific gene flow. This may be due to strong selection exerted by host-specific pollinators required to enter the ostiole (Starr et al., 2013; Yoder et al., 2013) and local adaptation to distinct environments (e.g., F. squamosa being restricted to riparian habitats; Fitzpatrick et al., 2015). Hybrids expressing only binary parental phenotypes could be considered a case of reduced phenotypic expression of genetic variance (de Casas et al., 2007).

4.5. Implications under predicted climate change scenarios

As climate change has been shown to shift hybridization zones (Britch et al., 2001; Walls, 2009; Chunco, 2014) it is likely that our findings mean we may expect altered and possibly augmented cases of hybridization events occurring at ecosystem scales. As such, as well as predicted extinction events under climate change models, future planetary warming will likely drive increased hybridization and introgression in some geographic regions which may even yield new species through a process of hybridization speciation (Grant and Grant, 1992).

5. Conclusions

We confirm hypothesized predictions that two closely related species sharing an extensive contact range only display signatures of hybridization at their climatic margins. Moreover, an examination of detailed natural history data surrounding the two focal species and their regional abiotic circumstances insinuates the mechanistic drivers of the observed patterns. We believe this is the first thorough geographic examination of this phenomenon and indicates that climatic margins may be hotspots for hybridization potentially leading to speciation events. This also has ramifications for predicting biodiversity fluxes in coming decades under current models of climate change. As climate change has been shown to shift hybrid zones we may expect altered and possibly augmented cases of hybridization events occurring at ecosystem scales.

Acknowledgments

We thank Da-Rong Yang, Jenjira Fungjanthuek, Pei Yang for helping with sample collection, Chun-Yang Xu and Yi-Xin Yan for DNA extraction, Gui-Xiang Liu for measuring the wasp size. This work was supported by the National Natural Science Foundation of China (31800313; 32261123001), Applied Basic Research Foundation of Yunnan Province (202301AT070378, 2019FB034), the "Light of West China" Program of the Chinese Academic of Sciences to J.-F. Huang, and Guangzhou Collaborative Innovation Center on Science-tech of Ecology and Landscape.

Author contributions

JFH and YQP designed this study and collected the samples; JFH generated and processed the data; JFH wrote the paper; CTD contributed to conceptual development of the manuscript and edited the paper.

Declaration of competing interest

The authors declare no conflict of interest.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2023.08.003.