1. Introduction

The mountains of Southwest China are one of the world's biodiversity hotspots (Myers et al., 2000). Understanding the phylogeographic patterns within species, and the evolutionary forces that produced them, is a major objective of evolutionary biologists (Crisp et al., 2011; Usinowicz et al., 2017; Mao et al., 2021). The complex topography and environmental gradients generated by the uplift of the Qinghai-Tibet Plateau (QTP) in Southwest China can restrict dispersal and lead to diversification and intraspecific differentiation of mountain plants. Quaternary climate oscillations, however, may have either facilitated intraspecific differentiation by strengthening existing geographic barriers or may have enhanced gene flow between differentiated populations by bringing them into secondary contact. Recent studies have suggested that four general phylogeographic breaks exist in this region, including the Mekong-Salween Divide, the Tanaka-Kaiyong Line, the Sichuan Basin, and the ca. 105°E line (Ye et al., 2017). These may have existed as long-time gene flow barriers between intraspecific lineages, resulting in their differentiation. However, the geographic location of phylogeographic breaks is not always consistent with known geographic barriers (Soltis et al., 1997; Bond et al., 2001; Puorto et al., 2001; Song et al., 2021), with species-specific biological characteristics likely to be the main confounding factor (Irwin et al., 2001). Therefore, external (e.g., geological and climate history) and intrinsic factors (e.g., biological traits) may work simultaneously and interact with each other to produce present-day distribution patterns of intraspecific genetic variation (Hewitt, 2000; Macqueen et al., 2011).

The Sichuan Basin, which was formed during the middle Miocene (~13 million years ago (Mya)), is one of the most intriguing geographic features of Southwest China (Liu and Xu, 2003). The Basin is a lowland plain surrounded by mountains – the Hengduan Mountains in the west, the Qinling Mountains in the north, the Daba Mountains in the northeast, the Wushan Mountains in the east, the Dalou Mountains in the southeast, and the Yunnan-Guizhou Plateau in the southwest. The complex topographical features of these mountains have created diverse habitats (He and Jiang, 2014) along steep ecological gradients (Liu et al., 2014; Favre et al., 2015) that harbor great species diversity and endemism. Additionally, they have provided potential refugia for plants during climatic extremes, such as the Quaternary glacial cycles (Wen et al., 2017). However, the great elevational differences cause considerable habitat heterogeneity between the basin and the surrounding mountains, resulting in a "soft" genetic barrier for the dispersal of some species and impeding gene flow between populations (Bennett and Provan, 2008; He et al., 2019).

Many previous studies have investigated the phylogeographic patterns of species distributed around the Sichuan Basin (Qiu et al., 2009; Qi et al., 2012; Xie et al., 2012; He and Jiang, 2014; Sun et al., 2014; Ma et al., 2015; Qiao et al., 2018; Wang et al., 2021). For example, phylogeographic analyses have suggested that populations of Davidia involucrata were divided into an eastern and western lineage by the Sichuan Basin, which also separates the Sino-Himalayan and the Sino-Japanese Forest subkingdoms (Ma et al., 2015). Similar phylogeographic breaks have also been identified in Primula ovalifolia (Xie et al., 2012) and Dysosma versipellis (Qiu et al., 2009; Guan et al., 2010). The Sichuan Basin also acts as geographic barrier for closely related species pairs, such as Dipelta floribunda and D. yunnanensis (Tian et al., 2020). Because the genetic differentiation of these species occurred significantly later than the formation of the Sichuan Basin, the phylogeographic structure of these plants may be related to the diverse topography of the Sichuan Basin as well as a changing climate, rather than to the formation of the Sichuan Basin. Furthermore, a few taxa with a ring-shaped distribution around the Sichuan Basin have been studied, such as Tetracentron sinense (Sun et al., 2014), Odorrana margaratae (Qiao et al., 2018), Apodemus draco (Wang et al., 2021), Phoebe zhennan (Xiao et al., 2020) and Cercidiphyllum japonicum (Qi et al., 2012).

In addition to the phylogeographic break of the Sichuan Basin, there are two other phylogeographic breaks, the Kaiyong Line and the 105°E line, which traverse the mountains around the Sichuan Basin (Ye et al., 2017). The "Tanaka-Kaiyong Line" (TKL) is a major phytogeographic boundary in Southwest China composed of the Tanaka line (Tanaka, 1954), located in Yunnan Province, and the Kaiyong line (Lang, 1994), located in Sichuan Province (Li and Li, 1997). It was hypothesized that the TKL is a result of the uplift of QTP and the action of different monsoon systems on areas located to the west and east of the line (Li and Li, 1997). Uplift in the Hengduan Mountains region, located on the southeastern margin of the QTP, is generally believed to have been recent, occurring mainly between the late Miocene and late Pliocene (Zhou et al., 2006; Royden et al., 2008; Wang et al., 2012, 2014; Meng et al., 2014; Deng and Ding, 2015; Fang et al., 2020; Mao et al., 2021). The Asia monsoon system, triggered by the uplift of QTP, led to dramatic climate change in China (Shi et al., 1999; Spicer et al., 2020), and the current monsoon regimes, which differ, were established on either side of the TKL after or from the Late Pleistocene (Wang, 1994; Mitsui et al., 2008; Fan et al., 2013). Strong genetic differentiation across the TKL has been reported previously (Fan et al., 2013; Wang et al., 2014; Tian et al., 2015; Zhao and Gong, 2015; Zhao and Zhang, 2015; Cheng et al., 2021). Some researchers have proposed that the TKL is a monsoon-driven barrier to present-day plant dispersal (Fan et al., 2013; Tian et al., 2015; Zhao and Zhang, 2015). However, in Cephalotaxus oliveri (Wang et al., 2014) and Apodemus draco (Wang et al., 2021), the intraspecific divergence of species was found to be related to the uplift of Yunnan-Guizhou Plateau and QTP.

The 105°E line coincides with the boundary of two floristic subkingdoms, the Sino-Japanese and the Sino-Himalayan Forest subkingdoms (Wu and Wu, 1996; Ye et al., 2017). The Sino-Japanese Forest subkingdom is located in the eastern, warm-to-cold temperate hilly and low-plain areas (0–500 m a.s.l.), a region presently under a Pacific monsoon climate (Qiu et al., 2011). The Sino-Himalayan Forest subkingdom, which has both Pacific and Indian monsoon regimes, consists of a mosaic of plateaus, mountains, basins, river valleys, and deep gorges with wide ranging elevations (500–4000 m a.s.l.) (Wu and Wu, 1996; Qiu et al., 2011). Many phylogeographic studies of plants have shown that a significant genetic barrier exists from 105 to 108°E (Yan et al., 2012; Guo et al., 2014; Sun et al., 2014; Xu et al., 2015). Despite their documented role as barriers to gene flow, whether these three phylogeographic breaks (i.e., the Sichuan Basin, the Kaiyong Line, and the 105°E line) produce common effects on different plant species and molecular markers from different genomes (nuclear, plastid and mitochondrial genomes) of the same species, and why, await further investigation.

Populus lasiocarpa Oliv. is a wind-pollinated and wind-dispersed deciduous tree species that occurs in temperate forests between 1700 and 3500 m around the Sichuan Basin (Luan et al., 2011; Tang et al., 2013). Its distribution, which is unique for poplar germplasm resources in southwest China, extends to Hubei, Shanxi, Guizhou, and Yunnan provinces (Jiang et al., 2015). The species has the largest leaf area in the genus Populus L., and possesses well-developed perianths (Hong et al., 1987), which make it an excellent candidate for exploring the relative contributions of extrinsic (e.g., geological and climate history) versus intrinsic (e.g., biological traits) factors in shaping distribution patterns of genetic variation within species.

In this study, we asked the following questions: (a) Does the distribution pattern of genetic variation in Populus lasiocarpa reflect known phylogeographic breaks around the Sichuan Basin? (b) How have the evolutionary history during the Quaternary and the biological traits of P. lasiocarpa affected its phylogeographic pattern? To answer these questions, we used three plastid fragments (ptDNA) and eight nuclear microsatellite loci markers (nSSRs) to examine the genetic structure of 265 individuals of P. lasiocarpa from 21 populations.

2. Materials and methods 2.1. Plant sampling collection

A total of 265 individuals from 21 populations between 816 and 3000 m a.s.l. was sampled, covering most of the distributional range of Populus lasiocarpa around the Sichuan Basin (Fig. 1; Table S1). Between 6 and 17 trees were sampled from each population, with all sampled individuals at least 100 m apart from one another. Fresh leaves were collected in silica gel and voucher specimens of each population were deposited at the Sichuan University Herbarium (SZ).

2.2. DNA extraction, amplification, and sequencing

Total genomic DNA was extracted using the CTAB method (Doyle and Doyle, 1987). We successfully sequenced three plastid markers that show high levels of intraspecific variability (matK, trnG-psbK and psbK-psbI; Table S2; Schroeder et al., 2012) in 193 individuals. A set of eight nuclear microsatellite loci were used to genotype 265 individuals (Table S2). Polymerase chain reaction (PCR) for all ptDNA fragments and nSSR loci was performed in 25 μL volumes, containing 50–100 ng DNA, 0.5 mM of each dNTP, 0.5 μL of each primer, and 1 × Taq buffer and 0.5 units of Taq polymerase (Beijing ComWin Biotech Co., Ltd., Chengdu, China). Reactions were run on a Master Cycler Pro thermal cycler (Eppendorf, Hamburg, Germany) using the following protocol: one cycle at 94 ℃ for 5 min, followed by 34 cycles at 95 ℃ for 30 s, 56 ℃ for 30 s, and 72 ℃ for 50 s, with a final extension at 72 ℃ for 5 min. The PCR products were checked on 1% agarose gels and then sent to Beijing Genomics Institute (Chengdu, China) for DNA sequencing and nSSR genotyping.

