Skip to main content

Potential distributions of seven sympatric sclerophyllous oak species in Southwest China depend on climatic, non-climatic, and independent spatial drivers


Key message

An ensemble modelling approach was performed to predict the distributions of seven sympatric sclerophyllous oak species in the Hengduan Mountains of Southwest China. Spatial eigenvector filters revealed missing factors in addition to commonly used environmental variables, thus effectively improved predictive accuracy for the montane oak species. This study identified a richness center of sclerophyllous oaks, which provides a reference for proper conservation and utilization of oak resources.


As key species and important trees for construction- and fuel-wood, montane sclerophyllous oaks (Quercus sect. Heterobalanus) in the Hengduan Mountains of Southwest China are threatened by climate change, habitat fragmentation, and human activities.


This study aims to simulate the potential distributions of seven sympatric sclerophyllous oak species with an emphasis on exploring the relative importance of climatic, non-climatic, and additional spatial factors.


We performed an ensemble modelling approach of six ecological niche models in combination with spatial eigenvector filters to predict the potential distributions of seven oak species.


The results elucidated that temperature seasonality, followed by land use/cover and the human influence index were the most critical variables controlling oak species distributions. Regardless of the selected algorithm, the best performing models for most oaks combined climatic and non-climatic factors as well as additional spatial filters.


It is necessary to strengthen the conservation of oak species at the junction of Sichuan and Yunnan Province where we found the richness center of the studied oaks. Our research provides essential insights for the rational conservation and management of sclerophyllous oak species, suggesting that spatial constraints might reflect limited ability of migration under future climate change.

1 Introduction

The Hengduan Mountains, located in the southeastern Tibetan Plateau, consist of a series of spectacular north–south ridges and valleys (Royden et al. 2008). This region contributes to the core of Asian freshwater’s supply by serving as an essential catchment zone for many great rivers in Asia (Favre et al. 2015). Moreover, it is recognized biogeographically as a land passage between the floras of Eurasia and the equatorial region (Huang et al. 2016) and a global biodiversity hotspot (López-Pujol et al. 2011). Benefiting from a high elevational variation and vertical climate gradient, the Hengduan Mountains provide intact habitats that range from tropical or subtropical zones to a cold temperate zone (Sun et al. 2017), sustaining over 12,000 vascular plants species of which 3500 are endemic (Xing and Ree 2017).

Quercus sect. Heterobalanus consists of several species, which are montane sclerophyllous oaks in the Hengduan Mountains; these trees have recently attracted much interest from physiological ecologist because of their specific adaptations to the harsh environment and unique traits (Yang et al. 2009; Zhou et al. 2015; Du et al. 2017). These oak species, relying on their high tolerance to low temperatures, drought, soil impoverishment, and strong ultraviolet radiation, can survive in extreme environments such as above the forest line and steep montane gorges where other broad-leaved trees are challenging to establish and grow (Zhang et al. 2005). Forests dominated by oak species of the section Heterobalanus are the second-largest forest type in the Hengduan Mountains after dark coniferous forests and play a pivotal role in biodiversity maintenance, soil and water conservation, and carbon storage (Yang 1990; Xu et al. 2016; Lu et al. 2018).

Currently, the oak species of the section Heterobalanus are roughly estimated to be distributed from the southeastern Tibet Plateau to Kashmir, with the Hengduan Mountains being their geographic distribution and differentiation center (Yang et al. 2009). Historically, these oaks have undergone considerable human disturbances, such as harvesting and burning, all of which continue today (Zhou 2001; Li et al. 2018). Despite their strong abilities to regenerate via sprouting (Zhu et al. 2012), frequent acute disturbances can lead to fragmented and isolated populations, or even the local extinction in some locations (Zhou 2001; Song et al. 2019). Therefore, a more comprehensive understanding of the oak species distributions across the Hengduan Mountains and the interactions of their driving factors is needed for the proper conservation and resource utilization measures.

Ecological niche models (ENMs, see Franklin 2009) quantify the correlation between species occurrences and various kinds of environmental variables to describe the realized ecological niche or habitat suitability of a species (Guisan and Zimmermann 2000). Various ENM algorithms have been implemented to successfully simulate the potential distributions of many tree species, such as keystone species (Liao et al. 2020a), rare species (Thurm et al. 2018), invasive species (Gallien et al. 2012), and economically important species (Zizka et al. 2020). To overcome the variability and to reduce the uncertainty caused by various algorithms, a recent improvement to ENMs has been to ensemble those forecast into a single, final prediction (Thuiller et al. 2009; Gritti et al. 2013; Breiner et al. 2018).

Generally, climate could determine a species distribution through some key variables, notably over a large scale (Vedel-Sørensen et al. 2013; Vessella et al. 2015). However, non-climatic factors, such as topography, geomorphic, edaphic, and land use variables, may be of greater importance for determining distribution ranges over smaller scales (Feng et al. 2016; Zuckerberg et al. 2016)—in some cases, even equaling or exceeding the importance of climatic factors (e.g., Estrada et al. 2016; Figueiredo et al. 2018). Furthermore, compared with natural factors, intense human activity could dramatically alter species’ habitat suitability and thus accelerate their distribution dynamics (Byg and Salick 2009; Li et al. 2018). Accordingly, some studies have argued that the results generated by climate-only ENMs are coarse predictions of actual distributions, and suggested integrating more comprehensive climatic and non-climatic predictors into the ENM process for improving its performance and accuracy (Austin and Van Niel 2011).

