1. Introduction

The Qinghai-Tibet Plateau (QTP) plays an important role in the evolution of mountain species due to its unique geological history and topography, complex climate, and diverse habitats (Liu et al., 2014; Wen et al., 2014; Favre et al., 2014; Wu et al., 2022). Comprising the Plateau Platform (also known as Qinghai-Tibet Plateau sensu stricto), the Himalayas and the Hengduan Mountains, the QTP contains more alpine plants than anywhere else outside the arctic (Li et al., 2014; Yu et al., 2018a; Mao et al., 2021). The rapid uplift of the QTP and climate oscillation during the Quaternary profoundly shaped the current distributions and genetic structures of many plant and animal species in this biodiversity hotspot (Lei et al., 2014; Liu et al., 2014; Muellner-Riehl, 2019). The impact of environmental changes since the Quaternary on speciation, genetic diversity, population genetic structure and demographic history can be studied using phylogeography and population genetics (Qiu et al., 2011; Liu et al., 2012; Wu et al., 2022). The last 20 years have seen an exponential growth of phylogeographic studies in this region by tracking the spatial distribution of genetic variation and lineage within specific species or between closely related species (Qiu et al., 2011; Liu et al., 2012, 2014; Yu et al., 2018b; Mao et al., 2021).

The most recent uplift of the QTP between the Late Miocene and Pliocene significantly altered the tectonic shape of the plateau, creating new mountains and rearranging river systems (Wang et al., 2014; Deng and Ding, 2015; Spicer et al., 2020). Many plant species exhibit a phylogeographic separation pattern consistent with the effects of these geological barriers that block or weaken gene flow (Fan et al., 2013; Zhao and Gong, 2015; Tian et al., 2020; Cheng et al., 2021). The Mekong-Salween Divide (MSD), which is located to the east of the Bershula and Gaoligong Mountains, roughly divides the southern and southeastern parts of the QTP into the eastern Himalayas and the Hengduan Mountains (Ward, 1921). The barrier, which was identified a century ago (Ward, 1921), has been shown to act as a phylogenetic break in some species by hindering gene flow between eastern and western populations (Sun et al., 2017), e.g., Taxus wallichiana Zucc., Sinopodophyllum hexandrum (Royle) T.S. Ying, and Marmoritis complanatum (Dunn) A.L. Budantzev (Gao et al., 2007; Li et al., 2011; Luo et al., 2017). The Tanaka-Kaiyong Line (TKL; Jiang et al., 2017; Qian et al., 2020), spanning Yunnan and Sichuan in southwestern China, is a dividing line between climate zones (Indian and East Asian monsoons to the west and east, respectively) and their associated biota (Li and Li, 1997; Ye et al., 2017). Distinct genetic lineages on either side of the TKL have been detected in various plant species, supporting its role as a phylogeographic dividing line that promotes differentiation (Fan et al., 2013; Tian et al., 2020). However, plant species with distinct traits, such as those related to dispersal ability, are known to occupy different ecological niches and geographic ranges (Papadopoulou and Knowles, 2016; Wu et al., 2018). Consequently, MSD and TKL may play diverse roles in shaping the phylogeographical patterns of different plant taxa.

Populus rotundifolia Griff (Salicaceae) is a deciduous tree native to southwestern China, naturally found on both sides of both the MSD and TKL. This species relies on wind for both pollination and seed dispersal (Fang et al., 1999). A recent study demonstrated that wind direction can systematically influence the phylogeographic patterns of species with wind-mediated dispersal and pollination (Kling and Ackerly, 2021). Similarly, our two published studies indicated that the phylogeographic breaks in two poplar species were significantly influenced by spring wind directions (Song et al., 2021; Li et al., 2023). Thus, the effects of a biogeographic divide may differ profoundly among wind-dispersed species, e.g., due to phenology. The March to May flowering and seed dispersal period of P. rotundifolia coincides with strong seasonal northwest winds in the Hengduan Mountains that contains both the MSD and TKL (Fang et al., 1999). Such strong winds might promote gene flow across biogeographical divides. The Hengduan Mountains, an area where the East Asian monsoon from the southeast meets the Indian monsoon from the southwest (Yun et al., 2014; Spicer 2017; Ha et al., 2018), is known as a global biodiversity hotspot (Mountains of South-Western China). Therefore, P. rotundifolia is an ideal plant model to test how these dividing lines apply to a wind-dispersed species, and whether they differ in their effects.

In this study, 14 microsatellite loci exhibiting low null allele frequencies and four chloroplast DNA fragments were used to study the genetic diversity and distribution patterns of 49 populations of Populus rotundifolia in southwestern China. We asked the following questions: (1) Are the Mekong-Salween Divide and Tanaka-Kaiyong Line partial or complete barriers to gene flow between populations of P. rotundifolia? (2) Does wind direction shape gene flow and genetic diversity patterns in P. rotundifolia, resulting in greater diversity downwind, and/or gene flow across geographic barriers? (3) What role did Quaternary climate history play during the phylogeographic history of P. rotundifolia? To answer these questions, we used Bayesian computation to infer gene flow and demographic history, and ecological niche modeling (ENM) to reconstruct the potential distribution of P. rotundifolia during the current time period, the last glacial maximum (LGM), the last interglacial (LIG), and the mid-Holocene (MH). Our findings will increase our understanding of the phytogeographic patterns around this biodiversity hotspot.

