The Multiscale Spatiotemporal Heterogeneity of Ecosystem Service Trade-Offs/Synergies and Bundles and Socioecological Drivers in the Yangtze River Delta Region of China

Abstract

A comprehensive exploration of the trade-offs/synergies and drivers of ecosystem services (ESs) is essential for formulating ecological plans. However, owing to the limited attention given to multiple scales, the relationship of ESs still needs to be further explored. Taking the Yangtze River Delta region of China as the study area, a multiscale data framework with a 1 km grid and 10 km grid and county was established, and six ESs were evaluated for 2000, 2010, and 2020. Then, the trade-offs and synergies between ESs were explored by Spearman’s correlation analysis and geographically weighted regression (GWR), and the ecosystem service bundles (ESBs) were identified by self-organizing maps (SOMs). Finally, the socioecological drivers of ESs were further analyzed via GeoDetector. The results showed that (1) the distribution of ESs exhibited spatial heterogeneity. (2) At the grid scale, there were very strong trade-off effects between crop production and the other ESs. The synergistic effects between ESs at the county level were further strengthened. (3) The ESBs identified at different temporal and spatial scales were different. (4) Land use had the strongest explanatory power for all the ESs. At the grid scale, climatic and biophysical factors had great impacts on ESs, whereas population density and night light remote sensing had significant impacts on crop production, carbon storage, and water yield at the county scale.

1. Introduction

Ecosystem services (ESs) are the benefits that humans directly or indirectly obtain from the ecosystem [1]. ESs are considered the bridge connecting the natural system and the socioeconomic system. ESs can be divided into four types of services: supply, regulation, support, and culture [2]. With global climate change and increased human activity, the sustainable supply of ESs is severely limited [3]. The spatiotemporal patterns and drivers of ecosystem service trade-offs and synergies (TOS) have become core scientific issues in the coordination of ecological protection and well-being among humans. With the increase in human activities (such as urbanization, agricultural activities, and industrial production), the interactions among ESs differ, and their complexity and scale dependence necessitate greater consideration when performing research. However, the interactions among ESs are affected by natural and social factors at multiple scales, and the laws governing the dynamic evolution of these interactions are still unknown. Cross-scale systematic research involving multiple methods is urgently needed to identify this law [4]. Understanding the spatial differentiation and drivers of these interactions is highly essential for regional ecological planning and the optimization of ecological security patterns.

The diversity and spatiotemporal heterogeneity of ESs and social economic policies bring about complex interactions among ESs [5,6], including trade-offs, synergies, and bundles [7,8]. An ecosystem service trade-off (EST) reflects that an increase in one ES corresponds to a decrease in other ESs, and ecosystem service synergy (ESS) indicates that pairs of or multiple ESs change in the same direction [9]. An ecosystem service bundle (ESB) refers to a series of ESs that occur repeatedly in space or time [10]. Analyzing TOS and identifying ESBs can help us to more systematically explore the interactions among ESs and more intuitively identify critical regions for ecosystem management [5]. Researchers have developed many methods to analyze the complex TOS among ESs. Correlation analysis is often used to explore the overall TOS among ESs. Because not all ESs follow a normal distribution, Spearman’s correlation analysis is widely used to explore correlations among them [11,12,13]. In addition, TOS may differ spatially due to differences in environmental conditions. Geographically weighted regression (GWR) considers the local effects of geographical elements, which is useful to identify the spatially explicit pattern of TOS between ES pairs [6,8]. Typical techniques for identifying ESBs include K-means clustering [11], principal component analysis [14], Gaussian mixture models [15], and self-organizing maps (SOMs) [8]. Among these techniques, SOMs are advantageous because they combine dimensionality reduction with clustering analysis, with a high fault tolerance, robustness, and multiple scales; thus, they represent a promising new method for identifying ESBs [6,8,16].

Although previous studies have explored the spatiotemporal heterogeneity of the relationships among ESs, there is still a lack of research considering multiple time series and multiple spatial scales. In terms of time, Mo et al. (2023) studied the TOS effect of ESs in the Dongjiang River Basin in 2020 at the township and grid scales [17], and Huang et al. (2023) studied the TOS of each ESB and the relevant drivers in the Wujiang River Basin in 2018 at the grid scale [18]. However, these studies focused on a single time point and therefore could not elucidate the complex time-dependent changes in interactions among ESs. In terms of space, Huang and Wu (2023) studied the temporal and spatial drivers of the TOS of ESs in the urban agglomeration of the Yangtze River Delta in China from 2005 to 2020 at the grid scale [12]. In addition, He et al. (2024) and Wang et al. (2023) studied the relationships among ESs and the relevant drivers at the county and town scales [11,19]. Unfortunately, as these studies focused on a single spatial scale, the effects of the interactions among ESs may vary significantly depending on the scale due to TOS among ESs and their drivers [7,8,20]. The interactions among ESs may undergo complex changes due to the dual effects of time and space [21]. Therefore, considering the different characteristics of ES interactions at different temporal and spatial scales, the spatiotemporal dynamics of these interactions need to be more comprehensively studied to more effectively support the implementation of spatial ecological planning [8].

The supply capacity of ESs is affected by many social and ecological factors, such as topography, climate change, biophysical conditions, landscape composition and socioeconomic conditions [6,8,11,12]. A variety of technical models have been developed and applied to identify the main drivers and analyze the driving paths of ESs. Xia et al. (2023) identified the main socioecological drivers of each ES in the Qiantang River Basin using GeoDetector [6]. Huang et al. (2023) investigated the socioecological drivers of ESs through a redundancy analysis [22]. Huang and Wu (2023) evaluated the impacts of climate, land use, and socioeconomic factors on ESs in the urban agglomeration of the Yangtze River Delta through a geographical and time-weighted regression analysis [12]. Liu et al. (2024) applied a partial least squares structural equation model to evaluate the direct and indirect effects of various driving factors on ESs in the Yellow River Basin [8]. When using GeoDetector, one does not need to assume that there is a linear relationship between variables, and GeoDetector can explore the interactions of multiple drivers, can quantify spatial differentiation, and has the ability to perform causal detection, making it a powerful method for exploring spatial heterogeneity and the impacts of multiple factors [23]. Notably, different driving factors may have different impacts on ESs at different temporal and spatial scales. Therefore, analyzing the impacts of driving factors on ESs and exploring their mechanisms at multiple temporal and spatial scales will hopefully improve the management of regional ecological planning.

The Yangtze River Delta Region (YRDR) is the most energetic region in terms of economic development in China and is a typical region with highly concentrated ecological and environmental problems [12]. Rapid urbanization has led to serious ecological problems (such as water quality deterioration, the loss of biodiversity, and ecosystem degradation). The sustainable development of the YRDR has emerged as a serious challenge. Research on the interactions among ESs and relevant drivers at multiple temporal and spatial scales in this region is still limited. We aimed to (1) analyze the patterns of spatiotemporal changes in ESs, (2) explore the multiple temporal and spatial effects of the complex relationships among ESs, and (3) determine the dominant socioecological drivers of ESs and analyze their driving mechanisms at multiple scales. This study can help policy makers fully understand the relationships among ESs and the relevant driving mechanisms, and provide some information for regional ecosystem management and ecological planning.

2. Materials and Methods

2.1. Study Area

The YRDR (114°52′48″–122°56′24″ E, 27°2′24″–35°7′48″ N) is located in eastern China and includes Shanghai, Jiangsu Province, Zhejiang Province, and Anhui Province (Figure 1), with a total area of 358,000 km². It is the only giant delta in China that covers both a megalopolis and multiple natural geographical units. The region has a subtropical monsoon climate, with an annual average temperature of 14–18 °C and an annual average precipitation of 1000–1600 mm. Rain and heat occur during the same period; however, the region is vulnerable to typhoons in summer and autumn, and floods occur frequently. Two-thirds of the land in the Northeast China Plain is covered by cropland. Rivers and lakes are densely distributed, and cities and towns are scattered. From 2000 to 2020, urban construction land expanded (Figure S1). Mountains and hills are concentrated in the southwest, and the area is rich in forests and biological resources. Coastal wetlands and estuarine beaches are rich in resources and are important post stations for migratory birds in East Asia and Australasia.

