Conservation of Chinese Theaceae species under future climate and land use changes

Climate and land use changes are two major pervasive and growing global causes of rapid changes in the distribution patterns of biodiversity, challenging the future effectiveness of protected areas (PAs), which were mainly designed based on a static view of biodiversity. Therefore, evaluating the effectiveness of protected areas for protecting the species threatened by climate and land use change is critical for future biodiversity conservation.


| INTRODUC TI ON
Earth's biodiversity is experiencing rapid loss globally, which itself can be considered as an important global change (Walker & Steffen, 1996). Human-induced land use change, such as agriculture and urban expansion, is perhaps the most prominent factor driving range contractions and species extinctions (Newbold et al., 2015;Sala et al., 2000). At the same time, anthropogenic climate change is having an increasingly negative effect on biodiversity (Garcia et al., 2014;Tang et al., 2022;Thomas et al., 2004) and may pose more critical threats to biodiversity than land use changes (Préau et al., 2022;Tang & Zhao, 2022a). Moreover, biodiversity in a specific area may be particularly affected by the combined effect of these two drivers, which could amplify the negative impact of each stressor in isolation (Newbold, 2018;Prestele et al., 2021;Tang & Zhao, 2022b). Owing to increasingly intensified global change, such combined effects will further threaten biodiversity in the coming decades (Newbold, 2018;Radinger et al., 2016;Zhang et al., 2017).
Therefore, to avoid dramatic loss of biodiversity and set effective conservation priorities, there is an urgent need to identify which species will be threatened under the combination of future climate change and land use change and which regions tend to have the most significant changes in the distribution patterns of threatened species (Peng et al., 2022).
Over the last decades, the vast majority of studies have attempted to assess how climate and land use change will affect the future distribution of species, mainly focusing on species' range shifts under future climate and/or land use conditions (Chen et al., 2011;Mantyka-pringle et al., 2012;Newbold, 2018). This is mainly because understanding the redistribution of species can help to inform future conservation planning Tang et al., 2022). Those species experiencing severe range contraction are usually identified as threatened species and have been widely used as one of the most effective surrogates for identifying priority conservation areas and for future protected areas (PAs) planning (Peng et al., 2022;Zhang et al., 2015). Specifically, species might shift their habitat range out of the PAs to track their suitable habitat due to future environmental change. As a result, species of high conservation concern might not be well protected by current PAs under future environmental conditions, thus leading to conservation gaps Hoffmann et al., 2019). Moreover, although some species might gain suitable habitat in some PAs under climate change, the conservation and management effectiveness of these PAs may also be reduced because existing PAs were designed to protect current biodiversity (Li, 2019). Therefore, it is crucial to evaluate the effectiveness of PAs for conservating species threatened by future climate and land use change, as this can significantly help to improve the future effectiveness of current PA networks and reduce the extinction risk of threatened species in a changing world (Araújo et al., 2004;Peng et al., 2022).
In China, Theaceae is an economically crucial angiosperm family, including about 12 genera and 274 species (Shen et al., 2016). These species are mostly found in southern China's tropical and subtropical moist lowlands and mountains (Grote & Dilcher, 1989;IUCN, 2017).
Despite their importance for ecosystem function and services (e.g. eating, drinking, ornamental, wood utilization, skin care, beauty, and fertilizer), over 33% of Theaceae species are currently threatened, primarily due to human-induced land-use changes (e.g. urbanization; IUCN, 2017). Moreover, previous studies have suggested that climate change likely exacerbates these threats to the Theaceae species (Zhang et al., 2020). Therefore, these Theaceae species will face more severe challenges in the near future (Rao et al., 2018;Tang & Zhao, 2022b). To arrest the decline of these Theaceae species and other higher plant and vertebrate species, China has established over 2700 protected areas in the past decades (Xu et al., 2017). Despite these efforts, few studies have assessed the effectiveness of current PAs in protecting threatened Theaceae species under future climate and land use change.
To fill these knowledge gaps, we aim to (a) identify Theaceae species that will be threatened by future climate and land use conditions,

| Species occurrence data
This study was conducted in China, which harbours ~75% (12 genera) of the total genera and ~46% of the total species (274 species) of Theaceae worldwide. The species occurrence records for China were obtained from Zhang et al. (2020), in which a total of 12,055 occurrence records for 200 Theaceae species were collected. To match our climate and land use data (see below), we removed duplicate records in each 1 km × 1 km gird cell across the study area for each species. Finally, 11,992 occurrence records for the 200 Theaceae species (occurrence records per species: mean = 60; median = 28; range = 5-574) were included for modelling the distribution of Theaceae species (Figure 1a).

| Climate, land use and protected area data
The current  19 bioclimatic variables were obtained from the CHELSA dataset at a 30 arc-second resolution (~1 km; Karger et al., 2017). The same 19 bioclimatic variables for the 2070s (average for 2061-2080) were obtained from four global circulation models (GCMs; CanESM5, CNRM-CM6-1, IPSL-CM6A-LR, MIROC-ES2L) and obtained from WorldClim dataset (Fick & Hijmans, 2017) with 30 arc-second resolution (~1 km). We selected these four GCMs because they can produce the greatest commonality among projections and therefore can allow us to explore the uncertainties of model projections (Fajardo et al., 2020;Thuiller et al., 2005). For each GCM, three representative concentration pathways (i.e. RCPs 2.6, 4.5 and 8.5), representing an optimistic, an intermediate and a pessimistic climate scenario, respectively, were considered. Finally, all climate variable layers were resampled to 1 × 1 km. Current (2010) and future (under the same three RCP scenarios by the 2070s) 10 land use variables were obtained from F I G U R E 1 The occurrence records (a) and current potential species richness (PSR; b) of the 135 Chinese Theaceae species used in this study.
the FROM-GLC datasets, with a spatial resolution of 1 × 1 km (Li et al., 2016). Land use variables included the proportion of area in each grid cell covered by 10 different land-use types: bare land, cropland, forest, grassland, impervious, shrubland, snow/ice, urban green spaces, water and wetland. To assess current and future habitat suitability, we selected 10 variables from the above 29 environmental variables. The selected 10 variables had low multi-collinearity (all variables with Pearson correlations |r| < .7) and had the greatest biological relevance to Theaceae species (annual mean temperature, Isothermality, annual temperature range, annual precipitation, precipitation seasonality, the proportion of area in each grid cell covered by cropland, forest, grassland, shrubland and waters).
We obtained spatial Polygon data for 1028 protected areas in China from Tang and Zhao (2022a). All polygons were converted to raster at a 10 km resolution by using the 'Polygon to Raster' tool in ArcGIS 10.2. By using the 'Extract Values To Points' tool, we then assigned each grid in the species distribution raster to the protected area (PA) and non-protected area (non-PA) category according to the protected status of each grid.

| Species distribution modelling
For 159 species with ≥10 unique occurrence records at a 1 × 1 km grid cell resolution, the maximum entropy model (MaxEnt ;Phillips et al., 2006) was used to project consensus distribution maps because it performs better with small sample sizes relative to other modelling methods (Elith et al., 2006;Kumar & Stohlgren, 2009;Pearson et al., 2007). To account for potential sampling bias in species distribution data, we selected 10,000 background points using target-group sampling (Phillips et al., 2009). Specifically, these background points were sampled from areas where other Theaceae species have been recorded. This approach is more objective than random sampling and thus can produce more accurate projections (Mateo et al., 2010). Furthermore, to avoid overfitting, we used the 'ENMevaluate' function in the R package 'ENMeval' (Version 2.0.3; Kass et al., 2021) to tune MaxEnt models and select the best parameters (i.e. regularization multipliers and feature class) for each species. Based on the above selected parameters for each species, we partitioned the above data into an 80%-20% split where 80% of the data were used for training model, and the remaining 20% of the data were used for assessing model performance using the area under the relative operating characteristic curve (AUC). We repeated this procedure ten times to create predictions independent of the training data. Accordingly, 10 current and 120 future (10 replicates × 4 GCMs × 3 RCPs) modelled distributions were produced for each species. For each species, we calculated the weighted mean habitat suitability maps of the 10 replicate models to produce a single robust ensemble forecast for each GCM and for each RCP scenario (Araújo & New, 2007). We then converted these probabilistic maps to binary suitable and unsuitable habitat maps using the threshold that maximized the sum of sensitivity and specificity . Finally, we used unweighted ensemble means across the four GCMs to generate consensus species distribution maps for each RCP scenario.
For 41 species with 5-9 unique occurrence records in a 1 km × 1 km grid, distribution was modelled using a range-bagging algorithm, which is a machine learning method that uses an ensemble of convex hulls created from a reduced set of niche parameters and refined via bootstrap aggregation, or bagging (Andrew et al., 2021;Drake, 2015). Specifically, the range of each species was estimated by the above ten environmental variables at occupied grids, dimension of ranges to bag and the number of votes and then projected them into environmental space under current and future periods. By using the same method above, for each species, we obtained the binary suitable and unsuitable habitat maps for each GCM and each RCP scenario and then obtained consensus species distribution maps for each RCP scenario.
To ensure the reliability of this study, we further removed species with AUC values <0.7. Finally, 135 Theaceae species were retained in the subsequent analyses.

