The spatial and temporal evolution of habitat quality and driving factors in nature reserves: a case study of 33 forest ecosystem reserves in Guizhou Province

Study area

According to the list of nature reserves in Guizhou Province published by the Guizhou Forestry Bureau, as of 2018, there are 11 national-level and 89 local-level nature reserves in Guizhou. The majority of these reserves are forest ecosystem types. Considering the integrity and stability of ecosystems, we excluded protected areas that were too small to sustain complete forest ecosystems. We prioritized larger protected areas that exhibited higher biodiversity, contained habitats for endangered species, or possessed significant ecological functions, this study focuses on 33 nature reserves of forest ecosystem in Guizhou Province (Table 1 and Fig. 1) to analyze changes in habitat quality from 2000 to 2020. Each reserve will be assigned a code to facilitate the analysis.

The selected 33 nature reserves in this study are primarily located in the northern and southeastern regions of Guizhou (Fig. 1). These areas face the simultaneous challenges of economic development and ecological protection, leading to a particularly pronounced conflict between environmental conservation and poverty alleviation efforts. Additionally, these regions have recently become major tourism hotspots. For example, the Xijiang Miao Village and the Fanjing Mountain scenic area are situated adjacent to the Leigong Mountain National Nature Reserve and the Fanjing Mountain National Nature Reserve, respectively, intensifying the conflict between economic development and ecological protection.

Data resources