Meanwhile, naturally occurring plant species usually show non-random aggregated spatial patchiness (Fahrig 2003). Apart from the current environmental factors, such patterns can also be related to historical effects like time-limited expansions from refugia (Gaston 2009; Dullinger et al. 2012; Nobis and Normand 2014) as well as intrinsic biotic interactions, such as competition and predation (Svenning et al. 2014). Such spatial constraints may prevent species from adapting their distributions to changing climate conditions; in other words, species are left unable to unobstructed expand, even if the changing environments are suitable to their survival.

In contrast to the realized distributions that species successfully occupy in geographical space, their potential distributions yielded by ENMs are often calculated without such spatial constraints. Consequently, the species’ potential distributions are often over-predicted compared with the actual distribution (Dubuis et al. 2011). To rectify this mismatch, eigenvectors (spatial filters) that effectively capture otherwise unconsidered constraints like limited migration are used as additional factors in ENMs to produce more accurate results (Blach-Overgaard et al. 2010; Vedel-Sørensen et al. 2013; Cardador et al. 2014). Nevertheless, the number of such studies remains limited, especially for trees (Blach-Overgaard et al. 2010; Vedel-Sørensen et al. 2013).

To our knowledge, rough distributions of oak species in the Hengduan Mountains have been investigated based on local field expedition or forest inventory data (Zhou et al. 1995). These pioneer findings, however, may insufficiently reflect the entire distributions because a large number of inaccessible areas were not considered. Furthermore, there has been little evaluation of the relative importance of multiple factors that may drive those oak species distributions in general. Therefore, examining the effect of spatial constraints in addition to climatic and non-climatic drivers on the potential distributions of tree species is a promising methodological approach, especially in the Hengduan Mountains, where the topographic and climatic conditions are complex.

To sum up, we hypothesized that the reason for shaping the current distribution patterns of oak species of the section Heterobalanus in the Hengduan Mountains might simply be explained not only by climatic factors but also by non-climatic and additional spatial constraints reflecting unconsidered factors. Our study aimed to figure out the following: (1) the current potential distribution of sclerophyllous oak species across the Hengduan Mountains; (2) the relative importance of climatic, non-climatic, and other spatial constraints for the distributions of these oak species; (3) the location of the core region where sclerophyllous oak species co-occur in the Hengduan Mountains. Therefore, we selected seven sympatric sclerophyllous oaks of the Hengduan Mountains to simulate their potential distributions based on ensemble ENMs, with an emphasis on exploring the relative importance of climatic, non-climatic, and other spatial constraints. Our research provides insights for the conservation and use of sclerophyllous oak species, thereby enhancing the prospect of the Hengduan Mountains’ nature conservation and sustainable resource management.

2 Materials and methods

2.1 Study area

The Hengduan Mountains region consists of a group of mountain ranges and valleys on the edge of the Tibetan Plateau in Southwest China, most of which run roughly north to south. The range of Hengduan Mountains is still inconclusive, but traditionally, it includes the western Sichuan Province, the northwestern portions of Yunnan Province, and the easternmost part of the Tibet Autonomous Region (Fan et al. 2009). In our study, since the oak species of the section Heterobalanus are widely distributed on the southeastern Tibetan Plateau, the range of the Hengduan Mountains follows its broad definition (see Chen 1984). Therefore, the Hengduan Mountains was identified as an area (ca. 5.5 × 105 km2) within the geographical coordinates of 24°–34° N and 95°–105° E (Fig. 1).

Fig. 1
figure 1

Geographical extent of the Hengduan Mountains. The locations of cities and rivers are provided for reference

2.2 Occurrence data