| Statistical analysis
To identify species that will be threatened under future climate and land use conditions, we calculated four types of range changes in each grid cell between the current and future periods for each species using the 'BIOMOD_RangeSize' function provided in the biomod2 package: To identify the priority areas for protection under future climate and land use changes, following Xie et al. (2022), we calculated a 95% confidence interval for the species richness change for each RCP scenario and the 5% of the grid cells outside the confidence interval were considered as areas of conservation priorities. Furthermore, we classified all grid cells in the above-identified conservation priority areas into two categories: (i) 'areas worth exploring', representing the grid cells with increased species richness, and (ii) 'areas needing attention', representing the grid cells with decreased species richness. To identify conservation gaps in current PAs, we overlaid the conservation priority areas with the current protected area boundaries in China and calculated percentage of the above two types of conservation priority areas covered by current PAs.

| Threaten status of Chinese Theaceae species
Under current environmental conditions, the mean potential suitable habitat area of the 135 Theaceae species in China was 784,593 ± 712,035 km 2 (Table S1). The Camelli parvimuricata had the largest range, covering an area of 3,999,348 km 2 , while the Ternstroemia simaoensis had the narrowest range, covering only 39,109 km 2 (Table S1). Our projections of future habitat suitability for the 135 Theaceae species suggest that future climate and land use change will threaten many Theaceae species (Table S2; Figure 2).
Specifically, by assuming unlimited dispersal, the percentage of the 135 Chinese Theaceae species that were predicted to be vulnerable or critically endangered was 36.30%, 38.52% and 40.00% for RCPs 2.6, 4.5 and 8.5 by the 2070s, respectively. As expected, the results under zero dispersal are more severe than those under unlimited dispersal, despite having a similar trend under different threatening status (Table S2; Figure 2).
In addition, small regions of other provinces of China are also predicted to harbour many threatened Theaceae species.
By the 2070s, the richness distribution patterns of the threatened Theaceae species were similar under the current conditions. However, under RCPs 2.6, 4.5 and 8.5, each grid cell of the entire study area was predicted to lose, on average, 2.9 ± 3.7, 3.1 ± 4.1 and 3.3 ± 4.4 species, respectively, whereas only gaining 0.2 ± 0.6, 0.2 ± 0.6 and 0.2 ± 0.6 new species, respectively. This is predicted to result in a net reduction in richness of threatened Theaceae species in China when compared to the current period. Specifically, mean species richness was predicted to decline by 2.7 species, 2.9 species and 3.1 species under RCPs 2.6, 4.5 and 8.5, respectively. The     (Figure 4). More alarmingly, no more than 8% of the 'areas needing attention' (7.913% under RCP2.6, 7.074% under RCP4.5 and 6.913% under RCP8.5) would be covered by the current PAs (Figure 4).