2. Materials and methods 2.1. Sampling and data collection

Our previous work determined that Populus rotundifolia is predominantly distributed in Tibet (Xizang), Yunnan, and Sichuan, whereas Populus davidiana Dode is distributed in northern China, such as Gansu and Shaanxi (Zheng et al., 2017). In this study, we employed a comprehensive data set that included 583 individuals from 49 populations (refer to Table 1; see Fig. 1), covering the main geographic distribution of P. rotundifolia as outlined by Zheng et al. (2017). Among these, 156 individuals from 32 populations were previously sampled (Zheng et al., 2017), while an additional 427 individuals were newly sampled for the current study (Table 1). Fresh leaves were sampled with > 50 m distances between any two sampled individuals, to avoid inbreeding or sampling clones. Geographic information for each population was recorded by the Extrex eTrex GIS unit (Table 1). Voucher specimens were deposited at Sichuan University. To determine the wind direction at each population and adjacent regions, we downloaded wind direction data for all available weather station across the distribution range of P. rotundifolia from the National Meteorological Science Data Center (https://data.cma.cn/site/index.html) for the months March, April and May from 1981 to 2010.

We examined 16 nuclear microsatellites (nSSR; Tuskan et al., 2004; Ma et al., 2013; Jiang et al., 2016) and four chloroplast DNA (cpDNA) fragments of 49 sampled populations. Total genomic DNA was extracted from silica gel-dried leaf tissue using a modified CTAB method (Porebski et al., 1997). Three of the cpDNA fragments, matK, trnG-psbK, and psbK-psbI, were amplified using primers described by Feng et al. (2013), whereas ndhC-trnV amplification followed Zheng et al. (2017). An individual of Populus adenopoda Maxim. was sampled as an outgroup, to root the cpDNA haplotype network. PCR was performed in a 20 μL reaction volume containing 50–100 ng DNA, 10 × PCR reaction buffer, 2.5 mM MgCl₂, 0.02 mM dNTP mix, 1 μL of each primer, and 0.1 μL of Taq DNA polymerase. Amplification was performed as follows: initial denaturation for 2 min at 98 ℃, followed by 34 cycles of 10 s at 98 ℃, 10 s at 55 ℃, 12 s at 72 ℃, and final extension for 2 min at 72 ℃. The PCR products were checked on 1% agarose gels to confirm gel bands of the expected fragment size. The PCR products were sent to Tsingke Biological Technology (Beijing, China) for DNA sequencing.

2.2. Statistical analyses of nSSR data

A total of 16 nSSR loci were genotyped and two steps were taken to minimize potential errors at each locus. Detailed data processing is shown in Appendix S1. Two loci that showed high null allele frequencies (F[Null] > 0.4) were excluded from all subsequent analyses (Dakin and Avise, 2004), leaving 14 loci for subsequent examination. A Bayesian clustering approach was implemented in STRUCTURE v.2.3.3 (Pritchard et al., 2000; Hubisz et al., 2009) based on the nuclear DNA sequences to examine population genetic structure. The optimal value of K, out of values from 1 to 8, was determined using the method of Evanno et al. (2005) and Pritchard et al. (2000). To cross-validate the results of STRUCTURE, we also conducted a Principal Coordinates Analysis (PCoA) on the nSSR data using GenAlEx v.6.5 (Peakall and Smouse, 2012).

In addition, Monmonier's maximum difference algorithm in BARRIER v.2.2 was used to identify biogeographic boundaries or regions exhibiting the largest genetic discontinuities between population pairs (Manni et al., 2004). The robustness of these boundaries was evaluated by running BARRIER on 100 bootstrap replicates of population pairwise difference matrices for cpDNA sequence data, and 1000 bootstrap replicates for nSSR loci. The bootstrap matrix for cpDNA sequence data was generated using SEQBOOT and ARLEQUIN v.3.5 (Excoffier and Lischer, 2010), and the bootstrap matrix for nSSR loci was generated using MICROSATELLITE ANALYSER (MSA) v.4.05 (Dieringer and Schlötterer, 2003).

Genetic diversity indices, i.e., the average number of alleles (A_a), the effective number of alleles (A_e), observed heterozygosity (H_o), expected heterozygosity (H_e), Shannon's information index (I), allelic richness based on five samples (A_r_5), fixation index (F), and the percentage of polymorphic loci (PPL; Table S5) were estimated using GenAlEx v.6.5 (Peakall and Smouse, 2012) and Microsatellite Analyser (MSA) v.4.05 (Dieringer and Schlötterer, 2003). Nei's expected heterozygosity and F-statistics were examined using analysis of molecular variance (AMOVA) as implemented in ARLEQUIN v.3.5 (Excoffier and Lischer, 2010). We divided all sampled populations into three regional groups (Western, Eastern, and Central) based on their geographical distribution relative to significant landmarks: west of the Mekong-Salween Divide (MSD), east of the Tanaka-Kaiyong Line (TKL), and the intervening area. Genetic variation was hierarchically partitioned into three levels: among groups, among populations within groups, and within populations. Finally, to test the significance of isolation by distance, we performed a Mantel test on the matrix of genetic distances and geographical distances (Somers and Jackson, 2022).

2.3. Statistical analyses of cpDNA sequence data

The cpDNA sequences were aligned and refined using MEGA v.11 (Tamura et al., 2021), while haplotypes and variable sites were detected using DNASP v.5.1 (Librado and Rozas, 2009). Insertions/deletions (indels, excluding mononucleotide repeats) were coded by the simple gap coding method as implemented in Gapcoder (Young and Healy, 2003). A median-joining haplotype network with Populus adenopoda as an outgroup was constructed using NETWORK v.4.6 under the maximum parsimony criterion for the 49 populations (Bandelt et al., 1999). The nucleotide diversity (π) and haplotype diversity (H_d) were calculated using DnaSP v.5.1 (Librado and Rozas, 2009). In addition, average gene diversity within populations (H_S), total gene diversity (H_T), coefficient of genetic variation overall populations (G_ST), and coefficient of genetic variation influenced by both haplotype frequencies and genetic distances between haplotypes (N_ST) were estimated by the program PERMUT v.3.0 (Pons and Petit, 1996). Analyses of molecular variance (AMOVA) were carried out in ARLEQUIN v.3.5 (Excoffier and Lischer, 2010), and the classification of grouping information and genetic variation is the same as in the previous section. To further examine if there are possible genetic barriers among P. rotundifolia populations, we also performed spatial analysis of molecular variance for cpDNA using SAMOVA v.2.0 (Dupanloup et al., 2002).

2.4. DIYABC and MIGRATE analyses

To decipher the demographic history of Populus rotundifolia, we also employed an approximate Bayesian computation framework in DIYABC v.1.0 (Cornuet et al., 2008). The detailed parameter settings are presented in Table S2 and Appendix S2. To transfer evolutionary timescale in generation numbers into absolute timescale in calendar years, the generation time of P. rotundifolia was set as 40 years, with a span of 20–60 years (Macaya-Sanz et al., 2012). Additionally, we employed the software MIGRATE v.3.3.1 (Beerli and Palczewski, 2010) to assess gene flow among different regional groups of P. rotundifolia, drawing on data from both nSSR and cpDNA data. This process involved calculating the population-scaled mutation rate (θ) and migration rates (m), enabling the estimation of the magnitude and direction of gene flow per generation from adjacent populations. For nSSR data, we adopted a mutation rate (μ) of 10⁻³ per gamete per generation, as per the methodology outlined by Lexer et al. (2005).

2.5. Ecological niche modeling

To investigate how Populus rotundifolia responded to Quaternary climate change history, we employed ecological niche modeling (ENM) to predict their potential distribution during the Last Interglacial (LIG; ca. 130 Ka), the Last Glacial Maximum (LGM; ca. 21 Ka), the mid-Holocene (MH; ca. 6 Ka), and the present. We obtained nineteen bioclimatic variables for these periods at 2.5 arc-minute resolutions from the WorldClim database (Fick and Hijmans, 2017). Data for LIG, LGM, and MH were sourced from the circulation model simulation of the Community Climate System Model (CCSM; Collins et al., 2006). To address multicollinearity, we assessed pairwise Pearson correlations (r) among the 19 bioclimatic variables (Bio1–19) and excluded one variable from each pair with r > 0.7. Subsequently, seven environmental variables were retained for subsequent modeling and analyses (Table S3).

The ENM was implemented using the maximum entropy method implemented in MAXENT v.3.3.3 (Pearson et al., 2007; Phillips et al., 2006; Phillips and Dudík, 2008). To ensure the accuracy of our records, only coordinates corresponding to P. rotundifolia population records obtained during our field surveys were used for the ENM, yielding a total of 49 records (Table 1). This approach was adopted to prevent the inclusion of potentially inaccurate data from public databases. See Appendix S3 for specific parameter settings.

To determine the ecological divergence among different groups of P. rotundifolia, we conducted a Principal Component Analysis (PCA) of the bioclimatic variables that were used for SDMs using two R packages: FactoMineR (Lê et al., 2008) and FactoExtra (Kassambara and Mundt, 2017).

3. Results 3.1. Variations and distribution pattern of cpDNA sequences

The total length of the four cpDNA regions sequenced was 2032 bp. A total of 62 chloroplast haplotypes (25 non-singleton haplotypes and 37 singleton haplotypes) were identified across all 49 populations (Table 1). The phylogenetic network of these 62 haplotypes revealed five distinct haplotype lineages (Ⅰ, Ⅱ, Ⅲ, Ⅳ, and Ⅴ; Fig. 1a–c). The most frequent haplotype (H3) was located near the center of the haplotype network, and private haplotypes were nearly all located in the eastern and central parts of the species range. Haplotype lineages Ⅰ and Ⅱ were mainly found on the western side of the Tanaka-Kaiyong Line (TKL), whereas lineages Ⅳ and Ⅴ were primarily distributed on the eastern side of the Mekong-Salween Divide (MSD). The presence of all five lineages within the area bounded by the TKL and MSD suggests the possibility of gene flow across these biogeographical barriers (Fig. 1a–c).

AMOVA based on cpDNA found that 64.09% of genetic variation was distributed within populations, 28.85% was portioned among populations within groups, and 7.06% between groups (Table 2). For the Western group, 28.90% of the variation was among populations, while the remaining 71.10% was within populations. For the Eastern group, these values were 21.67% and 78.33% respectively (Tables S7 and S8). The genetic differentiation index (F_ST) was 0.3591 (P < 0.05), indicating a pronounced level of genetic differentiation among populations (Table 2). The total genetic diversity (H_T = 0.88) was higher than the genetic diversity within populations (H_S = 0.65), indicating that the intraspecific differentiation within a population was larger than that among different populations. A permutation test (G_ST = 0.264; N_ST = 0.323) suggests that there exists a phylogeographical structure in Populus rotundifolia (N_ST > G_ST; P < 0.01; Table 3).

3.2. Genetic diversity, population structure, and group division based on nSSR markers

For all populations, the number of alleles per locus varied from 5 to 16, with an average of 10 (Table S4). The F_ST value averaged across all loci was 0.232, which indicated a pronounced level of genetic differentiation among populations (Table S4). The mean values for genetic indexes were as follows: number of alleles (N_a) was 2.872, the number of effective alleles (N_e) was 1.963, the Shannon index (I) was 0.666, the observed heterozygosity (H_o) was 0.387, the expected heterozygosity (H_e) was 0.372, the allelic richness based on five samples (Ar_5) was 2.354, fixation index (F) was −0.023 and the percentage of polymorphic loci (PPL) was 79.01% (Table S5). Moreover, AMOVA results indicated that 82.68% of genetic variation occurred within populations, whereas only 15.05% and 2.26% of the variation was partitioned among populations within groups and among groups, respectively, with a weak separation among populations in the Eastern, Central, and Western (F_CT = 0.0226, P = 0.02; Table 2). The Mantel test suggested that the correlation between geographic and genetic distances is weak and non-significant (R² = 0.0066, P = 0.09; Fig. S4).

Bayesian clustering analyses suggested that the optimal number of free mating meta-populations is two (Fig. S2). When K = 2, two genetic clusters were identified. One was predominantly distributed on the western side of the MSD, while the other occurred on the eastern side (Fig. 2a). This pattern suggests that the MSD may serve as a biogeographical barrier for Populus rotundifolia. When the analysis was extended to K = 3, some samples situated between the TKL and MSD exhibited a distinct genetic component (Fig. 2c), highlighting the complex genetic landscape of P. rotundifolia across these biogeographical divides.

BARRIER analyses identified strong biogeographical barriers with populations for nuclear cpDNA and nSSR data, but these biogeographical barriers are messy and difficult to demarcate, especially in the central regions, because support values for most of the biogeographical barriers were relatively low (< 80%; Fig. S7). The results derived from SAMOVA analysis revealed varying groups under different K settings (Table S13), indicating the presence of weak genetic barriers. The PCoA analysis revealed a distinct separation among populations west of the MSD and those east of the TKL, with populations located between the MSD and TKL showing admixture between them (Fig. 2b). Similar patterns were observed in the haplotype network of cpDNA data (Fig. 1).

3.3. Demographic history and migration rates of Populus rotundifolia

DIYABC analyses were used to simulate and compare seven plausible scenarios concerning effective population size changes for either the Western or Eastern group of P. rotundifolia (Fig. S1). Both direct estimate (DE) and logistic regression (LR) methods suggested the Western group of P. rotundifolia expanded abruptly (N_e increased from 1870 to 94,000) around 616 Ka (Fig. 3 and Table S11). The Eastern group is assumed to have expanded abruptly (Ne increased from 2080 to 101,000) around 598 Ka (Fig. 3 and Table S12). This scenario (scenario 2) is much more likely than the other six scenarios for each of the Western (P_DE > 0.70, P_LR > 0.97) and Eastern (P_DE > 0.64, P_LR > 0.92) lineages (Figs. S5 and S6). Type Ⅰ and Ⅱ errors for this scenario were 0.14 and 0.00, for both the Western and Eastern groups (Tables S9 and S10). PCA of the summary statistics for all simulated scenarios and posterior distribution for parameters of this scenario in both Western and Eastern groups are shown in Fig. S6.

Utilizing nSSR data, analyses conducted with MIGRATE revealed that the estimated effective population sizes for Populus rotundifolia within the Western (θ = 0.4557) and Central (θ = 0.4587) groups are significantly greater than that of the Eastern group (θ = 0.0187; Fig. 2d). The average estimated gene flow (4N_em) demonstrated more frequent movement from Western to Central, Eastern to Central, and Eastern to Western, in comparison to their reverse directions (501.20 vs 91.08, 20.92 vs 11.52, and 25.06 vs 6.56 respectively; Fig. 2d).

A parallel assessment using cpDNA via MIGRATE yielded congruent findings. According to cpDNA analysis, the Central group exhibited the largest effective population size (θ = 0.0136), followed by the Western group (θ = 0.00096), with the Eastern group displaying the smallest size (θ = 0.00008; Fig. 1d). Gene flow was substantially higher from the Western to Central group (N_em = 191.60) and from the Eastern to Central group (N_em = 124.32), compared to flows in the opposite directions (62.43 and 2.9; Fig. 1d). The gene flow from Western to Eastern group (N_em = 22.57) was slightly lower than in the reverse direction (N_em = 25.06; Fig. 1d).

3.4. Wind direction and ecological niche modeling

We found that the prevalent wind direction in March and May is northwest for most weather stations (Fig. 1b). However, northeast wind was found to be dominant in weather stations close to P41–P47, and southwest wind were found for weather station close to P1. Five weather stations were found to have experienced a shift of wind direction from March to May, and a shift to northwest wind direction was found in weather stations close to P1, P41, and P43/P44 (The China Meteorological Data Service Center; http://data.cma.cn/).

The Ecological Niche Modeling (ENM) based on the Current period resulted in an area under the receiver operating characteristic curve (AUC) value of 0.983 ± 0.010, indicating the potential distribution predicted for the Current period was in keeping with the actual current distribution of Populus rotundifolia (Fig. 4). The MH and Current simulation have similar distributions, although with a more restricted distribution for the LGM and LIG (Fig. 4). The ENM results indicated that Isothermality (Bio3) and Mean Temperature of Driest Quarter (Bio9) were the main factors influencing the distribution of P. rotundifolia (59.8% and 25.8%, respectively; Table S3).

PCA revealed niche differentiation between different groups. There was no niche overlap between Western and Eastern groups, slight niche overlap between Central and Eastern groups, yet obvious niche overlap between Western and Central groups (Fig. 5).

4. Discussion 4.1. Genetic diversity and genetic structure of Populus rotundifolia

In this paper, following the removal of loci with high null allele frequencies, 14 microsatellite loci and four chloroplast DNA fragments were used to study the genetic diversity and distribution patterns of 49 populations of P. rotundifolia in southwestern China. The results of microsatellite data analysis showed that the effective number of alleles (Ne), expected heterozygosity (He) and the total genetic diversity (H_T) of P. rotundifolia were all slightly higher than the average value for woody plants (Hamrick and Godt, 1996; Gong et al., 2016; Hahn et al., 2017; Fan et al., 2018a; Chung et al., 2020). In addition, AMOVA analysis showed that the genetic differentiation coefficient (F_ST) based on uniparentally inherited molecular markers, i.e., chloroplast DNA (cpDNA: 0.3591), is substantially higher than the F_ST for biparentally inherited molecular markers, i.e., nuclear microsatellites (nSSRs: 0.1732; Table 2). Wind dispersal of both pollen and seeds, which facilitates gene flow within and between populations, may have contributed to higher levels of genetic diversity in P. rotundifolia than it does in other woody plant species. The degree of genetic differentiation in P. rotundifolia is consistent with that of other Populus species, but higher than that of most other wind-pollinated species examined (Hamrick, 2004; Zheng et al., 2017; Fan et al., 2018b). For example, the genetic differentiation coefficients of P. davidiana based on cpDNA and microsatellites are 0.54 and 0.13, respectively (Song et al., 2021).

Inter-population differentiation within P. rotundifolia is significantly greater (P < 0.05) when both haplotype frequency and phylogenetic relationships among haplotypes are considered (N_ST = 0.32) than when only haplotype frequency is considered (G_ST = 0.26). Hence haplotypes tend to occur in geographical proximity to their close relatives (Table 3). Conversely, the contribution of geographical isolation to the genetic differentiation of P. rotundifolia at microsatellite loci was not obvious, because there was no significant correlation between genetic distance and geographical distance among populations (R² = 0.0066; P = 0.09; Fig. S4). The uneven distribution pattern of genetic diversity in P. rotundifolia across populations (Tables 1 and S5) may be related to the climate refugia and recolonization history of the Quaternary and population demographic history (Hewitt, 2000; Hickerson et al., 2010; Gavin et al., 2014).

4.2. Contributions of geographic and environmental factors to phylogeographic structure of Populus rotundifolia

Our analysis identified 62 chloroplast haplotypes of P. rotundifolia, which were further classified into five haplotype lineages. Lineages Ⅰ and Ⅱ have a widespread presence on the western side of the TKL and extend across the MSD. In contrast, lineages Ⅳ and Ⅴ were primarily found on the eastern side of the TKL, although there were exceptions with some samples located on the western side of the TKL. Lineage Ⅲ was exclusively observed in the region between the MSD and TKL, highlighting its specific geographic confinement. Indeed, P36, P41 and populations between them contained all the haplotype lineages (Fig. 1), suggesting that either they represent refugia or a region open to gene flow from other parts of the range. Our analysis of cpDNA data indicate that populations from either side of the MSD share lineages Ⅰ and Ⅱ; thus, MSD does not appear to act as a phylogenetic barrier. In contrast, we found a clear divide along the TKL that block lineages Ⅰ, Ⅱ and Ⅲ, albeit with evident gene flow from east to west across its southern end.

Although lineages Ⅳ and Ⅴ share a common distribution, they exhibit distinct evolutionary relationships in the haplotype network (Fig. 1c). Lineage Ⅳ is more closely related to lineage Ⅱ, whereas lineage Ⅴ appears more aligned with outgroup species, indicating divergent evolutionary paths. This divergence finds support in Zheng et al. (2017), where chloroplast haplotypes of some Populus rotundifolia populations on the eastern side of the TKL matched those of P. davidiana, hinting at gene flow between these species. Such evidence points to the possibility that the Eastern group of P. rotundifolia may have derived chloroplast genomes from multiple sources, including gene flow from Central or Western groups as well as from closely related species like P. davidiana. New haplotypes may have then emerged under various evolutionary forces, such as mutations, natural selection and genetic drift, with incomplete lineage sorting possibly contributing. Crucially, regardless of the origins and the complex evolutionary dynamics at play, the TKL has acted as a significant biogeographical barrier. This barrier has facilitated genetic divergence between the populations located on its eastern and western sides, underscoring the TKL's pivotal role in shaping the genetic landscape of P. rotundifolia through isolation and the promotion of divergent evolutionary trajectories.

Bayesian clustering analysis based on nuclear simple sequence repeat (nSSR) data delineated two principal genetic lineages in P. rotundifolia, highlighting a pronounced division where populations on the western side of the MSD were distinctly separated, in contrast to the eastern populations which were largely grouped together, as depicted in Fig. 1, Fig. 2a. Expanding the analysis to consider K = 3 revealed a more nuanced genetic structure, albeit less sharply defined, yet clearly indicating three genetic compositions aligned with the delineated geographic regions. This pattern underscores the influence of both the MSD and the TKL in driving genetic divergence within P. rotundifolia.

Notably, the subtler differentiation observed at K = 3 underscores the more significant impact of the MSD over the TKL on the nuclear genetic divergence of the species. This observation suggests that, despite the less conspicuous genetic structures at this level, the MSD acts as a critical geographic barrier, profoundly shaping the genetic diversity and structure of P. rotundifolia and thus contributing to its evolutionary complexities. This highlights the complex interplay between geography and genetics in the evolutionary dynamics of P. rotundifolia, with the MSD playing a key role in its nuclear genetic differentiation.

In strong contrast to cpDNA data, nuclear data indicates a divide across the MSD but only minor differences across the TKL. However, this divide is incomplete and concerns proportions of two genetic lineages, rather than their presence/absence (Fig. 1a–c). Hence, for P. rotundifolia, there is some gene flow across the MSD, in contrast to the strong dividing lines seen in Taxus wallichiana (Gao et al., 2007), Marmoritis complanatum (Luo et al., 2017), Biston falcata Warren (Cheng et al., 2019), the fungus Boletus reticuloceps (M. Zang, M.S. Yuan & M.Q. Gong) Q.B. Wang & Y.J. Yao (Feng et al., 2017), and Schizothoracinae fish, for which the MSD contributed to rapid radiation in the QTP region (Tang et al., 2019).

Previous studies have shown Tacca chantrieri André and Bombax ceiba L. exhibited only very limited sharing of cpDNA haplotypes across the TKL (Zhao and Zhang, 2015), whereas Sophora davidii (Franch.) Skeels more closely parallels P. rotundifolia in that both cpDNA and nuclear haplotypes are shared across the TKL (Fan et al., 2013). Hence the strength of this barrier clearly varies between taxa, and while the signatures of both the MSD and TKL are clearly visible for P. rotundifolia, both barriers appear weaker for this species than for many others examined. However, the two lines, the TKL and MSD, act together to create a higher degree of separation between the eastern and western parts of the distribution, indicating that the Hengduan Mountains themselves might be better regarded as the gene flow barrier for this species, despite the central group being present there. It is possible that much of the central range is lost during glacial maxima and then recolonized, with seed flow from the west complemented by greater pollen flow from the east, or northward expansion from the southern extreme of the central area (populations 36–40), whose high diversity indicates possible glacial refugia.

In short, we find a series of incomplete genetic barriers moving from east to west, all hindering gene flow but none preventing it entirely (Fig. S7). Two interpretations are possible. One is that the central group, running though the Hengduan Mountains between the MSD and TKL, in its entirety serves as a broad phylogenetic barrier. An alternative explanation is that it forms the middle of a cline from east to west, however, this does not fit well with the multiple barriers detected (Fig. S7), nor with the discontinuous distribution of genetic variation across the region.

We argue that the current patterns of genetic variation within Populus rotundifollia were shaped by several factors in concert, namely, the MSD, TKL, ecological niche differentiation, spring wind direction, and the complex topography of the Hengduan Mountains. The most likely ecological barrier between the Eastern and the Central + Western populations of P. rotundifolia is the TKL. On one hand, PCA detected clear niche differences between the Eastern and Western groups (Fig. 5). The climate of the Western group and the Eastern group is affected by the Indian monsoon and the East Asian monsoon, respectively, which generates clear differences in temperature and precipitation between the two regions (Ye et al., 2017). On the other hand, although both monsoons meet in the Hengduan Mountains, where Central group occurs, PCA detected strong niche overlap between Western and Central groups and obvious niche difference between the Western + Central and Eastern groups (Fig. 5). Although the prevailing winds in these regions facilitate the transportation of pollen or seeds across the TKL, natural selection ensures that only a small percentage of these can develop into fertile individuals. This did not happen for the MSD due to the high similarity of ecological niches between the Western and Central groups (Fig. 5).

Second, the Hengduan Mountains, including those around the MSD, are geographic barriers that hinder seed flow. Based on nSSRs, most of the genetic barriers detected in the BARRIER results were concentrated within the Central group, that is, the Hengduan Mountains (Fig. S7). The existence of these barriers, overlapping with geographic barriers such as high mountains, will eventually hinder the gene flow between populations in the Western and Eastern groups, which will have a certain impact on the phylogeographical pattern of P. rotundifolia. MIGRATE analysis indicated a great difference in effective population size between the eastern and western subpopulations of P. rotundifolia, from which a low level of gene flow between them can be inferred (Fig. 1, Fig. 2d). At the same time, we noted that there are a series of genetic barrier close to the MSD but not overlapping. Mountains around the MSD, such as the Boshulaling Mountains, could have acted as strong barrier for pollen flow, as well as a weak barrier for seed flow, between P16+P17+P18 and P21+P22+P27, and a barrier for both pollen and seed flow between P23 and P24+P26+P28 (Fig. S7).

Third, the existence of gene flow across both the MSD and TKL might be due to wind dispersal of both seed and pollen, permitting connectivity that other species might lack (Butcher et al., 2005; Fan et al., 2013; Tian et al., 2020). In wind-pollinated species, pollen grains can be transported over long distances via the wind. Furthermore, the relative contribution of pollen dispersal to gene pools has been shown to be much larger than that of seed dispersal (Sutherland et al., 2010). Therefore, dominant wind directions during the flowering season may play a critical role in genetic variation of populations (such as P. davidiana; Song et al., 2021). In flowering season, downwind populations usually have higher genetic diversity (Kling and Ackerly, 2021). The East Asian monsoon and the Indian monsoon, both originated well before the Quaternary, brought prevailing winds in the southeast and southwest directions, respectively, and they meet near the MSD and TKL (Wang, 1994). The different monsoons in this region, may have facilitated both seed and pollen flow across the TKL and MSD, as seen from shared genetic variation between groups (Fig. 1, Fig. 2a). MIGRATE analysis outcomes indicate significantly stronger gene flow from the Eastern group to the Central group across the TKL for both cpDNA and nSSR markers, compared to gene flow in the opposite direction, as depicted in Fig. 1, Fig. 2d. Similarly, between the Western and Central groups, analysis shows a pronounced asymmetry in gene flow, with seed and pollen transport from the Western to the Central groups across the MSD being substantially greater than from Central to Western. These patterns highlight the significant influence of the East Asian monsoon and Indian monsoon on gene flow of this species.

At the same time, the populations with the largest number of haplotypes and the highest level of genetic diversity are P29, P36, P38, P41 and P44, and these populations are in the downwind of the prevailing northwest wind in the spring (Figs. 1 and S3; Table S5). Our results confirm a previous hypothesis that downwind populations have higher genetic diversity (Kling and Ackerly, 2021).

4.3. The influence of Quaternary climate change history on phylogeographical patterns of Populus rotundifolia

The fluctuations between cold and warm climate in the Quaternary had an important impact on phylogeographical patterns of plants in southwestern China (Wan et al., 2016; Zhao et al., 2019). One major hypothesis is that, during glacial maxima, plants retreated to refugia at lower latitudes and altitudes, then moved back during interglacial periods, forming the existing distribution pattern (Qiu et al., 2011; Liu et al., 2012; Muellner-Riehl, 2019). Haplotype diversity tends to be concentrated in climate refugia, with even small refugia often containing population-specific haplotypes (Abbott et al., 2000; Comes and Kadereit, 1998; Opgenoorth et al., 2010). Here, the existence of multiple refugia might contribute to both the extant diversity and the separation between eastern, central and western groupings. The western, central and eastern material contain respectively five, nine and five populations that include private haplotypes (Fig. S3), many of which might represent a distinct refugium. Furthermore, these populations tend to be to the southeastern edge of the distribution range and occur at lower altitudes, and occupy roughly the same area that the species occupied during the LGM, according to ecological niche modeling (Fig. 4). Moreover, this corresponds roughly to the main refugial areas detected for other species in the Hengduan Mountains (Chen et al., 2020; Lu et al., 2016).

Outside of refugia, populations have been founded by range expansion after the LGM and may be subject to founder or leading-edge effects, since they tend to acquire only part of the diversity present in source populations (Hewitt, 1999). This might explain the low chloroplast haplotype diversity, expected heterozygosity and number of effective alleles of some populations in the western, central, and eastern groups of P. rotundifolia (Fig. S3 and Table S5). Other populations outside of refugial areas have much higher genetic diversity scores, likely reflecting mixing of material from multiple refugia, as seen in Alnus glutinosa (L.) Gaertn (Havrdová et al., 2015). and Spartina alterniflora Loisel (Qiao et al., 2019). Nevertheless, we cannot rule out the possibility that these populations were derived from in situ climate refugia (Muellner-Riehl, 2019).

DIYABC simulations (Figs. 3, S5 and S6) indicate that both western and eastern material increased their effective population size by a factor of > 50, sometime before 600,000 years ago, roughly coinciding with the end of the coldest and most extensive glacial period on the Qinghai-Tibet Plateau (Zheng et al., 2002). This indicates population expansion following a severe bottleneck for both groups. However, the predicted range for the last interglacial period is smaller than that from the LGM (Fig. 4), so the bottleneck could reflect an interglacial period, rather than a glacial maximum, that imposed challenging conditions for the species. Either way, the bottleneck was followed by expansion, then presumably a series of less severe cycles of range shifts, expansions and contractions as further glacial maxima came and went. Periods of contact between the two formed the Central grouping, with rare chloroplast haplotype lineage Ⅲ becoming unique to the central populations, either arising there or going extinct in the west, forming the existing chloroplast haplotype distribution pattern (Fig. 1a–c). STRUCTURE analysis based on microsatellite data supports the origin of the central group through mixing of the other two under the optimal K = 2 (Fig. 2a–c), but gives it a separate identity under the next best K = 3 (Fig. 2c), which, like the unique haplotype lineage Ⅲ, is consistent with it beginning to develop its own genetic identity. Therefore, the complicated Quaternary climate oscillations and the history of local expansion are other important causes of the existing phylogeographical pattern of P. rotundifolia.

5. Conclusions

Our cpDNA and nSSR analyses of Populus rotundifolia revealed differentiated genetic lineages in the western and eastern groups, each of which expanded from very limited ranges around 600,000 years ago. A central group was mostly likely formed by contact between them, and this group had the highest diversity for both markers, which might reflect both admixture effects and the presence of multiple refugia there during the LGM. Although the central group is roughly separated from the western and eastern groups by the MSD and TKL phylogeographical breaks respectively, neither serves as a strong barrier to gene flow. Instead, multiple obstacles within the Hengduan Mountains, where the central group occurs, might make the whole area a large phylogeographic break, restricting gene flow between the western and eastern groups. Hence phylogenetic breaks do not always form simple lines, especially where complex history and/or topography is involved. Here, the phylogeographic break within P. rotundifolia reflects the combined impacts of biogeographic history, climate changes, monsoon direction and biological traits. Our research has emphasized the importance of employing multiple lines of evidence to infer the Quaternary evolutionary history of montane plants in topographically complex areas.

Acknowledgements

This study was funded by the National Natural Science Foundation of China (grants 41571054 and 31622015), the National Basic Research Program of China (grant 2014CB954100) and Sichuan University (Fundamental Research Funds for the Central Universities, SCU2021D006 and SCU2022D003; Institutional Research Funds, 2021SCUNL102).

Data availability statement

The microsatellite and chloroplast DNA data generated in this study have been deposited in the Dryad database (https://doi.org/10.5061/dryad.6djh9w125).

CRediT authorship contribution statement

Jieshi Tang: Writing – original draft, Visualization, Software, Formal analysis. Xiaoyan Fan: Software, Formal analysis, Data curation, Writing – original draft, Methodology. Richard I. Milne: Writing – review & editing. Heng Yang: Visualization, Software. Wenjing Tao: Software, Formal analysis. Xinran Zhang: Software. Mengyun Guo: Software. Jialiang Li: Writing – review & editing, Methodology, Conceptualization, Formal analysis, Software, Visualization. Kangshan Mao: Writing – review & editing, Supervision, Resources, Project administration, Investigation, Funding acquisition, Methodology, Conceptualization.

Declaration of competing interest

The authors declared that they have no conflicts of interest to this work.

Appendix A. Supplementary data

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