2.3. nSSR data analyses

Microsatellite data were scored in GeneMarker (Softgenetics, Pennsylvania, USA) and then corrected by the FlexiBin Excel macro (Amos et al., 2007). Allele size was scored at each locus and checked for possible genotyping errors, allelic dropout, and null alleles using CERVUS v.3.0 (Kalnowskl et al., 2007). Two loci (peussr-56336, peussr-83115) with high-frequency null alleles (F [Null] > 0.4) were discarded, the remaining eight nSSR loci (Table S2) were used to estimate genetic diversity indices employing GenAlEx v.6.5 (Peakall and Smouse, 2012). The Bayesian clustering method implemented in STRUCTURE v.2.3.4 was used to investigate population structure in Populus lasiocarpa (Pritchard et al., 2000) with the following parameters: K from 1 to 10 with twenty replicates each, and 1, 500, 000 MCMC cycles after 500, 000 burn-in cycles. The optimal K value was determined using two methods, the likelihood estimate (Pritchard et al., 2000) and the ΔK statistic method (Evanno et al., 2005). STRUCTURE results were summarized and visualized using Structure Harvester (Earl and vonHoldt, 2012). To cross-validate the results of STRUCTURE, we also conducted a Principal Coordinates Analysis (PCoA) using GenAlEx v.6.5 (Peakall and Smouse, 2012). In addition, to assess the spatial population structure, we conducted a spatial analysis of molecular variance using SAMOVA v.2.0 (Dupanloup et al., 2002). We ran the SAMOVA program 100 times simulating the annealing process for K groups (K = 2–10).

2.4. ptDNA sequence analyses

Subsequently, ptDNA sequences were edited, aligned, manually checked, and concatenated with ClustalW in MEGA v.5.0 (Tamura et al., 2011). Insertions/deletions (indels, excluding mononucleotide repeats) were encoded by the software Gapcoder (Young and Healy, 2003), and the 0/1 characters (absence/presence, except '–' gaps after coding) were replaced manually by A/T to use indel information (Havrdová et al., 2015). Haplotype variant sites were detected using DNASP v.5.10 (Librado and Rozas, 2009). NETWORK v.4.6 was used to infer network relationships between ptDNA haplotypes based on sequence variation (Bandelt et al., 1999). In addition, to assess the spatial population structure, we conducted a spatial analysis of molecular variance using SAMOVA v.2.0 (Dupanloup et al., 2002). We ran the SAMOVA program 100 times simulating the annealing process for K groups (K = 2–10). The optimal number of groups (k) was determined based on the highest variance between groups (F_CT), incorporating information on haplotype divergence and geographical proximity. Finally, to investigate distribution patterns of genetic variation within the ptDNA data set, analyses of molecular variance (AMOVA) were carried out in ARLEQUIN v.3.5 (Excoffier and Lischer, 2010), with significance tests based on 1000 permutations. Genetic variation was hierarchically partitioned into three levels: among groups, among populations within group, and within populations.

To calculate divergence times, we adopted a two-step approach in BEAST v.1.7.5 (Drummond and Rambaut, 2007), using node ages from Zhang et al. (2018) and four species as outgroups: Populus euphratica, P. wilsonii, P. davidiana, and P. adenopoda (Zhang et al., 2018).

2.5. Genetic barrier analyses

Monmonier's maximum-difference algorithm in BARRIER v.2.2 (Manni et al., 2004) was used to identify biogeographic boundaries or areas exhibiting the largest genetic discontinuities between population pairs. We used nSSR and ptDNA data for genetic barrier analysis. The robustness of these boundaries was assessed by running BARRIER on 100 replicates of population average pairwise difference matrices. These matrices were generated by bootstrapping of haplotype sequences in SEQBOOT (Felsenstein, 2005) and subsequent analyses in ARLEQUIN v.3.5 (Excoffier and Lischer, 2010). One thousand Nei's genetic distance matrices were used to estimate the robustness of gene flow barriers. Nei's genetic distance matrices were produced by bootstrapping loci with microsatellite analyser v.4.05 (Dieringer and Schlötterer, 2003). The genetic groups were distinguished based on genetic discontinuities and were then compared with genetic clades/clusters identified by the phylogenetic tree and STRUCTURE analyses.

2.6. Geographical and environmental effects on genetic divergence