| DISCUSS ION
Recent scientific consensus suggests that the redistributions of species due to climate and land use change occur globally, which would potentially change the richness distribution patterns of threatened species and challenge the conservation effectiveness of current protected areas Hoffmann et al., 2019). In this study, we provide the first study that identifies the number and distributions of Theaceae species threatened by future climate and land use changes and assess the effectiveness of current protected areas for protecting these threatened Theaceae species. Our results suggest that at least 36% of Theaceae species would be threatened and suggest net declines in the richness of these threatened Theaceae species by the 2070s. Further, we found that future climate and land use change will reduce the effectiveness of the current PA network for protecting these threatened Theaceae species. These findings highlight that judicious management and expansion of the PA system are necessary to ensure the resilience of these threatened Theaceae species into the future.
Due to species-specific habitat requirements and physiological tolerances, species may respond to environmental change differently (Levinsky et al., 2007;Prestele et al., 2021;Thuiller et al., 2011). Similarly, our results showed considerable variation in projected changes in suitable habitat among the 135 Theaceae species under future climate and land use change, ranging from −99.7% to 738.0% under unlimited dispersal and from −99.7% to 0.0% under zero dispersal (Table S1). Nonetheless, our results also revealed that at least 49 Theaceae species are expected to be threatened under future climate and land use conditions and found that the number of threatened Theaceae species increased from RCP2.6 to RCP8.5 by 2070s. This finding supports previous studies suggesting that a wide range of taxa is predicted to face greater extinction risk and population decline as climate change intensifies (Tang et al., 2022;Wang et al., 2022). Further, we found that under the zero dispersal scenario, more species would be threatened when compared to the unlimited dispersal scenarios. Consistent with previous studies (Gouveia et al., 2016;Ma et al., 2021;Senior et al., 2019), this finding suggests that habitat connectivity and habitat corridors are increasingly essential to ensure that species can immigrate into newly suitable habitats if dispersal is sufficient (Barbet-Massin et al., 2012;Ma et al., 2021).
However, the unrealistic dispersal modes (i.e. unlimited and zero dispersal) we used to estimate the future distributions of the Theaceae TA B L E 2 The mean (±SD) values of the 10 selected climate and land use variables used for species distribution modelling of Theaceae species in the identified conservation priority areas. Abbreviations: ANA, 'Areas need attention'; AWE, 'Areas worth exploring'; BIO1, annual mean temperature; BIO12, total annual precipitation; BIO15, precipitation seasonality; BIO3, temperature seasonality; BIO7, temperature annual range; CL, the proportion of area covered by cropland; FL, the proportion of area covered by forest; GL, the proportion of area covered by grassland; SL, the proportion of area covered by shrubland; WL, the proportion of area covered by waters.

Variables
species could potentially underestimate (under unlimited dispersal) or overestimate (under zero dispersal) the extinction risks for these Theaceae species (Engler & Guisan, 2009;Zanatta et al., 2020), as the dispersal ability of these species will undoubtedly fall somewhere between zero and unlimited. In addition, the SDMs only capture the realized environmental niche of species, ignoring the acclimation and adaptation of species to future climate and land use change, which may lead to the overestimation of species' vulnerability to future global change (Bush et al., 2016(Bush et al., , 2018 which also found particularly significant diversity changes in mountainous areas (Levinsky et al., 2007).
Our projections also suggested that the major conservation priority areas for protecting threatened Theaceae species are concentrated in southern China, of which the 'areas needing attention' (with decreased species richness) are mostly located in Sichuan, Tibet, Yunnan, Guizhou, Yunnan, Fujian, Zhejiang and Taiwan provinces, while the 'areas worth exploring' (with increased species richness) are mostly located in Xinjiang, Shaanxi, Fujian, Jiangxi provinces. However, we found that most regions in the identified conservation priority areas are highly fragmented. Moreover, our conservation gap analysis indicated that these conservation priority areas are far from adequately covered by the current protected area network in China. Specifically, no more than 21% of the 'areas worth exploring' are covered by current protected areas, and worse still, no more than 8% of the 'areas needing attention' are situated within the current protected area network. These conservation gaps, together with habitat fragmentation, pose severe challenges for the persistence of these threatened Theaceae species. Therefore, it is urgently needed to not only establish new protected areas, but also construct regional ecological corridors among protected areas to increase habitat connectivity to conserve the threatened Theaceae species (Tucker et al., 2018).

ACK N O WLE D G E M ENTS
We thank Dr Francesco Maria Sabatini, Dr. Ronald R. Swaisgood and two anonymous reviewers for their comments on the manuscript. This study was supported by the National Natural Science Foundation of China (32100401).

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

PEER R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-revie w/10.1111/ ddi.13744.

DATA AVA I L A B I L I T Y S TAT E M E N T
All occurrence records data of the Theaceae species, climate and land use data that used for modelling are openly available in Dryad