Throwing light on dark diversity of vascular plants in China: predicting the distribution of dark and threatened species under global climate change

Species distribution data and environmental variables

This study covered the entire area of China. We retrieved vascular plants occurrence data in each province and municipality of China from the Flora of China (http://frps.eflora.cn/) and merged data of municipalities into their nearest provinces, so we got 27 records in total (Fig. 1). In comparison with the data applied to the study by Pärtel, Szava-Kovats & Zobel (2011), derived from a global map and at the scale of 100 x 100 km, although our data are at larger scale, the observed species lists are much more complete, especially for that of Northern China, which was considered difficult to access (Kier et al., 2005). We constructed an individual species dataset for them to show which plants were present or absent. Although the floristic data have been collected since 1959, we assume that they are still relevant because plants, particularly those in the temperate forest, respond relatively slowly to environmental changes (Franklin et al., 2016; Menéndez et al., 2006). We only included plant data at the species level, and we merged all subspecies. We also excluded all cultivated species that were not considered a natural part of the vegetation based on their life history and dispersal traits. Synonyms and conservation status such as threatened species were checked using The Plant List (http://www.theplantlist.org/) and The IUCN Red List of Threatened Species (http://www.iucnredlist.org/). This observed species richness (number of species) distribution shows a latitudinal trend, increasing from north to south in China (Fig. 1).

Climatic variables, such as precipitation and temperature have a significant impact on species’ distribution, particularly at large spatial scales and over long timeframes (Beer et al., 2010; Kelly & Goulden, 2008; Loarie et al., 2009; Williams, Jackson & Kutzbach, 2007). We extracted 19 bioclimatic variables (BIO1–BIO19; Table S1) in years of current, 2050 and 2070, respectively, from Worldclim (resolution, 2.5 s; www.Worldclim.org/). As for the variables in 2050 and 2070, we selected those under representative concentration pathway 4.5 (RCP4.5) and based on the CCSM4 global climate model (GLM). There are four RCPs, ranging from very high (RCP8.5, the temperature will increase to 3.7 ± 0.7 °C) through to deficient (RCP2.6, the temperature will rise to 1.0 ± 0.4 °C) future concentrations (Collins et al., 2013). Based on current status of global warming and the goal that hold average temperature increase to well below 2 °C (IPCC, 2018), in this study, we chose RCP4.5 scenario as the future climate scenario. In addition, according to the projection by Ying (2012), until the year of 2018, the CCSM4 has been showing as the more accurate model than the other 16 GLMs under RCP4.5 scenarios, so in this study, the future climate projections are based on the CCSM4 climate model. Current and future bioclimatic projection data were used in the MaxEnt model to predict habitat suitability of six threatened dark species that will be detailed in the following sections.

Estimation of dark diversity and community completeness

To obtain the best possible estimates of dark diversity in China, we used Beals probability index (Beals, 1984), as recommended by Lewis, Szava-Kovats & Pärtel (2016). This index can be used to estimate the probability of a species occurring in a particular region based on its co-occurrence within other regions and was calculated using the package “vegan” (Oksanen et al., 2016) in R ver. 3.3.1 (R Development Core Team, 2016).

The Beals index is defined as: $$P_{ij}=\frac{1}{S_i-I_{ij}}{\displaystyle \sum _{k\ne j}\frac{N_{jk}{I}_{ik}}{N_k}}$$ where P_ij is the probability that species j will occur at community i, S_i is the number of species at community i, I_ij is the incidence (0, 1) of species j at community i. N_jk is the number of joint occurrence of species j and k (k ≠ j) at community i, I_ik is the incidence (0, 1) of species k at community i, and N_k is the number of occurrences of species k (for further details, see Münzbergová & Herben (2004) and Lewis, Szava-Kovats & Pärtel (2016)).

The probability of occurrence varies among species, regions and depends on the frequencies of species in a particular assemblage. Each species was assigned an individual probability threshold, which was a quantile calculated from a user-defined probability. Based on the comparison of the dark diversity from different probabilities (1%, 5%, 10%) (Table S2), to explore the broader possible species pool in China, it was defined 1% in all regions here. For each area, a species was included in the dark diversity when it was absent from a target region and its occurrence probability was higher than its threshold value.

In addition, we calculated the community completeness index for each region using ln (observed richness/dark diversity) (Pärtel, Szava-Kovats & Zobel, 2013).

Assessing variables that affect dark diversity

We used RF to assess the contribution of bioclimatic variables toward the dark diversity. The RF technique estimates the importance of a predictive variable by evaluating the Out-of-bag (OOB) error increase (Mutanga, Adam & Cho, 2012). In other words, the decrease of prediction accuracy, represented as the percentage of increased mean square error (% IncMSE) when OOB data for the specific variable is switched while all the other variables remained (Mutanga, Adam & Cho, 2012). Here, we extracted a matrix of 3,844 spot records covering the whole area of China, with their 19 current bioclimatic data, and applied it to RF, where their dark diversity values at the regional level are treated as the response variable.

Random forest is a machine learning approach based on classification and regression trees (CART; Breiman, 2001). The RF model-building process is similar to that of CART; only it combines numerous independent trees to reach a final decision (for further details, see Wiesmeier et al. (2011) and Da Silva Chagas et al. (2016)). We estimated the relative importance of each variable in the RF model using OOB randomly selected data. The mean square error (MSE) was calculated as: $${\text{MSE}}_{\text{OBB}}=\frac{1}{N}{\displaystyle \sum _{i=1}^n{\left(Z_i-{\widehat{Z}}_i^{\text{OBB}}\right)}^2}$$ where Z_i is the measured value of the variable and ${\widehat{Z}}_i^{\text{OBB}}$ is the average of all OBB predictions. The MSE_OBB is normalized as it depends on the unit of the response variable.

Random forest modeling was implemented by R ver. 3.3.1 (R Development Core Team, 2016), using the “RF” package (Liaw & Wiener, 2002). The number of trees (ntree) was 500 and the number of randomly selected predictor variables at each node (mtry) was three.

Predicting the distribution of rare dark species with future climate change

We used MaxEnt modeling (Version 3.4.1) to predict changes in the potential areas of six rare dark plant species under global climate change. MaxEnt is a widely used species distribution algorithm, and many studies have compared it with other SDMs to confirm its predictive ability (Elith & Graham, 2009; Elith, Kearney & Phillips, 2010; Tognelli et al., 2009; Williams et al., 2009). It is used to predict the geographic distributions of species based on incomplete information of species (presence-only datasets) occurrence data and environmental variables by calculating the MaxEnt of species distribution (Phillips, Anderson & Schapire, 2006; Phillips & Dudik, 2008).

In this study, we chose the near threatened, vulnerable, and endemic species: Eucommia ulmoides, Liriodendron chinense, Phoebe bournei, Fagus longipetiolata, Amentotaxus argotaenia, and Cathaya argyrophylla from the threatened dark species based on available data and species’ distribution range: all species have a limited distribution and a small dataset of occurrence sites. In total, we found 187 occurrence sites for these species in relevant references, local flora, and the flora of China, which we included in the MaxEnt model (A. argotaenia: 30 sites, C. argyrophylla: 22 sites, E. ulmoides: 32 sites, L. chinense: 35 sites, F. longipetiolata: 29 sites, P. bournei: 39 sites) (Fig. S1; Table 1). According to the records of the references, geographic coordinate data are within an accuracy of about three km. The 19 bioclimatic variables were filtered out some less important ones by using RF model and removed a few high correlated variables by using correlation coefficient. The correlation coefficient was implemented in R ver. 3.3.1 (R Development Core Team, 2016), using the “raster” package (Hijmans, 2018), the cut-off threshold is 0.80.

In our model, we split the dataset into two parts: 75% of the occurrence data were used as training set, while the remaining 25% was used as a test set to evaluate the strength of the model. In this model, we set the 10,000 points as the max number of background sampling entire China. As for the small size of our data, we tuned the regularization multipliers (0.1, 0.25, 0.35, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, and 3) in MaxEnt model, to find out the best one with the smallest corrected Akaike information criterion (AICc) value (Morales, Fernández & Baca-González, 2017). The AICc was calculated by ENMTOOLS (version 1.4.4) (Warren & Seifert, 2011). The area under the receiver operating characteristic curve (AUC) and omission error (minimum training presence as the logistic threshold) was used to quantify the strength of the current scenario models. AUC values range from zero to one, with the values of ≤0.5 indicating that the model performs worse than a random model and values close to 1 indicating that the model performs better than a random model (Swets, 1988). In this study, we considered values of 0.5–0.6 to represent no discrimination, 0.6–0.7 as unaccepted, 0.7–0.8 as accepted, 0.8–0.9 as excellent, 0.9–1 as outstanding and 1 as perfect (Mischler et al., 2012; Phillips, Anderson & Schapire, 2006). Finally, we used the Jackknife test to evaluate which variable made the greatest contribution to the distribution of each of these six species under the current climatic scenario.

Dark diversity and community completeness in China

Different regions in China show various dark diversities, ranging from 0 to 215. In comparison with their observed local richness (from 1,170 to 11,760), their dark diversities are smaller. Meanwhile, the observed richness showed a latitudinal gradient change, increasing from north to south (Fig. 1), while dark diversity exhibited no latitudinal trend, and the values declined from east to west, except for Qinghai. The highest levels of dark diversity were concentrated in three regions: Qinghai (215), the southeast of the northern China plain (Shandong, 159; Henan, 155; Jiangsu including Shanghai, 143), and Heilongjiang (142) (Fig. 2A).

The community completeness indexes showed a similar distribution as their observed richness, almost fully complete in the west communities (except for Qinghai Province), and less complete in the east. In addition, the completeness indexes were lower in the northeast than those in the southeast of China (Fig. 2B).

We identified eight threatened species among the dark diversity, representing a range of threat levels (Table 1).

Explanatory variables for dark species

The relative importance of the 19 bioclimatic variables regarding dark diversity in China, as assessed using RF, is shown in Fig. 3. We found that variables associated with temperature showed more importance than those associated with precipitation (Fig. 3). However, together these variables only explained 71.72% of the variance, thus few unknown factors remain. The most important variables were Temperature annual range (BIO7) and Temperature seasonality (BIO4), Mean temperature of warmest quarter (BIO10), and Precipitation seasonality (BIO15). Other significant variables were Isothermality (BIO3), Mean temperature of wettest quarter (BIO8), Max temperature of warmest month (BIO5), Mean temperature of coldest quarter (BIO11), Mean temperature of driest quarter (BIO9), Precipitation of driest quarter (BIO17), Precipitation of coldest quarter (BIO19), Mean diurnal range (BIO2), and Precipitation of driest month (BIO14). These variables were filtered to use in MaxEnt projection firstly.

Then according to the results of cross-correlation (Fig. S2), we selected seven bioclimatic variables used in the MaxEnt model for current climate scenario (Fig. S2), and the contribution of them differed between six species (Table 2). For these six species, BIO14 in current bioclimatic scenarios was the most important variable when projected their distribution, BIO5 also showed its significance to project the species distribution (Table 2). For the 2050 and 2070 climate scenarios, we selected six bioclimatic variables used in MaxEnt (Fig. S2).

Distribution of dark threatened species in the future

All the models under current bioclimatic scenarios were considered outstanding according to their omission error, AUC and AICc values (Table 2), suggesting that they could accurately predict species distributions.

By overlaying the distribution results from MaxEnt with those from Dark Diversity, we could find out that most of these dark plant potential regions were also their climatically suitable areas from current to 2070 in MaxEnt modeling. However, Jiangsu province was estimated having dark species F. longipetiolata while MaxEnt models predicted it not climatically suitable for F. longipetiolata either under current or future climatic scenarios (Fig. 4A).

In addition, for these six species, the projected climatic suitability in their potential regions would change under future climatic scenarios. Our MaxEnt modeling results showed that most species’ suitable areas would firstly increase by 2050, and then contract by 2070 (Fig. 4). However, E. ulmoides suitable area in Jiangsu province would shrink by 2050 and then expand a little by 2070 (Fig. 4A). We believed these shifts probably related to BIO14 decline from current to 2050, and BIO19 increase from current to 2070, and similar associations have also been found in other species (Bueno et al., 2017; Deb et al., 2017; Noulekoun et al., 2017; Wu et al., 2017).