Multiple processes such as isolation-by-environment (IBE) and isolation-by-distance (IBD) may affect genetic variation of populations (Sexton et al., 2014). IBD investigates the correlation between spatial and genetic distance of populations (De Matthaeis et al., 2000; Wang, 2013; Pelletier and Carstens, 2018). IBE implies that genetic divergence increases with environmental heterogeneity (Wang and Bradburd, 2014). Here, three models were examined (Table 1): (a) IBDe, where the geographical distance between populations was represented by a Euclidian distance assuming that Populus lasiocarpa disperses across the basin (coordinates of sample sites were used to calculate "as the crow flies" distances with the pointDistance function in the raster R package; Hijmans and van Etten, 2012); (b) IBDr, the geographical distance between populations was represented by a ring distance (calculated using the ArcMap v.10.2.1), assuming no gene flow across the unoccupied areas across the Sichuan basin; (c) IBE, environmental difference was calculated based on Euclidean distances of the SDM environmental factors (see 2.8. Species distribution modelling (SDM) for details) at each population's locality.

We used the philentropy R package (Drost, 2018) to calculate the Euclidean distance between the different environments. Pairwise F_ST values between populations were calculated in ARLEQUIN v.3.5 (Excoffier and Lischer, 2010). Finally, Spearman's correlation between matrices of genetic, geographical, and environmental distance was assessed by a Mantel test using the Mantel function of the R package vegan (Oksanen et al., 2020) with 999 permutations.

2.7. Demographic history analyses

Mismatch distribution analysis examines the distribution of the observed number of differences between pairs of haplotypes. This distribution is usually multimodal in samples drawn from populations at demographic equilibrium, reflecting the highly stochastic shape of gene trees, whereas it is usually unimodal in populations having passed through a recent demographic or range expansion (Slatkin, 1991; Rogers and Harpending, 1992) with high levels of migration between neighboring demes (Ray et al., 2003). Mismatch distribution was calculated in ARLEQUIN v.3.5 (Excoffier and Lischer, 2010) using the sum of square deviations (SSD) between the observed and the expected mismatches as a test statistic. Arlequin assesses the probability of expansion using the two parameters SSD (Sum of square deviations) and r (Harpending's raggedness index), where the hypothesis of group expansion cannot be rejected if the statistical tests for these two parameters are not significant (p > 0.05).

We used BOTTLENECK v.1.2.02 (Piry et al., 1999) to detect past bottlenecks in three Populus lasiocarpa populations (the eastern, southern and western populations; Table 2) choosing a two-phase model (Di Rienzo et al., 1994) and a two-tailed Wilcoxon signed-rank test to test for significance.

We used the coalescent-based approximate Bayesian computation (ABC) implemented in DIY-ABC v.2.0 (Cornuet et al., 2014) to infer the demographic histories of Populus lasiocarpa. To alleviate the impact of inter-lineage gene flow on the testing of demographic history scenarios, we excluded all individuals that were potential hybrids (according to Fig. 3c, K = 2) from the southern populations and analyzed the demographic history of the eastern and western populations, respectively. Because mismatch distribution analysis and bottleneck analysis indicated a significant expansion and bottleneck in P. lasiocarpa (Fig. S3; Table 2), we tested six models of population size changes (Fig. 2): an old shrinkage followed by expansion (scenario 1: "shrinkage-expansion"); an old expansion followed by shrinkage (scenario 2: "expansion-shrinkage"); a recent expansion followed by shrinkage (scenario 3: "expansion-shrinkage"); expansion followed by shrinkage and a new expansion event (scenario 4: "expansion-shrinkage-expansion"); expansion followed by shrinkage with a small population size and a new expansion event (scenario 5: "expansion-shrinkage-expansion"); a recent shrinkage followed by expansion (scenario 6: "shrinkage-expansion"). We choose the demographic scenarios that best fit the data and estimated the parameters of interest (Table S3). Using a direct approach and logistic regression analyses, the posterior probabilities of all scenarios were calculated and compared. Following Macaya-Sanz et al. (2014) we set the generation time of P. lasiocarpa to 40 years, with a span of 20–60 years.

2.8. Species distribution modelling

Species distribution modelling (SDM) was conducted using the maximum entropy method implemented in MAXENT v.3.2.1 (Phillips et al., 2006) to predict the distribution of potentially suitable habitat for Populus lasiocarpa in five time periods: the Last Glacial Maximum (LGM; ca. 21-thousand years ago (Kya)), the Last Interglacial (LIG; ca. 130 Kya), the Mid-Holocene (MH; ca. 6 Kya), the present, and 2050. Nineteen bioclimatic variables at 2.5 arc-minute resolution were downloaded from the WorldClim database (Fick and Hijmans, 2017). Strong co-linearity between bioclimatic variables may affect the accuracy of the model. Therefore, a Pearson correlation test was performed on these bioclimatic variables across the 21 sampled populations. Seven bioclimatic variables between which all pairwise Pearson correlation coefficients r were > 0.70 (Table S4) were retained and used for subsequent SDM analysis and IBE analyses. The maximum entropy model was simulated for 20 replicates, randomly selecting 20% of the distribution data as the test data, whereas the remaining 80% of the distribution data were used for the training simulation. The maximum number of iterations was set to 5000, and the convergence limit was set to 0.01. The output format was set to "logistic" and the bioclimatic conditions for each grid cell were probabilistic in the range of 0–1. The performance of models was evaluated by checking the area under the Receiving Operating Characteristics (ROC) curve (AUC) value. DIVA-GIS v.7.5 (Hijmans et al., 2001) was used to map the distribution of habitat suitability.

3. Results 3.1. Nuclear genetic differentiation

Using eight nSSR loci, STRUCTURE yielded the highest likelihood for K = 2 (Fig. 3a and b), suggesting the existence of two genetic clusters (Fig. S2a), which occur mainly in the western (populations 1–8) and eastern groups (populations 9–21). For K = 3, there was an eastern group (populations 10 and 12–21), a southern group (populations 7–9 and 11) and a western group (populations 1–6) (Fig. 3c). The PCoA based on individual genetic distances showed a similar population genetic structure, with the southern populations (7, 8, 9, and 11) clustered between the eastern and western populations (Fig. S1b). Spatial genetic analysis indicated that F_CT increased to its maximum value of 0.231 at K = 3 (Table S9a). For K = 4 (Figs. 3c and S2b), three populations (2, 3 and 6) in the western group were genetically differentiated to some extent from the other three populations (1, 4 and 5). In the PCoA based on population genetic distance, the distribution of these populations was concordant with geography and generally followed a ring pattern around the barrier (Figs. S1a and S2b).

Averaged over all populations, the mean number of alleles was Aa = 2.73, the mean number of effective alleles was Ae = 1.83, the observed heterozygosity was Ho = 0.25, the expected heterozygosity was He = 0.28, and the inbreeding coefficient at the population level was F_IS = 0.14 (Tables S5 and S6). The genetic differentiation index was F_ST = 0.27 (Table S5).

3.2. Genetic variation of ptDNA sequences

We obtained 1688 bp of concatenated plastid sequences (matK 797 bp; trnG-psbK 476 bp; psbK-psbI 415 bp). A total of 12 haplotypes were detected (Tables S7 and S8), revealing two distinct haplotype clades (Fig. 4c). Clade I comprised all blue haplotypes (H2, H3, H5, H7, H8 and H10). Among these, H2 is the most common haplotype and is present in all but one population (Fig. 4a), indicating that it is the ancestral haplotype of Populus lasiocarpa. Clade II consists of yellow haplotypes (H1 and H4) and red haplotypes (H6, H9, H11 and H12). The yellow haplotypes are mainly represented in the southern populations and the red haplotypes in the western populations (Fig. 4a). Molecular dating suggests that the blue haplotypes diverged from all other clades at ca. 3.66 Ma (Fig. 4c), and that the red diverged from the yellow haplotypes at ca. 2.09 Ma (Fig. 4c). Spatial genetic analysis of ptDNA haplotypes indicated that F_CT increased to its maximum value of 0.631 at K = 3 (Table S9b). The grouping pattern of populations corresponding to K = 3 was as follows: 1) population 5; 2) populations 7–9, 11–12 and 17; 3) populations 1–4, 6, 10, 13–16 and 18–21. Based on the SAMOVA groups, AMOVA showed that 75.16% (p < 0.001) of the total variation occurred among groups (Table S10).

3.3. Genetic barriers

Our BARRIER analysis based on both ptDNA data and nSSR data showed a strong genetic barrier between the southern and western sampling areas (Fig. 5) both received high bootstrap values (100%). In addition, we observed one genetic barrier between populations 2 and 7 (Fig. 5) that had a high support value (100%). There are also other barriers between peripheral populations according to both ptDNA and nSSR data (Fig. 5).

3.4. Geographic and environmental effects on genetic divergence

The Mantle test based on nSSR data showed that both geographical distance and environmental variables played a considerable role in the differentiation of populations (Table 1). However, the fit was significant when the geographical variable was based on ring (r = 0.635, p = 0.001) rather than Euclidean distances (r = 0.462, p = 0.001), with a significant correlation between genetic variation and ring distance. In addition, compared to the environmental variables (r = 0.554, p = 0.001), the ring distance (r = 0.635, p = 0.001) delivered a stronger effect on the differentiation of populations.

3.5. Demographic history

The mismatch distribution of Populus lasiocarpa was unimodal (Fig. S3), fitting the expected distribution under the sudden expansion model, supporting a significant expansion of P. lasiocarpa (p [SSD] = 0.44, p [RAG] = 0.19). There was a significant heterozygosity excess (Table 2, p < 0.05) under all three models in the western, southern and eastern populations, indicating a putative bottleneck that might have occurred in those population.

