Enhanced and asymmetric signatures of hybridization at climatic margins: Evidence from closely related dioecious fig species

Jian-Feng Huang ,Clive T.Darwell ,Yan-Qiong Peng *

a CAS Key Laboratory of Tropical Forest Ecology,Xishuangbanna Tropical Botanical Garden,Chinese Academy of Sciences,Mengla 666303,China

b National Center for Genetic Engineering and Biotechnology (BIOTEC),113 Thailand Science Park,Phahonyothin Road,Khlong Nueng,Khlong Luang,Pathum Thani 12120,Thailand

Keywords: Hybridization Heterospecific visitation Fig Pollinator wasp Asymmetric gene flow

ABSTRACT Hybridization plays a significant role in biological evolution.However,it is not clear whether ecological contingency differentially influences likelihood of hybridization,particularly at ecological margins where parental species may exhibit reduced fitnesses.Moreover,it is unknown whether future ecosystem change will increase the prevalence of hybridization.Ficus heterostyla and F.squamosa are closely related species co-distributed from southern Thailand to southwest China where hybridization,yielding viable seeds,has been documented.As a robust test of ecological factors driving hybridization,we investigated spatial hybridization signatures based on nuclear microsatellites from extensive population sampling across a widespread contact range.Both species showed high population differentiation and strong patterns of isolation by distance.Admixture estimates exposed asymmetric interspecific gene flow.Signatures of hybridization increase significantly towards higher latitude zones,peaking at the northern climatic margins.Geographic variation in reproductive phenology combined with ecologically challenging marginal habitats may promote this phenomenon.Our work is a first systematic evaluation of such patterns in a comprehensive,latitudinally-based clinal context,and indicates that tendency to hybridize appears strongly influenced by environmental conditions.Moreover,that future climate change scenarios will likely alter and possibly augment cases of hybridization at ecosystem scales.

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” (sensuGoldschmidt,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′e et al.,2002;Hossaert-McKey et al.,2010;Souto-Vilar′os 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′erard 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 figpollinator 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′e 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′alez-Megías et al.,2005;Kawecki,2008;Liu et al.,2015a;Chen et al.,2018;Fr′ejaville 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,b;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,mostFicusmembers 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 heterostylaandF.squamosa(sectionSycocarpus)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 ofF.squamosafigs by the nominal pollinator ofF.heterostylahas 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.

Fig.1.Schematic distribution of two targeted species and genetic structure of 27 studied populations.The capital letters“H”and“S”in the population codes are the abbreviation for Ficus heterostyla and F.squamosa,respectively.Pie size is proportional to the number of individuals sampled per population.The size of each color block of a pie chart represents the overall proportion of a given population in each of three clusters (K=3) as revealed by STRUCTURE analyses.

2.Materials and methods

2.1.Study species

Ficus heterostylais 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.squamosais 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 typicalF.squamosapollinators,the nominal pollinator ofF.heterostyla,Ceratosolenspp.,enters receptive figs of sympatricF.squamosafacilitated 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,fourFicus heterostylapopulations (H-tko,H-tpb,H-tch,H-tur) were also sampled from southeast Thailand where noF.squamosaare 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 14F.heterostylapopulations and 143 individuals from 13F.squamosapopulations were sampled across southern China,Myanmar,Thailand,and Laos (Table 1;Fig.1).No morphological intermediates were found.

Table 1 Detailed information of sampled populations of Ficus heterostyla and F.squamosa.

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 (Na),observed heterozygosity (HO),and expected heterozygosity(HE),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 (FST/(1 -FST)).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 varyingKfromK=1 toK=10.For each value ofK,ten runs were performed.The STRUCTURE HARVESTER online program (Earl and Vonholdt,2012) was used to detect the optimalKvalue (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 ofqregressing 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 (pureFicus heterostyla,pureF.squamosa,F1hybrid,F2hybrid,backcrosses toF.heterostylaand backcrosses toF.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 usedPP≥0.50 to designate genotypic class.IfPP<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 withFST/(1 -FST) 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 squamosais significantly larger thanF.heterostylain 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 forF.heterostylaand population Scym forF.squamosain 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-samplettest 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 176Ficus heterostylaand 143F.squamosaindividuals,respectively.The genetic diversity parameters for each population(Table 1) indicate an observed species-levelHOof 0.624 andHEof 0.868 forF.heterostyla,andHOof 0.495 andHEof 0.707 forF.squamosa.At the population level,Naranged from 1.92 to 10.23(mean 6.26);HOranged from 0.451 to 0.840 (mean 0.643);HEranged from 0.423 to 0.802 (mean 0.641) amongF.heterostylapopulations,whileNaranged from 1.46 to 5.15 (mean 3.50);HOranged from 0.269 to 0.658(mean 0.489);HEranged from 0.173 to 0.591(0.442)amongF.squamosapopulations,respectively.WithinF.heterostyla,populations collected from southeast Thailand (Htko,H-tur,H-tpb,H-tch) exhibit lowerNa,HO,HEcompared to populations in other locations.The results show thatF.heterostylahas a higher level of genetic diversity thanF.squamosaat both the species and population levels.

3.2.Population structure

STRUCTURE results indicate an optimalKvalue of three using the deltaKcriterion(Fig.S2),rather than an expected value of two corresponding to the two studied species (Fig.2a;Table S1).Membership proportion of collected samples atK=2 is also shown in Fig.S3.AtK=3,geographical distributions (see Fig.1) indicate two genetically divergent and geographically structured clusters amongFicus heterostylapopulations: four populations from southeast Thailand formed cluster 1 while the ten remaining populations were assigned to cluster 2.A fewF.heterostylaindividuals presented genetic admixture between the two intraspecific clusters,suggesting a low level of inter-cluster gene flow.All 13F.squamosapopulations constituted cluster 3.After excluding the four southeast Thailand allopatricF.heterostylapopulations,an optimalKvalue of two was determined;the pattern of genetic admixture between sympatricF.heterostylaandF.squamosapopulations remained unchanged by this(Fig.2b).PCoA results further corroborated the finding of two identifiedF.heterostylaclusters and supported the lack of a clear boundary defining clusters 2 and 3(Fig.3).

Fig.2.Bayesian clustering analysis of 27 (a) and 23 studied populations (b) based on 13 nuclear microsatellite loci,respectively.The populations are arranged from high to low latitudes within cluster.Barplots indicate estimated assignment of each individual into (a) three or (b) two K clusters.Solid black lines define boundaries between populations.

Fig.3.Two-dimensional scatter diagram based on principal coordinate analysis of genetic variation in the 27 studied populations.

A significant correlation between genetic and geographical distances was observed in bothFicus heterostyla(r=0.321,p=0.015,Fig.4a),andF.squamosa(r=0.794,p<0.001,Fig.4b).The partial Mantel test also supported a significant pattern of IBD inF.heterostyla(r=0.282,p=0.032) andF.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 ofF.heterostylafrom cluster 2 (r=0.495,p=0.045,Fig.S4a) andF.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 inF.heterostyla(r=0.160,p=0.180,Fig.4c),but was significant inF.squamosa(r=0.734,p<0.001,Fig.4d).However,the significant pattern of IBE amongF.squamosapopulations 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.

Fig.4.Regression of paired FST/(1 - FST) vs the geographic distance.Analyses were significant for nuclear microsatellite data in both F.heterostyla (a) and F.squamosa (b).

An AMOVA test revealedFSTof 0.313 (p<0.001),and 68.7% of genetic variation contained within populations,with 31.3% of genetic variation observed amongF.heterostylapopulations.18.1% of genetic variation was contained between clusters 1 and 2(FCT=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 amongF.squamosapopulations(FST=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 betweenF.heterostylaandF.squamosa,14.0% of the genetic variation was accounted for by interspecific differences,lower than the differentiation between the twoF.heterostylaclusters,suggesting interspecific gene flow(Table 2).

Table 2 Analysis of molecular variance (AMOVA) based on 13 microsatellites for Ficus heterostyla and F.squamosa.

3.3.Genetic admixture analysis

Ficus heterostylawas posited as a relatively pure species but showing two distinctly intraspecific clusters(Figs.2 and 3),with an east-west geographic differentiation.The eastern cluster (i.e.cluster 1) is allopatric withF.squamosa,while western cluster (i.e.cluster 2)is sympatric withF.squamosa.In contrast,fiveF.squamosapopulations,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¨ah¨a 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 asF.heterostylain populations H-tla and H-mmk,while the remaining 45 inferred hybrids were all morphologically classified asF.squamosa,indicating asymmetrical bidirectional gene flow.STRUCTURE identified the two admixed individuals morphologically classified asF.heterostylawithqvalue of 0.17 hailing fromF.squamosa,while the 45 admixtured individuals amongF.squamosapopulations haveqvalues ranging from 0.23 to 0.63 hailing fromF.heterostyla.When members of cluster 1 were excluded from STRUCTURE analyses,the twoF.heterostylaindividuals 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 F2individuals (PP≥0.60).All otherF.heterostylasamples were assigned to parental genotypes.The five inferred hybrids of population S-tph by STRUCTURE were assigned to F2(PP>0.76),while the two inferred hybrids of population S-tme were assigned to backcrosses toF.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 F2hybrid when all 319 samples were included in the analyses.Other samples of these three populations were assigned to parentalF.squamosawithPPvalues 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 F2hybrids.We suspect that intraspecific differentiation at a level similar to or higher than interspecific differentiation (see the pairwiseFSTbetween 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 theFicus squamosacluster 3 andF.heterostylacluster 2 (Fig.5).PCoA scatterplots of individuals featuring only in clusters 2 and 3 are presented in Fig.S6.

Fig.5.Two-dimensional scatter diagram based on principal coordinate analysis of genetic variation in the 319 studied individuals.Green dots,Ficus heterostyla individuals nested in cluster 1;light blue dots,F.heterostyla individuals nested in cluster 2;red dots,F.squamosa individuals nested in cluster 3;purple dots,inferred hybrids nested in cluster 2;yellow dots,inferred hybrids nested in cluster 3.

3.4.Spatial pattern of interspecific hybridization

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

Fig.6.Regression of admixture proportions against latitude for populations of Ficus heterostyla and F.squamosa.(a) Introgression of F.heterostyla alleles into F.squamosa.(b)Introgression of F.squamosa alleles into F.heterostyla.Admixture proportions obtained for each individual from STRUCTURE are represented by black circles.Higher admixture values represent populations with greater proportions of alleles derived from another species.

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.

Fig.7.Regression of pairwise FST/(1 - FST) between interspecific populations of Ficus heterostyla and F.squamosa vs pairwise geographic distance.

3.5.Wasp body size

In total,150 wasps from 15 figs of three maleFicus heterostylatrees and 110 wasps from 11 figs of five maleF.squamosatrees 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 ofF.squamosaare significantly bigger than the wasps ofF.heterostylain both head width(p<0.001)and front femur length(p=0.001)(Table 3).

Table 3 Body size of wasps sampled from male figs of Ficus heterostyla and F.squamosa.

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.,Poritescorals;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,b).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.nishimuraeandF.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′erard et al.,2020).Thus,tendency towards hybridization may be exacerbated (Muhlfeld et al.,2014).

4.1.Restricted intraspecific gene flow

BothFicus heterostyla(FST=0.313)andF.squamosa(FST=0.341)showed high population differentiation,more so than that of the Asian dioecious species,F.hirta(FST=0.038,Yu and Nason,2013)andF.pumila(FST=0.123,Liu et al.,2015b),and slightly higher than that ofF.tikoua(FST=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 dioeciousFicusmediated 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 toF.tikoua,the figs ofF.heterostylapartially 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.heterostylatypically 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 southeastF.heterostylapopulations included in cluster 1 may result from the more fragmented habitats of Thailand they occupy.Ecological niche modeling ofF.heterostylaindicated 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 ofF.squamosaare 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 ofF.squamosawasps have been recorded among XishuangbannaF.squamosapopulations (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,c).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 betweenF.heterostylaandF.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 fromF.heterostylaintoF.squamosawas 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.squamosais significantly broader thanF.heterostylain both female and male receptive figs(see Table 1 in Liu et al.,2015a)while the present results show that wasps ofF.squamosaare significantly bigger than wasps ofF.heterostyla(Table 3).Thus,smallerF.heterostylawasps should more easily enter largerF.squamosafigs,while biggerF.squamosawasps should find access to smallerF.heterostylafigs 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.squamosais a shrub bearing figs positioned along branches close to the ground while the figs ofF.heterostylagrow along the ground with figs often partially covered(Liu et al.,2015a)(Fig.S1).It may be problematic forF.squamosawasps to locate partially coveredF.heterostylafigs 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 XishuangbannaF.squamosapopulations.Figs ofF.squamosaare 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 femaleF.heterostylatrees bear fruit year-round with a well-defined summer peak.Wasps emerge fromF.heterostylafigs throughout the period when receptiveF.squamosafigs are present(butF.squamosapollinators are absent).These partially disjunct fruiting patterns may have facilitated the colonization ofF.squamosaby the typical pollinator ofF.heterostyla(Liu et al.,2015a).This situation is different from that ofFicus auriculatacomplex 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.squamosais found growing abundantly though restricted to stream sides(Pothasin et al.,2014),whileF.heterostylaoccurs under the forest canopy or along roadsides at low density.The pollinator wasps ofF.heterostylamore likely encounter a local reduction in receptive figs and may be forced to enter figs of closely relatedF.squamosa,while the high local abundance ofF.squamosameans 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 betweenF.heterostylaandF.squamosaamong 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 coadaptation 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 bothF.heterostylaandF.squamosaexhibit marked seasonal variation in Xishuangbanna.A partial disruption in male and receptive phase fig production occurs amongF.squamosaXishuangbanna 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 targetedFicusspecies (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 marginalF.heterostylaandF.squamosapopulations,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 betweenF.pumilaandF.thunbergiiin 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.squamosabeing 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.

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 competing interests.

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 (31 800313;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.

Appendix A.Supplementary data

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