DOI: 10.5061/dryad.ns1rn8ptx
title: Conservation of woody species in China under future climate and land-cover changes

description: Climate and land-cover changes are major threats to biodiversity, and
 their impacts are expected to intensify in the future. Protected areas
 (PAs) are crucial for biodiversity conservation. However, their
 effectiveness under future climate and land-cover changes remains to be
 evaluated. Moreover, the impacts of climate and land-cover changes on
 multi-dimensions of biodiversity are rarely considered when expanding PAs.
 Using distributions of 8732 woody species in China and species
 distribution models, we identified species that will be threatened by
 future climate and land-cover changes (i.e. species with significant
 projected loss of suitable habitats by the 2070s) under different
 dispersal scenarios. We then estimated the geographical patterns in
 species richness (SR) and phylogenetic diversity (PD) of these threatened
 species, evaluated the effectiveness (i.e. the changes in SR and PD) of
 Chinese PAs, and identified conservation priorities for future PA
 expansion. Approximately 12-38% of woody species will be threatened under
 different scenarios. These species tend to be clustered in the tree of
 life, and their SR and PD show consistent spatial patterns, being highest
 at low latitudes. PAs currently protect 90% of these threatened species.
 However, their SR and PD of threatened species within PAs will decrease by
 30-40% by the 2070s, which reduces the PA effectiveness, especially for
 PAs at low elevations and those with low topographic heterogeneity and
 high natural vegetation loss. The conservation priorities identified from
 the SR and PD of the threatened species are mainly in mountains in
 southern China, the Yunnan-Guizhou Plateau, and Taiwan Island. PA
 expansion and ecological corridors in these regions are needed to conserve
 these threatened species. Synthesis and applications. We present a
 systematic study of the impacts of future climate and land-cover changes
 on the conservation status of woody species and PA effectiveness in China.
 Our results suggest that future climate and land-cover changes will reduce
 PA effectiveness, and the spatial prioritization of biodiversity
 conservation should consider the influences of future global changes on
 biodiversity. These results shed new light on the conservation priorities
 for the post-2020 expansion of PAs in China.
