The current study was restricted to the Longitudinal Range Gorge Region (LRGR). According to the boundary defined by Pan et al. (2012), the LRGR covers an area in southwest China that spans ca. 10.33 degrees of longitude (106.12–95.79°E) and ca. 9.45 degrees of latitude (30.60–21.15°N; Figure 1). Based on global 3 arc second (ca. 90 m) SRTM data (Farr et al., 2007), we estimated the average elevation of the LRGR to be 2369.06 m (above mean sea level), ranging from 76.4 to 6,740 m, and a total area of 410,251.64 km². It spans both the subtropical and tropical monsoon climatic zones, with some parts of the high-elevation LRGR being temperate. Hence, temperatures vary widely with latitude and elevation, with the area supporting a diverse altitudinal zonation including tropical rain forests, subtropical evergreen broadleaf forests, and alpine vegetation.

We retrieved published articles on population genetics and phylogeography from Google Scholar.¹ The literature search was conducted on 11/01/2019 using the terms “Hengduan Mountains” OR “southwest China” OR “Longitudinal Range Gorge Region” OR “Yunnan” OR “Three Parallel Rivers region” AND “plant” AND “phylogeography” AND “genetic structure” AND “cpDNA.” Among the retrieved articles, for consistency, we only considered studies that were based on chloroplast DNA (cpDNA) sequence markers; cpDNA sequences have been widely applied in population genetics and phylogeographic studies of plant species in the region; and hence, georeferenced population-level data are more easily accessible. Previously, cpDNA markers have been successfully applied to map genetic diversity patterns and clarify determinants thereof (Fan et al., 2018; Yu et al., 2019). Furthermore, studies in which the majority of the species range (>65%) fell outside the LRGR were excluded. Following these filtering steps, a total of 13 published studies, reporting genetic data for 15 plant species (Figures 1, 2; Table 1), met the selection criteria. For each of the populations of the 15 species, we gathered the geographic co-ordinate (longitude and latitude) from the articles and extracted their corresponding elevation data using the Extract Multivalues to Points toolbox in ArcGIS Pro v2.8.2.² The referenced haplotype sequences were then downloaded from GenBank,³ using the accession numbers provided in the source studies.

Based on the haplotype frequencies for each population and the mutation information in the haplotype sequences, we reconstructed the population and species sequence matrices in Geneious v2020.0.3⁴ using the MAFFT alignment algorithm (Katoh and Standley, 2013), and where necessary, alignments were trimmed to be of equal length. To ensure the validity of our population genetics analyses, only populations with a sample size of N ≥ 5 were retained for downstream analyses.

DnaSP v5.10 (Librado and Rozas, 2009) was used to define sequence sets, i.e., populations and to generate haplotype files for each of the 15 species. Subsequently, Arlequin v3.5.2 (Excoffier and Lischer, 2010) was used to calculate haplotype diversity H_D (Nei, 1987), nucleotide diversity π (Tajima, 1983; Nei, 1987), and population pairwise genetic differentiation F_ST (Tamura and Nei, 1993). As we detected a significant correlation (r = 0.728; p < 0.01) between H_D and π, we only used H_D in the subsequent analyses, since its calculation is insensitive to sample size, evolutionary forces, mode of reproduction, and ploidy level (Nei, 1973).

The population-level genetic diversity (H_D) and pairwise divergence (F_ST) values were used to generate genetic diversity and divergence landscapes using a modification of the Genetic Landscapes GIS Toolbox (Vandergast et al., 2011), which we implemented in ArcGIS Pro v2.8.2 with Python scripts. The landscape surfaces were generated at a resolution of 2.5 min (~5 km at the equator) based on the inverse distance weighting (IDW) interpolation method. Ordinarily, the interpolation method in the Genetic Landscapes Toolbox is based on the Minimum Bounding Geometry (convex hull type) for each species distribution point, which invariably leads to stretching of the boundary beyond the area where actual geometry exists. To resolve this challenge, we clipped the individual species landscapes to the boundary of the study region prior to combining them in the multi-species landscape, thus generating a multi-species landscape in which all the 15 species landscapes overlap throughout the study area. Subsequently, genetic diversity “hotspots” were defined as regions with standard deviations of genetic diversity >1.5 from the mean (Vandergast et al., 2008). Next, the geographic coverage of these hotspots was compared with that of protected areas, based on data obtained from the World Database on Protected Areas (WDPA),⁵ and validated with available information on China’s national nature reserves (You et al., 2018). To confirm whether the observed genetic differentiation was related to geographic distance, we performed isolation-by-distance (IBD) analysis for each species by testing correlation of the Euclidean geodistance with population pairwise genetic differentiation (F_ST) using Mantel tests in the vegan package in R (Oksanen et al., 2011).

A total of 66 explanatory variables (2 anthropogenic, 57 climatic, and 7 topographic variables) were extracted from various databases at various spatial resolutions (Supplementary Table S1). To avoid collinearity among the extracted variables within a given category, we omitted one variable for every pair that exhibited a high correlation coefficient (r > 0.7). After this, a total of 16 explanatory variables (10 climatic, 4 topographic, and 2 anthropogenic) were retained for subsequent analysis (Supplementary Table S2).

The response variable values (H_D) were extracted from the extrapolated plastid genetic diversity landscape at the same spatial resolution as the corresponding explanatory variables (2.5 min resolution; Supplementary Table S2). To avoid bias that would be caused by the uneven number of variables across the categories, we selected only the two best variables in each of the climatic and topographic categories using the importance function within the Random Forest model implemented in the randomForest R package (Liaw and Wiener, 2002), plus the sole two in the anthropogenic category. Here, importance was defined as the degree to which the inclusion of each predictor decreases the residual model variance. We then carried out regression analyses on the selected variables with the corresponding H_D values based on generalized additive models (GAM) with integrated smoothness estimation using the “gam” function in the mgcv R package (Wood, 2017). Subsequently, the relative contributions of climatic, topographic, and anthropogenic factors, as well as their overlaps, in explaining the observed patterns of genetic diversity, were determined through redundancy analysis (RDA) and variation partitioning (Borcard, 1992) using the vegan R package.

Population-level genetic diversity indices varied widely, with nearly half (49.7%) of the populations lacking genetic variation (i.e., comprising a single detected haplotype), and ~5% showing maximum haplotype diversity (H_D = 1, meaning every sampled individual had a different haplotype; Supplementary Table S2). After interpolation on the LRGR landscape surface, the average genetic diversity across the 15 species ranged from 0.25 to 0.47 (Figure 3A). Our results indicated that the northern area of the LRGR generally harbors higher levels of genetic diversity, with the maximum haplotype diversity residing in the Yulong Mountain region, near the first major meandering of the Yangtze River. Relatively lower genetic diversity was observed toward the eastern edges of the study region (Figure 3A).

The interpolated genetic divergence (F_ST) ranged from 0.39 to 0.59 and showed high divergence approximately along the 25°N latitudinal line, with notable peaks at the western and eastern ends of this line (Figure 3B). Conversely, relatively low genetic divergence between populations was observed for populations in the northern and southern peripheries of the study area. Based on IBD analysis, only 5 of the 15 species showed significantly positive correlation between geographic and genetic distances, i.e., Leucomeris decora (r = 0.3, p = 0.04), Ligularia vellerea (r = 0.47, p = 0.002), Quercus kerrii (r = 0.19, p = 0.01), Rosa soulieana (r = 0.39, p = 0.03), and Sophora davidii (r = 0.65, p = 0.001; Supplementary Table S3).

Based on the set threshold of 1.5 standard deviations (which corresponded to an H_D value of 0.38), we identified seven genetic diversity “hotspots” (labeled A-G in descending order of size) having a total area of 15,686.81 km² (Figure 4). Of these, the largest (A) was located in the Yulong Mountain region and the second largest (B) near Meili Mountains. The remaining five were much smaller, and of these, hotspot D was located in Nushan Mountain, whereas the other 4 were distributed at lower latitudes in the southern LRGR (Figure 4). We found that about 74.81% of the genetic diversity “hotspots” land area was located outside of the established protected areas (Figure 4).

The model selection analysis showed that cld/wet and lat/long were the best-fitting models in the climatic and topographic categories, respectively, and these were therefore combined with the two anthropogenic variables hfp and ahf for the RDA and variation partitioning. The selected variables showed weakly negative (long, cld, ahf, and hfp) and weakly positive (lat and wet) correlation with genetic diversity, though the relationships were significant in all cases (p < 0.001; Figure 5). Despite these general trends, we noted slight deviations in genetic diversity along the scales of some explanatory variables. For instance, although genetic diversity generally increased with increasing lat (R_adj.² = 0.19), there was a marked decrease in genetic diversity at ~22.5–25° and at 28.5°. A similar trend was observed for cld and hfp, in which genetic diversity increased slightly before decreasing.

In the subsequent RDA, the three categories of variables explained 78% of the total genetic variation, with topography alone accounting for the highest proportion (40%) of the variation. Climatic and anthropogenic factors explained 3 and 1% of the observed genetic variation, respectively (Figure 6). Notably, the variable combinations of climate/topography and topography/anthropogenic effects explained 14 and 11% of the total genetic variation, respectively.

We mapped patterns of plastid genetic diversity and divergence in the LRGR and determined the relative contributions of climatic, topographical, and anthropogenic factors in shaping the observed diversity patterns. Our results indicated that, within the LRGR, the Yulong Mountain region harbored the highest genetic diversity (Figures 3A, 4). The Yulong Mountain is located on the southeast edge of the Hengduan Mountains, a region that is known for its high species richness (Wu, 1987; Mittermeier et al., 2004) and glacial refugia (Liu et al., 2013; Jian et al., 2015). Indeed, the general distribution of genetic variation in the LRGR was reflective of the congruence between species and genetic diversity patterns. In our study, most populations were concentrated in the vicinity of the Yulong Mountain (approximately at the intersection of latitude 27°N and longitude 100°E). The high species richness of the northern LRGR is consistent with the view that orogeny (which was more intense toward the northern LRGR) was responsible for species diversification in the Himalaya/Hengduan Mountains region (Wambulwa et al., 2021). The observed congruence in species and genetic diversity patterns is a well-known phenomenon (Vellend and Geber, 2005; Vellend et al., 2014) and has been observed in both plant and animal communities (Fan et al., 2018; Hu et al., 2021).

The location of genetic diversity “hotspots” in the northern LRGR may also be explained by the fact that the Hengduan Mountains served as glacial refugia for many plant species in the region (Qiu et al., 2011 and references therein). Due to long isolation periods, populations located in refugia would be expected to harbor relatively high genetic diversity (Abbott et al., 2000; Hewitt, 2000). Since the distribution range of some of the species in our study extends northwards to the eastern Himalaya and Tibetan Plateau (e.g., Wang et al., 2011; Yang et al., 2012; Fan et al., 2013; Jian et al., 2015; Luo et al., 2018), we argue that the high genetic diversity areas observed in the current study represent the glacial refugia for some of the analyzed species (e.g., Qiu et al., 2011; Wang et al., 2011; Fan et al., 2013; Liu et al., 2013; Jian et al., 2015). Although the current study offers new insights into the distribution of genetic diversity in the region, a considerable proportion of the populations in our study lacked genetic variation, and this might have decreased the statistical power during the generation of genetic diversity landscapes and the subsequent extraction of the response variable values.

Moreover, the uneven distribution of populations in our study area might also have biased F_ST estimates. However, we found that the majority of the species did not show isolation-by-distance, suggesting the involvement of environment-driven differentiation in shaping the genetic structure of the majority of the species. The genetic divergence landscape in our study indicates the presence of gene flow barriers along a latitudinal range located approximately centrally in the LRGR, which is generally consistent with the findings of Fan et al. (2018), who mapped the genetic divergence of woody seed plants in subtropical China. Despite the insights obtained regarding the patterns of genetic diversity and divergence in the LRGR, our study may have been limited by the choice of molecular markers (cpDNA). The conserved nature of the cpDNA genome, coupled with its complex evolutionary patterns that arise from incomplete lineage sorting and interspecific introgressions, may decrease the capacity of cpDNA markers to discriminate disparate populations of a particular species.

Understanding the forces that shape the genetic structure of species across a given landscape is a fundamental goal in landscape genetics. Our analyses indicated that genetic diversity was weakly but significantly correlated with both latitude and longitude, tending to increase northwards and westwards (Figure 3A). The northern parts of the LRGR (located in the Hengduan Mountains), which harbored high levels of genetic diversity in the present study, are known for their high topographic heterogeneity. This heterogeneous landscape creates numerous and diverse within-population environmental niches that ultimately lead to strong patterns of isolation-by-environment (IBE; e.g., Liu et al., 2013). On the flip side, mountain ridges and valleys in such heterogeneous landscapes may also provide dispersal corridors among populations (Yu et al., 2017), as supported by the low levels of genetic divergence within the northern part of the LRGR in the current study (Figure 3B). The north-to-south latitudinal genetic cline is consistent with the general topology of the LRGR, i.e., the northern part of the region is characterized by extreme relief with deep river valleys while the southern part exhibits lower relief. Our results are also consistent with the idea that downwind populations are likely to harbor higher genetic diversity, as the rate of accumulation of genetic variation (caused by higher rates of net immigration) in such populations is higher than the variation that is lost due to selection or genetic drift (Kling and Ackerly, 2021). The South Asia Monsoon (Indian Summer Monsoon) blows over the LRGR in the general south to north direction; thus, a latitudinal gradient of genetic diversity would be expected in this region, although this effect may be limited to the wind-dispersed species (Q. kerrii and Q. schottkyana).

A similar trend in genetic diversity along a latitudinal gradient was recently found in subtropical evergreen oaks in southwest China (Xu et al., 2020). However, a meta-analysis of Tibetan plant species (Yu et al., 2019) did not find significant relationships between genetic diversity and either longitude or latitude, whereas Fan et al. (2018) found a positive correlation between genetic diversity and longitude in subtropical China. Although these previous studies do offer a fair axis of comparison with our results, plant communities in the two regions were subject to different historical and contemporary processes (e.g., glaciation patterns, orogeny, and climatic gradients; Qiu et al., 2011), and hence, are expected to show varying responses. Moreover, the inconsistencies could be attributed to the marked differences in spatial scale and the number of species between the two studies. In a similar scenario, a global-scale study recently found a relatively higher population genetic differentiation at lower latitudes (tropical/subtropical regions) as a result of restricted gene flow in these areas (Gamba and Muchhala, 2020). Though genetic patterns are often scale-dependent (Cushman and Landguth, 2010; Manel and Holderegger, 2013), thus rendering comparison of studies at different scales untenable, the patterns uncovered in the current study may offer useful insights into the partitioning of genetic diversity at larger spatial scales (Levy et al., 2014).

Genetic diversity had a weak but significant positive correlation with annual mean wet day frequency (wet) (Figure 5D). Such a link between water resources and plant genetic structure has been noted before (De Kort et al., 2021), though with precipitation being negatively related to plant genetic diversity (the correlation again weak but significant). This once again underscores the influence of spatial scale on the distribution of genetic diversity. Genetic diversity also had a weak but negative correlation with annual mean cloud cover percent (cld), thus reinforcing the link between genetic diversity and water. Although it is somewhat counterintuitive for wet and cld not to be correlated (they show discordant relationships with genetic diversity), the finding is not entirely unexpected, as cloud-precipitation anti-correlations have been noted particularly for low, optically thin cloud types (Jin et al., 2018). Nonetheless, our results support previous findings that climatic factors such as solar radiation, potential evaporation, and precipitation can affect the genetic structure of plant species (Bradbury et al., 2013; Kitamura et al., 2020). Although the mechanism by which such climatic factors influence genetic structure is as yet unclear, selective pressure for loci associated with photosynthesis and osmotic responses may have a role (Cruz et al., 2019). Our results hint at the importance of the interplay between energy- and water-related variables in shaping patterns of genetic diversity (e.g., Ramírez-Valiente et al., 2018; Li et al., 2019; De Kort et al., 2021), though this aspect may require a more thorough re-examination in the LRGR. As global warming is expected to cause large variations in energy and water-related factors, the evolutionary potential of most plant populations could be at risk owing to the uncertain effects of these climate changes on genetic variation. For this reason, more specific research inquiries into the role of climatic factors on plant genetic variation in the LRGR are urgently needed.

The LRGR is divided latitudinally into two climatic zones (the southern tropical monsoon and the northern subtropical monsoon climates), which might account for the relatively high proportion of genetic variation (14%) that was explained by the combination of climatic and topographic variables. Moreover, a far higher proportion (78%) of the observed genetic variation could be explained by the joint effect of climatic, topographic, and anthropogenic factors (Figure 6). This finding confirms the importance of these factors in determining the patterns of genetic structure in southwest China (e.g., Ju et al., 2018; Yang et al., 2019). Although our data show that topography is the main determinant of genetic structure in the LRGR, this relationship can be better understood in the context of climatic factors (e.g., Mosca et al., 2014; Li et al., 2019). The climatic dissimilarity and the ensuing IBE patterns in such a heterogeneous landscape (Bennie et al., 2008; Wang and Bradburd, 2014) might offer a plausible explanation for the considerably strong joint effect of climate and topography on the distribution of genetic diversity in the LRGR. The strong climate gradient in the LRGR might also explain the position of the break in genetic divergence (Figure 3B). This climatic dissimilarity, in the long run, may facilitate genetic divergence of populations along the established environmental gradients through the interaction of adaptive and non-adaptive mechanisms (Müller et al., 2017; Ramírez-Valiente et al., 2018; Galliart et al., 2019). Although longitude (long) and latitude (lat) were the best-fitting models in the topography category, other variables such as terrain roughness index (tri) and elevation (elev) might have a more practical and direct effect on genetic diversity patterns (Figure 5). We argue here that long and lat are simply proxies of these two variables (tri and elev).

Besides climate and topography, our analysis also found a weakly negative (but significant) relationship between the two anthropogenic variables and genetic diversity, suggesting the role of human activities in shaping the genetic structure of plant communities in the LRGR. However, this relationship appears to be driven partly by the huge sample size as well as other unknown outliers; thus, more anthropogenic variables may be required to clarify and validate the relationship. Though our data indicate a negative effect of human activities on the distribution of genetic diversity in the LRGR, some recent studies have alluded to the positive impact of ancient anthropogenic forces on biodiversity (e.g., Novák et al., 2019; Fajmonová et al., 2020), a trend that may also be true for genetic diversity, but more data are needed to confirm this relationship. A better understanding of the role of contemporary or recent human activities on the genetic structure of plant communities in the region would require characterization of habitats using more direct approaches. Nevertheless, the upsurge of human activities in the LRGR over recent decades is likely to exert an unprecedented pressure on natural ecosystems in the region (Liu et al., 2018), as demonstrated previously on natural reserves (including forest ecosystems) in Yunnan Province (Qiu et al., 2018).

In the wake of the current global environmental changes, it is important to preserve species adaptability and evolutionary potential, which will be achieved by considering not only species diversity but also intra-specific genetic diversity in conservation planning. Our analysis identified seven genetic diversity “hotspots” that showed an H_D value of >0.38 averaged across the study LRGR (Figure 4). As a significant proportion (~74.81%) of the identified “hotspots” land area was located outside of the protected areas, these “hotspots” should be prioritized in future conservation interventions, particularly in the context of climate change and increasing human activity. The larger “hotspots” (A, B, and D) were located in the northern part of the LRGR, where the relief is generally higher. The smaller “hotspots” (C, E, F, and G) geographically correspond to the isolated mountains in the southern part of the LRGR. However, more data from more species may be required to confirm the presence of the smaller hotspots since the southern part of the study area had relatively fewer representative species/populations. Conservation of the larger “hotspots” would require less effort since they are located closer to, and mostly overlap with, the major protected areas. However, conservation of the smaller “hotspots” may necessitate such interventions as the establishment of connections among the hotspots, as well as between the hotspots and the nearest protected areas, in order to limit gaps in the protected area network (Schoville et al., 2018). Moreover, establishment of protected areas close to cities, and corresponding to “hotspots” C, E, F, and G, would promote eco-tourism and mitigate the negative effects of urbanization on biodiversity (Conner, 2005; Puppim de Oliveira et al., 2011). To ensure sustainability and success in the conservation of these hotspot areas, it will be important for future investigations to consider population-level genetic dynamics (e.g., the effect of life history traits and the potential geophysical barriers on population-level genetic diversity) in the conservation decision-making process; conservation interventions should first recognize that the level of population genetic diversity (particularly in plants) is often a natural phenomenon that should not necessarily imply evolutionary capacity (see De Kort et al., 2021).