In the DIYABC analyses, the scenarios best supported based on direct estimate, logistic regression and PCA plots were "expansion-shrinkage-expansion" (scenario 5) for the eastern populations and "shrinkage-expansion" (scenario 6) for western populations (Figs. S4 and S5). These scenarios were also the least error prone (Table S11). The "expansion-shrinkage-expansion" (scenario 5) best-fit for the eastern populations indicates that Populus lasiocarpa initially expanded before the LGM period but not earlier than the LIG period, ca. 62 Kya (95% CI: 20.2–88.8 Kya), then experienced an obvious reduction in population size (ca. 11.6 Kya), and expanded again before the MH, ca. 6.7 Kya (95% CI: 1.5–11.3 Kya) (Table S12). By contrast, the "shrinkage-expansion" (scenario 6) best-fit for western populations indicates that a strong bottleneck occurred during the postglacial period, ca. 12.0 Kya (95% CI: 5.1–22.6 Kya), and expanded again before the MH, ca. 6.4 Kya (95% CI: 0.47–11.6 Kya) (Table S12).

3.6. Species distribution predicted by modelling

The predictions of potentially suitable distribution areas under the present and historical climate conditions are shown in Fig. 6. All models have high AUC test values (> 0.95), indicating the good performance of all predictions. Variable jackknife analyses suggested that Max Temperature of Warmest Month (56.4%) and Mean Monthly Temperature Range (22.2%) were the environmental variables that contributed the most to the distribution modelling for Populus lasiocarpa (Table S4).

The projected distribution of this species at present encompasses most of the sampling locations. During the LIG period, the range of P. lasiocarpa was narrow and scattered, relative to the current distribution, and it seemed to have been restricted mainly to the west and east of the Sichuan Basin, with a more fragmented distribution in the east (Fig. 6b). The range of P. lasiocarpa expanded from the LIG to the LGM (Fig. 6c) including extended distribution areas northwards and southwards of the Sichuan Basin. There was no obvious expansion or contraction for the distributions of P. lasiocarpa under present (Fig. 6e) and MH (Fig. 6d) climate conditions, yet only a very tiny distribution of the species was predicted for 2050 (Fig. 6f).

4. Discussion

In this study, we used a suite of phylogeographic analytical tools to investigate the evolutionary history of Populus lasiocarpa, a tree species that occurs in the mountains surrounding the Sichuan Basin. Bayesian clustering based on nSSR markers (Fig. 3) indicated that the populations of P. lasiocarpa were divided into three population groups (eastern, southern and western) by three phylogeographic breaks, the Sichuan Basin, the Kaiyong Line, and the 105°E line (Fig. 3d). Distribution patterns based on ptDNA haplotypes poorly matched these phylogeographic breaks (Fig. 4a). A Mantle test suggested that both geographical distance and environmental variables played a considerable role in the differentiation of populations, especially using the ring distance (Table 1). In addition, SDM analyses indicated a much more restricted distribution during the LIG than at present, and a slightly greater area at the LGM than that predicted under current climatic conditions (Fig. 6). Finally, the DIYABC analyses suggested a postglacial range contraction and a middle Holocene expansion for both western and eastern populations (Table S12).

4.1. The influence of multiple forces on the contemporary genetic structure

Results from STRUCTURE analysis using nSSR makers (Fig. 3) indicated that the Sichuan Basin acted as a phylogeographic boundary that divided Populus lasiocarpa populations into two genetic groups, one to the east and another to the west of the basin. In addition, our analyses showed that this geographical barrier led to a ring diversification pattern in P. lasiocarpa. The mountain ranges surrounding the Sichuan Basin in southwestern China are one of 100 top candidate regions for potential ring species (Monahan et al., 2012). Few polytypic taxa are in fact recognized by modern taxonomy as ring species (Irwin et al., 2001). Therefore, "ring diversification" (Monahan et al., 2012) is used to describe the situation in which populations diverge around a geographical barrier in a ring-like manner, regardless of whether reproductive isolation occurs at the terminus, such as for Hyla orientalis (Dufresnes et al., 2016), Odorrana margaratae (Qiao et al., 2018), Tetracentron sinense (Sun et al., 2014), and Apodemus draco (Wang et al., 2021). The populations of P. lasiocarpa are distributed in the mountains surrounding the Sichuan Basin in a ring like fashion (Figs. 3d, 4a, S1, and S2). The results of the SDM clearly show that the Sichuan Basin must have been an unsuitable habitat for P. lasiocarpa in historical and contemporary times, preventing gene flow across the basin (Fig. 6). The genetic divergence of P. lasiocarpa showed a higher correlation with the ring distance than with the Euclidean distance and environmental variables (Table 1), indicating limited dispersal across the barrier. Thus, the Sichuan Basin, which has been a strong geographical barrier since the middle Miocene (~13 Mya) (Liu and Xu, 2003), not only divided P. lasiocarpa into two genetic groups located in the western and eastern parts of the basin but caused the formation of a ring-shaped distribution.

