Spatial map of Chinese protected areas
The database of protected areas (PA) distribution in China and a digitized spatial map thereof were compiled from Zhao et al.40 and Zhang et al.41. In total, we obtained the information of 2622 protected areas in China, which also included marine reserves. In order to evaluate the representation of terrestrial protected areas, we excluded marine reserves from our analyses. We also excluded Taiwan because we did not have the spatial distribution data for nature reserves in Taiwan. Finally, we had the boundary information of 2572 protected areas covering about 15.2% land area in China.
Species’ range maps
Range maps of threatened vertebrates (birds, mammals, amphibians and reptiles) were obtained from the IUCN’s Red List42. Distribution data of threatened plants were compiled from Flora of China, Atlas of woody plants in China, provincial and local floras, checklists of nature reserves, various inventory reports across China and peer-reviewed papers. We obtained the information for critically endangered (CR), endangered (EN) and vulnerable (VU) species. The conservation status of vertebrates was obtained from IUCN Red List42, while that of plants was obtained from Qin et al.43. In total, we obtained the distribution information of 103 birds, 86 mammals, 134 amphibians, 50 reptiles, and 2983 plants in China (see Supplementary Data 1). We estimated the number of species in each PA by overlaying the map of PA with the species’ range maps in ArcGIS 10.2 (ESRI, Redlands, CA). In order to validate the distribution of species, we further verified the presence of species in respective PA by checking their inventory reports.
Human footprint data
In order to measure the extent of human pressure on the protected areas, we obtained the most comprehensive global map of human pressure i.e., human footprint (HFP) from https://wcshumanfootprint.org. The human footprint measures the cumulative impact of direct pressures on environment from human activities and is based on data from built environments, agricultural lands, pasture lands, human population density, night-time lights, railways, roads and navigable waterways44. It is one of the most complete and finest terrestrial datasets on cumulative human pressure on the environment. The human footprint maps of two time periods (1993 and 2009) are available at present. We downloaded the maps of both time periods at the spatial resolution of 1 km × 1 km to quantify the change in human pressure within Chinese protected areas over a 16-year period. It should, however, be noted that any point estimate of the change in HFP might include errors due to the resolution and reliability of the component layers. For example, one of the components of HFP is the night-time lights, which changed over time from incandescent to mercury vapor to light emitting diode. This means that the change in night light is due to more than development. As a result, the systemic bias in regional economy could likely cause low HFP in wealthy as compared to rural areas. While this issue does not invalidate our analyses, such comparisons should be applied with caution.
Climate data
In order to represent the climatic conditions of the past and the present, we obtained 20 years climate data comprising monthly minimum temperature, monthly maximum temperature and monthly precipitation for two time periods (past: 1961–1970 and present: 2010–2019). We used 10 years window for each time period to capture the variability in climatic conditions in order to prevent over- or underestimation of the past and present climate. In total, we obtained 12 months × 10 years × 3 variables × 2 time periods = 720 global raster layers from the Climate Research Unit (CRU TS v. 4.04) database (http://www.cru.uea.ac.uk/data) at the spatial resolution of 0.5° × 0.5°45. We then calculated the monthly mean values for the three variables for each time period separately. From these monthly mean values, we estimated mean annual temperature (MAT), mean temperature of warmest quarter (MTWQ), mean temperature of coldest quarter (MTCQ), mean annual precipitation (MAP), precipitation of driest quarter (PDQ) and precipitation of wettest quarter (PWQ) for each time period using biovars function in the R package ‘dismo’46. We then estimated the average value of climate variables in each PA using zonal.stats function in the R package ‘spatialEco’47. In order to reduce dimensionality and collinearity of the 6 climate variables, we performed principal component analysis (PCA) using prcomp function in R v4.0.248. Following Carroll et al.37, we used climate data for both past and present based on the first 2 PCA axes, which explained 89.4–89.8 % of the variance (Supplementary Tables 1-2). Additionally, we also estimated the change in mean annual temperature to identify PAs that have experienced climate warming higher than the Paris Agreement threshold of 1.5 °C. All the analyses were performed in R version 4.0.248.
Vulnerability mapping
We calculated three indices of vulnerability within each protected area: (i) species vulnerability, (ii) anthropogenic vulnerability, and (iii) climate vulnerability. In order to measure species vulnerability, we first assigned numerical value to each IUCN threat category using a geometric progression49,50. We gave scores of 2, 4, and 8 to species belonging to categories VU, EN, and CR, respectively. We, then, summed the score of all species in each PA and standardized the value to the range of 0–1 using minimum–maximum normalization. We performed these steps separately for birds, mammals, amphibians, reptiles and plants to calculate the vulnerability score of each group. We also calculated the cumulative score by combining the total scores of all five groups. Values close to 0 indicated low species vulnerability and values close to 1 indicated high species vulnerability. Although species diversity is highly correlated with the species vulnerability metric used herein (Pearson’s correlation coefficient = 0.94 − 0.99, p < 0.001), the species vulnerability metric is more useful from a conservation perspective because it offers trade-off between proactive (i.e., protecting high biodiversity) and reactive (i.e., protecting threatened species) conservation practices. While proactive approaches are useful to keep the ‘common species common’, reactive approaches are necessary to prevent the extinction of already threatened species51. Our species vulnerability metric captures both species diversity and the threat level of a species. Therefore, high species vulnerability may either mean high diversity of threatened species (CR, EN, VU) or higher number of highly threatened species, such as “Critically Endangered” or both.
In order to measure anthropogenic vulnerability, we first estimated the average HFP value in each PA using the zonal.stats function in the R package ‘spatialEco’47. We did this separately for 1993 and 2009. We, then, subtracted the average HFP value in 1993 with the average HFP value in 2009 to estimate the change in HFP (ΔHFP) in each PA. Finally, we standardized the ΔHFP value to the range of 0–1 using minimum–maximum normalization. Value close to 1 indicate high anthropogenic vulnerability (i.e., greater increase in human pressure since the past), while value close to 0 indicate low anthropogenic vulnerability.
In order to measure climate vulnerability, we estimated the change in climate in each PA by calculating the Euclidean distance between the past and present distributions in a 2-dimensional climate space (i.e., PCA1 and PCA2 axes). Larger distance between the past and present distributions means higher deviation from original climate and hence represent high climate vulnerability therein. We standardized the climate change value for each PA to the range of 0–1 using minimum–maximum normalization. Value close to 1 indicate high climate vulnerability (i.e., greater change in climate since the past), while value close to 0 indicate low climate vulnerability.
We computed the Pearson’s correlation coefficients between climate, anthropogenic and species vulnerabilities to evaluate concordance between these three measures. In order to account for spatial autocorrelation, we computed the significance levels of all the correlation coefficients using modified t-test52. We performed these analyses in Spatial Analysis in Macroecology (SAM) version 4.053.
Classification of protected areas
We classified Chinese PA based on three vulnerability indices explained above. We first identified (i) PA with significantly high species vulnerability (species vulnerability hotspots), (ii) PA with significantly high anthropogenic vulnerability (anthropogenic vulnerability hotspots), and (iii) PA with significantly high climate vulnerability (climate vulnerability hotspots). We used a null model approach to assess significance and identify the hotspots of each vulnerability category. We first randomized respective vulnerability scores of each PA. In doing so, we kept a track of where the observed value was greater than the randomized value. We repeated this process 1000 times. The significance of each PA was then assessed by evaluating the rank relative position of the original values against those of the random realizations. For each vulnerability category, we applied a one-tailed test to identify significantly high values. If the observed value was within the highest 5% of the null distribution for that PA, it was identified as significantly high (p < 0.05).
We overlaid the spatial maps of these three hotspots in ArcGIS 10.2 and counted the number of overlaps. The PAs with three overlaps, i.e., those with significantly high value for all three vulnerability indices, were classified as “Level 1” PA. The PAs with two overlaps, i.e., those with significantly high value for two out of three vulnerability indices, were classified as “Level 2” PA. The PAs with no overlap, i.e. those with significantly high value for one out of three vulnerability indices, were classified as “Level 3” PA. The protected areas that did not have significantly high values for any of these three vulnerability indices were classified as “Level 4” PA. Level 4 PAs have comparatively low scores for species, climate and anthropogenic vulnerabilities.
Intensity of climate and human footprint changes in species vulnerability hotspots and coldspots
In order to evaluate the intensity of climate change and human pressure change on the vulnerability hotspots of individual species groups (i.e. birds, mammals, amphibians, reptiles, and plants), we identified the species vulnerability hotspots of each group that, respectively, overlapped with the climate vulnerability hotspots and the anthropogenic vulnerability hotspots. We counted the total number of overlaps to identify (i) the mutual hotspots of species and climate vulnerabilities and (ii) the mutual hotspots of species and anthropogenic vulnerabilities.
We also compared the intensity of climate change and human footprint change (ΔHFP) in the species vulnerability hotspots and coldspots to evaluate if species vulnerability hotspots and coldspots have been under similar influence of climate and human pressure changes since the past. We defined species vulnerability hotspots as the PAs with significantly high species vulnerability and species vulnerability coldspots as the PAs with significantly low species vulnerability. We applied similar randomization approach (see above) to identify hotspots and coldspots of species vulnerability However, we applied a two-tailed test to recover both significantly high and low values. If the observed value was within the highest 2.5% of the null distribution for that PA, it was identified as significantly high (hotspots) and if the observed value was within the lowest 2.5% of the null distribution for that PA, it was identified as significantly low (coldspots).
We compared the differences in climate change and human footprint change between species vulnerability hotspots and coldspots by bootstrapping the samples 1000 times. We calculated the difference in means between hotspots and coldspots in each repetition. Finally, we used this bootstrapped means to create the null distribution and find statistically significant difference in climate change and human footprint change values between species vulnerability hotspots and coldspots.
In order to evaluate the magnitude of future climate warming in each PA, we obtained the future climate data from the worldclim database (http://www.worldclim.org/) at the spatial resolution of 30 arc seconds. We obtained the data simulated from three general circulation models: BCC-CM1-1, CCSM4, and MIROC-ESM and chose the lowest emission scenario (i.e., rcp 2.6) for the year 2070. We calculated the average future mean annual temperature (MAT) from these three data layers using the raster calculator tool in ArcGIS 10.2 and subsequently calculated the average MAT for each PA using the zonal.stats function in the R package ‘spatialEco’47. To enable comparison, we also obtained the data for the current climate (MAT) from the worldclim database at the spatial resolution of 30 arc seconds and similarly calculated the average value for each PA. We estimated the change in temperature in each PA by subtracting the future MAT with the current MAT.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
