Similar phylogeographic history in a fig species and its obligate pollinators forms parallel genetic structure

Abstract

Closely related and co-distributed species usually share a common phylogeographic history, but it remains unclear whether ecologically interacting species can respond synchronously to historical climate changes. Here, we focused on a fig–pollinator mutualism comprising Ficus pumila var. pumila and its obligate pollinators (morphospecies Wiebesia pumilae), and collected samples across most of their distribution ranges. We employed cytoplasmic DNA sequences and nuclear microsatellite loci to reveal the species composition within the pollinators and to test whether the two mutualists exhibited similar postglacial phylogeographic patterns. We identified three cryptic pollinator species, with two dominant cryptic species exhibiting parapatric distributions in the northern and southern parts of the plant’s range, respectively. Similar current spatial genetic structures were detected in the two dominant cryptic pollinator species and the host plant, with both showing eastern and western genetic clusters. Moreover, evidence for postglacial expansion was found for all three species, and their potential refugia during the Last Glacial Maximum were located in the eastern and western parts of their distribution ranges. These results suggest synchronous responses to historical climate changes. Our study demonstrates congruent phylogeographic patterns between obligate mutualists and highlights the role of biogeographic factors in shaping the current biodiversity across trophic levels.

摘要

薜荔与其传粉小蜂间相似的亲缘地理史形成平行的遗传结构

分布相同且紧密相关的物种通常具有相似的亲缘地理历史，然而相互作用物种是否对历史气候变化表现出同步响应仍不清晰。本研究聚焦于薜荔(Ficus pumila var. pumila)及其专性传粉小蜂(形态种：薜荔传粉小蜂Wiebesia pumilae)组成的互惠共生体系，从它们分布范围内广泛取样，采用胞质DNA序列和核微卫星标记分析传粉小蜂物种组成，并检验薜荔及其传粉小蜂是否呈现相似的冰期后亲缘地理格局。我们在传粉小蜂中发现了3个隐存种，其中2个优势隐存种分别位于薜荔的北部和南部分布区，且总体呈异域分布。这2个优势隐存种与其宿主植物具有相似的空间遗传结构，均可划分出东、西2个遗传类群。此外，薜荔与这2个传粉小蜂隐存种在分布范围的东部和西部均存在末次冰盛期的潜在避难所，并表现出冰期后扩张的信号。上述结果表明，薜荔与其专性传粉小蜂对历史气候变化的响应具有同步性。本研究揭示了专性互惠共生物种间一致的亲缘地理格局，强调了生物地理因素在塑造不同营养级现存生物多样性格局中的重要作用。

INTRODUCTION

Historical biogeographic factors can affect speciation, extinction and migration, and ecological interactions maintaining the coexistence of species (Bansal et al. 2022; Lu et al. 2018; Yin et al. 2023). In the subtropical and temperate zones of Eurasia, Pleistocene climatic oscillations, especially the Last Glacial Period, have profound effects on current species distribution (Hewitt 2000; Jiang et al. 2021, 2024; Tong et al. 2022). Many closely related and co-distributed species show concordant historical demographic dynamics and migration trajectories as they often have similar life history characteristics and might synchronously respond to past climate changes (Fu et al. 2024; Kreger et al. 2020; Marske et al. 2020). Such consistent phylogeographic patterns help to understand the formation of current biodiversity pattern at a single trophic level. However, ecosystems are composed of species linked by interspecific interactions across different trophic levels (Akesson et al. 2021; Lan et al. 2022; Thrall et al. 2007; Wang et al. 2023). The impacts of historical events on the current distribution of tightly interacting species with long coevolution history may be similar, causing consistent phylogeography (Garrick et al. 2013; Yu and Nason 2013). Nevertheless, this hypothesis has seldomly been tested.

There are only a few comparative phylogeographic studies involving interacting species at different trophic levels, but they exhibited distinct results. For example, the white campion (Silene latifolia) and its pollinating moth Hadena bicruris showed incongruent phylogeographic patterns (Magalhaes et al. 2011), but striking similarities in genetic structuring have been found in some parasites and their hosts (Bothma et al. 2021) and in some plants and their pollinating insects (Tian et al. 2015). Such an inconsistent phylogeographic patterns of associates may be explained by coevolutionary constraints, host shifts, disintegration of interactions during glacial periods and utilization of different refugia due to differences in adaptative strategies (Althoff et al. 2012; Maron et al. 2019). Taken together, these reasons reflect the importance of specificity in determining the degree of phylogeographic matching. When interacting species exhibit high specificity, pollinators have little chance to shift to other host plants, and thus they must follow their specific hosts and results in shared responses to historical climatic events. We thus hypothesize that specialized obligate mutualists are highly likely to share postglacial phylogeographic history.