The YRDR is one of the regions with the most active economic development and the fastest urbanization in China. The population increased by 19%, from 197.09 million in 2000 to 235.38 million in 2020. Meanwhile, the GDP increased nearly tenfold (from CNY 2265.6 billion to CNY 24,452.2 billion), and the urban and area increased by 80% (from 24,300 km² to 43,715 km²). In 2015, the urbanization level of the YRDR reached 68.2% [24]. In 2019, the regional integration development of the Yangtze River Delta became a national strategy. However, rapid economic development and urbanization also led to ecological degradation and resource constraints, such as surface water pollution, nongrain farmland and coastal eutrophication.

2.2. Data Collection

We evaluated the changes in key ESs on the basis of multisource data from 2000, 2010, and 2020. The raster data with various resolutions were resampled to a uniform spatial resolution of 1 km × 1 km and uniformly projected into WGS_1984_UTM_Zone_50N. More details about the datasets are shown in Table 1.

2.3. Quantification of Ecosystem Services

This study identified six ESs through the following methods: (1) selecting those in the four categories (including provisioning services, regulating services, supporting services, and cultural services) proposed by the Millennium Ecosystem Assessment; (2) comprehensively considering the resources, environment, and socioeconomic conditions of the study area and selecting the most representative ecosystem services; and (3) considering the feasibility of data collection. We finally selected six ESs (Table 2), namely, crop production (CP), carbon storage (CS), soil conservation (SC), water yield (WY), habitat quality (HQ), and leisure recreation (LR). Cropland is widely distributed in the plains of the northeast area and is an important grain production base in China; so CP was selected as the provisioning service. Combined with the relevant “carbon peaking” and “carbon neutralization” policies proposed by the Chinese government, CS was selected as the regulating service. Owing to the complex landform and intensive human activities, soil erosion has become the main ecological problem in the YRDR. WY is very important for ecology and domestic water, and can directly benefit residents. Mountains and forests are widely distributed in the southwestern area and act as important barriers preventing the maintenance of biodiversity and ecological security in the YRDR. Therefore, HQ was selected as the supporting service. Finally, the residents need to enjoy natural scenery and entertainment facilities, and LR can improve residents’ sense of happiness.

We selected three spatial scales to evaluate the ESs. The 1 km grid scale is suitable for fine and personalized ecosystem management, the 10 km grid scale is similar to the township scale for local landscape management [7], and the county scale is the basic scale for landscape planning and management [20]. We first quantified ES values for 2000, 2010, and 2020 at the 1 km grid scale and then used the regional statistical toolbox in ArcGIS 10.8 to calculate the average ES values at the 10 km grid and county scales. The detailed methods for assessing ecosystem services are provided in the Supplementary Material.

2.4. Data Analysis Methods

2.4.1. Data Preprocessing

Before the statistical analysis of ES values, it was essential to exclude abnormal values via a box plot. Then we tested the normal distribution for the data of each ES by constructing a histogram and assessed the linear correlations among the ES variables via a scatter diagram. The results revealed that the data were not normally distributed for all ESs, and some of the relationships between ES pairs were nonlinear (such as CP-CS, CP-SC, CP-WY, CP-HQ, CS-SC, CS-HQ, SC-HQ, and SC-LR). In addition, to eliminate the difference in amplitude among ES values, the data were normalized by the maximum–minimum value.

(1)${ES}_i={\displaystyle \frac{ES-{ES}_{min}}{{ES}_{max}-{ES}_{min}}}$

where ${ES}_i$ is the standardized ES value, and ${ES}_{max}$ and ${ES}_{min}$ are the maximum and minimum ES values, respectively.

2.4.2. Trade-Off and Synergy Analyses Between ESs

Analyses of trade-offs or synergies can guide the realization of a sustainable supply of ESs [26]. In this study, we used Spearman’s correlation analysis to analyze the overall TOS among the ESs at three spatial scales (1 km grid, 10 km grid, and county scales) from 2000 to 2020. A positive correlation indicates a synergistic effect between ESs, a negative correlation indicates a trade-off, and the absolute value of the correlation coefficient reflects the strength of the TOS. We used the “corrplot” package in R 4.3 for Spearman’s correlation analysis and significance testing [11].

The spatial patterns of TOS among the ESs were further elucidated by GWR. Similar to Spearman’s correlations, positive regression coefficients represent spatial synergy, whereas negative regression coefficients represent spatial trade-offs. The “GWR” package in R 4.2.3 was used to achieve geographically weighted regression [27].

2.4.3. Identifying Ecosystem Service Bundles

ESBs are an important tool for identifying the internal connections among ESs and displaying the spatial relationships of various ESs, which are very important for the division of landscape ecological function areas [11,28,29]. We used SOMs to identify the bundles of the six ESs in the YRDR in 2000, 2010, and 2020 at multiple spatial scales to understand the spatial distributions and structural characteristics of the different ESBs. To scientifically determine the number of clusters, we used the Davies–Bouldin index to quantify the optimal number of clusters [27]. SOMs were constructed via the “Kohonen” package in R 4.2.3 software.

2.4.4. Identification and Analysis of the Socioeconomic and Ecological Drivers of ESs

The natural and socioeconomic factors were used to analyze the driving mechanism of ESs. On the basis of the representative socioecological drivers used in previous studies, we divided the potential factors into four categories (climate, biophysical indicators, landscape composition, and socioeconomic indicators) and selected 16 indicators (Table 3). To explore the key socioecological drivers that affect the spatial heterogeneity of ESs, we used GeoDetector to determine the explanatory power of potential socioecological factors [6]. The “GD” package in R 4.3 was used to perform the GeoDetector method; ES was the dependent variable, and the natural and socioecological factors in Table 3 were the independent variables.

The GeoDetector statistic is a method used to measure the degree of spatial stratification heterogeneity and reveal the driving forces of various natural and socioeconomic processes [30], and the statistic is calculated as follows:

(2)$q=1-{\displaystyle \frac{{\sum }_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }_h^2}}=1-{\displaystyle \frac{SSW}{SST}}$

(3)$SSW={\sum }_{h=1}^L{N}_h{\sigma }_h^2$

(4)$SST=N{\sigma }^2$

where $q$ is the explanatory power of the driving factors; $h$ is the stratification of variables; $N$ and $N_h$ are the numbers of cells in the entire area and layer h, respectively; ${\sigma }^2$ and ${\sigma }_h^2$ are the variances in the entire area and layer h, respectively; and $SSW$ and $SST$ are the sums of intralayer variance and total variance, respectively. The range of the $q$ value is [0, 1], and the higher the value is, the stronger the explanatory power of the driving factors.

3. Results

3.1. Spatiotemporal Variation in Ecosystem Services

This study evaluated six ESs in the YRDR from 2000 to 2020. The results revealed that the ESs exhibited spatial heterogeneity, and the spatial distribution of each ES remained essentially stable over time (Figure 2). CS, SC, and HQ showed similar spatial distribution patterns [12]. Specifically, the areas with a high ES supply were distributed mainly in the mountains in the southwest, with less human disturbance and high biodiversity levels, whereas the areas with a low ES supply were distributed in the vast plains in the northeast, with relatively low naturality and habitat suitability. However, CP showed the opposite spatial distribution pattern, with the highest values concentrated in the vast plain area in the northeast, which benefits from appropriate farming conditions and intensive human activities. In addition, the areas with a high WY supply were mainly distributed along the Yangtze River and in the southern mountainous region. The areas with high LR were located mainly in southern mountains and in metropolitan areas.

The total supply of ESs fluctuated significantly over the past 20 years. CP, SC, and WY values tended to increase [12], rising by 19.5%, 59.0%, and 40.1%, respectively; CS, HQ, and LR values, however, tended to decrease [12], decreasing by 4.1%, 5.4%, and 3.4%, respectively (Figure 3).