Three groups showed strong isolation signals in the STRUCTURE analyses (K = 3) and PCoA (Fig. S1) based on the nSSR data. The SAMOVA groups based on ptDNA haplotypes and nSSR data also indicated that K = 3 was the best grouping pattern (Table S9). However, the groupings based on ptDNA poorly matched that of nSSRs as the same ptDNA haplotypes, H2 and H1, occur in two or more nSSR groups. The barrier effect of the Sichuan Basin is so strong that pollen and seeds rarely cross it. Hence, in addition to seed dispersal, shared ancestral polymorphism may also explain this pattern. Regardless, nSSR results suggested there are three geographical genetic groups, located at the western, southern, and eastern portions of the basin edge (Figs. 3d and S1a). These groups approximately reflect three subregions: (a) the Wushan and Daba Mountains at the eastern edge of the Sichuan Basin; (b) the Great Lou Mountains and the Great Liangshan Mountains at the southwestern edge of the basin; and (c) the Great Snowy Mountains, Qionglai Mountains and Minshan Mountains, located on the edge of the Hengduan Mountains and the QTP. Moreover, BARRIER analysis showed a strong genetic barrier between populations 2 and 7 based on both nSSR and ptDNA data (Fig. 5), which coincides with the Kaiyong Line. We observed one genetic barrier between populations 11 and 12 that coincided with the 105°E line (Fig. 2d), but it had a low support value (Fig. 5b). This low support value may be explained by the abnormality that occurs when genetic and geographic distances are coupled, as BARRIER analysis is less sensitive to barriers between long-distance populations. Interestingly, the geopositions of three genetic groups in Populus lasiocarpa correspond to three phylogeographical breaks, including the Sichuan Basin, the Kaiyong line and the 105°E line. Therefore, for P. lasiocarpa, the Sichuan Basin is a stronger geographical barrier, and there are also important genetic barriers at the Kaiyong line and the 105°E line.

Molecular dating in our study suggested that the intraspecific divergence time of haplotypes H1, H4 and haplotypes H9, H11, H12 on both sides of the Kaiyong line occurred approximately 2.09 Mya (Fig. 4c). This intraspecific divergence time coincided with uplift events of the QTP, which are known as the Qingzang Movement (~3.6–1.7 Mya) (Li and Fang, 1999). The intraspecific divergence time of A.podemus draco surrounding the Sichuan Basin (1.7 Mya) also coincided with the tectonic events of the QTP (Wang et al., 2021). The spatial and temporal concordance among these species suggests that tectonic events in the QTP were the same drivers in building their genetic differentiation. Thus, our results support the hypothesis that the historical tectonic events may be factors contributing to the formation of genetic barriers at the Kaiyong line. At the same time, we cannot rule out the possibility of monsoon driven hypothesis, because the latitude and the annual precipitation of western populations is higher than those of the southern populations; however, we found that a glacial driven hypothesis is less likely considering the ecological niche model of the LGM (Fig. 6c), where the potential distribution is even larger compared to the Current model (Fig. 6e).

Similar to Populus lasiocarpa, research on phylogeography of plants distributed around the Sichuan Basin has indicated genetic differentiation of some plants at the 105°E line, including Fagus hayatae (Ying et al., 2016), Liriodendron chinense (Yang et al., 2019), and Davidia involucrata (Ma et al., 2015). Notably, SDM analysis predicted low habitat suitability and gaps in the distribution of P. lasiocarpa around the 105°E line during three time periods (LGM, MH, and the present). The genetic differentiation of P. lasiocarpa at the 105°E line might be the result of different climatic and topographical processes on either side (Qiu et al., 2011; Ye et al., 2017).

Interestingly, our nSSRs data suggest the existence of a gene flow barrier at the 105°E line between populations 11 and 12 (Fig. 2d), yet a strong genetic barrier between populations 12 and 13 was supported by the ptDNA data (Fig. 5a). In plants, seed dispersal is often considerably less effective than pollen dispersal (Petit et al., 2005). Therefore, maternally inherited organelle markers (ptDNA or mtDNA) that are only dispersed by seeds should be more frequently introgressed (Petit and Excoffier, 2009). Moreover, due to the prevalence of sex-biased dispersal (Greenwood, 1980), different DNA regions are subject to varying levels of gene flow as a consequence of their mode of inheritance (biparental, maternal or paternal). Based on our nSSR data, we speculate that pollen and seeds of Populus lasiocarpa likely crossed this phylogeographic break at the 105°E line with the help of the monsoons from both the Indian and Pacific oceans. We suggest the idea that nuclear molecular markers, which have higher rates of gene flow, might be better indicators of phylogeographic breaks than organelle markers, which have lower rates of gene flow, in part because higher rates of intraspecific or local gene flow can prevent interspecific or between-lineage introgression (Petit and Excoffier, 2009). This is evident in P. lasiocarpa because wind-dispersed seeds probably facilitate introgression of ptDNA between geographical genetic groups compared to seeds dispersed by gravity, and large effective population size of trees favor the retention of these introgressed haplotypes. In addition to ancestral polymorphism and wind-dispersed seeds, hybridization-induced plastid capture events between closely related species may also have influenced the distribution pattern of ptDNA haplotypes. Wang et al. (2020) identified some gene flow and hybridization between P. lasiocarpa and P. wilsonii. The evolutionary history of P. lasiocarpa and its relative species can be studied in depth in the future using genomic-based approaches to better reveal hybridization and resulting introgression.