According to the Flora of China (Editorial Committee of Flora of China CAS, 1999), there are seven species of Quercus sect. HeterobalanusQuercus aquifolioides Rehder & Wilson, Q. guyavaefolia H.Lév., Q. pannosa Hand.-Mazz., Q. longispica Hand.-Mazz., Q. senescens Hand.-Mazz., Q. spinosa David and Q. monimotricha Hand.-Mazz.—that are mainly distributed in the Hengduan Mountains (Huang et al. 1999; Zhou et al. 2003). The scientific names of the seven oak species were verified as acceptable based on the Backbone Taxonomy (, These species have similar habitats, preferring areas situated between 2500 and 3600 m a.s.l. (see Appendix Table 3). We collected 1741 village-level occurrences of the seven oak species from the following five sources: (1) National Plant Specimen Resource Center (NPSRC,, (2) Chinese National Specimen Information Infrastructure (NSII,, (3) Global Biodiversity Information Facility (,, (4) private observation databases, and (5) published literature. Because some of these specimens lacked geographical coordinates, the collections were georeferenced using the available location descriptions found on the labels of those specimens. We excluded all cultivated occurrences according to observation notes and their locations, that is, occurrences in parks, botanical gardens, and nearby buildings. Due to the longevity (Plomion et al. 2018) and low dispersal ability (Moran and Clark 2012) of the oak species, we assumed that the mismatch between the period 1960–2015 (occurrence data) and the period 1960–1990 (current climatic data) would not bias the results following Dyderski et al. (2018). To avoid pseudo-replications, we used “GeoRes” option of “SSDM” package (Schmitt et al. 2017) in R (R Development Core Team, 2017) for thinning the occurrences within a grid of a resolution of 30 arc-seconds (~ 1 km at the equator). This step reduced sampling bias while retaining a high amount of information. Finally, the number of the remaining occurrences of each species were as follows: Q. aquifolioides (237 out of 452); Q. longispica (103/144); Q. monimotricha (148/204); Q. pannosa (159/222); Q. spinosa (275/538); Q. senescens (75/118); Q. guyavaefolia (54/63).

2.3 The selection of predictors

For climatic predictors, we initially selected 19 bioclimatic variables with 30 arc-seconds (~ 1 km at the equator) spatial resolution for the period of 1960–1990. These layers were obtained from WorldClim v1.4 (; Hijmans et al. 2005). We first ran the ensemble models (parameter settings and variable evaluation were consistent with the subsequent modelling as described below) to examine the relative importance of 19 bioclimatic variables for each species according to the output of ensemble models. As extreme temperature and precipitation are often considered the most critical limiting factors affecting tree growth in montane regions (Zhang et al. 2005), such nine variables (Appendix Table 4) with average contribution more than 5% were thus selected for the subsequent analysis. To remove multicollinearity among the remaining nine variables, we deleted highly correlated factors based on two criteria, that is (a) the variance inflation factor (VIF) value of each factor less than 10 (Kubota et al. 2015), and (b) the pairwise Pearson correlation coefficients |r| < 0.7 (Dormann et al. 2013). The above analyses were performed using the “vifstep” function of the R package “usdm” (Naimi 2017) and the “layerStats” function of R package “raster” (Hijmans 2019), respectively. Four independent bioclimatic variables were finally selected to represent the climatic predictors (CLIM), i.e., (1) BIO3, isothermality; (2) BIO4, temperature seasonality; (3) BIO15, precipitation seasonality; and (4) BIO19, precipitation of coldest quarter (Table 1; Appendix Fig. 7).

Table 1 The environmental predictor sets used in the different ensemble models

For the non-climatic predictors (NCLIM), we integrated topography and geomorphic, edaphic, land use/cover classes and human impact variables of the study area (Appendix Fig. 7). These variables have demonstrated their importance in previous studies of species distributions (Flojgaard et al. 2011; Zuckerberg et al. 2016; Figueiredo et al. 2018). For topography and geomorphic factor, we chose profile curvature (POC), which was downloaded from EarthEnv (, Amatulli et al. 2018). POC is a measurement of the rate of change of slope and is considered to be related to some ecological processes, such as soil moisture and nutrient cycling (Thomas et al. 2015). Classification of soil texture (CST) was used as an edaphic predictor and was obtained from the China Soil Map-Based Harmonized World Soil Database (HWSD) (, which classifies Chinese soil types into 11 different categories. The layer of land use/cover classes (LUCC) with 22 categories is provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) ( The human impacts were represented by a human influence index (HII) that measures human influence on the terrestrial ecosystem using several human activity effects factors (Sanderson et al. 2002). This layer was obtained from the Last of the Wild Project v2 ( The above non-climatic predictor layers were all at the same spatial resolution with climatic variables (Table 1; Appendix Fig. 7).

Spatial variables (SPAT) derived from spatial eigenvector mapping (hereafter called spatial filters) were calculated using SAM v3.0 software (Rangel et al. 2010). This method presupposes that the spatial arrangement of data (i.e., sample locations or, as in our study case, the regular grid cells of the whole study area) can be translated into a set of predictor variables that capture the spatial effects at different spatial scales. In SAM, spatial filters were computed by constructing a pairwise distance matrix among all grid cells by using their geographical coordinates, that is latitude and longitude. This matrix was submitted to principal coordinates of neighbor matrices (PCNM) to generate the spatial filters (Borcard and Legendre 2002). Because of the computational limitations of SAM, spatial filters were created using coordinates at a coarse resolution (~ 20 km) and then interpolated to the same resolution as ca. 1 km by using the inverse distance-weighted method (IDW) in ArcGIS v10.3 (Esri, Redlands, CA, USA,, following Blach-Overgaard et al. (2010). How to choose a suitable spatial filter is still controversial; here, the approach we followed is consistent with that of De Marco et al. (2008) and Blach-Overgaard et al. (2010), both of which used an equal number of spatial filters with other predictors. In our study, the number of spatial filters is eight based on the above four climatic and four non-climatic variables. The spatial filters represent relatively large- and medium-scale spatial patterns (Appendix Fig. 8). Finally, we calculated the VIF values of all 16 predictor layers and calculated pairwise Pearson correlation coefficients among all predictors to ensure all VIF values < 10 (Table 1) and |r| < 0.7 (Appendix Fig. 9).

2.4 Overlap of species niches

The environmental niches of the seven sympatric oak species were described and compared from occurrence and environmental data using the framework in Broennimann et al. (2012). This framework applies kernel smoothers to densities of species occurrence in gridded environment space to calculate metrics of niche overlap. We applied principal component analysis (Pearson 1901) to transform all environmental predictors plus the spatial filters (CLIM + NCLIM + SPAT) into two principal components (see Liao et al. 2020b for details). The environmental space was then defined by the axes of principal component analysis and was divided into a grid of 100 × 100 cells. The smoothed density of occurrence in each cell for each species was calculated by applying Gaussian kernel density function (Silverman 1986). Meanwhile, the occupancy of the environment was calculated and compared to estimate the niche overlap metric (Schoener 1970; reviewed in Warren et al. 2008).

2.5 Ensemble models

To simulate the potential distributions of seven oak species, we used an ensemble modelling approach encompassing six algorithms: generalized linear models (GLM), boosted regression trees (GBM), multivariate adaptive regression splines (MARS), maximum entropy (MAXENT), random forests (RF), and support vector machines (SVM). These algorithms were available in the “SSDM” package (Schmitt et al. 2017). For successful running each algorithm, 1000 pseudo-absence points were randomly selected from the background localities (as suggested by Barbet-Massin et al. 2012). To assess the average predictive performance of the ensemble model, we repeated a holdout method ten times and randomly selected a subset of 70% of presences and pseudo‐absences for calibration and the remaining 30% for validation each time. To reduce the uncertainty of each algorithm, we performed ten replications of each algorithm to calculate the mean area under the curve (AUC) of the cross-validation (Bahn and McGill 2013). Models with AUC values of > 0.7 were included in the final ensemble ENM, and the weighting of each algorithm prediction was based on its true skill statistic (TSS) following Liu et al. (2005). Further, we adopted the best default settings for each algorithm recommended by Schmitt et al. (2017).

In the present study, we performed three models with different combinations of the predictors for each oak species: namely (1) climatic predictors alone (CLIM); (2) environmental variables that included both climatic and non-climatic predictors (CLIM + NCLIM); (3) all environmental predictors plus the spatial filters (CLIM + NCLIM + SPAT), which was the full model. We used the maximum training sensitivity plus specificity (MTSS) as a threshold recommended by Liu et al. (2005), and divided the habitat suitability into two classes: not suitable (0–MTSS), suitable (MTSS–1). Model performance was evaluated by comparing the area under the curve (AUC) of the receiver operation curve (ROC) as well as the TSS (Schmitt et al. 2017) using the holdout data. Additionally, the differences in AUC values between paired models produced by different predictor sets were tested using Wilcoxon signed-rank tests (Blach-Overgaard et al. 2010).

2.6 Statistical analysis

To assess the unique importance of each predictor, we used the output of the SSDM package based on the Pearson’s coefficient between a full model and a model with each environmental variable omitted in turn (Schmitt et al. 2017). Furthermore, variation partitioning was performed to quantify the unique contribution of the three predictor sets to the variation of predicted habitat suitability (Borcard et al. 1992; Heikkinen et al. 2004). The variation in the predicted habitat suitabilities of all seven oaks were decomposed into the three groups of explanatory variables, i.e., CLIM, NCLIM, and SPAT variables, using a series of (partial) regression analyses with redundancy analysis (RDA), as implemented in the SAM v3.0 software (Rangel et al. 2010).

3 Results

3.1 Model validation

The AUC value of all 126 models (i.e., seven species × six algorithms × three predictor sets, see Liao et al. 2020b for details) ranged from 0.701 to 0.886 (mean 0.778 ± 0.036), while the TSS value ranged from 0.403 to 0.772 (mean 0.555 ± 0.072), indicating huge variability of model performance existed among species, predictor sets and algorithm (Fig. 2). For most species, except for Q. monimotricha and Q. senescens, the full model CLIM + NCLIM + SPAT showed superior predictive ability than the CLIM and CLIM + NCLIM models by its slightly but not significantly higher median AUC and TSS values (Fig. 2 a and b). Furthermore, predictive performances varied among the six modelling algorithms without considering the differences between species and predictors (Fig. 2 c and d). Three of the modelling algorithms, namely MAXENT, RF, and SVM, shared a better predictive capacity; these three algorithms produced higher median AUC and TSS values compared with the average predictive abilities (Fig. 2 c and d). In our case, the best performing algorithm was RF (median AUC = 0.793; median TSS = 0.586), while the lowest-ranked method was GLM (0.766; 0.531).

Fig. 2
figure 2

Predictive abilities of the ensemble models for three predictor sets (a, b), and of six modelling algorithms for seven species of Quercus section Heterobalanus (c, d). Predictive ability is expressed as the area under the curve of the receiver-operating curve (AUC) and the true skill statistic (TSS). The same letters or “ns” indicate no significant differences between models by Wilcoxon sign rank test comparisons. Dashed lines in red represent a cutoff threshold for AUC (0.7) used for selecting the algorithms and represent a good performance for TSS (0.5), while in black represent the average AUC and TSS value of all models

3.2 Niche overlap and potential distributions

The niche overlap between species varies from 0.44 to 0.82 (Appendix Table 5). The highest overlap occurred between Q. monimotricha and Q. spinosa, while the lowest occurred between Q. aquifolioides and Q. longispica (Appendix Fig. 10 and Table 5).

Overall, the potential distributions simulated by the CLIM and CLIM + NCLIM models were less different but more extensive than those produced by the CLIM + NCLIM + SPAT model (Fig. 3; Appendix Table 6). For example, only three oaks analyzed in our study truly grow in the eastern part of the Tibetan Autonomous Region; nevertheless, the western margin of Hengduan Mountains was deemed suitable habitat for all oaks by both CLIM and CLIM + NCLIM models. Among the oak species, Q. aquifolioides had the relatively larger range irrespective of the predicted model sets. The CLIM and CLIM + NCLIM models indicated suitable areas of Q. aquifolioides which fully covered the Hengduan Mountains, even extending northward to south Ganlan Prefecture and westward to east Bowo Prefecture. However, these areas were absent in the full (CLIM + NCLIM + SPAT) model, which also identified the significant differences in the southeast corner of Hengduan Mountains when compared with our partial (CLIM or CLIM + NCLIM) model (Fig. 4a). Q. longispica was found to have a relatively consistent range in the southern Hengduan Mountains by the CLIM and CLIM + NCLIM models, while the CLIM + NCLIM + SPAT model confined its range to southwestern Hengduan Mountains, specifically in Diqing and Lijiang Prefectures and southern Ganzi Prefecture (Fig. 4b).

Fig. 3
figure 3

The simulated potential distribution areas of ensemble models for all studied species based on three combinations of predictor variables

Fig. 4
figure 4

The occurrence records and the predicted distributions of seven oak species based on the ensemble models that used different variable sets. Binary maps were generated based on MTSS threshold, i.e., Q. aquifolioides (0.180), Q. longispica (0.234), Q. monimotricha (0.208), Q. pannosa (0.205), Q. spinosa (0.154), Q. senescens (0.168), Q. guyavaefolia (0.266)

For Q. monimotricha, Q. pannosa, and Q. spinosa, the CLIM + NCLIM + SPAT results showed these species primarily distributed in the southeastern part of Hengduan Mountains, i.e., from Nujiang-Aba Prefecture eastward to the edges of Hengduan Mountains. By contrast, their distributions predicted by the CLIM and CLIM + NCLIM models extended farther, up to eastern part of Linzhi and Changdu Prefectures in the west, along the valleys of the Lancangjiang, Nujiang, and Jinshajiang Rivers (Fig. 4 c, d, and e). Potential distributions for Q. senescens and Q. guyavaefolia were largely restricted to the southern part of Hengduan Mountains, covering the northern Yunnan Province and southwestern Sichuan Province, when predicted by the CLIM + NCLIM + SPAT model (Fig. 4 f and g).

When overlapping the highly suitable areas (suitability value > 0.6) of seven oak species simulated by their CLIM + NCLIM + SPAT models in our study area, a richness center of oak species was located within a triangular area situated among the Shangri-La, Eryuan, and Muli Prefectures, at the junction of Sichuan and Yunnan Provinces (Fig. 5).

Fig. 5
figure 5

A sketch of the presumed contemporary geographic distribution center for the seven oak species. The triangular region at the junction of Sichuan and Yunnan Provinces, China, is based on the habitat suitability (> 0.6) from the full ensemble models (CLIM + NCLIM + SPAT models)

3.3 The drivers of species potential distributions

As the CLIM + NCLIM + SPAT model better predicted the potential distributions of oak species across the Hengduan Mountains, we focused on the full models’ predictor variables when evaluating their relative importance. Generally, predictors differ in the importance of shaping the distribution of each species; however, temperature seasonality (BIO4) and land use/cover classes (LUCC) were consistently indispensable to all species’ distributions. Human influence index (HII), followed by isothermality (BIO3) and profile curvature (PC) together supplemented influences of their distributions greatly (Table 2). Besides, some spatial eigenvectors (e.g., S4 and S6) were ranked into the top-five factors in the distribution of specific species (Table 2). For all studied species, about 75 to 87% of the total variation in predicted habitat suitability was explained by CLIM, NCLIM, and spatial variables together, leaving corresponding 13 to 25% of the variation unexplained (Fig. 6). The main source of the variation was the rather dissimilar unique contributions from the three groups of predictors. The fractions of the pure effect of climate alone ranged from 21.5 to 31.0%; however, the fractions of the pure effect of non-climate alone were all less than 10%. Notably, the fraction for the spatial filters varied considerably, with four species less than 10% (Fig. 6 c, d, e and f), two close to 20% (Fig. 6 b and g), and one species over 41.3% (Fig. 6a).

Table 2 Mean relative importance of each variable to the final ensemble model of each species studied for all environmental predictors plus spatial filters model (CLIM + NCLIM + SPAT). The variables highlighted in italics are the five top-ranked importance factors for each species
Fig. 6
figure 6

Results of variation partitioning for the predicted habitat suitability simulated by our ensemble models in terms of fractions of variation explained. Variation of the occurrence probability is explained by three groups of explanatory variables: CLIM (climatic predictors), NCLIM (non-climatic predictors) and SPAT (spatial predictors). ac Unique effects of CLIM and NCLIM factors and SPAT variables, respectively. dg Fractions indicating their joint effects. Negative values are shown as “< 0%” in the figure. The font-size increases with the intensity of the effect. Note that for large datasets like our study, negative explained variances are possible when explanatory variables explain no more variation than random normal variables would (Legendre 2008)

4 Discussion

4.1 Potential distributions of the oak species

Accurately simulating the potential geographical distributions of plant species at broad scales in the Hengduan Mountains is a challenge because of data deficiencies, which mainly arise from many inaccessible high mountains and deep ravines (Zhang et al. 2014), but also from potentially incomplete sets of model predictors. Based on a comprehensive collection of available data, we used an ensemble modelling approach to analyze the potential distributions of seven oak species in a remote mountainous area, the Hengduan Mountains. This approach could better balance the different results produced by different algorithms when calibrating the same dataset (Thuiller 2004; Guisan et al. 2007), like in our study.

In the present study, jointly judged by the resulting AUC and TSS values, and predicted area, the full ensemble model, that is CLIM + NCLIM + SPAT model, was better than both the partial models (CLIM and CLIM + NCLIM models). In addition, the full models are more consistent with previous field surveys which indicated most of the oaks occurring between 25°–32° N and 98°–102° E (Liu et al. 2008). These results support Bucklin et al. (2015) who suggested that using spatial filters as additional species-related limiting factors and coupling them with climate factors in ENMs is more likely to approach the actual distributions in nature. The full models, however, has the highest number of predictors and might therefore be more adjusted to the observations with an increasing risk of overfitting. Nevertheless, the model evaluation was based on a cross-validation strategy using holdout data, which would have detected significant overfitting (Merow 2014).

Interestingly, the richness map showed that areas of high suitability were almost always located on high mountains and on both sides of river valleys across the Hengduan Mountains, yet they were missing from flat hilly plateaus, such as found at the junction of Sichuan and Qinghai Provinces, and from western and northern parts of Hengduan Mountains. In reality, the oak species may distinguishable as two life-history types with corresponding terrain requirements. Oak forests adapted to cold-dry, or sub-humid environments typically grow at high elevations, whereas the warm or hot, dry valleys, or extremely cold areas are inclined to develop an oak scrub-like community (Yang 1990; He et al. 1997). Therefore, our simulations matched this habitat preference of the oak species. Although the modelling results indicated differences between oak species, their distribution areas overlap largely. This characteristic of sympatry differs considerably from the evergreen sclerophyllous Quercus species of the Mediterranean that seem morphologically, and physiologically similar to the oak species of the section Heterobalanus (Manos et al. 2001) yet have disparate distributions (López-Tirado et al. 2018). As noted above, the oak habitats can range from deeply incised valleys to the high mountains, thus forming a clear vertical distributional gradient at the Hengduan Mountains (Yang 1990). Previous researchers found that all oak species of this section share the similar habitat requirements and have a similar tolerance of drought and cold (Yang et al. 2009; Tang 2015), exhibiting inconspicuous vicarious distribution in altitude (Yang 1990; He et al. 1997). In this way, these oak species are able to coexist on mountains to form a mosaic distribution (Yang 1990; He et al. 1997), which is also supported by the results of overlapping species niches (Fig. 10 and Appendix Table 5).

Based on field sampling records, He et al. (1997) proposed the region located at “the mid-point of the W-shape of Jinshajiang River (towards the north of Lijang)” as the most species-rich area for oak species of the section Heterobalanus. In this study, we proposed a “triangular”-shaped region consisting of Shangri-La, Muli, and Eryuan prefectures which locate at Jinshajiang River and Yalongjiang River Basins (Fig. 5). Although this region also includes areas without a single species, a supplemented richness hotspot close to Muli might also be of high importance with respect to genetic diversity (Zhang et al. 2014). More importantly, sporopollenin and fossils of oak species of the section Heterobalanus have been found at numerous locations in and near to this region (Huang et al. 2016), which indicates a suitable environment for oak since long, including interglacial-glacial cycles of late Pleistocene. However, under the rapid climate change, ecotypes at southern (or hottest/driest) range limits might be very important as well adapted proveniences for reforestation or (quasi-intraspecific) assisted migration (Liao et al. 2020a). Therefore, focusing only on the detected region of highest predicted species richness might not be appropriate under future climate change, although it is of great value to the biodiversity conservation of these oak species at present.

4.2 Climatic factors governing the oak species distribution across Hengduan Mountains

The results of variation partitioning showed that climatic variables have a stronger influence on all oak species, compared to the non-climatic variables in our study. This finding is consistent with the current perspective that climate is the primary determinant for the geographic distributions of plant species over large spatial scales. Evidence from the fossil record confirms that paleoclimate change has shaped the current oak species distribution patterns (Zhou et al. 2007; Meng et al. 2017). Earlier, Zhou et al. (1992) had proposed that ancestors of the oak species of our study were distributed within subtropical evergreen broad-leaved forests, as companion species, at the junction of China, Vietnam, and Myanmar during the Miocene. Because of successful adaptive evolution, they have gradually become prevalent in the eastern Himalayas and Hengduan Mountains under the development of an arid and cold climate caused by the uplift of the Tibetan Plateau. In particular, the ancestors of these oak species evolved specific ecophysiological traits, such as dense hairs, thick cuticles, lignified epidermal cell walls and cuticles, and a low stomatal density, which collectively increased their tolerance of dry-cold climates (Yan et al. 2019). Other studies have further indicated that the oak species of our study are favored by the seasonal alternation of dry and wet periods (Yang 1990; Tang 2015), similar to the Mediterranean and North America evergreen sclerophyllous Quercus spp. (Yang 1990). However, unlike water availability, which predominately determines the geographical distribution ranges of evergreen sclerophyllous Quercus species in the Mediterranean (Ogaya et al. 2003; Hidalgo et al. 2008), our results showed that effects of temperature-related variables exceeded those of water-related variables. The unique geography and climate conditions of the Hengduan Mountains probably explain this discrepancy. Firstly, the oak species of our study mainly inhabit the high mountains in southeastern Hengduan Mountains, where the annual precipitation is approximately > 900 mm due to the East Asian and Southwest monsoon, thus making it a humid region (Yang 1990). Those high mountains also always give rise to low-temperature environments, which more strongly govern montane plant growth than drought (Yang 1990; Zhou et al. 2015).

Notably, temperature seasonality was the most crucial variable for most oak species (Table 2). The Hengduan Mountains climate is characterized by hot, rainy summers and cold, dry winters; hence, the driest quarter is traditionally from February through April and accompanied by low temperatures (Zhang et al. 2005; Yang et al. 2009). A recent popular perspective was reported that the temperatures of the non-growing season (i.e., winter) are more closely linked to montane plant distributions than the annual or summer temperatures (Ladwig et al. 2016; Choler 2018). Extreme low winter temperatures could lead to freezing and frost injuries to plants that inhibit their growth in profound ways, such as via membrane damage and photosynthesis reduction (Ladwig et al. 2016). Alternatively, too high winter temperatures will reduce the soil moisture available to plants via increased evapotranspiration and disturb tree phenology and root activities (Choler 2018). Over the past several decades, both an increased winter temperature and the response of montane plants to winter warming on the Tibetan Plateau have been extensively observed (IPCC 2014). Hence, our finding also suggested these oak species would be at risk in the global context to future climate warming.

4.3 Effects of non-climatic variables on oak species distribution

Researchers have elucidated that when more non-climatic variables are added to ENMs—such as edaphic conditions, land cover, or human influence—along with climate, they could prove useful because they refine the potential distribution range (Hidalgo et al. 2008; Estrada et al. 2016; Figueiredo et al. 2018). In our case, the non-climatic variables of land use/cover classes, human influence index, and profile curvature were of high importance for many of the oak species distributions, while the classification of soil texture was opposite. However, the pure fraction of non-climatic variables explained less than 10% of the total variation in predicted habitat suitability, and for some species even less than 3%, e.g., Q. senescens and Q. guyavaefolia. Nevertheless, the importance of non-climatic factors showed an overlap with climatic factors as indicated by fraction d in Fig. 6, which strongly suggests some critical potential role of non-climatic variables for the distribution of these oak species.

In contrast to cultivated plants introduced elsewhere by humans, the oak species of our study maintain an almost natural distribution pattern in the Hengduan Mountains. Being the most important firewood in the eastern Tibetan Plateau, they were reported to be more disturbed by logging than by livestock grazing, planned burns, or infrastructure (Zhou 2001). Despite the aboveground biomass lost in oak forests destroyed by logging, these sclerophyllous oak species can quickly recover via their excellent resprouting ability to form scrub vegetation (Zhu et al. 2012). For this reason, the human impact on the distribution of oak species might in our study be lower than climate effects, which is consistent with findings of López-Tirado and Hidalgo (2014). The effect of land use/cover, however, seems to be one more noteworthy driver in comparison to other non-climatic variables. The transformation of land use through a number of disruptive measures (defragmentation, deforestation, fire, harvesting, etc.) would dramatically alter the distribution of species (Sala et al. 2000).

Even though human activity did not seem to be stronger related to the oak species ranges on a large scale, negative effects should not be dismissed. Hengduan Mountains regions, in particular, the “three parallel river-running areas,” are the richest in hydropower resources in China and the building of critical hydropower projects are underway. Reservoir water impounding and related infrastructure could damage or extirpate many oak populations, thereby affecting their distribution ranges at local or regional scales. Unfortunately, these oak species currently lack effective protected areas as well as relevant protection strategies. Therefore, it is necessary to strengthen the conservation of proveniences for oak species in the detected hotspot region, especially in the existing provincial nature reserve, that is Lashihai wetland, Haba Snow Mountain, and Jade Dragon Snow Mountain and to formulate effective conservation policies, such as the prohibition of excessive logging or unsustainable expansion of hydropower.

4.4 Spatial constraints and its implications

In our study, the CLIM + NCLIM + SPAT model took into account the spatial filters and exhibited highest AUC and TSS value for most species, suggesting that it was superior to models excluding this type of predictor. Spatial eigenvectors are sensitive to migration lags (Blach-Overgaard et al. 2010) and could reveal missing “classic” variables including climate predictors. For example, for Q. aquifolioides, pure spatial predictors explained 41.3% of the total variation, which even exceeded climatic variables in this case (Fig. 6a). These findings from the Hengduan Mountains are in accordance with other studies showing that spatial filters help to better predict palm distributions at continental (Blach-Overgaard et al. 2010), national scales (Tovaranonte et al. 2015), and regional scales (Vedel-Sorensen et al. 2013). We noticed that in the CLIM + NCLIM + SPAT model, its results removed some areas where no occurrence records have been found and the potential distribution ranges seem to resemble a large extent the actual distribution ranges (Fig. 4). This discrepancy means that these suitable habitats were more likely to have been over-predicted due to the lack of more accurate predictors. On the other side, spatial filters seemed to cause also some false negatives in regions with only few observations (Fig. 4). Therefore, we suggest that adding spatial filters in ENMs can be essential for obtaining more realistic potential distribution ranges, but at the same time, the risk of false negatives should be taken into account.

Although the mechanisms underpinning spatial range constraints are not yet thoroughly understood, many researchers believe that limited migration involving barriers of seed dispersal and/or time-lagged migrations from refugia play an important role. Meng et al. (2017) and Du et al. (2017) revealed the distributions of both Quercus sect. Heterobalanus and Q. aquifolioides were fairly stable during Last Glacial Maximum. Accordingly, we suppose that geographic features may chiefly impede the dispersal of these oak species. They are inferred to have originated in the northern Indo-China Peninsula, and its surrounding areas (Zhou 1992) and the Hengduan Mountains contains a large number of north–south mountains and rivers, which could provide channels for animal and water dispersing seeds of oak species to expand northwards along the river valley and mountain ridges (He and Jiang 2014). However, according to two paleogeographic studies (Qu et al. 2011; Yang et al. 2012), this special north–south direction topography is apparently a barrier for many species of Hengduan Mountains expanding westward to the interior Tibetan Plateau. This could partly explain why the oak species of our study insufficiently occupy the western part of the Hengduan Mountains which are suitable according to the CLIM model. Moreover, whether a species fills its entire potential range must depend on its intrinsic dispersal and migration ability. Apart from dispersal limitation, other researchers have posited that spatial constraints are associated with mechanisms of biotic interactions, for example, competition (Svenning et al. 2014). Early investigators, actually, partly attributed the extant distribution patterns of oak species of the section Heterobalanus to species competition (Yang 1990; Zhou 1992). They found most oak species were at a disadvantage when competing with broad-leaved species in the tropical and subtropical climate, with only Q. spinosa able to colonize eastern China from the Hengduan Mountains (Zhou 1992). Similarly, on the hilly plateau plane, in the western part of Hengduan Mountains, oak species could not outcompete alpine grassland plants that are better adapted to drought (Yang 1990). Therefore, the weak competitive ability could be another reason for explaining why these species were unable to expand westward through the Hengduan Mountains.

On the Tibetan Plateau, endemic species, including oaks, are faced with global climate changes, such as winter warming (IPCC 2014), which would strongly alter their extant distribution patterns. Interestingly, it is noteworthy that species with the broadest distribution (Q. aquifolioides) and those with smallest distributions (Q. longispica and Q. guyavaefolia) were most affected by spatial factors (Figs. 3 and 6a, b, and g). Therefore, we could speculate that factors behind these spatial constraints, especially for Q. guyavaefolia, could hinder their colonization of suitable climate habitats in the future that would amplify the negative effect of global climate change for these oak species. The decrease or disappearance of key species like the oak species of our study would affect species composition, ecosystem functions and local economy, e.g., local farmers’ income by affecting the production of the lucrative edible mushroom, Tricholoma matsutake. On the other hand, suitable but currently uncolonized areas detected by the climate-only models may serve as candidate areas for future range expansion and assisted migration. Although other uncertainty factors like biological interaction or dispersal ability (Thibaud et al. 2014) were not considered in this study, the results might fill data gaps in survey records and estimate priority regions for future surveying and monitoring efforts.

5 Conclusions and recommendations for future research

Our results revealed ensembles of ENMs in combination with spatial eigenvector filters as a useful tool to model the potential distributions of sclerophyllous oak species in the Hengduan Mountains. The temperature-related variables in climatic predictors were the most critical factors and non-climatic variables, including land use/cover and human impact, had a clearly visible influence on oak species’ distributions on a large scale. The “triangular”-shaped region identified by our models at the junction of Sichuan and Yunnan Provinces, is the current geographic distribution center for seven oak species in the Hengduan Mountains. Our predictions also highlight climatically suitable areas, which can be used for developing more sustainable reforestation strategies, although nature reserves currently protect only a small part of these areas. Hence, these tree species now need improved protection policies, such as reduced logging intensity and the reinforcing of the existing protected reserves in the above triangle. Furthermore, the availability of past, present, and future potential distributions of oak species would help to detect putative climate refugia. Potential changes in the distributions of oak species under future climate change should be further studied to better understand both the effects of global change and possible management options.

Data availability

The used online database of all species distribution data and environmental layers have been cited in the manuscript and indicated with accessible URL. The observations used and datasets generated and/or analyzed during the current study are available in the GitHub repository (Liao et al. 2020b) at


Download references


We appreciate the help of Ms. Zhang Fengying, Ms. Tan Xue, Ms. Huang Dan, and Mr. Tang Tianwen with the collection of the occurrence data and its verification. We thank Dr. Olusanya Abiodun Olatunji for his valuable suggestions while preparing our manuscript and for language polishing on it. We also thank the handling editor, Dr. Laurent Bergès, and the anonymous reviewers for their valuable comments.


This study was funded by the National Key Research and Development Program (Grant No. 2016YFC0502101) and Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (Grant No. 2019QZKK0303).

Author information

Authors and Affiliations


Corresponding authors

Correspondence to Kaiwen Pan or Lin Zhang.

Ethics declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Additional information

Communicated by Laurent Bergès

Handling Editor: Laurent Bergès



Table 3 Overview of studied oak species according to the description of the Flora of China
Table 4 Pearson correlation coefficients for 19 bioclimatic variables based on layers within the range of the Hengduan Mountains are shown. The nine variables highlighted in bold indicated that their average contributions to the initial ensemble ENMs for seven studied species exceeded 5%. Uncorrelated variables highlighted in italics were represented for climatic predictors and were ultimately selected to simulate the potential distributions of seven oak species
Table 5 The niche overlap between species. The metric varies between 0 (no overlap) and 1 (complete overlap)
Table 6 Predicted suitable areas based on three combinations of predictor variables: climatic predictors (CLIM), climatic and non-climatic predictors (CLIM + NCLIM), and all environmental predictors plus spatial filters (CLIM + NCLIM + SPAT)
Fig. 7
figure 7

The distributions of four climatic and four non-climatic variables used in the models: Isothermality (BIO2/BIO7) (× 100) (°C), temperature seasonality (standard deviation × 100) (°C), precipitation seasonality (coefficient of variation: mean/SD × 100), precipitation of coldest quarter (mm), land use/cover classes, profile curvature, classification of soil texture, and human influence index

Fig. 8
figure 8

The geographic patterns of the eight spatial filters used, showing the spatial relationship among cells (1 km × 1 km grids). Increasingly lighter (white) colors indicate larger numerical values of the eigenvectors

Fig. 9
figure 9

Pearson correlation coefficients for all 16 predictors, i.e., four climatic (bio3, bio4, bio15, and bio19), four non-climatic (CST, HII, LUCC, and POC) and eight spatial predictors (S1~S8). All coefficients are shown as percentages and are less than 0.7

Fig. 10
figure 10

Niches of the seven sympatric oak species along the two first axes from the principal component analysis (PCA) of all 16 predictors. The color that changes from red to blue shows the density of the occurrences of the species by cell

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liao, Z., Nobis, M.P., Xiong, Q. et al. Potential distributions of seven sympatric sclerophyllous oak species in Southwest China depend on climatic, non-climatic, and independent spatial drivers. Annals of Forest Science 78, 5 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: