Based on the LULC data, the landscape types of HTRNP were divided into five categories: cultivated land, forestland, grassland, water, and construction land. This study used the transfer matrix to analyze the landscape types transition in two adjacent time nodes (1980–1990, 1990–2000, 2000–2010, and 2010–2020) during 1980–2020.

Based on the actual situation of HTRNP and related study (Li L. et al., 2022), ten landscape indices were selected from Landscape and Class levels: percent of landscape (PLAND), number of patches (NP), patch density (PD), mean patch size (MPS), landscape shape index (LSI), largest patch index (LPI), Shannon’s diversity index (SHDI), Shannon’s evenness index (SHEI), contagion (CONTAG), and aggregation index (AI). These indices were calculated based on Fragstats v4.2 platform.

The biophysical value of WS was assessed using the AWY module. This module is based on the water balance principle and considers several factors, such as transpiration, evaporation, and precipitation. The calculation formula is as follows:

where Y_xj is the annual water yield on the pixel x of the landscape type j (mm); AET_xj is the actual annual evapotranspiration on the pixel x of the landscape type j (mm); P_x is the annual precipitation on the pixel x (mm) (Qi et al., 2019). The maximum root depth (Root_depth), plant evapotranspiration coefficient (K_c), and actual evapotranspiration value (LULC_veg) in the biophysical coefficient table (Supplementary Table 1) were set according to Zheng et al.’s (2019) research on the central mountainous region of Hainan Island. Zhang’s coefficients for each year were set by referring to the measured data of the Fucai hydrological station in HTRNP (Li A. et al., 2022) after several adjustments. The monetary value of WS was assessed using Alternative costing method. The calculation formula is as follows:

where V_ws is the monetary value of WS (CNY/a); Q_ws is the total annual water supply (m³/a); C is the engineering cost of reservoir constructing per unit capacity (CNY/m³). According to the “Yearbook of China water resources,” the average reservoir capacity cost is 2.17 CNY/m³, then C was calculated as 7.0547 CNY/m³ based on the 2009 price index (3.251) (Fang et al., 2013).

The biophysical value of WP was assessed using the NDR module, which inversely characterizes the WP service function by the total nitrogen (TN) and total phosphorus (TP) exports. The higher the TN and TP exports per unit area, the weaker the WP service function. The calculation formula is as follows:

where ALV_x is the adjusted loading value on the pixel x; pol_x is the export coefficient on the pixel x; HSS_x is the hydrologic sensitivity score on the pixel x, λ_x is the runoff index on the pixel x, λw¯ is the average runoff index in the watershed of interest (Qi et al., 2019). Nitrogen (N) and phosphorus (P) export coefficients and retention efficiency (Supplementary Table 2) for landscape types were set with reference to relevant studies (Zhe et al., 2013; Zheng et al., 2019). The monetary value of WP was assessed using Alternative costing method. The calculation formula is as follows:

where V_wp is the monetary value of WP (CNY/a); T_N and T_P are the total amount of TN and TP purification (t/a); P_N and P_P are the purification unit-prices of TN and TP (CNY/t), and the values are 1,750 and 2,800 CNY/t, respectively, according to the “Levy standard of pollution discharge fees and calculation methods” of the National Development and Reform Commission of China (Fan and Li, 2020).

The biophysical value of CS was assessed using the CSS module, which estimates CS based on each landscape type and its corresponding four primary carbon pools. The calculation formula is as follows:

where C_total is the total CS (t/hm²); C_above, C_below, C_soil, and C_dead are aboveground CS (t/hm²), belowground CS (t/hm²), soil CS (t/hm²), and dead organic CS (t/hm²), respectively. As the vegetation of HTRNP is mainly evergreen broad-leaved species of zonal forest type, its dead organic CS is tiny and difficult to estimate (Gong et al., 2022). This part of CS was not calculated in this study. Other three types of carbon density were set (Supplementary Table 3) with reference to Liu et al.’s (2022) study. The monetary value of CS was assessed using a combination of Alternative costing and Carbon tax methods. The calculation formula is as follows:

where V_cs is the monetary value of CS (CNY/a); C_t is the total CS (t/a). To obtain more accurate monetary value assessment results, this study referred to Wang et al.’s (2017) study, the unit-price of forest carbon sink in Alternative costing method used the arithmetic average of four unit-prices (251.40, 260.90, 273.30, and 305.00 CNY/t); the Carbon tax method used the international Swedish carbon tax rate (150 USD/t = 1,017.675 CNY/t, the exchange rate was calculated as 100 USD = 678.45 CNY on 29 January 2023).

The SDR module assesses SR’s biophysical value based on the Revised Universal Soil Loss Equation (RUSLE). The calculation formula is as follows:

where A_x is the SR amount on the pixel x (t⋅hm^–2⋅a^–1); R_x is the rainfall erosivity factor on the pixel x (MJ⋅mm⋅hm^–2⋅h^–1⋅a^–1); K_x is the soil erodibility factor on the pixel x (t⋅hm²⋅h⋅hm^–2⋅MJ^–1⋅mm^–1); LS_x is the topographical factor on the pixel x; C_x is the crop-management factor on the pixel x; P_x is the support practice factor (Qi et al., 2019). R was calculated using the rainfall erosion force model (Zhang and Fu, 2003), R_j is the rainfall erosion force in year j (MJ⋅mm⋅hm^–2⋅h^–1⋅a^–1), P_j is the rainfall in year j (mm). K was calculated using the Erosion-Productivity Impact Calculator (EPIC) model (Williams, 1995), the K_EPIC was corrected according to Zhang et al. (2008). The values of C and P were set (Supplementary Table 4) according to the relevant studies (Xiao, 1999; Zheng et al., 2019). The monetary value of SR (including the value of both sedimentation reduction and non-point source pollution reduction) assessed using the Alternative costing method. The calculation formula is as follows:

where V_sr is the monetary value of SR (CNY/a); V_sd is the value of sedimentation reduction (CNY/a); V_dpd is the value of non-point source pollution reduction (CNY/a); λ is the sedimentation coefficient, taking the value of 0.24 (Sheng et al., 2010); Q_sr is the SR amount (t/a); C is the cost of reservoir desilting project per unit area (CNY/m³), according to the “The building built water engineering budget norm” of the Ministry of Water Resources of the People’s Republic of China, the value of C is 17.63 CNY/m³ (Yu et al., 2020); ρ is the soil capacity (t/m³), according to the average of the measured data of each forest ecosystem type in the central mountainous region of Hainan island in 2008, the value of ρ is 1.28 t/m³ (Liu et al., 2009); C_i is the pure content of N (or P) in soil (%), and it was determined that the content of N and P in soil of China are 0.370 and 0.108%, respectively (Wang L. et al., 2017); P_i is the degradation cost of N (or P), according to the “Levy standard of pollution discharge fees and calculation methods,” the degradation costs of N and P are 1,750 and 2,800 CNY/t, respectively (Fan and Li, 2020).

The biophysical value of HQ was assessed using the HQ module. The value of the HQ index is in the range of [0,1], and the higher value indicates the higher level of biodiversity. The calculation formula is as follows:

where Q_xj is the HQ on the pixel x of the landscape type j; H_j is the habitat suitability of the landscape type j; D_xj is the total threat level on the pixel x of the landscape type j; K is the half-saturation coefficient, generally taking a value of 0.5; z is a scaling factor, generally taking a value of 2.5 (Yang et al., 2021). Referring to the studies on HTRNP and its neighboring regions (Lei et al., 2022; Yao et al., 2022), paddy fields, dry land, rural residential areas, other construction lands, and expressways were selected as threat factors in this study. The impact distance and weight of the threat factors (Supplementary Table 5) and the sensitivity of each landscape type (Supplementary Table 6) were set with reference to the model manual and the above studies. Referring to related research (Xiao et al., 2014; Sun et al., 2019), the monetary value of HQ was assessed based on the area of excellent HQ regions, and the ecological benefit of the excellent HQ region is about 191.28 × 10⁴ CNY/km². The excellent HQ regions were obtained by referring to Lei et al.’s (2022) study using the Natural Breaks method to classify the HQ spatial distribution into four classes (Poor: 0–0.3, Medium: 0.3–0.7, Good: 0.7–0.9, and Excellent: 0.9–1.0).

Testing whether there is a correlation between landscape indices and various types of ESVs is a critical a priori step in exploring the correlation between ESVs and landscape pattern changes, to test whether correlations exist among variables, whether the direction and magnitude of the correlations are as expected, and whether they apply to subsequent more complex multivariate analyses (Cen, 2016; Ge, 2020).

The ranking analysis is widely used in ecological studies to explain the response relationships between species and environmental variables, this method can effectively downscale the multi-dimensional information, and its analysis results are concise and intuitive (Legendre, 2008; Cen, 2016). According to the results of DCA (length of gradient <3), this study applied the Redundancy Analysis (RDA) based on the Canoco v5.0 platform to rank the correlation between ESVs and landscape patterns. The RDA is a ranking analysis method combining regression and Principal Component Analysis (PCA). It is an extension of multiple regression analysis and can be used to model multivariate response data. The RDA ranking chart can visualize the relationship between environment and response variables (Rao et al., 2016). The angle between their arrows reflects the correlation between them: when the angle <90°, it indicates a positive correlation between the variables, and the smaller the angle, the stronger the positive correlation; when the angle >90°, it indicates a negative correlation between the variables, and the larger the angle, the stronger the negative correlation; when the angle = 90°, it shows that the variables are not correlated with each other (Li C. et al., 2022).

From both the distribution (Figure 3) and proportion (Figure 4) of landscape types, forestland is the most dominant landscape type in the HTRNP, followed by grassland, and construction land occupies a minor proportion. During 1980–2020, the proportion of forestland fluctuated between 89.86 and 90.29%, with little overall change; the proportion of cultivated land, grassland, and water changed relatively significantly. The conversion between different landscape types became frequent after 2000, especially during 2000–2010. The total area conversion reached 13,446.00 hm² and then gradually eased after 2010.

(1) During 2000–2010, grassland, cultivated land, and forestland areas were transferred out more obviously, accounting for 52.06, 33.18, and 14.72% of the total area transferred out, respectively. Although forestland was one of the primary sources of transfer out area, its proportion showed positive growth because of the significant area input of grassland and cultivated land (Figure 4). The conversions of cultivated land in the northwest and grassland in the southeast to water and forestland, respectively, were evident (Figure 3). (2) The magnitude of changes in landscape types during 2010–2020 was less than in the previous period (Figure 4), and spatially (Figure 3), only the conversion of part of the northwestern water to cultivated land was slightly noticeable. Although the conversion between forestland and grassland was apparent (Figure 4), the main areas’ conversion was limited between these two types, and the number of areas transferred between them was similar. Hence, the changes in the proportion of these two types were relatively small.

During 1980–2020, at the Landscape level (Figure 5): NP and PD showed a trend of “slow decline – rapid decline – rapid rise,” with a slight overall decline, indicating that landscape fragmentation and heterogeneity of HTRNP tended to decline; SHDI and SHEI showed a trend of “slow decline – rapid decline – slow rise,” with an overall decline, indicating that the homogeneity of the HTRNP’s landscape decreased, the dominant landscape patch types tended to be prominent, and the distribution of landscape patch types in HTRNP tended to be uneven, which may be related to the shrinkage of both grassland and cultivated land and the expansion of both forestland and water. CONTAG and AI showed a trend of “slow rise – rapid rise – slow decline,” these two indices tended to increase overall, indicating that the connectivity of the dominant landscape in HTRNP was enhanced, and the degree of landscape agglomeration was increased.

At the Class level (Figure 5): (1) Cultivated land: NP and PD tended to increase, PLAND, LPI, and MPS tended to decrease, indicating that the cultivated land tended to be fragmented and its landscape dominance decreased; LSI tended to increase, indicating that the landscape shape tended to be complex. (2) Forestland: NP and PD tended to decrease, indicating an increase in forestland aggregation; the three indices of PLAND, LPI, and MPS of forestland were the highest among the five landscape types, and all tended to increase, indicating the landscape pattern of HTRNP were dominated by large patches of forestland, and the landscape dominance of forestland tended to increase; LSI tended to decrease, indicating its landscape shape become more regularized. (3) Grassland: grassland had the highest NP and PD among the five landscape types, indicating the highest degree of fragmentation in grassland landscape, but its fragmentation tended to decrease; LPI changed less, PLAND and MPS tended to decrease, indicating a decrease in its fragmentation was related to the conversion of some small scattered patches to other landscape patch types; LSI tended to decrease, indicating that the landscape shape tended to be regularized. (4) Water: NP and PD tended to increase but at a lower rate; PLAND, LPI, and MPS all tended to increase, indicating an increase in landscape dominance; LSI tended to decrease in general, indicating that the landscape shape tended to be regularized. (5) Construction land: NP, PD, PLAND, LPI, and MPS increased slightly, indicating its area was increasing in fragmented patches; LSI tended to increase, indicating its shape became complex. In summary, landscape indices at both levels changed significantly during 2000–2020.

All the five ESs of HTRNP showed varying degrees of trends during 1980–2020 (Figure 6).

(1) WS: WS’s biophysical value decreased continuously during 1980–2010; the decrease was more evident during 1990–2010, the total WY decreased by 14.16%, and the average WY depth decreased from 1,265.83 to 1,090.93 mm; it rebounded after 2010. There has been a declining trend in WS over the past 40 years, with total annual WY decreasing by 7.82%. (2) WP: the average TN and TP exports both decreased rapidly during 2000–2010, with the average TN export decreasing from 9.19 × 10^–2 kg/hm² to 9.09 × 10^–2 kg/hm² and the average TP export decreasing from 7.13 × 10^–3 kg/hm² to 6.91 × 10^–3 kg/hm². The total exports of both TN and TP had decreased slightly over the past 40 years, by 1.24 and 2.88%, respectively. (3) CS: the average CS declined rapidly during 2000–2010, from 12.97 t/hm² to 12.92 t/hm², and the total CS decreased by 0.53%; it rebounded slightly after 2010. There has been a slight overall downward trend in CS over the past 40 years, with a total reduction of 0.35%. (4) SR: the biophysical value of SR decreased significantly during 2010–2020, with a 25.81% decrease in total SR and a decrease in SR capacity per unit area from 908.68 t/hm² to 673.60 t/hm². SR decreased from 1980 to 2020, and the total SR decreased by 33.18%. (5) HQ: the biophysical value of HQ rose rapidly during 2000–2010, with the average HQ index peaking at 0.964; it dropped slightly to 0.963 after 2010. The overall HQ remained high and tended to increase during 1980–2020, the proportion of excellent HQ zone remained above 80%, and the average HQ index remained above 0.960.

The mean values of ESs at the sub-watershed scale were calculated based on the Zonal Statistics tool of ArcGIS v10.2 platform, and the ESs were normalized from 0 to 1 (Low to High) for each year to facilitate comparison across years (Xia et al., 2023). All five ESs of HTRNP showed different degrees of spatial heterogeneity during 1980–2020 (Figure 7).

(1) WS: the spatial heterogeneity of WS is “high in the east and low in the west.” This spatial heterogeneity tended to be significant during 2000–2010, the WS in the southwestern region (including JFL NNR) continued to decrease, and the low-value regions expanded toward the BWL and YGL NNRs in the central region. This trend decreased in 2020. The high-value regions were mainly concentrated in the southeast (including WZS and DLS NNRs). (2) WP: the spatial distribution of both TN and TP exports showed a gradient pattern increasing from the southwestern region (including JFL NNR) to the central region (including BWL and YGL NNRs) and then to the southeastern region (including WZS and DLS NNRs). The high-value regions in the northwest showed a transformation trend to the low-value regions. (3) CS: CS’s spatial heterogeneity did not vary significantly, the low-value regions were concentrated in the northwest, and the high-value regions were concentrated in the southwest (including JFL NNR). (4) SR: SR’s spatial heterogeneity was observed. The low-value regions in the southwest (including JFL NNR) showed a significant expansion during 2000–2010. This trend decreased in 2020. (5) HQ: the spatial heterogeneity of HQ was not observed and did not change significantly. Relatively, the HQ in the northwestern region was slightly lower.

The landscape indices, landscape areas, and ESVs of HTRNP during 1980–2020 were used as data sources to test the correlation between these variables. The correlation test was performed based on the IBM SPSS Statistics v26 platform, using the Spearman correlation coefficient as a measure and a two-tailed t-test for the significance of the correlation coefficient. The results showed (Supplementary Figure 1) that the variables had different degrees of correlation.

Based on the previous analysis of the HTRNP’s landscape pattern changes, it can be noted that the trends between some of the indices are very similar, indicating that there may be some degree of correlation between them. The strong correlation between indices will result in indices that do not satisfy the statistical properties of mutual independence, causing duplication and redundancy in the mathematical and theoretical significance of subsequent studies, thus affecting the accuracy of the results (Cen, 2016; Ge, 2020). Therefore, in this study, the landscape indices with low correlations were further screened based on the results of the correlation test for subsequent ranking analysis: (1) Landscape level: PD, SHDI, and CONTAG; (2) Class level: PD, LPI, and LSI.

The RDA ranking analysis (Figure 9) was conducted based on the Canoco v5.0 platform using ESVs as response variables and the landscape areas and indices as environmental variables. Based on the correlation test results between landscape indices, the environmental variable groups are as follows: (1) Landscape level: PD, SHDI, and CONTAG; (2) Class level: LPI, PD, and LSI.

(1) The changes of landscape types over the past 40 years have impacted ESVs. The explanatory degree of cultivated land (84.3%), forestland (85.8%), grassland (85.7%), and water (85.1%) areas are all higher to ESVs, indicating that changes in these four landscape types had significant effects on ESVs. In comparison, the explanatory degree of construction land (39.5%) area is lower, indicating that its changes had weaker effects on ESVs, which may be related to its small area. WP and CS have a strong positive correlation with the areas of grassland and cultivated land and a strong negative correlation with the areas of forestland and water; SR has a strong negative correlation with the area of construction land; HQ has a strong positive correlation with the areas of forestland and water and a strong negative correlation with the areas of grassland and cultivated land.

(2) At the Class level: landscape indices’ effects on ESVs across landscape types vary.

Cultivated land: the explanation degree of LPI to ESVs is the highest, reaching 85.7%, followed by PD; LSI has the lowest explanatory degree and did not explain changes in ESVs significantly (p > 0.10), and also had a weak correlation with ESVs. WP-LPI and CS-LPI have strong positive correlations, while HQ-LPI and SR-PD have strong negative correlations.

Forestland: the explanation degree of LPI to ESVs is the highest, reaching 85.8%, followed by LSI; PD has the lowest explanatory degree and did not explain changes in ESVs significantly (p > 0.10). WP-LSI, CS-LSI, and HQ-LPI have strong positive correlations, while WP-LPI, CS-LPI, and HQ-LSI have strong negative correlations.

Grassland: the explanation degree of PD to ESVs is the highest, reaching 83.7%, followed by LPI; LSI has the lowest explanatory degree and did not explain changes in ESVs significantly (p > 0.10). WP-PD and CS-PD have strong positive correlations, while HQ-PD has strong negative correlations.

Water: the explanation degree of LPI to ESVs is the highest, reaching 84.9%, followed by LSI; PD has the lowest explanatory degree. HQ-LPI has strong positive correlations, while WP-LPI and CS-LPI have strong negative correlations.

Construction land: compared with other landscape types, the explanatory degrees of landscape indices of construction land to ESVs are lower. The index with the highest explanatory degree is LSI (explanatory degree only reached 56.1%), followed by LPI; PD has the lowest explanatory degree and did not explain changes in ESVs significantly (p > 0.10). SR-LSI and SR-LPI have strong negative correlations.

(3) At the Landscape level: the explanation degree of CONTAG to ESVs is the highest, reaching 82.0%, followed by PD; SHDI has the lowest explanatory degree and did not explain changes in ESVs significantly (p > 0.10). HQ-CONTAG has strong positive correlations, while WP-CONTAG and CS-CONTAG have strong negative correlations.