Nevertheless, different current genetic structure has also been recorded in some obligate mutualisms (Deng et al. 2020; Liu et al. 2013). This is likely caused by distinct dispersal ability between mutualists, leading to asynchronous expansions from refugia (Liu et al. 2015; Yu et al. 2019). Obligatory pollination mutualisms are one of the most appropriate systems (Liu et al. 2013; Yu et al. 2019), but even in these mutualisms, the migration of pollinators only represent the pollen movements of their host plants. Furthermore, the causes of genetic structure may be more complicated when multiple obligate pollinators pollinate a single host plant species. Such complicated genetic structure has been documented in several systems, including figs and fig wasps (Deng et al. 2023; Yang et al. 2015; Yu et al. 2019), yuccas and yucca moths (Godsoe et al. 2008), globeflowers and globeflower flies (Després et al. 2002; Espíndola et al. 2014), and Silene and Hadena (Kephart et al. 2006). In this case, it is crucial to describe the composition of pollinators and their phylogenetic relationships before comparing phylogeographic patterns.

Fig trees (Ficus, Moraceae) and their pollinating wasps (Agaonidae, Hymenoptera) constitute a classic example of obligate mutualism (Weiblen 2002). In a fig–pollinator mutualism, the fig species produces enclosed inflorescences called figs which are pollinated by a unique or only a few pollinating wasp species, whose larvae can only grow in the galled ovules inside host figs (Wang et al. 2019; Yang et al. 2023). Therefore, these highly specialized mutualisms provide an ideal system for testing our hypothesis. Previous comparative phylogeographic studies on several fig–pollinator mutualisms have revealed conflicting results in the matching degree of phylogeographic patterns (Rodriguez et al. 2017; Tian et al. 2015; Yu et al. 2019), probably due to the differences in time scales and the complicated historical background in tropical zones. The recent phylogeographic history of fig–pollinator mutualisms distributed in subtropical and temperate zones has been poorly explored, despite the well-documented influences of the Last Glacial Period in these regions.

In this study, we chose Ficus pumila L. and its pollinating wasps as our study system. Ficus pumila is a functionally dioecious creeping fig tree, which is widely distributed in subtropical areas of East Asia (Liu et al. 2015). Wiebesia pumilae is the host-specific pollinator of F. pumila, but molecular studies have evidenced three cryptic species within this morphospecies (Chen et al. 2012; Zhang et al. 2022), and additional cryptic species may exist as these previous studies only conducted sampling within a limited area in the natural distribution range of F. pumila. Moreover, though genetic structure and demography of F. pumila and W. pumilae had been separately studied (Chen et al. 2012; Liu et al. 2015), a rigorous comparison (based on pairwise sampling in each location) of genetic structure and population dynamics between these two mutualists across their main natural distribution range has not yet been performed. In this study, we expanded our sampling range to cover most of their distribution areas in China and conducted a comparative phylogeographic study using both nuclear DNA markers and cytoplasmic DNA makers, to answer: 1) How many cryptic pollinator species of F. pumila are present in China and what is their distribution pattern? 2) Do F. pumila and the pollinators present congruent current spatial genetic structure? and 3) Do F. pumila and the pollinators share similar signatures of postglacial expansion?

MATERIALS AND METHODS

Study species and sample collection

Ficus pumila contains two varieties, F. pumila var. pumila and F. pumila var. awkeotsang, while the former has a wide distribution range, the latter occurs part distribution of the species in the southeastern coasts of mainland and Taiwan Island (Liu et al. 2015; Wang et al. 2013). The three cryptic pollinator species of F. pumila (Wiebesia sp. 1, 2 and 3) have been reported to have a largely parapatric distribution in southeast China (Chen et al. 2012). Wiebesia sp. 1 and Wiebesia sp. 2 are distributed in south of the Yangtze River in China, with the Nanling and Wuyi Mountains as the main boundary, and Wiebesia sp. 3 is an insular species mainly found on some islands nearby the southeastern coast of China (Tong et al. 2021). Sympatric occurrence of two cryptic species has been observed in the contact zones (Chen et al. 2012).