description: <b>Species distribution
 data</b> The county-level
 distributions of woody species in China were compiled from the
 <i>Atlas of Woody Plants in China: Distribution and
 Climate</i>, which contains a total of 11405 native woody species
 (Fang et al., 2011; Wang et al., 2011). We substantially updated the
 database using more than 40 volumes of national and regional floras
 published from 2009-2018 (<i>see</i> Appendix
 <i>S1</i>), and available specimen records in China
 (http://www.nsii.org.cn/). To further improve the accuracy of species
 distributions, the elevational range and habitat types of each species
 were compiled and used to refine their distributions within each county.
 To eliminate the bias of unequal areas on the estimates of species
 diversity, we transformed species distributions into 20 × 20
 km<sup>2</sup> grid cells with the Albers equal-area
 projection. After removing the incomplete grid cells with area less than
 200 km<sup>2</sup> within national borders and coasts, a total
 of 23718 grid cells (ca. 94%) remained. Previous studies suggested that
 the performance of SDMs depends on the sample size of the species
 distributions. Therefore, we excluded species with fewer than 20
 distribution records from our analyses (Thuiller et al., 2011). Finally,
 8732 species were included for the subsequent analyses. <b>Climate and soil
 data</b> The current
 climatic data during the period 1960-1990 with a spatial resolution of
 2.5' (ca. 5×5 km at the equator) were obtained from the
 <i>WorldClim</i> (ver. 1.4; http://www.worldclim.org). In
 total, 19 variables were obtained, which were divided into three groups:
 1) temperature 2) precipitation and 3) climate seasonality. We calculated
 the Pearson’s correlation coefficients among all variables, and they were
 highly correlated with each other within groups (r &gt; 0.7), but not
 among groups. For each group, we chose those with high ecological
 relevance to woody plant distributions following previous studies
 (Thuiller et al., 2011; Zhang et al., 2017). Finally, 6 climate variables
 were chosen: mean temperature of the warmest and the coldest quarters,
 temperature seasonality, precipitation of the wettest and the driest
 quarters, and precipitation seasonality. The climatic data of each grid
 cell were calculated as the average of all data points within
 it. Climate data for
 the 2070s under three climate scenarios (i.e. RCP 2.6, RCP 6.0 and RCP
 8.5) were also obtained from <i>WorldClim</i>. We chose 2070
 as a critical point since a global assessment show that nearly 1/3 species
 will go extinct by 2070 (Román-Palacios & Wiens, 2020). Here, we
 used the data projected by the BCC-CSM1-1.1 model, which could better
 reflect the magnitude of future climate change in China compared with
 other models (Xin et al., 2013). We also included
 four commonly-used soil variables reflecting physical (i.e. sand-clay
 ratio, bulk density) and chemical properties (soil water pH and organic
 carbon concentration). Soil data were obtained from the SoilGrids 250m
 database (https://www.soilgrids.org), and the value for each 20×20 km grid cell was calculated as the average of all data points within it. We used data for the top soil (0-15cm depth), and assumed soil indices were constant over time.    <b>Land cover data</b> Current and future land-cover data with a spatial resolution of 1×1 km<sup>2</sup> were obtained from Li et al. (2017). This dataset includes the following layers: forests, shrublands and grasslands (including savanna-like grasslands), barren, farmland, urban and water. The first three layers reflect the coverage of natural land-cover types, while farmland and urban reflect human-related ones. Decreases in natural land-cover types and increases in human-related ones could significantly influence biodiversity patterns (de Chazal & Rounsevell, 2009). We calculated the proportion of each land-cover type within each grid cell for the current and future periods, and adopted each as an independent predictor in the SDMs. We excluded water because they represent habitats not suitable for terrestrial plants. It is noteworthy that the data of future land cover were projected under earlier emission scenarios (i.e. A1B, A2, B1 and B2). We matched three of these emission scenarios, i.e. B1, A1B and A2, with RCP2.6, RCP6.0 and RCP8.5, respectively (IPCC, 2014). Finally, six bioclimatic variables, four soil variables and five land-cover variables were selected for SDM analyses. <b>Protected area data</b> The PA data in China were obtained from Zhang et al. (2015). We only incorporated national and provincial PAs as most of prefectural and county-level PAs are inefficiently managed due to inadequate funding, staff and equipment (Quan et al., 2009). We further excluded another 30 PAs located in the "mangrove biome", which are normally considered as marine PAs. Finally, our database included 1161 PAs. We rasterized these PA polygons at the same resolution as the species distribution data via layer intersection. A grid cell was identified as a "protected cell" if the percentage of its area covered by PA polygons was higher than a threshold. Here, three thresholds (10%, 30% and 50%) were used, and preliminary analyses suggest that the results based on these thresholds were highly consistent. Therefore, we only presented the results based on the threshold of 30% following the post-2020 conservation target (i.e. a global 30% coverage). To explore the drivers of the effectiveness of PAs under climate and land-cover changes, we evaluated the impacts of five PA characteristics: 1) PA area; 2) topographic heterogeneity; 3) elevation; 4) human footprint and 5) the percentage of natural vegetation change. The area of each PA was obtained from Zhang et al. (2015). Human footprint represents the cumulative impacts of human pressures on environments, and was calculated using eight variables: human population density, buildings, electric infrastructure, roads, railways, navigable waterways, cropland, and pasture (Venter et al., 2016). Elevational range within each PA was estimated as the difference between the highest and lowest elevation and was used to represent the topographic heterogeneity (Stein & Kreft, 2015). Elevation data with a spatial resolution of 30″were obtained from the U.S. Geological Survey (https://topotools.cr.usgs.gov/gmted_viewer/viewer.htm). The natural vegetation here included forests, shrublands and grasslands (including savanna-like grasslands) (Li et al., 2017). <b>Species distribution models (SDMs)</b> Five SDM algorithms were used, including classification tree analysis (CTA), the generalized linear model (GLM), the generalized boosting model (GBM), random forest (RF) and maximum entropy (MaxEnt). These SDMs are widely used in previous studies (e.g. Zhang et al. 2017). For each species, the distribution data were randomly partitioned into training (80%) and testing (20%) data, which were used to calibrate and test the SDMs, respectively (Franklin & Miller, 2010). We repeated this process 10 times for each species, and each SDM, and evaluated the model performance based on true skill statistics (TSS) (Allouche et al., 2006). Models with TSS &gt; 0.5 were used to project species distributions under current and future conditions (Thuiller et al., 2009). As different SDM algorithms have different performance, an ensemble method was used to generate the predicted suitability of each species. Specifically, we used the median of all model predictions. We then further converted the predicted habitat suitability of species into binary maps with the threshold maximizing the TSS (Zhang et al., 2017). SDMs often tend to overestimate the current and future potential ranges of species since dispersal limitations are not considered. Following previous studies (e.g. Kremen et al., 2008), we used a 200 km buffer minimum convex polygon (MCP) around the current distribution range of a species to clip its predicted ranges. Dispersal affects species' ability to track favorable habitats under global changes, and should be considered in the evaluations of global change impacts on species distributions (Acevedo et al., 2020). Here, we assumed three dispersal scenarios: 1) no dispersal; 2) full dispersal; and 3) 20 km/decade dispersal beyond predicted current ranges (Chen et al., 2011). The third scenario was conducted similarly using a 200 km (i.e. 10 decade * 20km/decade) buffer MCP around the predicted current distribution range of each species. All analyses were performed using the biomod2 package in R (R Core Team, 2019). <b>Identification of threatened species </b> We identified species that would be threatened in the future under climate and dispersal scenarios based on the relative loss of their suitable areas (LSA). Following similar logic to criterion A3 of the IUCN Red List (IUCN, 2012), we classified all species into five categories: 1) extinct (EX, LSA = 100%); 2) critically endangered (CE, 80%≤LSA &lt; 100%); 3) endangered (EN, 50% &lt;= LSA &lt; 80%); 4) vulnerable (VU, 30% &lt;= LSA &lt; 50%) and 5) least concerned (LC, LSA &lt; 30%). Species with LSA≥30% (i.e. EX+CE+EN+VU) were recognized as threatened species. <b>Phylogenetic signal and PD of threatened species </b> A dated genus-level phylogeny of seed plants in China was obtained from Lu et al. (2018). To obtain the species-level phylogeny covering all 8732 woody species included in our study, we first inserted all the studied species into the genus-level phylogeny as basal polytomies of their corresponding genera, and then randomly resolved the polytomies within each genus using the birth-death bifurcation model. We repeated this process 100 times, and generated 100 random phylogenetic trees for the subsequent analyses. To evaluate whether threatened species were randomly distributed across the phylogeny, we tested the phylogenetic signal in the conservation status of species using <i>D</i>-statistic calculated by the "phylo.d" function of the <i>caper</i> package (Fritz & Purvis, 2010). Here, threatened species were coded as "1" and non-threatened as "0". <i>D</i>-statistic normally ranges from 0 under the Brownian evolution to 1 under random distribution of threatened species across the phylogeny. The significance of <i>D</i>-statistic was tested by comparing the observed <i>D-</i>statistic with the frequency distributions of the simulated values under 1) phylogenetic randomness and 2) Brownian evolution, which was generated from 999 randomizations for each phylogeny. The current and future geographical patterns in SR and PD of threatened species were estimated. Here, Faith’s PD was used (Faith, 1992). PD was calculated using each of the 100 randomly-resolved phylogenies, and the average across these phylogenies was finally used. Since Faith's PD is normally positively correlated with SR, we conducted a generalized additive model (GAM) between PD and SR, and extracted its residuals (hereafter termed residual PD) to represent the net contribution of phylogenetic relationships to PD after SR is accounted for. Positive (or negative) residual PD suggests higher (or lower) PD than expected given the level of SR (Thuiller et al., 2011). <b>Protected area network analyses</b> To assess the effectiveness of the PAs under future climate and land-cover changes, we first estimated the total number and the percentage of PD of all threatened species protected by the current PA network. Second, we calculated the changes in SR (∆SR) and PD (∆PD) of the threatened species within each PA, respectively. The frequency distributions of ∆SR or ∆PD were separately estimated for PAs with ∆SR or ∆PD<i>&gt;0</i> (group I) and ∆SR or ∆PD<i>&lt;0</i> (group II). The PAs with ∆SR or ∆PD<i> </i>&gt; the 50% quantile of group I was considered to be "effective", while those with ∆SR or ∆PD&lt; the 50% quantile of group II was considered to be "ineffective". The remaining PAs of both groups, together with those with ∆SR or ∆PD<i>=0</i> were considered to be "stable". The five PA characteristics were compared between the effective and ineffective PAs using the Kruskal-Wallis test. Moreover, to evaluate the relative importance of the five PA characteristics in explaining PA effectiveness, we conducted Random Forest analysis with PD or SR changes as the response variable, and the five PA variables as the predictors. To identify the conservation priorities under future climate and land-cover changes, we identified the grid cells 1) not covered by the current PAs, and 2) with the highest SR or PD in the future. Following the "top 5% richness algorithm" (Grenyer et al., 2006; Pollock et al., 2017), the top 5% of the grid cells with the highest SR or PD were considered as areas of conservation priorities (Priority<sub>5%</sub>). As additional PAs covering 15% of land area in China are required to achieve the post-2020 conservation target, we also used a 15% cut-off. To facilitate the conservation planning, we further overlaid the Priority<sub>5%</sub><sub> </sub>with the county boundaries in China. A county should be given high conservation priorities in the future if: 1) more than 50% of its area was covered by Priority<sub>5%</sub>; or 2) the sum of the overlapping area was larger than the median area of all counties (i.e. 2081 km<sup>2</sup>).