The data used in this study primarily includes GIS spatial data as inputs for the model. The threat factors consist of two types: cropland and construction land. Moreover, the land use data for the years 2000, 2010, and 2020 were obtained from the China Land Cover Annual Dataset (CLCD) (https://doi.org/10.5281/zenodo.5816591), which provides raster data with a spatial resolution of 30 m. Furthermore, in this study, the land use types in the study area were categorized into seven types: forest, farmland, shrub, grassland, water bodies, impervious surfaces, and wasteland.

Additional data, including the digital elevation model (DEM), average annual temperature, average annual precipitation, population density, and gross domestic product (GDP), were obtained from the Data Center for Resource and Environmental Sciences of the Chinese Academy of Sciences (http://www.resdc.cn). Slope and elevation values were extracted from the DEM using ArcGIS 10.8. All data were sourced from official platforms as detailed in Table 2. Moreover, Fig. 2 presents a schematic representation of our research methodology.

The landscape pattern index

Markov models serve as the foundation for the Land use transfer matrix (Jawarneh et al., 2024). In a recent study, changes in land-use types between the start and end of the study period were determined, and the original sources of each land-use type at the end of the study period were identified using the Markov model (Feudjio Fogang et al., 2023). The Markov model is, therefore, suitable for the assessment of land-use changes (Ait El Haj, Ouadif & Akhssas, 2023; Congalton, 1991; Dietzel & Clarke, 2006). The formula for land-use change is as foll (1)${\displaystyle S_{ij}=\left\{\begin{array}{cccc}\hfill S_{11}\hfill & \hfill S_{12}\hfill & \hfill \cdots \hfill & \hfill S_{1n}\hfill \\ \hfill S_{21}\hfill & \hfill S_{22}\hfill & \hfill \cdots \hfill & \hfill S_{2n}\hfill \\ \hfill S_{n1}\hfill & \hfill S_{n2}\hfill & \hfill \cdots \hfill & \hfill S_{nn}\hfill \end{array}\right\}}$ where S_ij represents the area of land initially classified as type i and converted to type j by the end of the study, and n denotes the total number of land-cover types.

The dynamics of integrated land use reflect the overall speed and magnitude of land-use changes within the study area. A higher value of integrated land-use dynamics indicates a faster rate of integrated land-use change in the region. (2)${\displaystyle L_c=\frac{\sum _{i=1}^n\Delta {LU}_{i-j}}{2\sum _{i=1}^n{LU}_I}\times \frac{1}{T}\times 100\text{\%}}$ where L_C represents the integrated land-use dynamics, LU_i denotes the area of category i land-use type at the beginning of the study, and ΔLU_i−j denotes the absolute value of the area of the i land-use type that was converted to a non-i land-use type during the study period.

The InVEST HQ model and spatial autocorrelation analysis

Spatial autocorrelation analysis

Spatial autocorrelation was analyzed using the global Moran’s I index to assess whether the distribution of habitat quality in protected areas exhibited a clustered pattern (Jinrui et al., 2019). Analyses of hotspots and coldspots aim to identify and quantify the degree of clustering and the statistical significance in geospatial data. Hotspots and coldspots represent clusters of high and low values, respectively, which are visualized by this method. Spatial distribution patterns and aggregation characteristics of the data can thereby be revealed (Song et al., 2020b; Xiao-Han & Xiu-Juan, 2024). The formulae used are as follows: (5)${\displaystyle I=\frac{n\sum _{i=1}^n\sum _{j=1}^n{w}_{i\mathrm{j}}\left(x_{\mathrm{i}i}\bar{X}\right)\left(x_j\bar{X}\right)}{\sum _{i=1}^n{\left(x_i\bar{X}\right)}^2\left(\sum _{i=1}^{\mathrm{n}}\sum _{j=1}^{\mathrm{n}}{w}_{ij}\right)}}$ (6)${\displaystyle G_i^{\ast }=\frac{\sum _{j=1}^n{w}_{ij}{x}_j-\sum _{j=1}^n{w}_{ij}}{S\sqrt{\frac{\sum _{j=1}^n{w}_{ij}{x}_j-{\left(\sum _{j=1}^n{w}_{ij}\right)}^2}{n=1}}}}$ (7)${\displaystyle \bar{X}=\frac{\sum _{j=1}^n{x}_j}{n}}$ (8)${\displaystyle S=\sqrt{\frac{\sum _{j=1}^n{x}_j^2}{n}}-{\overline{X}}^2}$

where I refers to Moran’s index; $G_i^{\ast }$ is the hotspot index; w_ij represents the spatial weight between the i and j spatial cells; x_i and x_jdenote the values of the i and j cells, respectively; $\bar{X}$ represents the mean value of the cells; and n is the total number of cells in the study area. Moran’s I index ranges between −1 and 1, with larger values indicating a stronger spatial correlation, values less than 0 indicating a negative correlation, and a value of 0 indicating a random distribution.

Optimal parameters-based geographical detector

In this study, seven indicators were selected from both natural environmental and socioeconomic dimensions to investigate the main factors influencing the habitat quality of forest ecosystem nature reserves in Guizhou Province, and the evolutionary mechanisms, taking into account the specific characteristics of these reserves (Table 5). Topography, a fundamental natural environmental factor, forms the basis for the development and evolution of geographic environments (Wu, Sun & Fan, 2021). In this study, elevation and slope were used to represent topography. Precipitation and temperature, which influence plant growth and development, were represented by average annual precipitation and annual temperature. Land-use intensity reflects effects of human activity and was categorized based on the degree of environmental disturbance: forests, shrubs, grasslands, and water bodies were classified as Category 1 and represent a relatively natural state; agricultural land was classified as Category 2; and impervious surfaces were classified as Category 3. Agricultural land and impervious surfaces typically suggest more intensive human activities, often manifesting as land modification, vegetation destruction, and habitat degradation, all of which can severely affect ecosystems. In contrast, Category 1 areas, such as forests and water bodies, may offer better ecosystem protection. Gross Domestic Product (GDP) serves as a key indicator of the economic status and developmental level of a region. Population density provides information on the distribution of people within a unit area. By selecting GDP and population density, both the degree of human activity and economic transactions (e.g., tourism) affecting nature reserves can be captured. Higher levels of disturbance exert greater pressure on the protected areas, leading to ecosystem degradation, which in turn negatively affects species and habitat quality within these reserves.

Traditional geographic detectors require manually setting up the parameters when discretizing continuous variables, leading to issues of subjectivity and suboptimal discretization. To overcome these limitations, the Optimal Parameters-based Geographic Detector (OPGD) model was used to analyze the drivers of spatial and temporal changes in habitat quality of the forest ecosystem nature reserves in Guizhou Province (Dan-Dan & Yong-liang, 2023; Song et al., 2020b).

First, the explanatory power (q value) of each continuous-type driver is calculated under different classification methods and varying numbers of classifications. The q value ranges from 0 to 1, with a higher q value indicating a stronger explanatory power and a lower q value indicating a weaker explanatory power of the driver for spatial distribution of disaster sites. The OPGD model selects the optimal discretization method and number of breaks for each variable based on the maximum q value (Cen, Zhou & Qiu, 2024), thereby ensuring the highest possible explanatory power. The two-way interaction detector evaluates whether the combined effect of two factors enhances or weakens the explanatory power of the dependent variable Y (Yang et al., 2024). This analysis was conducted using the “GD” package in RStudio.

The interaction detector evaluates whether the combined effect of two influencing factors enhances or reduces the explanatory power of the dependent variable Y. The relationship between the two factors is presented in Table 6.

Land-use change analysis

To investigate the land use type transitions within the nature reserves of forest ecosystem in Guizhou Province over the past two decades, this study conducted a comparative analysis of land use data from two periods: 2000 to 2010 and 2010 to 2020. Corresponding land use transition matrices were constructed (Table S1). The analysis of these continuous time period matrices revealed the patterns of land use evolution in the Guizhou reserves.

The land use transition and chord diagrams (Fig. 3) indicate that forests dominate the study area, covering approximately 90% of the total area. A ranking of land use transitions by area across the two time periods reveals that over the past 20 years, land use transitions in the study area predominantly involved conversions from forest to agricultural land, forest to impervious surfaces, and forest to shrubland. Specifically, from 2000 to 2010, a total of 5,518.89 hectares of land experienced transition, characterized by the predominant conversion of forest to agricultural land and shrubland, with minimal transition to grassland. The increase in impervious surfaces primarily resulted from the conversion of agricultural land. In the subsequent period from 2010 to 2020, shrubland transitioned mainly back to forest, with a smaller portion converting to agricultural land. Forests continued to primarily transition to agricultural land and shrubland, but the outflow was less than the inflow. Most of the impervious surface increase originated from agricultural land, with a smaller contribution from forests, and the rate of increase was higher than during the 2000 to 2010 period.

The overall land use dynamics index rose from 0.06 in 2000 to 0.1 in 2020, indicating an acceleration in the rate of land use change and greater overall land use activity (Table 7). However, the rates of change varied among different land use types, with agricultural land, impervious surfaces, and water bodies exhibiting increased dynamics, while forests, grasslands, and shrublands showed decreased dynamics. The dynamic level of forests rose during the 2000 to 2010 period but declined in the 2010 to 2020 period, with a slowdown in the rate of area reduction. Meanwhile, the dynamic level of impervious surfaces increased from 2000 to 2010 but declined from 2010 to 2020, indicating a reduced rate of area increase. This suggests that during this period, efforts to protect and restore forests were strengthened, which aligns with the aforementioned land use analysis results.

Distributional characteristics of habitat quality in nature reserves

Spatiotemporal variations in habitat quality within nature reserves

Based on land use and threat factor data, this study employed the habitat quality module of the InVEST model to assess habitat quality across 33 nature reserves of forest ecosystem from 2000 to 2020. Spatial data on habitat quality were obtained for the years 2000, 2010, and 2020, with habitat quality values of 0.7089, 0.5400, and 0.4911, respectively, with standard deviations of 0.1539, 0.1635, and 0.1657 (Fig. 4). These results indicate a downward trend in habitat quality, with the rate of decline gradually decreasing over time. This suggests that the spatial heterogeneity of habitat quality within the nature reserves has progressively increased, consistent with the results of the prior land-use dynamic analysis.

The results of the equal interval classification method in ArcGIS 10.8 showed that, habitat quality was categorized into four grades, and a table presenting the proportion of area in each grade was generated (Table S2). Integrating the results from Table S2, it was found that in 2000, habitat quality within the study area was primarily categorized as Highest Level and Higher Level. By 2010, the area classified as Highest Level had decreased substantially, while Higher Level remained nearly unchanged, and Lower Level showed a marked increase. In 2020, the area with low-quality habitat continued to expand, though the rate of decline in Highest Level slowed significantly compared to the previous decade. Additionally, throughout the study period, habitat quality in national-level nature reserves generally remained high, with Higher Level habitat quality consistently dominant. In contrast, habitat quality in local-level nature reserves was relatively average, with a more pronounced rate of decrease in Highest Level habitat area. Overall, the habitat quality in the study area was moderate. However, the distribution of high-, medium-, and low-habitat-quality areas within different levels of the nature reserves was uneven. For example, in 2010, high-quality habitats accounted for only 10.36% of the total area, with the majority (10.13%) concentrated in national-level nature reserves. The only local-level reserves with high-quality habitats were Congjiang County Moon Mountain Nature Reserve (R1), Rongjiang County Moon Mountain Nature Reserve (R2), and Baiqing-Huanglian Nature Reserve in Tongzi County (R6), accounting for 0.23% collectively. The proportion of good-quality habitat was 55.70%, with national and local reserves contributing 28.01% and 27.69%, respectively. Poor and lowest-quality areas accounted for a combined 33.95%, with national and local reserves comprising 10.32% and 23.63%, respectively. From 2000 to 2020, the high-quality-habitat area showed a declining trend, while relatively-high-quality habitats generally increased. Simultaneously, the low- and relatively-low-quality habitats increased annually. This pattern of decreasing high-quality-habitat areas and increasing relatively-high-, relatively-low-, and low-quality-habitat areas aligns with the previously observed declining trend in the habitat quality index. From 2000 to 2020, the habitat quality of the 33 forest ecosystem nature reserves exhibited an overall declining trend, with average decreases of 0.1689 and 0.0489, respectively, and a deceleration in the rate of decline.

Analysis of hotspot and coldspot ratios

A spatial autocorrelation analysis (Moran’s I) was conducted using ArcGIS 10.8, to examine whether changes in habitat quality across the 33 studied nature reserves in Guizhou Province exhibited spatial aggregation from 2000 to 2020. Sixteen nature reserves passed the significance test in 2000, 2010, and 2020, with global Moran’s I Z-scores exceeding 1.65 and p-values less than 0.1. Among these, eight reserves had Z-scores greater than 2.58 and p-values less than 0.01 (Table S3), suggesting a clear pattern of spatial aggregation. Additionally, the habitat quality of some reserves, such as Chishui Nature Reserve and Houshui River Nature Reserve in Suiyang County, transitioned from random distribution to significant spatial aggregation from 2000 to 2020, while others, including Mafulin Nature Reserve in Wuchuan County, Dashahe, and Kuankuoshui, exhibited a shift from significant spatial aggregation to random distribution.

The Getis-Ord Gi* statistical index in ArcGIS was the basis for a hotspot analysis method that was used to find habitat quality (Hotspots) and low value clustering (cold spots) in nature reserves. The interaction zones between hotspots and coldspots in nature reserves may result from developmental activities near major transportation routes and tourist facilities, leading to localized declines in habitat quality. For example, in Fanjingshan National Nature Reserve, habitat-quality hotspots areas were primarily concentrated in the core of the reserve, while coldspots areas were mainly distributed along the periphery, reflecting pressures from human activity at the periphery. In Rongjiang Moon Mountain Nature Reserve, habitat-quality hotspots showed a fragmented pattern in 2000 but evolved into a continuous distribution centered around coldspots by 2020, suggesting better habitat protection within the reserve resulting from reduced human activities (Figs. S1 to S3). Overall, from 2000 to 2020, most of the habitat-quality changes in the 33 nature reserves showed a trend of spatial aggregation that initially weakened and then strengthened. Hotspots were mainly concentrated in the core areas of the reserves, while coldspots were predominantly observed on built-up land and certain cultivated areas at the reserve boundaries.

Factors influencing the spatiotemporal evolution of habitat quality: optimal parameters-based geographical detector

Changes in habitat quality are likely driven by complex interactions of multiple factors, which, collectively, may have either a promoting or constraining effect. This study uses the bi-factor interaction detection module of the optimal parameter geographical detector model to examine how interactions between any two driving forces affect changes in habitat quality within protected areas. This method makes it possible to evaluate the interaction effects between any two elements by analyzing the interaction impacts of different influencing factors on the evolution of habitat quality. Bi-factor enhancement and non-linear enhancement are the main characteristics of interactions between any two driving factors, with no indication of separate interaction effects, according to Tables S4 and S5, which only display the top four interaction effect values.

National nature reserves

In 2000, the habitat quality of most nature reserves was predominantly influenced by natural factors, with a particularly significant impact observed in certain reserves (e.g., Fanjingshan and Guizhou Chishui Alsophila National Nature Reserve). However, the influence of natural factors in these areas gradually diminished by 2010 and 2020. This trend may indicate that the natural ecological conditions within these reserves have gradually stabilized or have benefited from effective protection measures. Meanwhile, socioeconomic factors have shown a steadily increasing influence across multiple reserves, especially in 2020 (e.g., in the Maolan and Fodingshan reserves, where the impact of land use factors has increased significantly). This trend is likely associated with socioeconomic development, population growth, and intensified human activities in areas surrounding the reserves, leading to a gradual rise in the impact of socioeconomic factors on habitat quality.

In some reserves (e.g., Leigongshan and Dashanhe), natural factors remain the dominant determinants of habitat quality, for example, elevation is a large influence among natural factors. However, in other reserves (e.g., Xishui and Kuankuoshui), the impact of socioeconomic factors has either approached or exceeded that of natural factors. This variation in influencing factors between reserves underscores the distinct environmental pressures and management challenges each reserve faces.

Although the factors influencing habitat quality vary among reserves, overall, natural factors exert a more pronounced effect on habitat quality evolution, especially geographical factors, whose explanatory power surpasses that of socioeconomic factors. Natural geographical factors provide the foundation for habitat quality evolution in national nature reserves in Guizhou Province, playing a pivotal role in either promoting or constraining the evolution of habitat quality (Fig. 5). In general, from 2000 to 2020, the influence of various factors on habitat quality within reserves has tended to stabilize or decrease.