3.2. TOS Between ESs

3.2.1. Overall TOS Between ES Pairs

Using Spearman’s correlation analysis, we analyzed the overall TOS among the six ESs at three spatial scales from 2000 to 2020 (Figure 4). The significance of TOS between ES pairs was determined using the p value (p < 0.05), and most ES pairs showed significance. At the 1 km grid scale, CP had a very strong trade-off effect with CS, SC, HQ, and LR in all years (p < 0.001), and CP–LR had the greatest trade-off effect (r < −0.58). CS, SC, HQ, and LR had significant synergy in all years, and synergy was highest for the CS–SC pair (r > 0.74), followed by the HQ–LR pair (r > 0.66). At the 10 km grid scale, CP also had a very strong trade-off effect with CS, SC, HQ, and LR in all years (p < 0.001), and the trade-off effects were greater for the CP–LR and CP–HQ pairs (r < −0.53). CS, SC, HQ, and LR had significant synergy in all years, and the CS–SC pair had the greatest synergy (r > 0.81), followed by the HQ–LR pair (r > 0.69). WY had a trade-off effect with CP and HQ in 2010 and 2020 and showed synergy with other ESs. At the county scale, the CP–LR pair had a very strong trade-off effect in all years (p < 0.001). CP, CS, SC, and WY had very strong synergy (p < 0.001); the synergies of the CS–WY pair and the CS–SC pair were the greatest (r > 0.93, r > 0.83). HQ had very strong synergy with CS, SC, WY, and LR (p < 0.001), and LR had significant synergy with SC and HQ (p < 0.001). In general, synergy existed mainly between regulating services, supporting services, and cultural services, and trade-offs were related mainly to provisioning services (especially crop production) [13]. Among them, the CS–SC pair had the greatest synergistic effect, whereas the CP–LR pair had the greatest trade-off effect.

3.2.2. Spatiotemporal Characteristics of the TOS Between ES Pairs

Correlation analysis was used to determine the overall TOS among the ESs, but could not capture the spatial heterogeneity of TOS among the ESs. We used GWR to further reveal the spatial distributions and intensities of the TOS among the ESs at three spatial scales. As shown in Figure 5, there was spatial heterogeneity in the TOS between ES pairs. There were three spatial distribution patterns of EST/ESS at the 1 km grid scale. First, for the CP–CS, CP–HQ, and CP–LR pairs, the areas of high trade-off were concentrated in the northern plain agricultural area, and the areas of high synergy were concentrated in the eastern Taihu Lake Basin and urban built-up areas, showing a trend consistent with the urban–rural gradient. Second, the areas of high trade-off for the CS–SC, CS–WY, SC–HQ, SC–LR, WY–HQ, and WY–LR pairs were concentrated along the Yangtze River, and the trade-off intensity gradually decreased in South China and North China, which was related to urbanization and complex human activities along the Yangtze River. Third, synergistic effects of the CS–HQ, CS–LR, SC–WY, and HQ–LR pairs were observed in most regions. The areas of high synergy for the CS–HQ, CS–LR, and HQ–LR pairs were concentrated in the vast plains in the northeast, the vast mountains in the southwest, and the intermountain basin, respectively. At the 10 km grid and county scales, the three EST/ESS spatial distribution patterns were similar to those at the 1 km grid scale. In particular, the spatial distribution pattern of EST/ESS was smoother at the 10 km grid scale and tended to converge at the county scale.

3.3. Identification and Spatiotemporal Dynamic Patterns of ESBs

SOMs were used to identify ESBs at different spatiotemporal scales. According to the ES characteristics of the landscape functional area, the spatiotemporal dynamic patterns of ESBs at three spatial scales were as follows (Figure 6):

At the 1 km grid scale, two ESBs were identified in all years. (1) The crop production bundle (CPB) covered 70% of the YRDR and was mainly distributed in the northern plain and southern inland basin, and the land use type was cropland. The CPB had a high CP supply capacity, and the supply capacities of the other ESs were low. The region covered by the CPB is an important grain production base, and the ecological environment in this region has been seriously threatened by human activities. (2) The integrated ecology bundle (IEB) covered 30% of the YRDR and was mainly distributed in the western Dabie Mountains and southern mountainous and hilly regions, and the land use type was forest. Except for CP, the supply capacity of the other ESs in this region was high. IEB was the important ecological barrier in the YRDR.

Three, two, and four ESBs were identified at the 10 km grid scale in 2000, 2010, and 2020, respectively. (1) The lake coastal habitat bundle (LCHB) covered 14.9% and 13.8% of the YRDR, respectively, in 2000 and 2020 and was mainly distributed in lakes and eastern coastal areas. The supply capacity of HQ was high, whereas the supply capacities of CS, SC, and WY were relatively low. (2) The water yield bundle (WYB) covered 31.3% of the YRDR in 2020, was mainly distributed along the Yangtze River and Hangzhou Bay, and was sporadically distributed in the northern plain urbanization region and the southern intermontane basin. The supply capacities of WY, HQ, and other ESs were the highest, lowest, and moderate, respectively. (3) CPB was similar to that for the 1 km grid scale in 2010, differentiated into the LCHB in 2000, and further differentiated into the WYB in 2020. (4) IEB was similar to that for the 1 km grid scale in all years.

At the county level, three ESBs were identified in all years. (1) CPB covered 43.4%, 41.0%, and 39.1% of the YRDR in 2000, 2010, and 2020, respectively, and was mainly distributed in the northern plain area. CP had the highest supply capacity; CS and WY had moderate supply capacities; and SC, HQ, and LR had the lowest supply capacities. (2) IEB covered 31.7%, 30.0%, and 29.4% of the study area in 2000, 2010, and 2020, respectively, and was mainly distributed in the western Dabie Mountains and southern hilly areas. Except for CP, the supply capacity of the other ESs was very high, reflecting the comprehensive supply capacity of regulating, supporting, and cultural services. (3) The leisure recreation bundle (LRB) covered 24.9%, 29.0%, and 31.5% of the YRDR in 2000, 2010, and 2020, respectively, and was mainly distributed in cities and surrounding areas. The supply capacity of all ESs was very low, and that of LR was relatively high, which was related to the many tourist attractions and leisure places in urbanized areas.

3.4. Socioecological Drivers of ESs

In this work, GeoDetector was used to calculate the explanatory power of the natural socioeconomic drivers of each ES to determine the main driving factors and analyze how they changed with time and scale (Table 4, Table 5 and Table 6). The study revealed the following:

(1). CP was driven the most by P_C at all scales and years, and its impact continued to increase over time. In addition, CP was driven mainly by SOC and PRE at the grid scale and by NTL and PD at the county scale.

(2). CS was driven mainly by P_F, DEM, PRE, SOC, and Sand at the 1 km grid scale; by PRE, DEM, NDVI, and Sand at the 10 km grid scale; and by P_B, PD, NTL, P_F, and DEM at the county scale. At the 1 km grid scale, P_F had strong explanatory power for CS, whereas at the 10 km grid and county scales, the explanatory power of P_F was greatly weakened. The explanatory power of DEM for CS decreased with increasing scale. The driving forces of PD and NTL to CS gradually increased with time at the county scale.

(3). SC was driven mainly by DEM, P_F, PRE, and SOC at all scales and years. The q value of DEM was between 0.68 and 0.85, the q value of P_F was between 0.66 and 0.78, the q value of PRE was between 0.49 and 0.71, and the q value of SOC was between 0.44 and 0.55. At the county scale, SC was also driven by P_B (0.50 < q < 0.68).

(4). WY was driven mainly by PRE (0.44 < q < 0.70) in 2000 and 2010 at the 1 km and 10 km grid scales and was driven by TSR in 2010. WY was driven by P_B, PD, and NTL in all years at the county level.

(5). HQ was driven mainly by P_C, P_B, and P_F in all years at the 1 km and 10 km grid scales and was also driven by PD, DEM, and SOC in all years at the 10 km grid scale. At the county scale, HQ in all years was driven mainly by P_B, P_F, PD, DEM, and AET. The explanatory power of P_F and P_B for HQ increased with increasing scale, whereas that of P_C first increased but then decreased.

(6). In all years, LR was driven mainly by P_F, SOC, DEM, and P_C at the 1 km and 10 km grid scales and by SOC, P_C, PRE, and P_F at the county scale. The explanatory power of SOC for LR increased with increasing scale.

In general, the main driving factors of the different ESs exhibited significant differences at the different scales and years. The main drivers of some ESs were robust in terms of scale and time (e.g., CP and SC), some changed over time (e.g., WY), and some differed at different scales.

4. Discussion

4.1. Complexity of TOS Among ESs

Exploring the complex TOS among ESs is the basis for robust ecosystem planning and management [6]. The interactions among ESs often reflect conflicts and contradictions related to land use in a region [11]. For example, the supply of other ESs in regions with high crop production is often low, and forests and grasslands provide more regulating and supporting services than other types of land do [31]. In addition, owing to the low evapotranspiration of impervious surfaces, the supply of WY in urbanized areas is relatively high. In this study, there were obvious trade-offs between CP and other ESs (Figure 4) because cultivated land tends to yield more crops, but its ability to provide regulating, supporting and cultural services is weak [32,33]. Therefore, regional ecological planning managers need to pay attention to the trade-offs between CP and other ESs and coordinate and optimize the layout and configuration of farmland and other natural ecosystems to ensure the sustainability of human well-being in the future [6].

This study revealed that the CS–HQ, CS–LR, SC–WY, and HQ–LR pairs had synergistic effects in most regions. In general, areas with good habitat quality tend to have high carbon storage, and good ecological environments are also destinations for vacations and leisure. Areas with good soil conservation have increased water yields due to the reduction in soil erosion. These synergies among ESs are consistent with previous research results [6,34]. However, along the Yangtze River, there were high trade-offs among CS, SC, WY, HQ, and LR, which may be attributed to the high urbanization rate along the Yangtze River and the frequent and complex effects of human activities on the ecosystem, resulting in landscape fragmentation and thus damaging the structure and function of the ecosystem. Therefore, during the process of regional urbanization, we should consider protecting the ecological environment and ensure that humans and nature coexist harmoniously by optimizing the layout of cities and towns and the construction of ecological cities.

The interrelationships among ESs often differ according to scale. In this study, the CP–CS and CP–SC pairs showed trade-offs at the 1 km and 10 km grid scales and synergies at the county scale; however, while the CP–HQ pair showed trade-offs at the grid scale, CP and HQ were not correlated at the county scale. Moreover, the CS–LR and SC–LR pairs showed synergies at the grid scale, while CS and LR and SC and LR were weakly correlated at the county scale. In addition, the synergies of the CS–WY and SC–WY pairs became increasingly stronger as the scale increased, which was supported by Zheng et al. (2020) [35]. The differences in the direction and intensity of ES interactions at different scales highlight the complexity of TOS among ESs [6]. The effect of scale on the TOS of ESs varies with the specific environment and is related mainly to regional natural and socioeconomic conditions.

4.2. Driving Mechanism of ESs

Identifying the impacts of socioecological factors on ESs can aid in understanding the driving mechanism of ES changes, help formulate corresponding ecological planning and management strategies, and provide support for improving the versatility of regional ESs and promoting regional sustainable development [19]. First, land use was intrinsically linked to all ESs. In this study, P_C had a significant effect on CP and HQ; P_F had a significant positive effect on CS, SC, HQ, and LR; and P_B had a significant effect on CS, SC, WY, and HQ at the county scale. This shows that land use and landscape composition are the basis for the formation of ecosystems and their services. Second, climate factors (especially rainfall) had significant effects on CP, CS, SC, and WY at the grid scale, whereas Liu et al. (2024) reported that the climate in the Yellow River Basin had the strongest impact on SC, CP, and WY at the subbasin scale [8]. This may be because the climate in the Yellow River Basin was dry and had greater spatial homogeneity, whereas the climate in the YRDR was warm and humid and had more obvious spatial heterogeneity. In addition, biophysical conditions also had a great impact on the ESs. DEM significantly affected CS and SC, and SOC also had a strong effect on CP, SC, and LR. Finally, PD and NTL, which represent the degree of urbanization, affected mainly CP, CS, and WY at the county scale.

From the perspective of multiscale effects, the main drivers of ESs varied with scale. Except for WY, all the ESs were driven mainly by landscape composition (especially P_C, P_F, and P_B). They were also driven mainly by DEM, SOC, and PRE at the grid scale and by PD and NTL at the county scale. Notably, WY was driven mainly by PRE and TSR at the 1 km and 10 km grid scales, whereas WY was driven mainly by P_B, PD, and NTL at the county scale. A change in spatial scale may alter the interactions between ESs and driving factors and the intensity of these interactions [6,36]. Therefore, policy and planning managers should focus on improving soil and allocating water and heat resources reasonably at the grid scale, pay attention to population distribution and economic layout at the county scale, and plan land use reasonably at all scales. Research on the driving mechanisms of ESs at various scales is very important for the formulation of cross-scale differentiated sustainable ES management strategies.

4.3. Management Insights Based on ESBs in the YRDR

Ecological function zoning management based on ESBs can improve and optimize regional ecosystem service functions and maximize ecological benefits. As shown in Figure 6, the CPB represents an ecologically fragile area due to the strong effect of agricultural production activities [17]. While maintaining the red lines (also known as “ecological conservation red lines”, which refer to spatial boundaries and environmental management limits that must be strictly observed in ecological conservation services, environmental quality and safety controls, and natural resource utilization [37]) of cultivated land, actively creating high standard farmland and ensuring regional food security, land consolidation should be implemented to reduce farmland fragmentation and develop efficient and intensive ecological agriculture. The IEB has a complete ecological structure, complete ecological functions, and excellent habitat quality, which play important roles in maintaining regional ecological security. This approach is based mainly on the forest ecological service function to support and regulate the quality of the entire region’s ecological environment. As hard ecological red lines, development and construction are strictly prohibited. Ecological public welfare forests, water source protection areas, and biodiversity protection networks can be built to further protect the ecological environment. In addition, ecological compensation will help to change the local primitive agricultural production mode and develop ecological forestry and leisure tourism according to the local conditions. The LCHB is an ecologically sensitive area of lakes and coasts and is characterized mainly by the synergy of HQ and LR. Managers should demarcate ecological protection red lines and establish water source protection areas to ensure the safety of drinking water sources. For the LRB, while leisure tourism is developing, the construction of environmental protection infrastructure should also be considered. Through the construction of greenways, urban parks, green buildings, and other landscape facilities in a city, a green city with moderate living space can be built. For the WYB, we suggest performing water quality monitoring and water ecological restoration in the region and strengthening the protection of water sources [16]. Taking appropriate ecological restoration measures in areas with poor WY, such as the use of nature-based solutions, can improve water resource conservation and utilization capacity [5]. In particular, the CPB and IEB were identified at three scales, which constitute the basic ecological zoning patterns of the region. Integrating the spatiotemporal dynamic characteristics of ESBs into regional planning and management will deepen the understanding of the multiple functions of ESBs in regional planning and effectively prevent conflicts among stakeholders.

4.4. Limitations and Prospects of Research

