1. Introduction

Mountains are defined by their prominence, steep sides, and significant height above the surrounding regions (Fang et al., 2004). They harbor exceptionally high levels of biodiversity and a great number of endemic species worldwide, which are also extremely sensitive to the ongoing global climate change (Zu and Wang, 2022). However, the underlying mechanisms that generate and maintain mountain biodiversity are not well understood (Körner et al., 2016; Rahbek et al., 2019).

Numerous hypotheses have been proposed to explain the large-scale distribution of species diversity (Wang et al., 2009). The contemporary climate hypothesis – e.g., environmental energy (Currie, 1991), water-energy dynamics (O'Brien, 1998; Hawkins et al., 2003) and climate seasonality (Connell and Orias, 1964; Stevens, 1989) – propose that climate (i.e., temperature, moisture, and seasonality) limit the distribution of species by affecting physiological processes. The historical climate change hypothesis posits that a greater number of species survive during periods with stable climates than during period with unstable climates, e.g., the Quaternary period (Liu et al., 2017; Zou et al., 2019). In contrast, the habitat heterogeneity hypothesis states that heterogeneous environments promote species diversity by providing more niches, geographical isolation, and climate refuges (Kerr and Packer, 1997; Stein et al., 2014). Despite extensive theoretical work, it remains unclear whether these hypotheses explain observed diversity in mountain plant species.

An additional theory that may help understand mountain species biodiversity is island biogeography theory (IBT). This theory explains that species diversity on islands is mediated by variation in island area and distance to the mainland (i.e., isolation degree) (MacArthur and Wilson, 1967). Specifically, island area is inversely correlated with extinction rate, whereas isolation degree is inversely correlated with immigration rate (Wu, 1990; Han, 1994). In addition, studies have confirmed that isolation degree is negatively correlated with species similarity (Sklenář et al., 2014; Ménde-Castro et al., 2018). Most studies of IBT focus on ‘true islands’, i.e., ocean islands (Han, 1994; Costanzi and Steifetten, 2019), which can be described as suitable habitats for populations from the mainland, surrounded by unsuitable matrix habitats (e.g., ocean) (Haila, 2002). However, IBT has been extended to various habitat islands on land that may not have a clear boundary with matrix habitats (Han, 1994; Kreft et al., 2008; Matthews, 2021), such as forests (Lovei et al., 2006), grasslands (Lindgren and Cousins, 2017), and mountains (Sklenář et al., 2014). Mountains are surrounded and separated by lowland, just as ‘true islands’ are surrounded by water, and hence they have been referred to as ‘islands on land’ or ‘sky islands’ (Howell, 1947; Heald, 1967; Carlquist, 1974; Han, 1994; Ding and Zheng, 1996; Fang et al., 2004; Quammen, 2004; Sklenář et al., 2014). However, few studies have determined whether the same factors drive and maintain biodiversity on mountains and ‘true islands’ (Costanzi and Steifetten, 2019; Mendez-Castro et al., 2021).

Rhododendrons are an excellent model for examining the underlying factors that regulate mountain biodiversity (Shrestha et al., 2018b). Rhododendron is a large genus of more than 1000 species that occur in every possible habitat in Asia, Europe, North America and Oceania (Fang and Ming, 1995; Gibbs et al., 2011). Rhododendron is also one of the most important genera of mountainous regions (Kumar, 2012; Yu et al., 2017). It is the dominant species in the understory of coniferous and broad-leaved mixed forests in subalpine coniferous forests, and plays a very important role in maintaining the stability of the structure and function of mountain ecosystems (Gibbs et al., 2011; Kumar, 2012). Previous studies have tried to explain the large-scale patterns of Rhododendron diversity with various environmental hypotheses (Shrestha et al., 2018b; Xia et al., 2022), these studies relied on evolutionary and genetic explanations of patterns of Rhododendron diversity. However, we try to explore the explanatory power of IBT in explaining Rhododendron diversity pattern, and the applicability and extensibility of IBT when regarding mountains as islands.

In this study, we determined which of six hypotheses best explains the formation and maintenance of Rhododendron diversity. These hypotheses included contemporary climate parameters (environmental energy, water availability, climate seasonality), Quaternary climate change, habitat heterogeneity, and IBT. After our analysis indicated that Rhododendron diversity can be largely explained by the IBT, we examined whether the main tenants of this theory are all valid for ‘mountain islands’, including the relationships between diversity and island area, diversity and distance from the mainland, as well the number of shared species and distance from the ‘mainland’.

2. Materials and methods 2.1. Species data

The global distribution records of Rhododendron species were collected from the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/). We identified and deleted records with no coordinates, outlier coordinates, zero coordinates, identical latitudes/longitudes, and invalid coordinates based on the R package ‘CoordinateCleaner’ (Zizka et al., 2019). We supplemented these data with distribution records from herbaria and botanical museums at Kunming Institute of Botany (KUN), Chinese Academy of Sciences (CAS); South China Botanical Garden (SCBG), CAS; Wuhan Botanical Garden (WBG), CAS; Sichuan University of Botany (SUB); Sichuan Forest School (SFS); and Lushan Botanical Garden (LBG). The scientific names of Rhododendron were checked by the R package ‘plantlist’ (Zhang et al., 2022), and only the correct records were kept. Ultimately, a total of 262,496 occurrence records, including 898 species of Rhododendron were retained.

To eliminate the influence of area on diversity estimation and statistical analysis, the distribution data were transferred into 1° × 1° grid cells. Grids located at the land boundary or coastal regions with a land area of less than 50% of the grid were removed. When calculating species diversity, the number of Rhododendron species in each grid was taken as the species diversity of the grid.

2.2. Predictors used to test environmental hypotheses

To explore the role of different processes in explaining species diversity, we calculated 17 predictors that represented measurable properties of the hypotheses associated with contemporary climate, Quaternary climate change, habitat heterogeneity and island biogeography theory (see Table 1 for details). The contemporary climate group was divided into three categories: environmental energy, water availability, and climate seasonality. All climate data were obtained from WorldClim (https://www.worldclim.org/) and the CGIAR Consortium for Spatial Information (https://cgiarcsi.community/2019/01/24/global-aridity-index-and-potential-evapotranspiration-climate-database-v3/) with a spatial resolution of 30 arc seconds. The climate data from the Last Glacial Maximum (LGM) of the Quaternary were calculated from the average of three commonly used global climate models (GCMs): CCSM4 (Gent and Danabasoglu, 2011), MIROC-ESM (Watanabe et al., 2011), and MPI-ESM-P (Giorgetta et al., 2013) at a spatial resolution of 2.5 arc minutes. Elevation data were downloaded from the GMTED2010 digital elevation model (https://www.usgs.gov/centers/eros/science/usgs-eros-archive-digital-elevation-global-multi-resolution-terrain-elevation) at a spatial resolution of 30 arc seconds. Global mountain polygon data (the ‘Basic’ map unit in the version of Global Mountain Biodiversity Assessment (GMBA) Mountain Inventory v.2_standard) were downloaded on EarthEnv (https://www.earthenv.org/mountains). All environmental predictors were aggregated to the 1° × 1° grid. The mountain area of each grid was estimated as the sum of planimetric area of mountain polygons on each grid, and the mountain area values less than half of the grid cell area were not included in the statistics to reduce error in the results. The distance to the mainland was calculated as the centroid distance between each grid and the ‘mainland’ of Rhododendrons (see section 2.3 for more information about ‘mainland’).

2.3. Testing the application of the island biogeography theory in mountain islands

In classic IBT, the mainland plays important roles as it is the species pool or colonist source of islands, and promotes potential immigration events (MacArthur and Wilson, 1967; Han, 1994). However, it is much more complicated to identify a ‘mainland’ in terrestrial habitat islands than it is for ‘true islands’. Studies have tried to identify the ‘mainland’ as regions with the longest species history (Allen et al., 2009) or with the most numerous species (Sanchez and Parmenter, 2002) or the largest and most continuous habitats in the study areas (Costanzi and Steifetten, 2019; Yan et al., 2023). The Hengduan Mountains to Himalayas region in southwestern China harbor a large and continuous distribution of mountains suitable for Rhododendrons, and it has been recognized as the diversity center of world's Rhododendrons (Yu et al., 2017) for most of extant Rhododendron species originated from evolutionary radiations occurred here in history (Shrestha et al., 2018b; Xia et al., 2022). These conditions make this region a good choice for being the ‘mainland’ of Rhododendrons.

To further test the ‘island area-diversity’ relationship, ‘distance-diversity’ relationship and ‘distance-shared species’ relationship of IBT in mountains, we regarded each mountain as an ‘island’ and the Hengduan Mountains region as the ‘mainland’ (Fig. 1). A total of 6717 mountains were initially included in global mountain polygon data downloaded from GMBA Mountain inventory v.2.

We calculated the planimetric area of each mountain as the ‘area of islands’ (i.e., the projected area of the mountains on the horizontal plane). Because mountains are three-dimensional, we also calculated the rough volume of each mountain, which was calculated as: the planimetric area × the highest elevation of the mountain. We then used both the planimetric area and the volume as indicators of the area of ‘mountain islands’, and counted the number of Rhododendron species in each ‘mountain island’. To explore the ‘distance-shared species’ relationship of IBT in ‘mountain islands’, we calculated the centroid distance between each ‘mountain island’ and the ‘mainland’, and then counted the number of same species in each pair of ‘mountain island’ and ‘mainland’.

2.4. Statistical analysis

Because species diversity per grid cell normally fits a Poisson distribution (Wang et al., 2011), we first established generalized linear models (GLMs) with Poisson errors to estimate the relationship between species diversity and 17 explanatory variables. For each predictor, we established a GLM and the explanatory power of the predictor was estimated by the adjusted R², which was calculated as follows: R_adj²(%) = 100–100*(residual deviance/null deviance). To control the influence of spatial autocorrelation in the associations between species diversity and variables (Kissling and Carl, 2007), we further applied simultaneous autoregressive (SAR) models and used the Akaike information criterion (AIC) as an evaluation.

To explore the relative importance of environmental variable groups represented by six hypotheses in explaining Rhododendron diversity patterns, we performed a hierarchical partitioning analysis that estimated individual contributions of the six groups, and further partitioned their unique contribution, sharing the contribution by variance partitioning analysis. Since each environmental group contained highly correlated variables, this method effectively addresses the collinearity problems between variables (Lai et al., 2022). To maintain the same number of variables for each group, we selected the most explanatory predictor (bold predictors in Table 2) in each hypothesis as their representative.

Spearman correlation analyses were conducted to assess the relationships between species diversity and the planimetric area and the volume of ‘mountain islands’, species diversity and the number of shared species and the distance from ‘mountain islands’ to the ‘mainland’. We also performed a test of significance by ‘cocor’ package (Diedenhofen and Musch, 2015) for the difference in the correlation coefficients between the species diversity and the planimetric area and the volume of ‘mountain islands’.

All statistical analyses were performed in R 4.1.2 (R Core Team, 2021). Spatial autoregressions were conducted using the R package ‘spatialreg’ (Bivand et al., 2021). The hierarchical partition and variance partition analysis was conducted using the ‘rdacca.hp’ package (Lai et al., 2022).

3. Results 3.1. Geographical patterns of Rhododendron diversity

Rhododendron hotspots (regions with over 32 species) are mainly distributed in the Hengduan Mountains of southwestern China and the southern regions of the Himalayas (Fig. 2). Regions with species numbers ranging from 5 to 31 are scattered throughout Southeast Asia, including southeastern China, the Japanese archipelago, and the Pacific islands between Asia and Oceania. In contrast, Rhododendron diversity is low (regions with species less than 5) in Northeast Asia, South Asia, Europe, and North America.

3.2. The relative explanatory power of six environmental hypotheses on Rhododendron diversity patterns

The most influential predictor of global patterns of Rhododendron diversity was MA (R_GLM² = 35.07%, P < 0.001; AIC_SAR = 19,901, P < 0.001), a proxy for the IBT (Table 2 and Fig. 3). Specifically, MA was positively correlated with Rhododendron diversity. The second-best predictors of Rhododendron diversity patterns were proxies for habitat heterogeneity, i.e., ELER (R_GLM² = 29.69%, P < 0.001; AIC_SAR = 19,834, P < 0.001) and MATR (R_GLM² = 27.91%, P < 0.001; AIC_SAR = 19,853, P < 0.001). These proxies were also positively correlated with species diversity (Table 2 and Fig. 3). Proxies of the paleo-climate change hypothesis, MAT_vel, explained 21.08% of the total variation in Rhododendron diversity (P < 0.001). Variables with relatively low explanatory power included MD, TSN, ART, MTCQ, MAT, and MAP, with contributions ranging from 12.52% to 17.57% (P < 0.001). The contributions of the remaining predictors were each lower than 10%. The results from GLMs and corresponding SAR models were similar for all predictors.

Hierarchical partitioning analysis also revealed that IBT explains the largest variance (10.17%) in Rhododendron diversity patterns (Fig. 4). The individual contribution of habitat heterogeneity explained 7.31% of the total variance. The individual effects of several hypotheses explained similar, low levels of variance, i.e., water availability (4.72%), environmental energy (4.49%), paleo-climate change (3.68%), and climate seasonality (3.35%). All hypotheses combined explained 33.7% of the variance of Rhododendron diversity. The decomposed variance of individual effects in the fractions showed the shared effect between IBT and habitat heterogeneity was the highest, accounting for 6.66%.

3.3. The application of island biogeography theory in explaining mountainous Rhododendron diversity

Rhododendron diversity was positively correlated with the planimetric area of mountains (R = 0.39, P < 0.001) (Fig. 5a). This correlation was significantly stronger between species diversity and volume of mountains (R = 0.41, P < 0.001) (Fig. 5b). In addition, Rhododendron diversity of each ‘mountain island’ was negatively correlated with isolation degree, i.e., distance from the ‘mainland’ (R = −0.28, P < 0.001) (Fig. 5c). The number of Rhododendron species shared between ‘mountain islands’ and the ‘mainland’ were also negatively correlated (R = −0.52, P < 0.001) (Fig. 5d).

4. Discussion 4.1. Patterns of Rhododendron diversity can be explained by island biogeography theory and habitat heterogeneity

Our results demonstrate that Rhododendron diversity patterns are best explained by the IBT, then by habitat heterogeneity. Regions that contain large mountain areas tend to have higher species diversity. For example, the Himalayas-Hengduan Mountains (HHM) (Myers et al., 2000) and the Andes Mountains (Hoorn et al., 2013) have both been identified as biodiversity hotspots (Körner and Spehn, 2019). IBT holds that large islands have lower extinction rates and higher species diversity because they provide more available resources, ecological niches, and habitat types (MacArthur and Wilson, 1967; Wu, 1990; Han, 1994). Studies on the growth habits of Rhododendrons seem to confirm that this theory also applies to mountains regions (Gibbs et al., 2011), as mountains with large areas provide sufficient habitats. Furthermore, mountains with smaller areas have lower species diversity (Si et al., 2017), most likely due to increased competition between species. This relationship between island area and species diversity has been shown to exist in other terrestrial habitat islands, such as grassland patches in the agro-pastoral ecotone (Zhang et al., 2021), forest islands (Lovei et al., 2006), alpine plants (Sklenář et al., 2014).

Habitat heterogeneity also helps explain Rhododendron diversity. This finding is consistent with previous studies on Rhododendron evolution and genomics worldwide (Shrestha et al., 2018b; Xia et al., 2022). The main reason that habitat heterogeneity is thought to impact Rhododendron diversity is because environments with complex topography and diverse climates can provide a rich variety of habitats, creating more niches to accommodate more species (López-Pujol et al., 2011; Stein et al., 2014). For example, vegetation zones along elevation gradients in mountain regions provide various habitat conditions for the growth of Rhododendron species (Cox and Cox, 1997; MacKay and Gardiner, 2016). Furthermore, the complex topography of mountains has provided refuges for species under extreme climatic conditions in the past, e.g., during the large temperature changes of the Quaternary glacial and interglacial periods (Fjeldså et al., 2012). These refuges may also reduce extinction rates by minimizing the need to have strong dispersal ability to track climatic fluctuations (Sandel et al., 2011). Finally, complex topography creates geographical isolation, facilitating new species through allopatric speciation and adaption to various environmental conditions (Steinbauer et al., 2016; Xing and Ree, 2017).

The shared contribution of IBT and habitat heterogeneity on Rhododendron diversity patterns explained the largest part of decomposed variance (Fig. 4), suggesting that IBT and habitat heterogeneity are strongly associated. Several studies have confirmed that the impact of island area on species diversity may include the mediating effect of diverse habitats (Han, 1994; Keppel et al., 2016; MacDonald et al., 2018). In addition, mountains have strong climatic and topographic heterogeneity (Perrigo et al., 2019). Hence, large ‘mountain islands’ increase inherently high habitat heterogeneity, further promoting species diversity. The HHM contains a dense array of mountains and an extremely complex environment (Kumar, 2012). This unique environment creates many habitats, ecological niche spaces, and favorable conditions for Rhododendron (Shrestha et al., 2018a; Namgay and Sridith, 2020), which may provide one explanation for why it is global center of Rhododendron diversity (Shrestha et al., 2018b). Rhododendron is an ancient plant that colonized the southeast of the Himalayas around 30–40 Ma (Shrestha et al., 2018b). Previous studies showed that Rhododendron and several other seed plants may have experienced rapid diversification during the Miocene (Wen et al., 2014; Xia et al., 2022), which coincided with the formation of the HHM (Harrison et al., 1992). This further indicates that the emergence of Hengduan and surrounding mountains and the development of environmental heterogeneity have formed rich ecological niches and geographical isolation, which jointly enhanced the diversity of Rhododendron and provided an important ‘cradle’ for plants in East Asia (Xing and Ree, 2017).

Although the relative contributions of the contemporary climate and paleoclimate change hypotheses was not high, together these hypotheses accounted for nearly half of the total explanatory variance (Fig. 4). Rhododendron generally prefer environments with high rainfall, high humidity, and mild climates (Gibbs et al., 2011; Kumar, 2012). The climate fluctuation and glacial movement since the Quaternary may have affected the modern distribution of Rhododendron (Fang and Ming, 1995; Dynesius and Jansson, 2000). Furthermore, several refuges concentrated in the HHM may have ensured the survival of many plants and animals during the glacial period, including Rhododendron species (Zhan et al., 2011; Srinivasan et al., 2014).

4.2. The application of island biogeography theory in mountain islands

This study confirms that IBT is applicable ‘mountain islands’. This claim is supported by our findings that Rhododendron diversity conform to predictions of IBT. Specifically, we found that Rhododendron diversity is significantly correlated with the ‘area of mountain islands’, the ‘volume of mountain islands’, ‘distance to the mainland’ (isolation degree), and that the number of shared species is correlated with ‘distance to the mainland’. These findings are consistent with previous studies that have found mountains share attributes with ‘true islands’, including limited size, isolation, and dispersal limitation (Ding and Zheng, 1996; Fang et al., 2004; Hughes and Atchison, 2015).

Most studies of classical IBT have used the planimetric area of islands as the measurement of ‘island species-area’ relationships (Lindgren and Cousins, 2017; Zhou et al., 2019). Here, we showed that species diversity is more highly correlated with the volume of ‘mountain islands’ than with their planimetric area. Previous studies have found a similar relationship between species diversity and the volume of islands. For example, tree and beetle species diversity are positively correlated with the volume of some forest vegetation (LaRue et al., 2023); similarly, specialist herbivore species diversity has been positively correlated with the volume of desert shrub islands (Sanchez and Parmenter, 2002). The stronger correlation between species diversity and ‘mountain island’ volume can be explained by the increased niche space and environmental heterogeneity along the elevational plane of mountains (Coverdale and Davies, 2023), especially for Rhododendron, which occupy habitats across vast elevational ranges (Fang and Ming, 1995).

Our findings also confirmed that both species diversity and the number of shared species are influenced by the distance between ‘mountain islands’ and the ‘mainland’. Specifically, Rhododendron diversity is lower in ‘mountain islands’ that are greatly isolated from the ‘mainland’. In addition, fewer species are shared between ‘mountain islands’ and the ‘mainland’ when they are separated by large distances. Rhododendron have small, light, and winged seeds, which enable them to be spread over long distances by wind and animals (Stephenson et al., 2007; Wang et al., 2014; Ma et al., 2022). However, long distances between ‘mountain islands’ and the ‘mainland’ increase dispersal resistance and provide fewer opportunities for Rhododendron species to migrate, decreasing the number of island species and creating large differences between populations. Conversely, short distances promote high species diversity and species similarity. The ‘mountain islands’ surrounding the ‘mainland’ may have more similar natural environments than those farther away, which may explain their species similarity (Ménde-Castro et al., 2018).

Although studies on island species diversity have largely focused on the balance between immigration/colonization and extinction, speciation (e.g., isolation and divergence, hybridization and polyploidy in plant) is also thought to play a significant role (Presgraves and Glor, 2010; Flantua et al., 2020), either through adaptive changes after colonization or by acting as another form of colonization (Heaney, 2001; Alzate et al., 2018). Several studies have confirmed that the formation rate of island endemic species through in situ speciation are positively correlated with island area and isolation (Whittaker et al., 2008; Valente et al., 2020). Local speciation may play a key role in present diversity patterns of Rhododendron. Phylogenetic analyses of global Rhododendron (Shrestha et al., 2018b; Xia et al., 2022) revealed that the evolutionary radiations of this genus occurred in the HHM and Malay Archipelago in the Miocene, resulting in higher speciation rate of Rhododendron in Southeast Asia than in other regions of the world (e.g., Europe, North America). However, it remains unclear how speciation, immigration/colonization, and extinction interact to shape regional Rhododendron diversity.

4.3. Diversity patterns and conservation implications for Rhododendron

Understanding the spatial patterns of species diversity plays a critical role in the conservation of species and diversity (Yu et al., 2017; Shrestha and Wang, 2018). Mountains have long been biodiversity hotspots and refuges that, play an important role in biodiversity conservation planning (Körner et al., 2016; Rahbek et al., 2019; Li et al., 2022). Our study confirms that mountain regions are hotspots of Rhododendron, especially in the HHM region. Therefore, strengthening mountain conservation is a priority for protecting Rhododendron diversity. However, several studies have shown that a large proportion of mountain regions with Rhododendron hotspots are not protected by nature reserves (Mao, 2010; Shrestha and Wang, 2018; Zhao et al., 2023). Furthermore, global climate change has been shown to disproportionally affect species with wide elevational distributions (Ma et al., 2014; Elsen et al., 2018; Vincent et al., 2019). We therefore strongly recommend the establishment of protected areas for Rhododendrons that inhabit mountain regions, and strengthen the elevation gradient protection of Rhododendron species.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (NO. 41901060). We are grateful to anonymous reviewers for their thoughtful comments on the manuscript.

CRediT authorship contribution statement

Yanwei Guan: Writing – review & editing, Writing – original draft, Visualization, Formal analysis, Data curation. Yongru Wu: Writing – review & editing, Data curation. Zheng Cao: Writing – review & editing, Supervision. Zhifeng Wu: Writing – review & editing, Supervision. Fangyuan Yu: Writing – review & editing, Supervision, Funding acquisition, Data curation, Conceptualization. Haibin Yu: Writing – review & editing, Supervision. Tiejun Wang: Conceptualization.

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.