Regional systematic conservation planning needs to consider current and future biodiversity distribution patterns under a changing environment (Schloss et al., 2011; Heller et al., 2015; Jones et al., 2016; Tukiainen et al., 2017; Reside and Adams, 2018). Many studies have revealed that ecosystems and species’ ranges will shift with environmental changes (Hoffmann and Sgrò 2011; Aguilée et al., 2016; Levine et al., 2016), which may prevent existing reserve networks from effectively conserving biodiversity in the long-term (Schloss et al., 2011; Scriven et al., 2015; Regos et al., 2016); this is of particular concern in high mountains and plateaus (Ackerly et al., 2010; Zhang et al., 2012; Zomer et al., 2015; Lehikoinen et al., 2019). More critically, it is difficult to accurately predict future distribution patterns of biodiversity, considering the uncertainty around how environments will change and what impacts those changes will have on ecosystems and species distributions (Schloss et al., 2011; Aguilée et al., 2016; Jones et al., 2016).

To deal with this dilemma, some researchers have proposed dealing with systematic conservation planning by coupling environmental variables (e.g., topography, climate, soil) with biodiversity surrogates (Hufford et al., 2014; Heller et al., 2015; Tukiainen et al., 2017). The essence of idea is that conserving environmentally heterogeneous landscapes supports diverse species and communities in a region, which is consistent with niche theory (MacArthur 1970; Ricklefs 1977; Stein et al., 2014). Coupling environmentally heterogeneous landscapes into existing reserve networks on a regional scale would help conserve current and future biodiversity in a changing environment (Tingley et al., 2014; Van Schalkwyk et al., 2017).

The above studies have been largely based on positive regional heterogeneity-diversity relationships (HDRs) (Veech and Crist 2007; Stein et al., 2014; 2015). However, there is no universal method for understanding regional HDRs because of environmental variables and biodiversity and therefore the methods used to measure them vary among regions (Bar-Massada and Wood, 2014; Chocron et al., 2015; Stein et al., 2015). Regional HDRs could be positive, negative, unimodal, or without significant correlation (Lundholm 2009; Gazol et al., 2013; Laanisto et al., 2013; Bar-Massada and Wood, 2014). Hence, it is necessary to study HDRs in specific regions before selecting suitable regional environmental variables for systematic conservation planning.

Additionally, with the rapid development of human society and economies, inevitably, the biodiversity would be disturbed by human disturbance almost all over the world (Lehikoinen et al., 2019). While regions with a low degree of human disturbance and a high degree of naturalness (natural areas) have low conservation costs but high potential conservation value (Theobald et al., 2012; Triviño et al., 2018). And ecosystem structure in natural areas is primarily intact and largely unaffected by human influence; thus, it can effectively protect biodiversity and ecological processes (Kormos et al., 2016), thereby providing the ideal option for regional biodiversity conservation. Furthermore, only natural areas can present original HDRs meaningful for regional systematic conservation planning. Besides, Seiferling et al., (2014) performed a comprehensive analysis and found that the HDRs in natural areas showed complex and equivocal relationships. Therefore, the HDRs of natural areas would provide important guidance for the optimization of regional protected areas system.

Hence, this study selected all 20 national and 12 of the 38 provincial nature reserves in Yunnan, China to represent the relatively natural areas and to examine the relationship between environmental heterogeneity (EH) and biodiversity. Although only 32 nature reserves were selected for this study, these selected nature reserves represent over 80% of Yunnan’s geographic elements (soil types, topographic units, climatic units), 90% of its ecosystem types, and 90% of its national key protected wild animal and plant species. More importantly, nature reserves with biodiversity information have been collected as far as possible. For each selected nature reserve, we calculated three types of biodiversity indices and 17 EH parameters, and explored three questions: 1) how well does EH parameter correlate with the biodiversity index? 2) What are the primary and important EH parameters for the biodiversity index? and 3) to what extent can biodiversity index be explained by relevant primary EH parameter(s)? We hope to find explicit HDRs that can be used to optimize the Yunnan reserve network and a methodology that can be used in other regional studies.

Yunnan is located in the southwest border of China (97°3′E–106°12′E, 21°08′N–29°15′N) (Fig. 1), the southeast edge of the Qinghai-Tibet plateau, where there are strongly topographic and climatic gradient changes and significant environmental heterogeneity. As a well-known global biodiversity hotspot (Myers et al., 2000; Yang et al., 2004; Zhang et al., 2014; Zomer et al., 2015), although Yunnan has established 20 national, 38 provincial, 56 municipal or prefectural, and 46 county-level (160 in total) nature reserves, covering 7.4% of Yunnan Province (Forestry Department of Yunnan Province, 2017), the impact of environmental change has reduced the conservation effectiveness of nature reserves and increased the likelihood of species extinction (Zhang et al., 2014; Zomer et al., 2015). Therefore, optimizing nature reserve network of Yunnan for better conserving its biodiversity is facing environmental changes (Zhang et al., 2014; Zomer et al., 2015; Wang et al., 2018). However, knowledge is limited to the relationship between Yunnan’s biodiversity distribution and environmental heterogeneity.

Figure 1.  Location of the study area, showing 32 nature reserves, variation in terrain elevation within Yunnan, China

Each nature reserve had detailed scientific survey reports and other study materials (Appendix S1) for its biodiversity status (ecosystem types, plant and animal species). These reports or study materials were collected to gather data on vegetation formations and plant and animal species in the reserves.

The quantitative method (Li et al., 2011) was used to calculate three types of biodiversity indices: ecosystem diversity index (EI), plant diversity index (PI), and animal diversity index (AI). EI_i (i represents each nature reserve from 1–32, hereafter the same) was measured directly by the number of vegetation formations derived from the vegetation map of each nature reserve. PI_i and AI_i were measured by the numbers of plant and animal species recorded in nature reserve i using equations (1) and (2), respectively. Like a previous study (Li et al., 2011; Song et al., 2016), we assigned weights of 100 and 50 to national I- and II-level protected wild species, respectively, which we used to evaluate regional biodiversity conservation values (Song et al., 2016).

where VP_i refers to the number of vascular plant species recorded in nature reserve i, NPP_Ii and NPP_IIi are the numbers of national I- and II-level protected wild plant species, respectively, in nature reserve i (Yang et al., 2016).

where VS_i and IS_i refer to the numbers of vertebrate species and insects, respectively, recorded in nature reserve i, and NPA_Ii and NPA_IIi are the numbers of national I- and II-level protected wild animal species, respectively, in nature reserve i (MFPRC and MAPRC, 1988; Yang et al., 2016).

Previous studies have demonstrated that topographic, climatic, and soil heterogeneity show strong correlations with plant and animal diversity (Irl et al., 2015; Stein et al., 2014; 2015). Topographic heterogeneity plays a more important role in shaping the species distribution and the pattern of species diversity than elevation itself (Tukiainen et al., 2017). Meanwhile, several researches have shown that climatic heterogeneity determinate species diversity pattern of terrestrial vertebrates and vascular plants, especially variables associate with water and energy availability (Veech and Crist, 2007; Stein et al., 2014; 2015; Tukiainen et al., 2017). Additionally, Zhang et al (2012) have found that several climate factors, such as annual mean temperature, temperature annual range, annual precipitation, precipitation of driest month, and precipitation seasonality, are crucial to predict the distribution of plant diversity in Yunnan. Thus, it is very meaningful to explore the relationship between climate heterogeneity and biodiversity in Yunnan, China. Furthermore, edaphic heterogeneity is critical to driving the diversity pattern (Hufford et al., 2014; Hulshof and Spasojevic, 2020).

Collectively, combined with the accessibility of environmental data we derived three subject areas: topography, climate and soil, and 17 environmental variables in each nature reserve (Table 1). All variables were produced at the same resolutions: 30 arc-seconds (except soil data). To better measure the environmental heterogeneity of 32 nature reserves, the study divided the soil into 144 types following the Second National Soil Survey (National Soil Survey Office of Yunnan Province, 1996). With regard to topographic variables: altitude, slope, and aspect, we reclassified the altitude into 55 classes by 100m intervals. The slope was reclassified into six categories: flat (0–5°), gentle (5°–15°), pitched (15°–25°), steep (25°–35°), hard (35°–45°), and extreme (≥ 45°). The aspect was reclassified into north, northeast, east, southeast, south, southwest, west, northwest, and no slope, a total of 9 categories. According to Zomer et al., (2015), 13 climate variables selected for nature reserve in Yunnan can be reclassified into 33 classes. All those 17 reclassified environmental variables were used to calculate 17 EH measures at each selected nature reserve.

Studies have shown that the Shannon-Wiener index of environmental variables is an effective measure of EH (Stein et al., 2015). The 17 EH parameters (Table 1, Appendix S2) were calculated:

where $ {H}_{ij}^{{'}} $ the heterogeneity of environmental variable j (1–17) of the selected nature reserve i (1–32), k represents the reclassified type or class of each environmental variable, among which soil variable has 144 types, topographic variables (altitude, slope, and aspect) have 55, 6 and 9 classes, respectively, and each climatic variable has 33 classes. S_ijk is the area occupied by the type or class k of environmental variable j in nature reserve i. S_i is the area of nature reserve i.

We used Spearman rank correlation function in R package Corrplot (Dormann et al., 2013; Li and Wang 2013; Wei and Simko, 2016) to analyze the collinearity of EH parameters (Fig. 2) and performed a single predictor ordinary least squares (OLS) regression to examine how well EH parameters explain biodiversity index (by determination coefficient R²_adj). If the Spearman correlation coefficient was higher than 0.7, the EH parameter with higher explanation power was identified as the primary EH parameter. We used multi-predictor OLS regression to examine how well the subject explained the biodiversity index. Because the previous study found that correlations were observed among the areas, EH parameters, and biodiversity indices, hence, the area factor was also considered a key subject area in this study (Li and Wang 2013).

Figure 2.  The correlation coefficient matrix among the 17 environmental heterogeneity parameters and area of selected 32 nature reserves in Yunnan, China. Abbreviation of environmental heterogeneity parameters can be seen in Table 1; A, area (hm²)

Whereafter, we identified the primary environmental variables of EI were area and the EH measures of soil, slope, aspect, mean temperature of driest quarter, temperature annual range, annual precipitation, precipitation of coldest quarter and precipitation seasonality; the primary environmental variables of PI were area and the EH measures soil, slope, aspect, mean temperature of coldest quarter, temperature annual range, annual precipitation, precipitation of coldest quarter and precipitation seasonality; and the primary environmental variables of AI were area and the EH measures of soil, altitude, slope, aspect, temperature annual range, precipitation of wettest quarter, precipitation of coldest quarter and precipitation seasonality. Considering these environmental variables, the optimal EH measure interpretation model for each biodiversity index was constructed using the Akaike information criterion (AIC) (Quinn and Keough, 2002). The standard regression coefficient of each EH parameter in the optimal model reflected the degree of importance for the corresponding biodiversity index. All statistical analyses were performed in R version 3.4.4 (R Core Team, 2016).

The spatial distribution of biodiversity in the natural areas of Yunnan varied greatly. The three nature reserves with the highest ecosystem diversity were Xishuangbanna, Nangun River, and Gaoligong Mountains reserves. Gaoligong Mountain, Dawei Mountain, and Wenshan reserves showed the highest plant diversity, while the richest diversity of animal species was found in Xishuangbanna, Nangun River , and Tongbiguan reserves. Overall, Xishuangbanna, Gaoligong Mountains, and Tongbiguan reserves were the most diverse nature reserves (Table 2).

Based on the spatial distribution of biodiversity (Fig. 3), the distributions of plant and animal diversity were negatively correlated with latitude (R² = 0.251, P = 0.002; R² = 0.252, P = 0.002, respectively, (Figs. 3d, 3g). The PI and AI showed a decreasing trend from south to north, whereas no significant correlation was detected between ecosystem diversity and latitude (R² = −0.017, P = 0.491) (Fig. 3a). From west to east, ecosystem and animal diversities were significantly negatively correlated with longitude (R² = 0.123, P = 0.028; R² = 0.232, P = 0.003, respectively) (Figs. 3b, 3h), whereas this was not observed for plant diversity. Concerning elevation, ecosystem diversity was irregularly distributed, and PI and AI decreased with elevation (R² = 0.167, P = 0.011; R² = 0.113, P = 0.033, respectively) (Figs. 3c, 3f, 3i).

Figure 3.  The relationship between biodiversity indices and latitude, longitude, altitude within selected nature reserves of Yunnan, China. a, b and c refer to the ecosystem diversity (EI); d, e and f refer to the plant diversity (PI); g, h and i refer to the animal diversity (AI)

The results of single predictor OLS models indicated that temperature annual range heterogeneity was the strongest predictor of EI (R² = 0.469, P = 0.001) (Table 3), and precipitation of coldest quarter heterogeneity was the best predictor of the diversity pattern in PI (R² = 0.405, P = 0.001), followed by annual precipitation, precipitation seasonality, temperature annual range, and mean temperature of coldest quarter, which all had good explanatory power for the PI (R² = 0.252, P = 0.002; R² = 0.245, P = 0.003; R² = 0.199, P = 0.004; and R² = 0.159, P = 0.023, respectively). For the AI, only the temperature annual range heterogeneity was able to sufficiently explain the variation (R² = 0.281, P = 0.01).

The multi-predictor OLS models (Table 3) showed that for EI, PI, and AI, climatic heterogeneity had stronger explanatory power for biodiversity indices (R ²= 0.463, P = 0.001; R² = 0.423, P = 0.001; R² = 0.263, P = 0.014, respectively) than soil or topographic heterogeneity. Topographic heterogeneity almost had no explanatory power for three biodiversity indices. Soil heterogeneity could only influence EI (R² = 0.193, P = 0.007), and had very weak explanatory power for the other biodiversity indices. The area had a limited ability to interpret EI and PI (R² = 0.301, P = 0.001; R²= 0.103, P = 0.034, respectively).

The AIC optimal model (Table 4) constructed using the soil, mean temperature of driest quarter, temperature annual range, annual precipitation heterogeneity, and area effectively explained the spatial variation in EI (R² = 0.569, P = 0.001). The optimal model formed by the annual precipitation and precipitation of coldest quarter heterogeneity effectively explained the spatial variation in PI (R² = 0.545, P = 0.001). The model consisting of the temperature annual range and precipitation seasonality heterogeneity explained the spatial variation in AI (R² = 0.296, P = 0.002). According to the standard regression coefficients, the temperature annual range heterogeneity was the most important environmental factor affecting EI and AI (0.647, P = 0.001; 0.498, P = 0.01, respectively), and precipitation of coldest quarter heterogeneity was the most important environmental variable for the PI (0.525, P = 0.001). Hence, we concluded that the temperature annual range and precipitation of coldest quarter heterogeneity were not only positively correlated with biodiversity, but also the primary driving forces of the biodiversity patterns in the natural areas of Yunnan.

Our analysis revealed that soil heterogeneity best explained the EI (R² = 0.193, P = 0.007, Table 3), which is consistent with the thesis that different combinations of soil types could provide a variety of nutritional levels (habitat selection) for biological communities (Hufford et al., 2014). Given the AIC optimal model for EI, which couples soil, climatic (temperature, precipitation) heterogeneity, and area factors; this further confirms the basic conditions (complex physical environment) required to form a variety of ecosystems (Lapin and Barnes, 1995). Soil heterogeneity cannot effectively explain species diversity (AI, PI), perhaps because the heterogeneity of soil structure and organic matter composition may relate more to species diversity at the finer scale (Cramer and Verboom, 2017).

Himalayan orogeny is generally believed to drive environmental factors to rapidly compress in space, and extreme changes in topography may be one of the reasons for the formation of high species diversity in Yunnan (Xing and Ree, 2017). However, only the altitude heterogeneity presented explanatory power for the PI ( P = 0.018) in this study. Topographic heterogeneity most likely had weak explanatory power because its direct effects on biodiversity were not significant at the non-uniform scale of ‘nature reserves’. More explicitly measured environmental heterogeneity explained biodiversity patterns better than did crude topographic measures such as mean slope, altitude range, and mean aspect (Bailey et al., 2017). Additionally, topographic heterogeneity may have explained biodiversity so well because it is an excellent proxy for several sources of climatic and soil heterogeneity (Bailey et al., 2017). For example, the highest plant diversity in this study was in the Gaoligong Mountains (Table 2), which are distributed among the Hengduan Mountains, one of the world's biodiversity hotspots (Xing and Ree, 2017). Complex regional topographic conditions in these areas can increase the climatic heterogeneity (Wang et al., 2015), which may have a greater effect on the distribution of species diversity in topographic heterogeneous regions (Irl et al., 2015). Perhaps we can assume that topographic heterogeneity might indirectly affect the distribution of biodiversity by affecting climatic heterogeneity in Yunnan. Hence, the synergistic effect of topographic and climatic heterogeneity on biodiversity should be explored when HDRs are studied in topographically complex areas.

We found that climatic heterogeneity sufficiently explained biodiversity (EI, PI, AI), and its variables influenced biodiversity the most, particularly temperature annual range and precipitation of coldest quarter heterogeneity (Table 3). Yunnan is located at low latitudes, where the environmental tolerances of species are weaker because where the annual temperature range is lower than in high latitude regions. According to the climate stability hypothesis, Stevens (1989) suggested that regions with a stable climate are more likely to promote the formation of narrow niches for species. Additionally, Klopfer and MacArthur (1960, 1961) proposed that a smaller annual range of climatic conditions reduces niche overlap and supports species with narrower niches. This is consistent with a large number of species with narrow ranges in Yunnan and the fact that Yunnan is an important global centre for endemic species (Li, 1994; Wang and Zhang, 1994; Huang et al., 2012; 2016). In this study, the water-related variable of precipitation of coldest quarter heterogeneity had the greatest influence on the PI. The result that precipitation of coldest/driest quarter is collinear (Fig. 2) indicates that precipitation heterogeneity is the main driving factor for the PI when species or communities face limited resources. Moreover, these results indicate that the water-energy dynamics hypothesis (Veech and Crist, 2007; Stein et al., 2014) could explain the plant diversity in the natural areas of Yunnan.

The spatial scale (spatial extent or cell size) represents an unavoidable problem in regional HDR research. In this study, the area size affected the HDRs, particularly on the ecosystem diversity (R² = 0.301, P = 0.001). Through climate-vegetation models, it is generally known that climate type is strongly correlated with vegetation distribution (Kaplan and New, 2006). Regarding species diversity, the results were consistent with a previous study concluding that the spatial scale does not affect the overall trend in HDRs (Seiferling et al., 2014). Additionally, the negative correlation in HDRs is primarily due to an increased degree of heterogeneity, reduces the effective habitat area of each species and increase the probability of random extinction (Laanisto et al., 2013; Chocron et al., 2015). Finer-scale environmental heterogeneity is likely to intensify habitat fragmentation and therefore threaten regional biodiversity (Stein et al., 2015). According to species-area curves, however, only areas of a certain size (such as a nature reserve) can effectively protect biodiversity. More importantly, larger areas can more effectively regulate the introduction of exotic species and the rate of species renewal (Stein et al., 2014; Bailey et al., 2017). Therefore, these discussions confirm that regional HDR research is more suitable on a macro scale and that HDRs could help generate more rational and effective conservation planning for biodiversity.

Despite data limitation, this study did not explore thoroughly the discrepancies of environmental heterogeneity on the drivers of biodiversity in different nature reserves, and could not elucidate the mechanisms of the formation of biodiversity patterns in different spatially distributed nature reserves. However, we provide a paradigm for studying regional HDRs. With sufficiently accurate geographic distribution information of species collected, conservationists can study micro-scale natural area HDRs in-depth and devising appropriate conservation strategies for a local area in future. More importantly, micro-scale HDRs can include critical thermal maximum and minimum physiological limitation factors, particularly, the role of environmental heterogeneity on the driving forces of animal diversity in natural areas of Yunnan, China.

Climatic heterogeneity had the best explanatory power for the biodiversity distribution patterns in relatively natural areas of Yunnan. The region’s diverse hydrothermal conditions produce selection pressures for species and promote interspecific diversification or even the formation of new species (Hua and Wiens 2013, Irl et al., 2015). Further studies have shown that more species could coexist within a climatic heterogeneous region, through improving the fitness of a species increases the probability that it exhibits phenotypic plasticity (Gianoli and Valladares, 2012, Lázaro-Nogal et al., 2015). Therefore, we argue that climatic heterogeneity is the primary driving force for the species diversification and biodiversity patterns in the natural areas of Yunnan. Additionally, climate conditions are more stable in climatic heterogeneous regions (Ackerly et al., 2010), which could allow species to migrate over less distance to locate suitable habitats and reduce species extinction rates under future climate change. Climatic heterogeneity could help mitigate the effects of climate change on biodiversity, and therefore we emphasize that climatically heterogeneous regions have large conservation significance for biodiversity in Yunnan under climate change scenarios.

The effective conservation of ecosystem diversity involves maintaining a region’s important ecological processes and ecological stability, particularly in the context of climate change (Levine et al., 2016). However, conservation in Yunnan currently only takes into account the conservation value of specific vegetation types (Zhang et al., 2013). Although studies showing that vegetation community heterogeneity is most likely a direct driver of species diversity (Stein et al., 2014; Levine et al., 2016), the significance of regional ecosystem diversity is still ignored. In the face of conservation gaps in ecosystem diversity in the province, the region’s complex soil composition and climatic heterogeneity must be prioritized.

Heterogeneity-based priority conservation areas represent a novel approach that could assist in effectively protecting species diversity under climate change, based on understanding regional HDRs (Heller et al., 2015; Paudel and Heinen, 2015). According to the AIC optimal models (Table 4), regions with heterogeneity of annual precipitation or preccipitation of coldest quater or temperature annual range or precipitation seasonality heterogeneous indicate areas with rich plant and animal diversity. Moreover, the results also demonstrate that water-related variables more effectively explained the spatial distribution of the PI, whereas temperature-related factors better explained the AI. In other words, the PI and the AI have different environmental drivers, indicating that one of these indices cannot replace the other. Consistently no single biodiversity surrogate can fully reflect regional biodiversity (Di Minin and Moilanen, 2014, Yang et al., 2016), and this is universal, even in areas with environmental heterogeneity. Therefore, coupling the environmentally heterogeneous regions, which have multiple dimensions of biodiversity, will help increase the effectiveness of biodiversity conservation within priority conservation areas based on environmental heterogeneity.

The study of heterogeneity-diversity relationships in natural areas will help conservationists and decision-makers to have a more explicit recognition of the process of shaping regional diversity patterns and provide scientific support for coupling environmentally heterogeneous areas in future systematic conservation planning. Our research explored the relationships between biodiversity and soil, topographic, and climatic heterogeneity in natural areas of Yunnan. We demonstrated that water-related and temperature-related factors are the most important environmental driver for plant and animal diversity, respectively. In general, climatic heterogeneity holds the most important role in the AIC optimal models and also have appreciable explanatory power to ecosystem (56.9%), plant (54.5%), and animal (29.6%) diversity.

Although this study has some limitation with its biodiversity data and spatial scales, and nor are there further studies of the synergies effects of environmental heterogeneity on biodiversity in different subject areas. Nevertheless, our study of natural area heterogeneity-diversity relationships indicates that climatically heterogeneous areas are maybe pivotal for coupling environmental heterogeneity in systematic conservation planning and optimizing existing protected areas of Yunnan Province in future. To achieve effectively protect the biodiversity of Yunnan under environmental changes, we have some suggestions that we should investigate: 1) the effect of spatial scales on regional HDRs and 2) the conservation effectiveness of coupled climatically heterogeneous regions into systematic conservation planning for biodiversity conservation in Yunnan, China under climate change.

We are grateful to HUA Chaolang at the Yunnan Institute of Forest Inventory and Planning for providing the Yunnan Nature Reserve information.