This study was subject to several limitations. First, we quantified six key ESs (CP, CS, SC, WY, HQ, and LR) only for the main ecological risks of the YRDR at present, which may have concealed other important ESs (such as flood control, crop pollination, water purification, and localism). In the future, other ESs can be included to enhance the comprehensiveness of the assessment of ESs. Second, although the InVEST model is a popular method for evaluating ESs at present, it uses limited parameters to simulate complex ecological functions, which may lead to discrepancies between the quantified results of the model and the actual ecosystem services [11]. Owing to the lack of measurement data, verifying the results of ESs was difficult, which may have led to some uncertainty in the results [19]. In future research, measurement data can be collected from field surveys to further optimize the parameters of the InVEST model to improve the accuracy of the model estimation results. Third, we studied only the supply capacity of ESs, but it may change greatly due to different human needs [12]. The supply/demand matching and flow mechanism of ESs among different regions may have greater impacts on the relationship between ESs and stakeholders. In future research, we will further explore the interactions among ESs and the related driving mechanisms from the aspects of supply/demand and flow to obtain more comprehensive and valuable insights for optimizing regional ecological management and ensuring human well-being. Finally, although we revealed the overall impacts of various drivers on ESs through GeoDetector, the driving mechanism of spatial heterogeneity and complex driving paths is still unclear. In the future, we should use a combination of Bayesian networks and structural equation models with spatial prediction models to further explore the spatial effects and mechanisms of various driving factors.

5. Conclusions

This study explored the multiscale spatiotemporal dynamics and socioecological drivers of the interactions among the ecosystem services (ESs). First, the total supply of ESs changed significantly over the past 20 years, with crop production, soil conservation, and water yield values showing increasing trends and carbon storage, habitat quality, and leisure recreation values showing decreasing trends. The distribution of ESs exhibited spatial heterogeneity but remained essentially stable over time. Second, the overall trade-off effects were mainly between crop production and carbon storage, soil conservation, habitat quality, and leisure recreation, while the spatial distribution of trade-offs/synergies among ES pairs exhibited different patterns, reflecting the complex interactions among cultivated land, human activities, and natural ecosystems in the study area. Third, the ecosystem service bundles identified at different temporal and spatial scales differed, reflecting the dynamic evolution and spatial differences in ecosystem functions and services. In addition, the main socioecological drivers of ESs differed at different scales and across years. Overall, land use had the strongest explanatory power for all ESs, climate factors and biophysical conditions had great impacts on ESs at the grid scale, and population density and night light remote sensing also had significant impacts on crop production, carbon storage, and water yield at the county scale. This research reveals that different spatial planning and ecological management strategies need to be implemented at different spatiotemporal scales. This study provides support for maintaining regional ecological functions, improving human well-being, and achieving sustainable development.

Author Contributions

Conceptualization, Z.Z. and Y.C.; methodology, Z.Z.; software, C.Y.; validation, Z.Z. and Y.C.; formal analysis, Z.Z. and Y.C.; investigation, Z.Z.; resources, Z.Z.; data curation, Y.C.; writing—original draft preparation, Z.Z.; writing—review and editing, Y.C.; visualization, C.Y.; supervision, Z.Z.; project administration, Z.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author.

Acknowledgments

The authors express their sincere appreciation to the anonymous reviewers for their insightful suggestions during the peer review phase, and acknowledge the editorial team’s professional coordination in advancing the manuscript evaluation workflow.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Study area.

Generate image description

Figure 2 Spatial distributions of the six ecosystem services in the YRDR from 2000 to 2020.

Generate image description

Figure 3 Annual changes in the total value of each ES in the YRDR from 2000 to 2020.

Generate image description

Figure 4 Correlations among ES pairs (** p < 0.01, and *** p < 0.001) at three spatial scales.

Generate image description

Figure 5 Spatial synergies and trade-offs of ES pairs at three spatial scales.

Generate image description

Figure 6 Spatiotemporal patterns of ESBs and composition of ESs in ESBs at three scales. Longer segments represent higher ES supplies (CPB, crop production bundle; IEB, integrated ecology bundle; LCHB, lake coastal habitat bundle; WYB, water yield bundle; LRB, leisure recreation bundle).

Generate image description

Table 1

Data types and sources.