From 2006 to 2012, we collected leaves of the two varieties of F. pumila and their adult pollinating wasps. Our sampling locations covered most distribution range of F. pumila in China (Fig. 1a and b). Healthy leaves were collected from each individual of F. pumila, and mature females of the pollinator were collected from mature figs. In total, F. pumila and the pollinators were sampled from 82 and 76 locations within their natural distribution ranges (Supplementary Table S1). Note that F. pumila var. awkeotsang samples were only used for constructing phylogenetic trees and were excluded from all analyses of population genetics and phylogeography due to its limited distribution.

DNA extraction, microsatellite genotyping and cytoplasmic DNA segment sequencing

For F. pumila, the total genomic DNA was extracted using the Plant Genomic DNA Kit (Tiangen, Beijing, China). To study its genetic structure, eight nuclear microsatellite loci (FP9, FP38, FP102, FP134, FP213, FP540, FP556 and FP601) were amplified and genotyped as described by Zhang et al. (2011). In addition, three non-coding regions of chloroplast DNA (cpDNA) (trnS–trnG, atpF–atpH and trnC–ycf6) were amplified and bidirectionally sequenced as described by Liu et al. (2015).

For the pollinators, the total genomic DNA was extracted from each sampled adult female, and genotyping of 10 nuclear microsatellites was conducted as described by Liu et al. (2013). Mitochondrial DNA (mtDNA) gene COI was amplified using the primer pair Jerry⁄Pat bidirectionally sequenced as described by Chen et al. (2012).

Acquisition of genetic variations and identification of cryptic pollinator species

The sequences of cpDNA and mtDNA fragments were aligned separately using CLUSTALW implemented in MEGA v. 6.06 (Tamura et al. 2013). Then, the sequences from different loci were merged together for the following analyses. DnaSP v. 5.10 (Librado and Rozas 2009) was adopted to calculate the number of haplotypes (k), haplotype diversity (h) and nucleotide diversity (π). Haplotype network was analyzed using a statistical parsimony approach implemented in TCS 1.21 (Clement et al. 2000). All nodes were connected with a >95% probability cut-off, and K2P (Kimura-2-parameter) distances of sequences were calculated using MEGA v. 6.06.

To assign all pollinating wasps to different cryptic species, we constructed the phylogenetic trees using a series of COI sequences of fig wasps as outgroups (Supplementary Table S2). Two methods were used for reconstructing phylogenetic relationships. First, the maximum likelihood (ML) approach was conducted by PHYML v. 3.0 (Guindon et al. 2010). The nucleotide substitution model selected was TVM+I + G, with the starting tree constructed using the maximum parsimony method. The tree topology was estimated using the Subtree Pruning and Regrafting option, and node support was estimated with 100 bootstrap replicates. Second, Bayesian phylogenetic trees were constructed using software BEAST v. 1.8.2 (Drummond et al. 2012), using the nucleotide substitution model chosen by JMODELTEST with COI data partitioned into two parts (1 + 2, 3) based on codon positions. BEAST was run under a strict clock, with tree prior as Yule process and MCMC steps as 60 000 000 (sampled every 1000 generations). Log files were checked using TRACER v. 1.7 (Rambaut et al. 2018), ensuring the effective sample size for each parameter exceeded 200. The phylogenetic tree was summarized by TreeAnnotator v. 1.8.2 followed by the removal of 10% burn-in. Node ages of the Bayesian tree were estimated and calibrated by the divergence time between genera Pegoscapus and Tetrapus, with prior age distribution specified by lognormal distribution (mean = 0.9; stev = 0.9; offset = 16; 95% highest posterior density [HPD] = 16.42‒30.35 Ma).

For microsatellite loci, linkage disequilibrium (LD) among loci per population and deviation from the Hardy–Weinberg equilibrium (HWE) of each locus in each population were tested using GENEPOP v. 4.0 (Rousset 2008), with false discovery rate correction (Benjamini and Hochberg 1995). The presence of null alleles was detected with MICRO-CHECKER (van Oosterhout et al. 2004). We used the following indices calculated by FSTAT v. 2.9.3 (Goudet 1995) and GenALEx v. 6.5 (Peakall and Smouse 2012) to characterize within-population genetic diversity: mean number of alleles per locus (N_A), allelic richness per locus (A_R), and observed (H_O) and unbiased expected heterozygosities (H_E). We adopted linear regression to test the correlation of genetic diversity between each cryptic pollinator species and its corresponding host plant populations. Populations with fewer than five samples were excluded from all analyses at the population level.

Genetic differentiation and isolation by distance

To test whether differentiation between F. pumila var. pumila populations correlated with that of the pollinators, we calculated the genetic differentiation based on microsatellite loci. We estimated both the overall F_ST value and those between populations and calculated unbiased F_ST (F_ST(ENA)) using the ENA (excluding null alleles) method in FREENA (Chapuis and Estoup 2007). The significance of overall genetic differentiation values was tested using 1000 permutations, and P values were adjusted using the Bonferroni correction method.

To evaluate the relative role of pollen dispersal in shaping the genetic differentiation among F. pumila var. pumila populations, we calculated the pollen-to-seed migration ratio (r) for all populations and those pollinated only by Wiebesia sp. 1 or Wiebesia sp. 2, using estimations of F_ST based on both microsatellite loci and cpDNA sequences as proposed by Ennos (1994). We did not include the populations pollinated by Wiebesia sp. 3 because it appeared in very few locations (Supplementary Table S1).

Then, isolation by distance was evaluated for both plant and pollinator populations using Mantel tests, with F_ST/(1 − F_ST) or F_ST(ENA)/(1 − F_ST(ENA)) as the estimates of genetic distances, in R package ‘vegan’ (Oksanen et al. 2019), setting 10 000 permutations.

Genetic clustering

STRUCTURE v. 2.3.1 (Pritchard et al. 2000) was used to infer genetic structure for each cryptic pollinator species and F. pumila var. pumila based on microsatellite loci. We ran the admixture model with correlated frequencies, with 10 independent runs for each K (from 1 to 10) performed with 1 000 000 MCMC repetitions and a burn-in of 300 000. We used both the mean maximum estimated logarithm of posterior probability LnP(D) (Pritchard et al. 2000) and ΔK values (Evanno et al. 2005) to determine the optimal K. Assignment probabilities were averaged across all 10 replicate runs for the optimal K.

To uncover the genetic structure of F. pumila var. pumila populations based on cpDNA sequences, we constructed ML and Bayesian phylogenetic trees using the same processes as mentioned above. Node ages of the Bayesian tree were estimated based on two secondary calibration points setting normal distribution (Cruaud et al. 2012; Xu et al. 2011): the divergence time between F. pumila and other Ficus species (mean = 12.0 Ma; stev = 0.6; 95% HPD = 10.82‒13.18 Ma) and that between F. pumila and Ficus erecta (mean = 23 Ma; stev = 0.5; 95% HPD = 22.02‒23.98 Ma).

Population demographic history

The historical effective population size of each cryptic pollinator species was estimated using the Bayesian skyline plot (BSP) method implemented in BEAST. Nucleotide substitution model for each species was specified by JMODELTEST, with uncorrelated exponential relaxed clocks, tree prior as Bayesian skyline and MCMC of 25 000 000 (sampled every 25 000 generations). To convert the estimates scaled by mutation rate to calendar years, we applied the substitution rate of 1.9% per million years as proposed by Lin et al. (2008). We did not conduct this analysis for F. pumila due to very few haplotypes in its cpDNA sequences (see Results).

Mismatch distributions and neutrality tests (Slatkin and Hudson 1991) implemented in ARLEQUIN v. 3.5 (Excoffier and Lischer 2010) were used to infer the historical population expansion of all cryptic pollinator species and F. pumila var. pumila. Both sudden expansion and spatial expansion models (Excoffier 2004) were used. The goodness-of-fit of models was tested based on the sum of squared deviations (SSD) with 10 000 bootstrap replicates. For neutrality tests, Tajima’s D (Tajima 1989) and Fu’s FS (Fu 1997) were calculated and tested, using 10 000 simulations.

Potential distribution at the Last Glacial Maximum of the pollinators and F. pumila var. pumila

Species distribution modeling was explored to reconstruct the potential distribution of the pollinators and F. pumila var. pumila at the Last Glacial Maximum (LGM). For the two major pollinators, population locations from our field investigation were used as the presence records. The present occurrence information of F. pumila var. pumila were obtained from our field investigation, Chinese Virtual Herbarium (https://www.cvh.ac.cn/) and Global Biodiversity Information Facility (http://www.gbif.org/). To avoid spatial autocorrelation of records and artificial spatial clusters of observations, occurrence records were filtered by OccurrenceThinner (Verbruggen 2012). There were 42, 87 and 209 locations for Wiebesia sp. 1, Wiebesia sp. 2 and F. pumila var. pumila, respectively.

Current global climate data were downloaded from the WorldClim website (http://www.worldclim.org) at 2.5′ spatial resolution. Data from both Community Climate System Model (CCSM) (Collins et al. 2006) and Model for Interdisciplinary Research on Climate (MIROC) (Hasumi and Emori 2004) were used for estimating the distributions of F. pumila var. pumila and its pollinating wasps at LGM. Climate attributes of 19 bioclimatic variables at occurrence records were extracted using ArcGIS v.10. The Variance Inflation Factor (VIF) analysis was used to detect multicollinearity between bioclimatic variables using the vifstep function in R package usdm (Naimi et al. 2014), and bioclimatic variables with the VIF >10 were excluded. Eight bioclimatic variables (BIO2, BIO3, BIO8, BIO10, BIO13, BIO15, BIO18 and BIO19) were retained to model the current and LGM distribution of the three species.

The software MAXENT 3.3.3k (Phillips and Dudik 2008) was used to generate the potential distribution of the three species using default settings. Ten independent replicates of subsampling procedures with the random test percentage of 20% and the median probabilities were projected onto the current and LGM climate. The area under the receiver operating characteristic curve (AUC) was used as an indicator of the accuracy of model estimation.

RESULTS

Identification of cryptic pollinator species

For the pollinators, 665 mature females from 76 locations were sequenced for the COI gene. The lengths of COI gene sequences were 906–909 bp, and 123 haplotypes were detected. The results of phylogenetic analysis showed that all samples were assigned to three major clades (Fig. 1c), each of which contained the haplotypes of a single cryptic species from Chen et al. (2012). In addition, the mean K2P distances within each clade were 0.70%, 0.62%, and 0.71%, which were far lower than the K2P distances among clades (ranging from 6.55% to 12.13%). These results confirm the presence of three cryptic pollinator species, with 379, 230 and 56 samples belonging to Wiebesia sp. 1, sp. 2 and sp. 3, respectively (Supplementary Tables S3–S5). The estimated node ages of Bayesian tree showed that Wiebesia sp. 3 diverged from the common ancestors of Wiebesia sp. 1 and sp. 2 at about 7.56 Ma (95% HPD: 5.32–10.99 Ma), and the divergence time between Wiebesia sp. 1 and sp. 2 was estimated to be c. 4.50 Ma (95% HPD: 2.92–6.71 Ma) (Fig. 1c).

The three cryptic pollinators are mainly parapatric (Fig. 1a), with the Nanling–Wuyi Mountains as the boundary separating Wiebesia sp. 1 from Wiebesia sp. 2, and Wiebesia sp. 3 is mainly distributed in Zhoushan Archipelago, where it is sympatric with Wiebesia sp. 1. Wiebesia sp. 3 was also found in Ningde and Taiwan Island, where it pollinates F. pumila var. awkeotsang (Fig. 1c).

Genetic variations

For the pollinators, we genotyped a total of 1471 mature females using microsatellite loci, with no significant signals of LD. The clustering analysis allocated all samples into three distinct clusters, perfectly matching the distribution of the three cryptic pollinator species inferred from COI sequences (Supplementary Fig. S1a). We therefore identified the 802, 507 and 162 samples as Wiebesia sp. 1, sp. 2 and sp. 3, respectively. For F. pumila var. pumila, a total of 1682 samples were genotyped using microsatellite loci, with no significant signals of LD.

In all three cryptic pollinator species and F. pumila var. pumila, significant deviations from the HWE were only detected in a few locations for each microsatellite locus, and no evidence of the presence of null alleles were detected in all loci. Values of genetic diversity indices are presented in Supplementary Tables S1 and S3–S5. For the pollinators, there were no significant correlations between genetic diversity and latitude or longitude, and all indices (N_A, A_R, H_E, H_O) of F. pumila var. pumila decreased significantly with increasing latitude and longitude (Supplementary Fig. S2).

Using the combination of three cpDNA fragments (a total length of 2413 bp), we obtained a total of 26 haplotypes in F. pumila from 407 samples, with two haplotypes specifically associated with F. pumila var. awkeotsang (Fig. 1b; Supplementary Table S1). Haplotype diversity (h) and nucleotide diversity (π) for all locations were 0.667 and 0.00122, respectively, and most locations had only one haplotype (Fig. 1b; Supplementary Table S1).

Genetic differentiation among populations

For the pollinators, the overall genetic differentiation (using microsatellites) of Wiebesia sp. 1, sp. 2 and sp. 3 were 0.118 (F_ST) and 0.113 (F_ST(ENA)), 0.078 (F_ST) and 0.077 (F_ST(ENA)), and 0.089 (F_ST) and 0.085 (F_ST(ENA)) (all P values < 0.05). Significant positive correlations between pairwise genetic differentiation and geographic distance were detected in Wiebesia sp. 1 and Wiebesia sp. 2 (Supplementary Fig. S3a–f).

For F. pumila var. pumila, the overall genetic differentiation (using microsatellites) was 0.160 (F_ST) and 0.156 (F_ST(ENA)) (both P values < 0.05). Mantel tests revealed a significantly positive relationship between genetic differentiation and geographic distance (Supplementary Fig. S3g and h). The genetic differentiation based on cpDNA sequences was 0.833, and the pollen-to-seed dispersal ratio (r) was 24.19. In addition, the genetic differentiation (F_ST) among populations of F. pumila var. pumila pollinated only by Wiebesia sp. 1 was 0.169 (microsatellites) and 0.809 (cpDNA), leading to an r value of 18.83. The F_ST values among the plant populations pollinated by only Wiebesia sp. 2 were 0.140 (microsatellites) and 0.799 (cpDNA), with an r value of 22.42.

Genetic clustering

For the pollinators, we detected obvious genetic structure in Wiebesia sp. 1 and Wiebesia sp. 2. All samples of these two cryptic species were assigned into two genetic clusters, i.e. the western cluster containing most samples from inland and some coastal locations and the eastern cluster mainly comprising samples from some coastal locations and nearby islands (Fig. 1a; Supplementary Fig. S1b and c). There was no clear genetic clustering pattern for Wiebesia sp. 3 even for K = 2 (Supplementary Fig. S1d). For F. pumila var. pumila, using microsatellites, two genetic clusters were also detected, with the eastern cluster containing most samples from locations in the northeast of our sampling range and the western cluster containing samples from the rest locations (Supplementary Fig. S1b and e).

Interestingly, the combination of the two eastern clusters and that of the two western clusters of the two major cryptic pollinator species roughly matched the ranges of the two clusters of F. pumila var. pumila, though the boundary of the pollinators slightly shifted eastwards compared with that of the plant (Fig. 1a and b). The phylogenetic tree and the haplotype network based on cpDNA haplotypes also showed the presence of the eastern and the western cluster, but the dividing line moved westwards to the Wuyi Mountains (Fig. 1b and c).

Historical demography of the pollinators and F. pumila var. pumila

Results from BSP analysis showed that the effective population sizes of Wiebesia sp. 1 and Wiebesia sp. 2 started to increase since about 20 000 years BP and 30 000 years BP but the effective population size of Wiebesia sp. 3 began to decrease from about 10 000 BP (Fig. 2a–c). Moreover, significant sudden expansion was found in Wiebesia sp. 1 and Wiebesia sp. 2 but not in Wiebesia sp. 3 (Table 1). We observed different mismatch distributions among the three cryptic pollinator species. Wiebesia sp. 1 and Wiebesia sp. 2 exhibited a flat and unimodal peak distribution curve (Fig. 2d and e). For Wiebesia sp. 3, a sharp decrease in frequency was detected when pairwise differences increased (Fig. 2f), indicating a recent bottleneck. In addition, the neutrality tests in Wiebesia sp. 1 and sp. 2 also showed evidence for expansions, with both Tajima’s D and Fu’s Fs having negative values though significant results were only detected for Fu’s Fs (P < 0.001, Table 1).

For F. pumila var. pumila, mismatch distributions using cpDNA sequences under both sudden expansion and spatial expansion models did not reject significant demographic and spatial expansions in the two genetic clusters (Table 1). Besides, the observed pairwise differences of these two clusters rapidly decreased to very low values (Fig. 2g and h), also suggesting recent expansions. Additionally, negative values were found in Tajima’s D and Fu’s Fs in both clusters, a sign of recent population expansions (Table 1).

Potential distribution at LGM of the pollinators and F. pumila var. pumila

For Wiebesia sp. 1, Wiebesia sp. 2 and F. pumila var. pumila, the cross-validation of the climate envelope models revealed high mean model fitness for both LGM-CCSM and LGM-MIROC, and their current potential distribution ranges predicted by MAXENT were consistent with the records of current distribution (Fig. 3).

The highest climate suitability supported by both CCSM and MIROC showed two potential refugia for Wiebesia sp. 1, with the eastern refugium including the east coast and some islands in the Ryukyu Islands and the western refugium covering some areas in the mainland (Fig. 3b and c). The overlapped results from both models also detected that the potential refugia for Wiebesia sp. 2 comprised some areas in the south of mainland (the eastern refugium) and some island in South China Sea and Ryukyu Islands (the western refugium) (Fig. 3e and f). For F. pumila var. pumila, two potential refugia were found to cover most area of the potential refugia for the two cryptic pollinator species by both models: the eastern one contained the east coast area and some islands in the Ryukyu Islands and the western one covered a large area in the mainland (Fig. 3h and i).

DISCUSSION

Both biotic and abiotic factors participate in shaping the phylogeographic patterns of organisms (Althoff et al. 2012; Hewitt 2000). Here, we sampled F. pumila var. pumila and its host-specific pollinators across their distribution ranges in China to reveal their current genetic structure and postglacial phylogeographic history. Consistent with Chen et al. (2012), we found three cryptic pollinator species within our study area, with the dominant two species distributing in the northern and southern part of the distribution range of F. pumila var. pumila (Fig. 1a). The two dominant cryptic pollinator species and the plant displayed parallel current spatial genetic structure, with all species containing an eastern and a western genetic cluster (Fig. 1a and b). Furthermore, these three species had two potential glacial refugia, with the western refugium located in the mainland and the eastern refugium covering the coastal area and some nearby islands (Fig. 3), followed by expansion after the LGM (Fig. 2). These results suggest that these three species share similar phylogeographic history, contributing to the current parallel genetic structure. Therefore, our study provides evidence for our hypothesis that similar phylogeographic patterns appeared in the Last Glacial Period in these tightly interacting species.

The presence of three cryptic pollinator species and the absence of corresponding genetic divergences within F. pumila var. pumila indicate that the pollinators evolved at a faster rate than their host plant. As the key coevolved traits (e.g. attractive chemicals and body size) did not change among cryptic pollinator species (Wang et al. 2021), speciation within W. pumilae is possibly driven by abiotic factors rather than coevolution with F. pumila. Both the current distribution and the potential refugia at LGM of Wiebesia sp. 1 are in the north of those of Wiebesia sp. 2, probably reflecting that Wiebesia sp. 1 adapts to the colder climate. Wiebesia sp. 3 is the only species pollinating F. pumila var. awkeotsang but only existed in a few F. pumila var. pumila populations. In addition, we did not find signals of postglacial expansion for this species. These suggest that it is a weaker competitor than other cryptic pollinator species, probably due to its later emergence during the flowering period of F. pumila (Liu et al. 2014). Nevertheless, source-sink dynamics has been found to assist the long-term maintenance of its island populations (Tong et al. 2021).

Consistent with Tian et al. (2015), which found identical genetic structure between two mutualists in another fig–pollinator mutualism, our results showed roughly matched phylogeographic patterns between the pollinators and F. pumila var. pumila. This confirms that obligate pollinators should have synchronous responses with their host plants to changes in biogeographic factors. Nevertheless, slight differences were still observed between the pollinators and the plant, possibly due to the dispersal syndrome of F. pumila var. pumila.

Genetic structure of plant species is shaped by both pollen and seed dispersal (Wang et al. 2011), and the range expansion of the plants mainly relies on seed dispersal (Chen et al. 2008). Though pollen movements may play the predominant role in shaping the genetic structure of F. pumila var. pumila (reflected by the high pollen-to-seed dispersal ratio), frequent long-distance dispersal of F. pumila var. pumila seeds is likely to occur because the common cpDNA haplotype (F1) in the eastern cluster also dominated some western populations (Fig. 1b). As an important food resource, mature figs of F. pumila can attract many long-distance seed dispersers (e.g. bats and large birds) (Corlett 2006). The presence of long-distance seed dispersal is likely to prevent genetic differentiation among F. pumila. var. pumila populations pollinated by different cryptic pollinator species, resulting in a single eastern/western genetic cluster. Nevertheless, the Nanling–Wuyi Mountains may serve as a natural barrier for seed dispersal based on the cpDNA results.

Restricted movements of W. pumilae has been reported by Wang et al. (2013), and our study also detected isolation by distance patterns in Wiebesia sp. 1 and Wiebesia sp. 2 (Supplementary Fig. S3a–d), probably because of their extremely short lifespan at adult stage (1–3 days) (Liu et al. 2014). Thus, expansions of the pollinators from the coastal and insular refugia may not exactly keep pace with the colonization of F. pumila var. pumila, possibly resulting in a smaller range occupied by the eastern clusters of the pollinators than that occupied by the eastern cluster of F. pumila var. pumila. Moreover, the absence of geographic barriers between the two genetic clusters for either cryptic pollinator species indicates that the current genetic structure is not at equilibrium. The western cluster may be more adapted to mainland environments while coastal populations in the eastern clusters can get recruitment from island populations.

The presence of eastern and western potential refugia at LGM has been found in many species in the Sino–Japanese region, where the emergence of the land bridge of the East China Sea often acted as a corridor or a refugium during the Last Glacial Period (Jiang et al. 2021). However, this theory has seldom been tested for pairs of interacting species from different tropic levels. Our results provide evidence and support that tightly interacting insects can closely follow the migration trajectories of their host plants, offering an excellent case for understanding the historical dynamics of vertical diversity in East Asia. Moreover, consistent with other plant species in the Sino–Japanese region (e.g. Wang et al. 2015), the refugia in southwest mainland China harbored most genetic diversity of F. pumila var. pumila, as population genetic diversity of F. pumila var. pumila decreased with increasing latitude and longitude and a high haplotype diversity was found in the southwest of the plant’s distribution range.

CONCLUSIONS

Species in an ecosystem are inevitably linked by interspecific interactions, and the ongoing global changes are likely to distort these interactions and thus threaten the maintenance of biodiversity (Segar et al. 2020; Urban et al. 2016). Our results show that tightly interacting species responded similarly to historical climate changes, providing insights into the processes forming the current distribution of vertical diversity across different trophic levels.

Supplementary Material

Supplementary material is available at Journal of Plant Ecology online.

Table S1: Population information of F. pumila.

Table S2: Outgroups used in the phylogenetic tree of fig wasps.

Table S3: Population information of Wiebesia sp. 1.

Table S4: Population information of Wiebesia sp. 2.

Table S5: Population information of Wiebesia sp. 3.

Figure S1: Results of genetic clustering for Wiebesia pumilae (Wiebesia spp.) (a) and its three cryptic species in W. pumilae (b–d) and F. pumila (e) using STUCTURE.

Figure S2: Relationships between genetic diversity indices (N_A, A_R, H_E, H_O) in F. pumila populations and latitude (a–d) and longitude (e–h) using linear regressions.

Figure S3: Results of Mantel tests of the genetic divergence (represented by both F_ST/(1 − F_ST) and F_ST(ENA)/(1 − F_ST(ENA)) among populations of the three cryptic species in W. pumilae (a–f) and F. pumila (g, h) and the geographic distances (log-transformed) among populations.

Authors' Contributions

Min Liu (Data curation, Formal analysis, Investigation, Methodology, Validation, Writing—original draft, Writing—review & editing), Man-Juan Huang (Formal analysis, Methodology, Validation, Visualization, Writing—review & editing), Finn Kjellberg (Conceptualization, Methodology, Writing—original draft, Writing—review & editing), Yan Chen (Data curation, Investigation, Methodology, Writing—review & editing), Jian Zhang (Data curation, Investigation, Methodology, Writing—review & editing), Yuan-Yuan Ding (Data curation, Investigation, Methodology, Writing—review & editing), Yang Yang (Data curation, Investigation, Methodology, Writing—review & editing), Jun-Yin Deng (Methodology, Writing—review & editing), Kai Jiang (Investigation, Methodology, Writing—review & editing), Yuan-Yuan Li (Investigation, Methodology, Writing—review & editing), Xin Tong (Investigation, Methodology, Writing—review & editing), Tong Luo (Visualization, Writing—review & editing), Rong Wang (Conceptualization, Formal analysis, Investigation, Project administration, Supervision, Writing—review & editing), and Xiaoyong Chen (Conceptualization, Formal analysis, Funding acquisition, Investigation, Project administration, Resources, Supervision, Validation, Writing—original draft, Writing—review & editing)

Funding

This study was supported by the National Natural Science Foundation of China (32171609, 32261123001) to X.-Y.C.

Acknowledgements

Many local people provided support in field sample collection.

Conflict of interest statement. The authors declare that they have no conflict of interest.

References

Author notes