4.2. Infrequent demographic histories of Populus lasiocarpa

Our SDM analyses indicated a much more restricted distribution of Populus lasiocarpa during the LIG than at present, and a slightly greater area at the LGM than that predicted under current climatic conditions (Fig. 6). The mismatch distribution supported a significant expansion (Fig. S3), and the bottleneck analysis indicated that P. lasiocarpa experienced a recent genetic bottleneck (Table 2). Consistent with this, the DIYABC analysis suggested that both western and eastern groups of P. lasiocarpa experienced a bottleneck 11, 600 to 120, 000 years ago, which matches the Younger Dryas (Cheng et al., 2020), an abrupt increase of paleotemperature during postglacial periods; a demographic expansion was also detected for both western and eastern groups around the middle Holocene, 6440–6720 years ago, which is consistent with previous observations that paleotemperature started to decline across much of the globe around 5000–7000 years ago (Steig, 1999).

These findings contradict the scenario where most plants experienced expansion during the LIG and then contraction during the LGM in southern China (Qu et al., 2009; Hofmann, 2012; Zhou et al., 2013; Wang et al., 2017, 2018). However, the climatic changes in the Sichuan Basin were not as pronounced as in other regions (such as Europe and North America) (Jablonski, 1993; Kose, 2014), and some species even had a constant distributional range during this period (Chen et al., 2015; Cao et al., 2016; Wang et al., 2021). As refugia, the mountains of Southwest China contributed to a more stable distribution pattern by providing diverse micro-environments that prevented over-heating or over-cooling (Wen et al., 2017). Paleoclimatic data have indicated that during the LIG the temperature was at least 5 ℃ higher than at present (Shi et al., 1999; Li and Pan, 2002; Li and Kang, 2006). For some cold-adapted or moisture-loving plants, climate warming such as that during the LIG or postglacial period could have led to range contraction instead. Taxus wallichiana and Picea likiangensis in the adjacent areas of the QTP have already been affected in this way in that their current predicted distributions have been reported to be smaller than those predicted at the LGM (Li et al., 2013; Liu et al., 2013). Biological characteristics of Populus lasiocarpa indicate its sensitivity to high temperatures: P. lasiocarpa has the largest leaf area in poplars, about 15–30 cm long and 12–20 cm wide (Hong et al., 1987), which makes it intolerant of drought and high temperatures due to strong transpiration and respiration (Rozendaal et al., 2006; Picotte et al., 2007; Li et al., 2011); in addition, variable jackknife analyses suggest that Max Temperature of Warmest Month (56.4%) was the environmental variable that contributed most to potential distribution modelling for P. lasiocarpa (Table S4). Therefore, P. lasiocarpa may have experienced a bottleneck event during the postglacial period, and most likely also during the LIG (Fig. 6b), when the paleoclimate was warmer, in contrast to the more stable and milder environment of the mountains surrounding the Sichuan Basin during the LGM that P. lasiocarpa populations might have favored. The demographic expansion of western and eastern groups of P. lasiocarpa following paleotemperature decline around the middle Holocene also suggests that this species favors mild environments. Hence, considering that climate warming during the LIG had a marked effect on the predicted distribution of P. lasiocarpa, coupled with a predicted range contraction of the distribution of this species in 2050 (Fig. 6f), population decline and range contraction of this species in the future cannot be ruled out. Taken together, our findings indicate that the conservation and management of this remarkable and beautiful poplar is an imperative task.

Acknowledgements

We thank Jia-Liang Li and Wen-Jing Tao for their assistance in sample collection. This project was supported by National Natural Science Foundation of China (grants 31971567 and 31622015) and Fundamental Research Funds for the Central Universities (YJ201936, SCU2020D003, SCU2021D006, SCU2022D003).

Author contributions

K.-S.M. and J.W. designed the study. H.-Y.Z., X.-Y.F, X.L. and L.Z. performed the field work; X.L. and H.-Y.Z. analyzed the data; X.L. and Y.W. wrote the initial draft of the manuscript. M.R. contributed new reagent to the manuscript. K.-S.M, X.L., M.R. and J.W. finalized the manuscript, and all authors contributed to interpretation/discussions of the data and revisions of the manuscript. All authors approved the final manuscript for publication.

Data availability statement

All sequences have been deposited in GenBank (accession nos. OL878409–OL878601 for psbK-psbI, OL878602–OL878794 for trnG-psbK, OL878795–OL878987 for matK). Microsatellite data and ptDNA sequences data have been deposited in the Science Data Bank and are available at http://doi.org/10.57760/sciencedb.j00143.00020.

Declaration of competing interest

The authors have no conflict of interest.

Appendix A. Supplementary data

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