Data Type	Format/Spatial Resolution	Data Source	Application
Land use/land cover	GeoTIFF/30 m	The 30 m annual land cover dataset and its dynamics in China from 1990 to 2019 [25]	CP, CS, SC, WY, HQ, LR, GeoDetector
Digital elevation model (DEM)	GeoTIFF/30 m	Geospatial Data Cloud (https://www.gscloud.cn/)	SC, GeoDetector
Potential evapotranspiration	GeoTIFF/1000 m	National Qinghai Tibet Plateau Science Data Center (https://data.tpdc.ac.cn)	WY
Annual total evapotranspiration	GeoTIFF/500 m	MODIS (https://modis.gsfc.nasa.gov/)	WY
Annual total precipitation	GeoTIFF/1000 m	National Qinghai Tibet Plateau Science Data Center (https://data.tpdc.ac.cn)	GeoDetector
Annual mean temperature	GeoTIFF/1000 m	National Qinghai Tibet Plateau Science Data Center (https://data.tpdc.ac.cn)	GeoDetector
Annual total solar radiation	GeoTIFF/1000 m	Geographic Data Sharing Infrastructure, global resources data cloud (www.gis5g.com)	GeoDetector
Rainfall erosivity index	GeoTIFF/30 m	Geographic remote sensing ecological network platform (www.gisrs.cn)	SC
Soil erodibility index	GeoTIFF/300 m	Geographic Data Sharing Infrastructure, global resources data cloud (www.gis5g.com)	SC
Soil organic carbon content	GeoTIFF/1000 m	ISRIC—World Soil Information (https://www.isric.org)	GeoDetector
Clay content in soil	GeoTIFF/1000 m	ISRIC—World Soil Information (https://www.isric.org)	GeoDetector
Silt content in soil	GeoTIFF/1000 m	ISRIC—World Soil Information (https://www.isric.org)	GeoDetector
Sand content in soil	GeoTIFF/1000 m	ISRIC—World Soil Information (https://www.isric.org)	GeoDetector
Volumetric coarse fragment content in soil	GeoTIFF/1000 m	ISRIC—World Soil Information (https://www.isric.org)	GeoDetector
Volumetric water content in soil	GeoTIFF/1000 m	ISRIC—World Soil Information (https://www.isric.org)	GeoDetector
Soil depth	GeoTIFF/250 m	ISRIC—World Soil Information (https://www.isric.org)	WY
Available water content in plants	GeoTIFF/250 m	ISRIC—World Soil Information (https://www.isric.org)	WY
Ecosystem carbon density	.xls file	National Ecosystem Science Data Center (https://nesdc.org.cn)	CS
Net primary productivity (NPP)	GeoTIFF/500 m	The Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov)	CP
Normalized vegetation index (NDVI)	GeoTIFF/1000 m	Geographic Data Sharing Infrastructure, global resources data cloud (www.gis5g.com)	GeoDetector
Agricultural output	.xls file	Statistical Yearbook	CP
Gross domestic product	GeoTIFF/1000 m	Resource and Environmental Science Data Platform (https://www.resdc.cn)	GeoDetector
Population density	GeoTIFF/1000 m	Oak Ridge National Laboratory (https://landscan.ornl.gov)	LR, GeoDetector
Night light remote sensing	GeoTIFF/500 m	National Earth System Science Data Center (https://www.geodata.cn)	GeoDetector

Note: CP, crop production; CS, carbon storage; SC, soil conservation; WY, water yield; HQ, habitat quality; LR, leisure recreation.

Table 2

Overview of the ecosystem services quantified in this study.

Category	Ecosystem Service	Abbreviation	Unit	Method(s)	Description
Provisioning services	Crop production	CP	ton	Measurable proxies	Yield of staple food crops
Regulating services	Carbon storage	CS	ton	Carbon storage and sequestration from the InVEST model	Amount of carbon stored by terrestrial ecosystems
	Soil conservation	SC	ton	Sediment delivery ratio from the InVEST model	Capacity of vegetation cover to retain soil
	Water yield	WY	mm	Annual water yield from the InVEST model	Annual yield of water
Supporting services	Habitat quality	HQ	Index (dimensionless)	Habitat quality from the InVEST model	Ability to provide conditions suitable for the persistence of individuals and populations by ecosystems (ranges from 0 to 1)
Cultural services	Leisure recreation	LR	Index (dimensionless)	Recreation opportunity spectrum model	Index compounded by naturalness, tourist attraction density, and population density (ranges from 1 to 100)

Table 3

Selected natural and socioeconomic factors.

Category	Factor	Abbreviation
Climate	Annual total precipitation	PRE
	Rainfall erosivity	RE
	Annual mean temperature	TEM
	Annual total solar radiation	TSR
	Annual total evapotranspiration	AET
Biophysical indicators	Elevation	DEM
	Soil organic carbon content	SOC
	Sand content in soil	Sand
	Annual maximum normalized vegetation index	NDVI
Landscape composition	Percentage of cropland	P_C
	Percentage of forest	P_F
	Percentage of water	P_W
	Percentage of built-up land	P_B
Socioeconomic factors	Gross domestic product	GDP
	Population density	PD
	Night light remote sensing	NTL

Table 4

Identification of drivers of ESs at the 1 km grid scale over time.

Factors	CP2000	CP2010	CP2020	CS2000	CS2010	CS2020	SC2000	SC2010	SC2020	WY2000	WY2010	WY2020	HQ2000	HQ2010	HQ2020	LR2000	LR2010	LR2020
P_C	0.567	0.630	0.614	0.194	0.186	0.156	0.418	0.469	0.454	0.031	0.048	0.092	0.497	0.477	0.442	0.518	0.500	0.444
P_F	0.339	0.352	0.332	0.725	0.732	0.728	0.739	0.737	0.735	0.141	0.316	0.068	0.423	0.431	0.410	0.579	0.598	0.604
P_W	0.063	0.061	0.048	0.237	0.232	0.196	0.109	0.146	0.139	0.150	0.102	0.150	0.100	0.100	0.098	0.055	0.067	0.067
P_B	0.302	0.340	0.343	0.273	0.328	0.383	0.307	0.409	0.424	0.015	0.054	0.085	0.421	0.471	0.480	0.432	0.455	0.486
PRE	0.288	0.415	0.391	0.562	0.543	0.453	0.659	0.517	0.497	0.450	0.698	0.334	0.234	0.242	0.300	0.434	0.440	0.348
RE	0.120	0.156	0.214	0.248	0.229	0.216	0.205	0.198	0.171	0.193	0.337	0.117	0.100	0.094	0.090	0.205	0.199	0.190
TEM	0.070	0.118	0.184	0.187	0.160	0.148	0.107	0.051	0.077	0.295	0.313	0.204	0.051	0.037	0.038	0.048	0.046	0.069
TSR	0.189	0.273	0.166	0.420	0.374	0.227	0.340	0.305	0.149	0.178	0.450	0.180	0.164	0.139	0.078	0.373	0.341	0.170
AET	0.116	0.163	0.141	0.366	0.426	0.273	0.193	0.218	0.141	0.055	0.124	0.040	0.136	0.139	0.147	0.288	0.238	0.183
DEM	0.320	0.321	0.303	0.712	0.711	0.701	0.689	0.738	0.753	0.134	0.283	0.109	0.394	0.388	0.375	0.508	0.527	0.530
SOC	0.264	0.365	0.409	0.520	0.521	0.520	0.527	0.501	0.519	0.183	0.324	0.173	0.282	0.269	0.258	0.459	0.458	0.451
Sand	0.207	0.233	0.245	0.476	0.492	0.490	0.383	0.466	0.495	0.092	0.210	0.084	0.246	0.242	0.234	0.236	0.246	0.251
NDVI	0.081	0.091	0.099	0.372	0.354	0.499	0.086	0.112	0.304	0.186	0.108	0.222	0.061	0.067	0.193	0.019	0.027	0.120
GDP	0.134	0.135	0.160	0.261	0.278	0.281	0.195	0.195	0.248	0.140	0.142	0.106	0.183	0.195	0.199	0.078	0.077	0.087
PD	0.330	0.365	0.249	0.260	0.205	0.195	0.245	0.317	0.298	0.132	0.046	0.119	0.344	0.373	0.369	0.291	0.346	0.282
NTL	0.010	0.017	0.047	0.044	0.097	0.155	0.009	0.041	0.102	0.032	0.009	0.090	0.030	0.081	0.167	0.003	0.010	0.046

Note: (1) Yellow represents landscape composition factors, blue represents climate factors, green represents biophysical factors, and purple represents socioeconomic factors. (2) The bold q values in the table represent the main influencing factors, whereas the red q values represent the most explanatory factors. (3) Table 3 shows abbreviations for socioecological drivers. (4) For ‘CP2000’, CP represents the ES, and 2000 represents the year. The same applies for the other abbreviations.

Table 5

Identification of drivers of ESs at the 10 km grid scale over time.

Factors	CP2000	CP2010	CP2020	CS2000	CS2010	CS2020	SC2000	SC2010	SC2020	WY2000	WY2010	WY2020	HQ2000	HQ2010	HQ2020	LR2000	LR2010	LR2020
P_C	0.610	0.676	0.676	0.155	0.173	0.254	0.493	0.540	0.553	0.025	0.081	0.157	0.547	0.515	0.487	0.474	0.499	0.511
P_F	0.396	0.435	0.441	0.394	0.435	0.438	0.665	0.731	0.780	0.080	0.296	0.065	0.440	0.445	0.448	0.556	0.577	0.583
P_W	0.164	0.137	0.117	0.254	0.240	0.211	0.241	0.289	0.317	0.168	0.115	0.116	0.185	0.194	0.192	0.045	0.063	0.077
P_B	0.373	0.364	0.386	0.186	0.226	0.257	0.408	0.482	0.534	0.029	0.123	0.038	0.462	0.471	0.488	0.316	0.304	0.325
PRE	0.406	0.533	0.434	0.540	0.555	0.414	0.714	0.650	0.552	0.504	0.696	0.358	0.282	0.317	0.344	0.475	0.482	0.422
RE	0.208	0.245	0.313	0.245	0.234	0.229	0.212	0.218	0.189	0.209	0.272	0.231	0.165	0.157	0.151	0.281	0.279	0.271
TEM	0.118	0.168	0.230	0.312	0.290	0.155	0.089	0.070	0.063	0.349	0.313	0.193	0.040	0.029	0.029	0.062	0.064	0.133
TSR	0.257	0.335	0.196	0.408	0.405	0.315	0.372	0.329	0.199	0.301	0.424	0.286	0.188	0.152	0.104	0.318	0.332	0.247
AET	0.271	0.386	0.245	0.343	0.420	0.383	0.320	0.338	0.343	0.151	0.266	0.276	0.248	0.322	0.279	0.366	0.377	0.310
DEM	0.361	0.346	0.324	0.503	0.523	0.531	0.744	0.805	0.816	0.135	0.314	0.099	0.418	0.417	0.406	0.460	0.492	0.503
SOC	0.419	0.497	0.532	0.375	0.364	0.362	0.541	0.553	0.532	0.272	0.326	0.356	0.424	0.409	0.397	0.564	0.568	0.565
Sand	0.291	0.256	0.246	0.415	0.417	0.417	0.444	0.513	0.560	0.233	0.281	0.268	0.301	0.303	0.297	0.240	0.257	0.268
NDVI	0.123	0.148	0.180	0.448	0.444	0.502	0.153	0.179	0.470	0.289	0.180	0.255	0.063	0.082	0.266	0.050	0.053	0.162
GDP	0.251	0.238	0.258	0.379	0.364	0.296	0.161	0.166	0.225	0.351	0.198	0.341	0.180	0.225	0.254	0.050	0.052	0.090
PD	0.432	0.387	0.345	0.333	0.232	0.298	0.301	0.351	0.443	0.212	0.101	0.205	0.438	0.442	0.446	0.256	0.243	0.285
NTL	0.038	0.060	0.208	0.016	0.045	0.096	0.029	0.097	0.226	0.052	0.006	0.051	0.074	0.155	0.273	0.004	0.006	0.061

Note: (1) Yellow represents landscape composition factors, blue represents climate factors, green represents biophysical factors, and purple represents socioeconomic factors. (2) The bold q values in the table represent the main influencing factors, whereas the red q values represent the most explanatory factors. (3) Table 3 shows abbreviations for socioecological drivers. (4) For ‘CP2000’, CP represents the ES, and 2000 represents the year. The same applies for the other abbreviations.

Table 6

Identification of drivers of ESs at the county scale over time.

Factors	CP2000	CP2010	CP2020	CS2000	CS2010	CS2020	SC2000	SC2010	SC2020	WY2000	WY2010	WY2020	HQ2000	HQ2010	HQ2020	LR2000	LR2010	LR2020
P_C	0.365	0.492	0.578	0.256	0.244	0.209	0.554	0.506	0.446	0.195	0.242	0.113	0.345	0.281	0.334	0.578	0.596	0.592
P_F	0.185	0.222	0.219	0.453	0.455	0.464	0.696	0.726	0.753	0.329	0.499	0.186	0.652	0.651	0.650	0.414	0.435	0.462
P_W	0.018	0.018	0.033	0.195	0.197	0.191	0.239	0.231	0.214	0.146	0.149	0.063	0.133	0.144	0.149	0.028	0.050	0.071
P_B	0.239	0.274	0.330	0.539	0.595	0.629	0.509	0.639	0.679	0.442	0.620	0.528	0.650	0.677	0.668	0.263	0.252	0.271
PRE	0.182	0.285	0.282	0.183	0.269	0.297	0.620	0.583	0.535	0.253	0.468	0.225	0.446	0.505	0.496	0.482	0.503	0.394
RE	0.105	0.126	0.175	0.083	0.086	0.087	0.228	0.222	0.176	0.106	0.183	0.059	0.290	0.276	0.272	0.277	0.278	0.285
TEM	0.151	0.215	0.303	0.085	0.068	0.063	0.028	0.014	0.008	0.071	0.026	0.058	0.029	0.024	0.029	0.145	0.164	0.242
TSR	0.140	0.212	0.105	0.189	0.163	0.140	0.347	0.291	0.184	0.200	0.281	0.135	0.341	0.326	0.270	0.314	0.414	0.228
AET	0.124	0.155	0.164	0.208	0.318	0.356	0.332	0.393	0.424	0.198	0.370	0.231	0.521	0.583	0.577	0.397	0.344	0.341
DEM	0.104	0.106	0.117	0.427	0.434	0.436	0.848	0.840	0.743	0.377	0.542	0.201	0.579	0.583	0.578	0.382	0.425	0.455
SOC	0.213	0.284	0.295	0.187	0.188	0.190	0.517	0.445	0.361	0.165	0.246	0.052	0.530	0.498	0.486	0.639	0.665	0.665
Sand	0.087	0.085	0.092	0.316	0.322	0.324	0.466	0.509	0.534	0.259	0.351	0.216	0.431	0.435	0.427	0.199	0.238	0.247
NDVI	0.152	0.241	0.254	0.351	0.427	0.557	0.187	0.211	0.484	0.317	0.310	0.371	0.190	0.211	0.420	0.115	0.142	0.182
GDP	0.083	0.095	0.130	0.109	0.024	0.122	0.154	0.102	0.208	0.083	0.017	0.064	0.225	0.108	0.314	0.106	0.089	0.108
PD	0.249	0.307	0.310	0.545	0.592	0.622	0.404	0.579	0.662	0.409	0.536	0.510	0.508	0.605	0.595	0.199	0.295	0.261
NTL	0.316	0.218	0.311	0.475	0.530	0.603	0.237	0.327	0.412	0.419	0.438	0.529	0.373	0.409	0.429	0.147	0.109	0.134

Note: (1) Yellow represents landscape composition factors, blue represents climate factors, green represents biophysical factors, and purple represents socioeconomic factors. (2) The bold q values in the table represent the main influencing factors, whereas the red q values represent the most explanatory factors. (3) Table 3 shows abbreviations for socioecological drivers. (4) For ‘CP2000’, CP represents the ES, and 2000 represents the year. The same applies for the other abbreviations.

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/su17167200/s1: Figure S1. Land use/land cover of the YRDR in 2000, 2010, and 2020; ecosystem service (ES) assessments.

References

1. Costanza, R.; d’Arge, R.; De Groot, R.; Farber, S.; Grasso, M.; Hannon, B.; Limburg, K.; Naeem, S.; O’neill, R.V.; Paruelo, J. The value of the world’s ecosystem services and natural capital. Nature; 1997; 387, pp. 253-260. [DOI: https://dx.doi.org/10.1038/387253a0]

2. Millennium Ecosystem Assessment (MA). Ecosystems and Human Well-Being; 1st ed. Island Press: Washington, DC, USA, 2005.

3. Wang, J.; Zhai, T.; Lin, Y.; Kong, X.; He, T. Spatial imbalance and changes in supply and demand of ecosystem services in China. Sci. Total Environ.; 2019; 657, pp. 781-791. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.12.080]

4. Tong, L.; Liang, X.; Zhang, J.; Geng, Y.; Geng, T.; Shi, J. Ecosystem service trade-off and synergy relationship and its driving factor analysis based on Bayesian belief network: A case study of the Loess Plateau in northern Shaanxi Province. Acta Ecol. Sin.; 2023; 43, pp. 6758-6771. (In Chinese)

5. Xu, J.; Wang, S.; Xiao, Y.; Xie, G.; Lei, G. Mapping the spatiotemporal heterogeneity of ecosystem service relationships and bundles in Ningxia. J. Clean. Prod.; 2021; 294, 126216. [DOI: https://dx.doi.org/10.1016/j.jclepro.2021.126216]

6. Xia, H.; Yuan, S.; Prishchepov, A.V. Spatial-temporal heterogeneity of ecosystem service interactions and their social-ecological drivers: Implications for spatial planning and management. Resour. Conserv. Recycl.; 2023; 189, 106767. [DOI: https://dx.doi.org/10.1016/j.resconrec.2022.106767]

7. Shen, J.; Li, S.; Liu, L.; Liang, Z.; Wang, Y.; Wang, H.; Wu, S. Uncovering the relationships between ecosystem services and social-ecological drivers at different spatial scales in the Beijing-Tianjin-Hebei region. J. Clean. Prod.; 2021; 290, 125193. [DOI: https://dx.doi.org/10.1016/j.jclepro.2020.125193]

8. Liu, Q.; Qiao, J.; Li, M.; Huang, M. Spatiotemporal heterogeneity of ecosystem service interactions and their drivers at different spatial scales in the Yellow River Basin. Sci. Total Environ.; 2024; 908, 168486. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2023.168486]

9. Cord, A.F.; Bartkowski, B.; Beckmann, M.; Dittrich, A.; Hermans-Neumann, K.; Kaim, A.; Lienhoop, N.; Locher-Krause, K.; Priess, J.; Schröter-Schlaack, C. . Towards systematic analyses of ecosystem service trade-offs and synergies: Main concepts, methods and the road ahead. Ecosyst. Serv.; 2017; 28, pp. 264-272. [DOI: https://dx.doi.org/10.1016/j.ecoser.2017.07.012]

10. Kareiva, P.; Watts, S.; McDonald, R.; Boucher, T. Domesticated nature: Shaping landscapes and ecosystems for human welfare. Science; 2007; 316, pp. 1866-1869. [DOI: https://dx.doi.org/10.1126/science.1140170]

11. Wang, K.; Gao, J.; Liu, C.; Zhang, Y.; Wang, C. Understanding the effects of socio-ecological factors on trade-offs and synergies among ecosystem services to support urban sustainable management: A case study of Beijing, China. Sustain. Cities Soc.; 2023; 100, 105024. [DOI: https://dx.doi.org/10.1016/j.scs.2023.105024]

12. Huang, Y.; Wu, J. Spatial and temporal driving mechanisms of ecosystem service trade-off/synergy in national key urban agglomerations: A case study of the Yangtze River Delta urban agglomeration in China. Ecol. Indic.; 2023; 154, 110800. [DOI: https://dx.doi.org/10.1016/j.ecolind.2023.110800]

13. Zhang, T.; Zhang, S.; Cao, Q.; Wang, H.; Li, Y. The spatiotemporal dynamics of ecosystem services bundles and the social-economic-ecological drivers in the Yellow River Delta region. Ecol. Indic.; 2022; 135, 108573. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.108573]

14. Baró, F.; Gómez-Baggethun, E.; Haase, D. Ecosystem service bundles along the urban-rural gradient: Insights for landscape planning and management. Ecosyst. Serv.; 2017; 24, pp. 147-159. [DOI: https://dx.doi.org/10.1016/j.ecoser.2017.02.021]

15. Liu, D.; Chen, H.; Li, T.; Zhang, H.; Geng, Y. Spatiotemporal differentiation of village ecosystem service bundles in the loess hilly and gully region and terrain gradient analysis. Prog. Geogr.; 2022; 41, pp. 670-681. [DOI: https://dx.doi.org/10.18306/dlkxjz.2022.04.011]

16. Yan, X.; Li, X.; Liu, C.; Li, J.; Zhong, J. Scales and Historical Evolution: Methods to Reveal the Relationships between Ecosystem Service Bundles and Socio-Ecological Drivers—A Case Study of Dalian City, China. Int. J. Environ. Res. Public Health; 2022; 19, 11766. [DOI: https://dx.doi.org/10.3390/ijerph191811766]

17. Mo, W.; Zhao, Y.; Yang, N.; Xu, Z. Ecological function zoning based on ecosystem service bundles and trade-offs: A study of Dongjiang Lake Basin, China. Environ. Sci. Pollut. Res.; 2023; 30, pp. 40388-40404. [DOI: https://dx.doi.org/10.1007/s11356-022-24782-z]

18. Huang, F.; Zuo, L.; Gao, J.; Jiang, Y.; Du, F.; Zhang, Y. Exploring the driving factors of trade-offs and synergies among ecological functional zones based on ecosystem service bundles. Ecol. Indic.; 2023; 146, 109827. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.109827]

19. He, L.; Xie, Z.; Wu, H.; Liu, Z.; Zheng, B.; Wan, W. Exploring the interrelations and driving factors among typical ecosystem services in the Yangtze River economic belt, China. J. Environ. Manag.; 2024; 351, 119794. [DOI: https://dx.doi.org/10.1016/j.jenvman.2023.119794]

20. Zhang, X.; Han, R.; Yang, S.; Yang, Y.; Tang, X.; Qu, W. Identification of bundles and driving factors of ecosystem services at multiple scales in the eastern China region. Ecol. Indic.; 2024; 158, 111378. [DOI: https://dx.doi.org/10.1016/j.ecolind.2023.111378]

21. Zhang, J.; Guo, W.; Cheng, C.; Tang, Z.; Qi, L. Trade-offs and driving factors of multiple ecosystem services and bundles under spatiotemporal changes in the Danjiangkou Basin. China. Ecol. Indic.; 2022; 144, 109550. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.109550]

22. Huang, J.; Zheng, F.; Dong, X.; Wang, X. Exploring the complex trade-offs and synergies among ecosystem services in the Tibet autonomous region. J. Clean. Prod.; 2023; 384, 135483. [DOI: https://dx.doi.org/10.1016/j.jclepro.2022.135483]

23. Wang, J.; Xu, C. Geodetector: Principle and prospective. Acta Geogr. Sin.; 2017; 72, pp. 116-134. [DOI: https://dx.doi.org/10.11821/dlxb201701010]

24. National Bureau of Statistics. Available online: https://data.stats.gov.cn/index.htm (accessed on 3 July 2025).

25. Yang, J.; Huang, X. The 30 m annual land cover dataset and its dynamics in China from 1990 to 2019. Earth Syst. Sci. Data; 2021; 13, pp. 3907-3925. [DOI: https://dx.doi.org/10.5194/essd-13-3907-2021]

26. Chen, T.; Feng, Z.; Zhao, H.; Wu, K. Identification of ecosystem service bundles and driving factors in Beijing and its surrounding areas. Sci. Total Environ.; 2020; 711, 134687. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2019.134687]

27. Chang, B.; Chen, B.; Chen, W.; Xu, S.; He, X.; Yao, J.; Huang, Y. Analysis of tradeoff and synergy of ecosystem services and driving forces in urban agglomerations in Northern China. Ecol. Indic.; 2024; 165, 112210. [DOI: https://dx.doi.org/10.1016/j.ecolind.2024.112210]

28. Dittrich, A.; Seppelt, R.; Vaclavik, T.; Cord, A.F. Integrating ecosystem service bundles and socio-environmental conditions—A national scale analysis from Germany. Ecosyst. Serv.; 2017; 28, pp. 273-282. [DOI: https://dx.doi.org/10.1016/j.ecoser.2017.08.007]

29. Lyu, R.; Clarke, K.C.; Zhang, J.; Feng, J.; Jia, X.; Li, J. Dynamics of spatial relationships among ecosystem services and their determinants: Implications for land use system reform in Northwestern China. Land Use Policy; 2021; 102, 105231. [DOI: https://dx.doi.org/10.1016/j.landusepol.2020.105231]

30. Wang, J.; Zhang, T.; Fu, B. A measure of spatial stratified heterogeneity. Ecol. Indic.; 2016; 67, pp. 250-256. [DOI: https://dx.doi.org/10.1016/j.ecolind.2016.02.052]

31. Wang, S.; Hu, M.; Wang, Y.; Xia, B. Dynamics of ecosystem services in response to urbanization across temporal and spatial scales in a mega metropolitan area. Sustain. Cities Soc.; 2022; 77, 103561. [DOI: https://dx.doi.org/10.1016/j.scs.2021.103561]

32. Boeraeve, F.; Dendoncker, N.; Cornelis, J.-T.; Degrune, F.; Dufrene, M. Contribution of agroecological farming systems to the delivery of ecosystem services. J. Environ. Manage.; 2020; 260, 109576. [DOI: https://dx.doi.org/10.1016/j.jenvman.2019.109576]

33. Cao, Y.; Li, G.; Tian, Y.; Fang, X.; Li, Y.; Tan, Y. Linking ecosystem services trade-offs, bundles and hotspot identification with cropland management in the coastal Hangzhou Bay area of China. Land Use Policy; 2020; 97, 104689. [DOI: https://dx.doi.org/10.1016/j.landusepol.2020.104689]

34. Yang, G.; Ge, Y.; Xue, H.; Yang, W.; Shi, Y.; Peng, C.; Du, Y.; Fan, X.; Ren, Y.; Chang, J. Using ecosystem service bundles to detect trade-offs and synergies across urban–rural complexes. Landsc. Urban Plan.; 2015; 136, pp. 110-121. [DOI: https://dx.doi.org/10.1016/j.landurbplan.2014.12.006]

35. Zheng, D.; Wang, Y.; Hao, S.; Xu, W.; Lv, L.; Yu, S. Spatial-temporal variation and tradeoffs/synergies analysis on multiple ecosystem services: A case study in the Three-River Headwaters region of China. Ecol. Indic.; 2020; 116, 106494. [DOI: https://dx.doi.org/10.1016/j.ecolind.2020.106494]

36. Sun, X.; Wu, J.; Tang, H.; Yang, P. An urban hierarchy-based approach integrating ecosystem services into multiscale sustainable land use planning: The case of China. Resour. Conserv. Recycl.; 2022; 178, 106097. [DOI: https://dx.doi.org/10.1016/j.resconrec.2021.106097]

37. China Key Words. Available online: http://keywords.china.org.cn/2015-09/07/content_36521200.html (accessed on 25 July 2025).