1. Introduction

Climate change caused by human activities poses a profound threat to global biodiversity (Bellard et al., 2012). Alterations in temperature and precipitation are expected to affect species distribution ranges (Pecl et al., 2017), modify regional species composition and structure, and lead to spatial and temporal variations in biodiversity patterns (Thuiller et al., 2011; Scheffers et al., 2016). These changes may compromise the effectiveness of conservation efforts in existing protected areas and exacerbate future conservation challenges (Araújo et al., 2011; Hoffmann et al., 2019). For example, the number of threatened endemic woody species in China is projected to increase by at least 50% compared to the current Red List, with over half of these newly identified species lacking effective protection (Peng et al., 2023). Further, by the 2070s, between 36.3% and 51.85% of Chinese Theaceae species are expected to face threats (Zhao et al., 2023). In southwest China, the overlap between existing nature reserves and potential future conservation hotspots for 2050 is only 13.71%, with most hotspots remaining inadequately protected (Wu et al., 2023). Consequently, predicting changes in species diversity and identifying potential conservation gaps within the current protected areas (PAs) are essential for guiding future conservation strategies and ensuring ecosystem integrity in a rapidly changing world.

Previous studies have shown that species across different regions and elevations respond variably to climate change. Generally, lowland species migrate upward, while tropical and boreal species tend to shift latitudinally towards temperate-dominated environments and polar communities in the future (Scheffers et al., 2016). For montane species, the risk of extinction may increase during upward migration due to the reduction or loss in available habitats (Dullinger et al., 2012). However, mountain range topographies vary across regions, displaying different patterns of change with increasing elevation, rather than following a simple monotonic decrease. Species partially distributed in lower montane zones may benefit from the increased surface area, potentially leading to range expansion as they move uphill (Elsen and Tingley, 2015). These differences in ecological conditions and geographic distributions may complicate species’ responses to climate change and present additional challenges for conservation planning. Accordingly, it is essential to identify the most vulnerable biological groups and to formulate targeted conservation strategies tailored to regional and habitat-specific contexts to avoid future biodiversity loss.

Currently, species distribution models (SDMs) are extensively applied to simulate potential distribution patterns and guide protected area planning (Liang et al., 2018; Li et al., 2021a; Peng et al., 2022). SDMs calculate ecological niches from species occurrence data and environmental variables using specific algorithms, projecting the results as probabilities across various spatial and temporal scales to predict species potential distribution (Peterson et al., 2011). Compared to single models, ensemble models (EMs) integrate predictions from multiple algorithms to circumvent the limitations of individual models and enhance predictive accuracy (Araújo and New, 2007; Thuiller et al., 2009). For species with limited occurrence records, ensembles of small models (ESMs) construct multiple bivariate models to increase the ratio of occurrence data to predictor variables, further improving predictive reliability (Breiner et al., 2015, 2018). Additionally, future projections from various global circulation models (GCMs) and climate change scenarios (i.e., shared socioeconomic pathways, SSPs) often show regional discrepancies (Buisson et al., 2010; Thuiller et al., 2019; Lawrence et al., 2021). Given the inherent uncertainties in species distribution modelling, combining ensembled SDMs (EMs and ESMs) with multiple GCMs and SSPs is critical for producing more robust predictions under future climate scenarios (Song et al., 2023).

Rhododendron is the largest genus of woody plants native to the Northern Hemisphere, comprising more than 1000 species worldwide (Chamberlain et al., 1996). China harbors the highest diversity of Rhododendron species, with over 600 species recognized to date (excluding intraspecies ranks), of which approximately 409 species are endemic (Fang et al., 2005). Alarmingly, 32% of these species are classified as endangered (Gibbs et al., 2011). The genus is primarily distributed in southern and southwestern China, with the Himalaya-Hengduan Mountains (HHM) serving as key center of its distribution and diversification (Shrestha et al., 2018). The subtropical evergreen broadleaf forests in southern China provide important habitats for Rhododendron. While in southwestern China, species occur not only in forests, but also above the treeline forming typical brush landscapes and alpine biomes (Zou et al., 2021). This adaptability of Rhododendron to diverse geographical regions and habitats underscores its ecological significance in alpine and subalpine ecosystems. Nevertheless, the extent to which these regional and habitat-specific variations shape its responses to anthropogenic climate change remains a critical knowledge gap that warrants further investigation.

In recent years, there has been a growing body of research examining the responses of Rhododendron species to climate change in China (Yu et al., 2017a, 2017b, 2019, 2021; Li et al., 2023). Notably, species with narrow geographic ranges and those inhabiting high elevations are particularly susceptible to the adverse impacts of climate change (Yu et al., 2019; Li et al., 2023). Wu et al. (2025) identified five threatened regions based on the changes in richness, endemism and phylogenetic diversity of 189 Rhododendron species. However, previous studies have seldom linked species-level changes with future alterations in community composition (i.e., temporal beta diversity), and identified priority protected areas have often overlooked potential refugia that could allow species to survive and adapt to the ecological consequences of climate change. Furthermore, research on the degree of habitat fragmentation under climate change as well as the relationship between range change and fragmentation remains limited. Therefore, the present study utilized detailed distribution data of 314 Rhododendron species (including 17 subspecies) and environmental variables to explore the spatial and temporal changes in their distribution and composition under climate change using ensembled SDMs. Specifically, we aimed to: (1) examine the impacts of climate change on species richness, weighted endemism, and temporal beta diversity of Rhododendron species; (2) identify priority conservation areas in the context of climate change, including climate refugia that preserve Rhododendron diversity (areas needing attention) and areas vulnerable to climate change (areas needing restoration and exploration); (3) compare the response of species from different regions and biomes in terms of similarities and differences in range change, elevation shift, latitude shift, threat levels and habitat fragmentation. By assessing the dynamic changes of Rhododendron diversity at both overall and local scales in relation to climate change, our goal is to enhance understanding of species’ ecological adaptation mechanisms in the context of global warming and provide recommendations for future conservation strategies.

2. Materials and methods 2.1. Distribution data

The distribution data of Rhododendron species analyzed in this study included over 70,000 records for 603 species, covering China and adjacent countries. These data were obtained from online databases, including the Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn), the Global Biodiversity Information Facility (GBIF, http://www.gbif.org), the Biodiversity of the Hengduan Mountains (BHM, http://hengduan.huh.harvard.edu/fieldnotes/), the herbarium of the Kunming Institute of Botany, Chinese Academy of Sciences (KUN, http://groups.kib.cas.cn/kun/), as well as from field surveys conducted by several research groups and collection teams of the Germplasm Bank of Wild Species (GBOWS) over the past two decades. We used the Flora of China (FOC, http://www.iplant.cn/foc) and Plants of the World Online (POWO, https://powo.science.kew.org) for species name standardization and correction. To ensure data validity and accuracy, incomplete geographic information, duplicates, cultivated plants, and specimen records from urban botanical gardens, parks, and museums were removed. Specimens with odd distribution points and conflicting identifications (e.g., identical collector and collection number but different species names) were re-identified. Only one record per grid cell was retained for each species to avoid uneven sampling. Species occurring in more than six distinct grid cells were selected to avoid poor model fit caused by inadequate occurrence points. For subspecies of certain species with clearly non-overlapping geographic distributions, we treated them as separate taxa due to their possible ecological niche differences. Different varieties of the same species were not distinguished and analyzed as one taxon. Ultimately, a total of 10,306 occurrence records for 314 Rhododendron species (including 17 subspecies) were retained for species distribution modelling analyses (Fig. S1) (See supporting information: Supplementary Methods for more details).

2.2. Environmental data

We employed three types of predictor variables (climatic, soil and topographic data) to model species distributions. The bioclimatic and elevation variables were downloaded from the Worldclim database (https://worldclim.org/, v.2.1, Fick and Hijmans, 2017) at 2.5 arc-minutes resolution (≈5 km). The data of soil physical and chemical properties (0–15 cm depth) were obtained from the SoilGrids database (https://maps.isric.org/). Aspect and slope were extracted from the elevation variable using the spatial analysis function in ArcGIS ver. 10.6 (ESRI, 2018).

Climate data utilized in this study comprised both current (1970–2000) and future (2061–2080 (the 2070s)) scenarios. Future climate data were derived from the SSPs outlined in the Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC). We calculated the mean values of seven GCMs (CNRMCM6-1, CNRM-ESM2-1, CanESM5, IPSL-CM6ALR, MIROC-ES2L, MRI-ESM2-0, MIROC6) under three climate scenarios: SSP126 (optimistic), SSP245 (intermediate) and SSP585 (pessimistic), which served as the basis for future predictions. SSP126 represented a sustainable pathway with low emissions (RCP2.6), achieved through global cooperation and green technology. SSP245 reflected an intermediate scenario (RCP4.5), with trends following historical patterns. SSP585 represented a fossil-fueled pathway with high emissions (RCP8.5), driven by rapid economic growth and limited climate mitigation efforts (O'Neill et al., 2017; Leão et al., 2021).

All variable layers were resampled to a 5 × 5 km resolution and projected to the WGS_1984_UTM_Zone_48N coordinate system (EPSG: 32648). To mitigate the risk of overfitting due to multicollinearity among variables, we calculated the variance inflation factor (VIF) and excluded variables with a VIF exceeding 10 using the package ‘usdm’ in R v.4.1.3 (Naimi and Araújo, 2016). Finally, 16 variables were retained for model development (Table S1).

2.3. Species distribution modelling

Before constructing models, we delineated the region for model formulation by creating a 200 km buffer around a minimum convex polygon encompassing all occurrence records (Mota et al., 2022). This buffer was included to account for potential species dispersal, thereby enhancing prediction accuracy and reducing overestimation of suitable habitat ranges (Anderson and Raza, 2010; Barve et al., 2011).

Two modelling approaches were applied based on the number of occurrence records: ESMs for species with 7–24 occurrences (194 species), and EMs for species with 25 or more occurrences (120 species) (Breiner et al., 2015). Both EMs and ESMs were constructed using five basic algorithms (i.e., classification tree analysis (CTA), generalized boosting model (GBM), generalized linear models (GLM), maximum entropy (MaxEnt), and random forest (RF)), implemented through the packages ‘biomod2’ (Thuiller et al., 2009) and ‘ecospat’ (Di Cola et al., 2017) in R v.4.0.5. For algorithms requiring absence data, 10,000 pseudo-absence points were generated by randomly sampling using the BIOMOD_FormatingData function in the package ‘biomod2’ (Phillips et al., 2009). To evaluate and calibrate the models, we conducted a split sample test, with 75% of the data for training and the remaining 25% for testing, and repeated this process 10 times. Each basic algorithm was run 10 times to generate 50 base models for EMs (5 algorithms × 10 runs) and 6000 bivariate models for ESMs (120 two-variable combinations from 16 environmental variables × 5 algorithms × 10 runs). Models were calibrated using true skill statistics (TSS) (Allouche et al., 2006), and those with TSS > 0.5 were selected for integration to improve the predictive performance (Peng et al., 2023). The performance of individual and aggregate models in both EM and ESM was assessed by comparing the TSS and area under the curve (AUC) values using the Kruskal–Wallis test (Kruskal and Wallis, 1952). Ultimately, we converted the species occurrence probabilities into binary maps using Max SSS (maximizing the sum of sensitivity and specificity) as a threshold (Liu et al., 2013). The default settings were applied for other parameters in both R packages.

2.4. Spatial pattern analyses

Species richness (SR) was defined as the total number of species present within each grid cell, while weighted endemism (WE) was determined as the sum of the inverse of the range size for each species (Williams et al., 1994; Crisp et al., 2001; Alves-Ferreira et al., 2024). The Sørensen dissimilarity index was utilized to assess compositional variations between grid cells under current and future scenarios (βt: temporal beta diversity), and further divided into replacement (β−3) and richness difference (βrich) (Schmera and Podani, 2011; Carvalho et al., 2012). The ratio of β−3/βt was also calculated to evaluate the relative contributions of species replacement versus richness difference under future scenarios. For all species, we calculated diversity patterns (SR, WE, and βt) using the R packages ‘phyloraster’ (Alves-Ferreira et al., 2024) and ‘divraster’ (Mota et al., 2023). Subsequently, standard deviational ellipses (SDE) and gravity center transfer were used to explore spatial processes and directional trends in species diversity over time (Lefever, 1926).

We identified areas vulnerable to climate change and areas capable of maintaining Rhododendron diversity under climate change as priority conservation areas (Li et al., 2024). The intersection of binary maps for individual species across current and future scenarios was extracted, which potentially served as long-term stable zones for species. Consistent grid cells ranking in the top 1% for both species richness and weighted endemism were identified as areas needing attention (ANA). Grid cells in the top 1% with the greatest loss (negative value) in species richness and weighted endemism were identified as areas needing restoration (ANR), while those with the highest gain (positive value) were designated as areas needing exploration (ANE) (Wu et al., 2025). Considering that priority conservation areas may differ across climate scenarios, we extracted the consistent priority areas across three scenarios and performed sensitivity analysis using threshold values of 5% and 10%. Additionally, we calculated the proportion of these priority conservation areas currently protected by PAs. All maps were generated using ArcGIS.

2.5. Statistical analyses

Variation partitioning analysis (VPA) and random forest methods were employed to identify the dominant factors influencing current potential distribution patterns, using the R packages ‘varpart’ (Dixon, 2003) and ‘randomForest’ (Breiman, 2001). The 314 modelled species were categorized into two groups based on their current distributions, following the classification of Chinese Floristic regions (Wu and Wu, 1996). The boundary, roughly delineated by the western edge of the Sichuan Basin and the eastern edge of the Yunnan Plateau, separated species into western and eastern regions (Fig. S1 and Table S2). The western region, located in southwestern China, primarily included the HHM, the Yunnan Plateau, and adjacent portions of the Qinghai-Tibet Plateau, comprising 215 species. The eastern region included hills and plains in the subtropical and tropical zones south of the Qinling-Huaihe line, encompassing 76 species. The other 23 species were not assigned to designated regions because their distributions either straddled the two regions or were located in the temperate zone of northeast China, far from both regions. In addition, species were assigned into three biome types according to their distributions relative to the treeline: alpine (above the treeline, 49 species), non-alpine (primarily below the treeline, 210 species), and both habitats (across the treeline in both forest and alpine habitats, 55 species) (Table S2). The assignment was based on habitat descriptions from the Flora of China and verified by specimen records.

To investigate species-specific differences in responses to climate change across regions and biomes, we selected the following indicators: (1) species range change (%); (2) change in the number of habitat patches (%); (3) change in the mean area of habitat patches (%); (4) elevation shift (m); (5) latitude shift (km) and (6) threat levels. Species range change was calculated as follows: ((future suitable habitat – current suitable habitat)/current suitable habitat) × 100%. Elevation shifts were quantified as the difference in mean elevation between future and current suitable habitats, while latitude shifts were computed based on the centroid displacement. The number of patches represented the total number of all individual patches within the suitable area of each species (where each patch was not connected to others in eight directions (including diagonals)), while the mean patch area was calculated as the average area of these separate patches (Hesselbarth et al., 2019; Cadena et al., 2023). These two metrics indicated the degree of habitat fragmentation. The change in the number of habitat patches was calculated as ((number of future habitat patches – number of current habitat patches)/number of current habitat patches) × 100%, and the change in the mean area of habitat patches was calculated as ((mean area of future habitat patches – mean area of current habitat patches)/mean area of current habitat patches) × 100%. Wilcoxon test was used to compare differences between species from different regions and biomes. To mitigate the impact of extreme values on significance assessments, we excluded species projected to lack suitable habitats in the future from the analysis of elevation and latitude shift. Furthermore, we also used standardized major axis regression to evaluate differences in the relationship between species range change and patch changes across regions. Species were classified into threat levels based on range change, following the IUCN Red List Categories and Criteria (A3c Indicators) (IUCN, 2012), which included Extinct (EX), Critically Endangered (CR), Endangered (EN), Vulnerable (VU), and Least Concern (LC). All statistical analyses were conducted using the packages ‘rstatix’ (Kassambara, 2023), ‘landscapemetrics’ (Hesselbarth et al., 2019) and ‘smatr’ (Warton et al., 2012) in R v.4.1.3.

3. Results 3.1. Model performance and key environmental variables

The performance evaluation of the two models (ESMs and EMs) showed good fits, as indicated by the high values of AUC and TSS (Fig. S2). For the EMs, the average AUC was 0.997 (range: 0.982–1), and the average TSS was 0.977 (range: 0.864–1). ESMs showed slightly lower average values, with an average AUC of 0.987 (range: 0.93–1), and an average TSS of 0.967 (range: 0.791–1). Our analyses indicated that EMs significantly outperformed the five basic algorithms across both metrics. While ESMs did not show a significant overall performance advantage, they, along with GLM, outperformed other baseline algorithms.

The variation partitioning results showed that the three predictor groups (climatic, soil and topographic variables) collectively explained 33.3% of the variation in Rhododendron distribution (Fig. 1a). Among these, climatic variables independently explained the largest portion of variation (10.8%). Specifically, temperature seasonality (BIO4) and precipitation seasonality (BIO15) were identified as the top important environmental variables affecting species richness and endemism, respectively (Fig. 1b).

3.2. Changes in patterns of species diversity

The current potential distribution patterns indicated that Rhododendron species were primarily distributed in the mountainous regions of southern and southwestern China, with the highest diversity observed in the Hengduan Mountains (Fig. 2a and b). Forward projections for the 2070s suggested that while overall distribution patterns remained similar to those observed currently, the core distribution area contracted and became more concentrated, decreasing by 10.85%, 10.33%, and 16% under the SSP126, SSP245, and SSP585 scenarios, respectively. Besides, under these three scenarios, the average species richness per grid cell was expected to decline by 0.484, 0.794 and 1.211, while the average weighted endemism remained relatively stable (Table S3). Concretely, areas with increased species richness were anticipated in the south and east of the Qinghai-Tibet Plateau and the northern Hengduan Mountains, with additional increases in parts of eastern and central China (such as Hunan, Jiangxi, Fujian and Zhejiang provinces). Conversely, declines in species richness were projected in southwestern China (such as East Himalaya, southern Hengduan Mountains and Yunnan Plateau) as well as southern China (notably Guangxi and Guangdong provinces) (Fig. 2c–e and g). Significant alterations in weighted endemism were mostly expected in the mountainous regions including the HHM, the Nanling Mountains in southern China and the Changbai Mountains in northeast China (Fig. 2d–f and h).

Climate change was also expected to alter the composition of Rhododendron species, with these changes intensifying under more extreme scenarios (Fig. 3). The average beta diversity per grid cell was projected to increase from 0.48 under SSP126 to 0.502 under SSP245, and to 0.568 under SSP585. The most pronounced changes in beta diversity were anticipated mainly in areas characterized by low species richness, especially at the edges of the distribution range, where richness differences were most prominent. These areas included northeastern China and parts of southwestern China (such as southeastern Qinghai, central and southern Xizang and southeastern Yunnan). In contrast, areas with higher species richness, where species remained stable (e.g., the southern Hengduan Mountains), predominantly exhibited the replacement component of beta diversity. Despite fluctuations in species composition due to migration, a higher proportion of similar components in these regions resulted in lower total beta diversity.

3.3. Priority protected areas and conservation gaps

The priority protected areas identified here remained substantially consistent across the three thresholds (1%, 5%, and 10%), with their spatial extent expanding at higher thresholds (Fig. 4). Given our focus on regions experiencing the most pronounced changes or maintaining higher diversity under climate change, we emphasized the results based on the 1% threshold.

Under the richness indicator, the northern Hengduan Mountains was classified as ANE, while the ANR was primarily located in the Himalaya and the southern Hengduan Mountains (Fig. 4a). Compared to the richness indicator, the ANE and ANR derived from the endemism indicator were co-distributed (Fig. 4b). The identified priority protected areas encompassed several other mountainous regions, including the Ailao Mountain, Wumeng Mountain – Daliang Mountain, Nanling Mountains, Taiwan Central Mountains, and Changbai Mountains. Additionally, ANA (Fig. S3), primarily concentrated in the Hengduan Mountains, may serve as a crucial refuge against climate change.

When the threshold was expanded to 5%, in addition to the previously identified regions, attention should also be directed to the Taibai Mountains in Shaanxi and the eastern part of Liaoning, which experienced substantial declines in endemism (Fig. 4d). At the 10% threshold, areas in eastern Chongqing and northern Guizhou, which showed significant declines in endemism, were identified as ANR (Fig. 4f). While eastern Zhejiang, showing a notable increases in richness, was classified as ANE (Fig. 4e).

The conservation gap analyses indicated that the existing network of PAs would be inadequate to effectively cover these priority protected areas under future climate change scenarios. The protection ratio declined significantly with increasing thresholds (Table S4).

3.4. Differences and consistencies in response to climate change between species in western and eastern regions

Shifts in regional distribution patterns, key environmental variables, and changes in species distribution collectively revealed consistencies and differences in species responses across regions. Significant disparities in species richness were observed between the western and eastern regions. The western region harbored more than twice the number of Rhododendron species compared to the eastern region, with a higher proportion of narrowly distributed species (Figs. 5d and S4a). In terms of environmental variables, BIO4 emerged as the most influential factor, accounting for 56.69% of all species and 65.58% of species in the western region (Table 1). Conversely, in the eastern region, BIO14, BIO4 and BIO15 were the top three variables, with proportions of 40.79%, 27.63%, and 11.84%, respectively.

From the regional distribution patterns, the main distribution area showed a trend of contraction in the western region, whereas it exhibited an expansion in the eastern region (Fig. S4 and Table S5). Furthermore, the distribution centers in both regions showed a northwestward shift, with migration distances in the western region being greater than those in the eastern region (69.06 km > 18.26 km, 92.43 km > 37.09 km, 128.9 km > 51.43 km) under the three SSP scenarios.

Species distribution changes indicated that species in the western region exhibited significantly greater elevation migration compared to those in the eastern region, while the latter showed significantly larger latitudinal migration distances (Fig. S5a and b). Both regions displayed similar trends in species range loss and habitat fragmentation in the future (Fig. S5c–e and Table S6), with no significant differences detected between them. However, changes in the number of habitat patches and mean patch area were significantly and positively correlated with changes in species distribution range (Fig. 5a and b). In the eastern region, the reduction in species distribution ranges was accompanied by a more pronounced decline in average patch area, while no regional differences were observed in the linear regression of changes in the number of habitat patches and species distribution ranges. Further analysis revealed that the contraction ratio (the ratio of current habitat loss in the future) significantly decreased as the current range size increased (Fig. 5c). Below the critical thresholds—where species from both regions exhibited identical contraction ratios—species in the eastern region showed higher contraction ratios. Conversely, above the thresholds, species in the western region showed higher contraction ratios.

The assessment of species threat status revealed that the western region had a higher overall proportion of threatened species than the eastern region under the three SSP scenarios (42.79% > 32.89%, 51.16% > 35.53%, 52.56% > 47.37%, respectively) (Fig. S5f).

3.5. Different responses of species from different biomes to climate change

The species-level distribution changes under all three SSP scenarios revealed that alpine species experienced significantly higher habitat loss and had the greatest proportions of threatened species compared to non-alpine species and those distributed across both alpine and non-alpine biomes (both habitats species) (Fig. S6c and f). Regarding the mean area of habitat patches changes, the alpine species showed significantly lower values than both habitats species. Additionally, the latitudinal migration distance for both habitats species was significantly shorter than that of non-alpine species. However, no significant differences were observed among the three categories concerning migration distance along elevation gradients or changes in the number of habitat patches. Notably, under the SSP585 scenario, with the most severe warming, alpine species showed a significantly greater reduction in average habitat patch area compared to the other two categories (Fig. S6e).

4. Discussion 4.1. Key environmental variables affecting species distribution of Rhododendron

At the macroscale, climate change primarily influences plant species diversity patterns (Punyasena et al., 2008), with hydrothermal conditions playing a critical role in shaping geographical distributions and regulating physiological processes (Woodward, 1987). Our findings, based on over half of the Rhododendron species occurring in China, confirmed that climate variables were the predominant factors influencing overall diversity patterns and the distribution of individual species (Fig. 1 and Table 1). Further analysis revealed that the observed differences in species richness between northern and southern China may be associated with temperature seasonality (BIO4). Areas with higher richness in the south were characterized by lower temperature seasonality, indicative of greater climate stability. According to the climate stability hypothesis, regions with stable climatic conditions can accommodate a broader ecological niche space, promoting species evolution and adaptation, and ultimately fostering higher species richness (Klopfer, 1959; Klopfer and Macarthur, 1961; Luo et al., 2012).

Our study also found that distinct environmental factors influenced the distribution of Rhododendron species in the western and eastern regions (Table 1). Alpine species located in the western region exhibited greater sensitivity to temperature fluctuations, corroborating previous studies on mountain flora (Dakhil et al., 2021; Liao et al., 2021; Li et al., 2023). In contrast, the distribution of Rhododendron in the eastern region was significantly affected by the precipitation of the driest month (BIO14), which reflected the species’ water requirements and tolerance to arid environments. Global warming is likely to exacerbate aridity in subtropical regions, where hygrophilous plants tend to be more responsive to changes in rainfall than temperature, particularly due to the limited precipitation during the early growth season (Li et al., 2021b).

Apart from climatic variables, soil also played a crucial role in the distribution of certain species in southwestern China (Table 1). Higher organic carbon content and suitable soil texture enhance water and nutrient absorption, thereby supporting plant growth (Martre et al., 2002; Raynaud and Leadley, 2004; Bornman et al., 2008). The complexity of soil conditions not only influences the distribution of individual species, but also contributes to the formation and maintenance of regional plant diversity (Fig. 1; Cai, 2022).

4.2. Diversity pattern change of Rhododendron species under climate change

Our projections indicated that climate change significantly affected the patterns of species diversity, with these impacts escalating as climate change progressed (Fig. 2, Fig. 3). While species richness changes were expected to occur across the entire distribution range, shifts in endemism were likely to be more spatially localized, particularly in mountainous regions. Mountains serve as natural barriers that isolate species populations and restrict their migration, resulting in more fragmented distributions and potential range contractions under climate change (La Sorte and Jetz, 2010). For Rhododendron, such contractions lead to a more concentrated distribution. The primary drivers of range contraction are the shrinking “warm edges” of most species’ ranges (manifested as a northward shift of the southern boundary and a westward shift of the eastern boundary for most species), coupled with the limited expansion of “cold edges” (which includes a northward shift of the northern boundary and a westward shift of the western boundary for a few species) (Fig. S7; Sun et al., 2020). The shifting boundaries imply that marginal populations are under greater environmental stress from climate change compared to those in core ranges, facing challenges such as disrupted biological interactions and reduced genetic diversity (Vilà-Cabrera et al., 2019). These factors contribute to higher mortality rates in marginal populations (Geppert et al., 2020).

Moreover, future variation in Rhododendron species composition would be more pronounced, especially in areas experiencing greater rates of species richness loss or acquisition. In areas currently characterized by low species richness, even minor fluctuations in species numbers can result in substantial changes in acquisition and loss ratios (Duan et al., 2016). In extreme cases, range reductions may lead to future communities being subsets of current ones, while species migrating into new habitats may integrate existing communities into future ones, thereby amplifying differences in richness gradients and increasing beta diversity (De Lima et al., 2019; Mota et al., 2022).

4.3. Spatial and biome heterogeneity of Rhododendron species in response to climate change

Under future climate scenarios, many Rhododendron species would experience habitat loss and fragmentation, migrating to higher elevations or latitudes (Fig. S5), which is consistent with previous studies (Yu et al., 2021; Li et al., 2023). Specifically, species in the western region were anticipated to move upward along elevational gradients, while species in the eastern region were more likely to shift towards higher latitudes. Consequently, a distinct richness pattern may emerge in southwestern China, characterized by increases at higher elevations and declines at lower elevations. In the eastern region (i.e., southern China), notable latitudinal changes may be observed. These patterns are attributable to the distinctive topographical features of China, characterized by higher elevations in the west and lower elevations in the east, resulting in a stepped landscape that declines significantly from west to east. Numerous studies across various species have corroborated these trends (He et al., 2019; Gan et al., 2024). Meanwhile, the ability of species to migrate to new habitats is influenced by factors such as dispersal capability, seedling establishment and habitat constraints (Lawlor et al., 2024). For instance, small, lightweight Rhododendron seeds are capable of long-distance dispersal facilitated by wind and animals (Stephenson et al., 2007; Wang et al., 2014). However, the complex topography and geographic barriers in the western region may prevent them from reaching climatically suitable zones, resulting in a contraction of the core distribution range (Freeman et al., 2018).

Spatial variation in species responses to climate change is also associated with their current area and elevational range (Yu et al., 2017a). In the western region, relatively narrow-range species with wider elevational ranges may exhibit greater resilience to climate variability (e.g., Rhododendron arizelum and R. setosum), leading to less contraction of suitable habitat compared to species in the eastern region (Valladares et al., 2014; Yu et al., 2019). Broad-range species, with expansive ecological niches and greater dispersal capacities, typically experience less habitat loss than narrow-range species (Williams et al., 2006; Lawlor et al., 2024). However, rapid elevational changes may hinder large-scale migration and reduce habitat connectivity, potentially increasing habitat fragmentation risks and resulting in higher contraction ratio for broad-range mountain species in the west (Trew and Maclean, 2021). Not surprisingly, species with narrow elevational ranges and restricted geographic distributions, such as R. aureum, R. pronum, and R. ziyuanense, which are confined to specific habitats, warrant particular conservation efforts (Yu et al., 2019).

Alpine Rhododendron species are particularly vulnerable to climate change, as indicated by significant declines in suitable area and mean habitat patch area compared to the other two categories (Figs. S6c, e and f). These species inhabit higher elevations and are at risk of “mountain-top” extinction (La Sorte and Jetz, 2010; Dullinger et al., 2012). On the other hand, species that inhabit both the alpine and non-alpine biomes demonstrate a high capacity to adapt to diverse habitats (Geppert et al., 2020). However, the limited space available in high elevation environments during upward migration makes these species moderately threatened.

4.4. Mountains as priority conservation areas for Rhododendron species

In the context of climate change, our research on Rhododendron species indicated that mountainous regions demonstrated contrasting characteristics: they served as critical refugia, providing stable habitats that buffer against climate variability (ANA), while concurrently being highly sensitive to climate change (ANE and ANR). The Hengduan Mountains exemplify this phenomenon. The complex landscape heterogeneity and micro-climate conditions provide diverse habitats and micro refugia for various species (Sun et al., 2020; Trew and Maclean, 2021; Guan et al., 2024). However, the accelerated warming rates in high-mountain regions (Pepin et al., 2015), along with the limited habitat tolerances of certain mountain species (Thuiller et al., 2005), render these areas vulnerable. Although significant protections are afforded by national nature reserves, such as those in the Gaoligong Mountain and the Baima Snow Mountain, connectivity between the reserves is inadequate. Therefore, key conservation strategies should prioritize enhancing the connectivity of existing nature reserves and protecting marginal populations, potentially through the establishment of temporary nature reserves (Lawlor et al., 2024).

Beyond the Hengduan Mountains, the ANE and ANR encompass additional mountainous areas experiencing significant changes in species diversity, including the Himalaya, Ailao Mountain, Wumeng Mountain–Daliang Mountain, Nanling Mountains, Taiwan Central Mountains and Changbai Mountains (Fig. 4). These dispersed and fragmented high-risk areas pose challenges to the existing conservation network. Effective strategies, such as the establishment of ecological corridors, are urgently needed in these regions to mitigate the impacts of landscape fragmentation on gene flow and to expand population ranges (Luo et al., 2024). Our research, using Rhododendron species as a case study, underscores the need to incorporate mountainous regions into priority protected areas to support the achievement of the ‘30 × 30’ target outlined in the Kunming-Montreal Global Biodiversity Framework.

4.5. Limitations of the study

This study primarily focused on the effects of environmental variables on the distribution of Rhododendron species. However, other factors, including resource competition, predation, anthropogenic activities, land use and biotic–abiotic interactions, also influence species distribution (Hampe and Petit, 2005; Dullinger et al., 2012). Additionally, we employed a 5-km resolution for the species distribution models, which may be inadequate for capturing the highly heterogeneous habitats of Rhododendron species, particularly in mountain regions like the HHM. Furthermore, different global circulation models contribute to model uncertainty, and relying solely on their mean without accounting for the full range of GCMs outputs may reduce prediction reliability. These uncertainties probably introduce bias into the projections (Pearson and Dawson, 2003). Despite these uncertainties, our findings on the dynamic changes in Rhododendron diversity patterns under climate change provide an important foundation for developing effective species- or genus-specific conservation strategies, underscoring the magnitude of the challenge that climate change presents to conservation efforts.

5. Conclusion

This study employed ensemble models (EMs) and ensembles of small models (ESMs) to predict the potential distribution of 314 Rhododendron species in China and evaluate how climate change-driven species redistribution may alter diversity patterns. Our findings revealed spatiotemporal variations in species’ responses to climate change across different regions: in the western mountainous regions, species exhibited upward migration and a contraction of their core distribution areas. While in the eastern hills and plains, species shifted northward, accompanied by an expansion of their core distribution areas. In addition, the projected future distribution changes identified priority conservation areas, primarily in mountainous regions, which may potentially experience significant changes or consistently maintain high species diversity. However, existing static nature reserves are insufficient to accommodate species with shifting distributions, and conservation agencies lack the necessary resources to protect newly suitable habitats. This underscores the need for more effective conservation strategies and substantial financial investments to mitigate the potential consequences of climate change on biodiversity.

Acknowledgements

We thank Professors Renlin Liu, Yongpeng Ma, and Fei Wang for providing the field collection data of Rhododendron species, Mr. Guangfu Zhu and Dr. Dan Xie for their suggestions on the modelling selection and data analysis, and Dr. Hantao Qin and Ms. Yajie Cui for their valuable comments on the manuscript revision. We also gratefully acknowledge the data support provided by online databases, including National Plant Specimen Resource Center, GBIF, BHM, KUN, E and GBOWS. This study was supported by the National Key Research and Development Program of China (2023YFF0805800), the Science and Technology Basic Resources Investigation Program (2021FY100200), the Key Basic Research program of Yunnan Province, China (202101BC070003), and the Chinese Academy of Sciences President's International Fellowship Initiative (2024PVA0087). The Royal Botanic Garden Edinburgh is supported by the Rural and Environment Science and Analytical Services Division (RESAS) of the Scottish Government.

CRediT authorship contribution statement

Ming-Shu Zhu: Writing – original draft, Software, Methodology, Formal analysis, Data curation. Zhi-Qiong Mo: Software, Methodology. Michael Möller: Writing – review & editing. Ting Zhang: Resources, Data curation. Chao-Nan Fu: Software, Methodology. Jie Cai: Resources, Data curation. Wei Zheng: Resources, Data curation. Ya-Huang Luo: Software, Methodology. De-Zhu Li: Conceptualization. Lian-Ming Gao: Writing – review & editing, Funding acquisition, Conceptualization.

Data availability statement

The species occurrence records used in this study will be openly accessible on ScienceDB upon publication (https://doi.org/10.57760/sciencedb.17856). Environmental data were sourced from online databases, with details provided in the methods. Full data on indicators of changes in the distribution of Rhododendron species in China under climate change are available at https://doi.org/10.57760/sciencedb.17893.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

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