Analysis of Ecosystem Service Contribution and Identiﬁcation of Trade-Off/Synergy Relationship for Ecosystem Regulation in the Dabie Mountains of Western Anhui Province, China

: The study of tradeoffs/synergies among ecosystem services (ESs) is highly signiﬁcant for land-use planning and regional ecosystem optimization. Land-use change and topographic factors have important implications for ESs. Strengthening the comparative analysis of the capacity of ESs provided by different land-use types in speciﬁc regions, studying the topographic gradient effects of ecosystem service trade-offs/synergies with slope changes, and identifying the dominate trade-off/synergy relationship among ESs will help us to carry out ecosystem regulation according to local conditions through land-use layout optimization at a ﬁne scale. Our research site was located in the Dabie Mountains of western Anhui Province, China (DBM), where, based on the InVEST software, R language, self-organizing maps (SOM), and GeoDA, the temporal and spatial variations of ﬁve typical ESs, including food supply, soil retention, water yield, carbon storage, and biodiversity maintenance from 2005 to 2020, were analyzed, and spatial distributions of the different ESs clusters were also recognized by using the SOM method. Moreover, the impacts of land-use type and slope on ESs, and the characteristics of trade-offs/synergies among the ﬁve ESs, were discussed. Results showed, ﬁrstly, that the total values of ESs showed a changing trend of “three increases and two decreases” from 2005 to 2020. Among the ESs, food supply, soil retention and water yield showed upward trends, with annual growth rates of 2.83%, 6.50% and 2.98%, respectively, whereas carbon storage and biodiversity maintenance showed downward trends, with annual decline rates of 0.03% and 0.07%, respectively. Second, the results showed that the Moran’s I index of the total ESs was 0.3995 in 2005 and 0.4305 in 2020, respectively, indicating that they had signiﬁcant spatial clustering characteristics. The Low-Low clustering regions with reduced changes were mainly in the central and northern parts of the study area, whereas the High-High clustering regions with increased changes were found


Introduction
Ecosystem services (ESs) are life-support products and services that can be obtained directly or indirectly through the structure, processes, and functions of the ecosystem to maintain the natural environmental conditions and utilities on which human beings depend [1]. Land-use change is an important driver of ecosystem service transformation [2]. Currently, nearly 60% of ESs have been partially or fully diminished due to irrational land use [3]. Rational land-use planning utilizes the spatial relationship between ESs towards synergistic optimization [4,5]. The scientific assessment of ESs has gradually become an important prerequisite for the formulation of ecological protection policy objectives, the determination of ecological compensation models, and the precise regulation and control of ecosystems. Studies have shown that ESs have provided a scientific basis for ecosystem management at the national and regional scales [6,7]. At present, research on ESs mainly focuses on the calculation of ES value, the assessment of ecosystem service supply and demand [8], ecosystem service flows [9], and the collaborative research on the relationship between ESs and human wellbeing [10,11]. ES assessment and trade-offs/synergies research provide comprehensive and practical methods for addressing the environmental challenges faced in development [12]. They also offer an important practical field for research on the core proposition of the relationship between man and land in geography [13,14].
The diversity of ecosystem service types, the spatial heterogeneity of ESs, and human management lead to nonlinear and complex interactions among ESs [15,16]. Thus, uncovering the complex interactions and identifying the dominating types of trade-offs and synergies among ESs constitute the basis and premise for sustainable ecosystem management [17][18][19]. Presently, many methods, such as correlation analysis [20] and Bayesian belief networks [21], have been applied to uncover the pairs of trade-offs/synergies relationships among ESs. Meanwhile, spatial autocorrelation analysis [22], the spatial mapping of ESs evaluation [23], and the geographically weighted regression (GWR) [24] are methods commonly used to study the spatial distribution characteristics of ESs. Furthermore, a growing number of recent studies have used the structural equation model (SEM) [25] and self-organizing maps (SOM) [26] to analyze ES bundles, and some new indicators, such as [27] the trade-off/synergy criteria (TSC) and trade-off/synergy index (TSI), have been designed to measure the strength of trade-offs and synergies between ESs. The method of multiple scenarios [28] has also been used to simulate the impact of land-use change on ESs. However, some of these recent studies were still limited to static assessments at a single time or spatial scale [26]. Realizing the trade-off/synergy effect of ESs is fundamental to achieving efficient ecosystem management and improved human wellbeing. Yet current research addressing the driving mechanisms behind the formation of ecosystem service relationships remains limited [29]. Strengthening the study of ecosystem service bundles can effectively reflect the interdependence of ESs, which is important for ecosystem regulation decisions [30]. ES model and method integration and algorithm design are important ways to research trade-off mechanisms and multi-dimensional and complex ESs [31,32]. Therefore, it is necessary to research the spatiotemporal variation and multi-scale characterization of ecosystem service interactions. For regions with complex terrain, it is valuable to pay attention to the sensitivity of the interaction between ESs to terrain differences. Identifying the dominant relationship between ecosystem service trade-offs/synergies is also key to the precise regulation of regional ecosystem functions. In addition, the analysis of ecosystem service contribution rates for different land types can clarify the composition of ecosystem service contributions for different land types, thereby contributing to land-use management and optimization, and ultimately effectively improving the value of regional ESs. The DBM is famous as a region with history and culture and rich tourism resources which enjoys a winning reputation internationally. In 2018, UNESCO approved the Dabie Mountain Geopark as a world geopark. We all know that targeted poverty alleviation is a great livelihood project in China and in the world, and the DBM is a core area for taking targeted measures in poverty alleviation which has received wide attention from scholars because of its unique location, geological landscape, and history and culture in recent years [33,34]. Therefore, it is important to strengthen research on ecosystem service assessment, contribution sources, and trade-off/synergy interactions in the DBM for ecosystem regulation and the enhancement of ecosystem service values here. With the complex topography of the DBM, the optimization of ecosystem function is important for providing ecological security construction and resident wellbeing in the old revolutionary areas. Research on ecological security in the DBM is extremely important for the optimization of regional ecosystem function, the sustainable development of the ecological economy, and the construction of ecological civilization [35,36]. In this study, multi-source data such as statistical data and spatial raster data were mainly used to quantitatively calculate five typical ESs in the study area by using the NDVI matching method and InVEST model. The spatiotemporal changes of ESs, the contribution of different land-use types to ESs, the trade-offs/synergies among ESs, the topographic effect, and the ES bundles were then analyzed. The contribution of different land-use types to ESs was revealed, the dominant pairs of ESs were identified, and the different ES bundles were also recognized, providing an important basis for the quantitative assessment of ESs and ecosystem regulation in the DBM. Overall, this study aimed to (1) quantitatively evaluate five typical ESs and analyze their temporal and spatial characteristics, (2) analyze the impact characteristics of land use and topography on ESs, and also discuss the spatial pattern and trade-offs/synergies among ESs, and (3) recognize ES bundles based on the five typical ESs by using the SOM method. Finally, this study revealed the spatiotemporal heterogeneity of ESs and the slop effect, identified the dominant types of trade-offs/synergies among ESs, and recognized the ES bundles and their spatial distribution in the study area. The ultimate goal was to provide a decision-making basis for ecosystem function optimization and ecosystem restoration and management in the DBM to a certain extent.

Study Area
The Dabie Mountains region (DBM) is located at the intersection of Hubei, Henan, and Anhui Provinces. This location is an important ecological function area in Central China and a significant ecological barrier in the middle and lower reaches of the Yangtze River [37,38]. The DBM of western Anhui Province is located in the northeast of the DBM region and the west of Anhui Province, China. It belongs to the core of the DBM, which is rich in tourism resources and is also one of the famous red scenic spots in China, as shown in Figure 1. The study area mainly includes five counties and two districts, namely Huoqiu County, Shou County, Jinzhai County, Huoshan County, Shucheng County, Yu'an District, and Jin'an District, covering 18,430 km 2 , being dominated by middle and low mountains and hills, and belonging to the north subtropical warm and humid monsoon climate region. The annual average temperature ranges from 14-16 • C, the annual average rainfall from 1100-1450 mm, and the annual average sunshine hours from 2000-2200 h. The area is in the watershed of the Yangtze River and Huaihe River Basins in Anhui Province. Based on the Second National Remote Sensing Survey of Soil and Water Loss, ecological problems, such as soil and water loss in counties and districts, are prominent, making this a typical ecologically fragile area [39]. The DBM of western Anhui Province plays a core part of the key ecological function within the mountain range, and it is also the main water source for the Pishihang Irrigation area. As a famous reservoir region in China, it contains six large reservoirs, including Foziling, Mozitan, Bailianya, Meishan, Xianghongdian, and Longhekou, as well as the head water control project of the Pishihang Irrigation area [40].

Data Source and Processing
The types of data sources mainly included grid, vector, and text data in this study: land-use data with a spatial resolution of 30 × 30 m in 2005 and 2020, a digital elevation model (DEM) with a spatial resolution of 30 m, soil data in this study generated from a Chinese soil type distribution map and the Chinese soil dataset HWSD_ China_ Subset_ V1.1, precipitation of meteorological data in 2005 and 2020 collected from 44 meteorological stations in the study area and surrounding areas, annual actual evapotranspiration data from a geographic remote-sensing ecological network platform, normalized vegetation index (NDVI) data extracted from the remote sensing images (these were downloaded from geospatial data cloud website), and socioeconomic data mainly collected from relevant statistical yearbooks from 2006 to 2021. The descriptions of the data are shown in Table 1. The five typical ESs and their formulas are shown below, including food supply, soil retention, water yield, carbon storage, and biodiversity maintenance in the list of ESs.
(a) Food supply (fs): Food supply service is an important service, in an agricultural ecosystem, which plays a vital role in human survival and regional development [41]. The type and quantity of food supply selected for this study were mainly derived from the food types and corresponding outputs provided in relevant statistical yearbooks. In terms of calculation methods, we allocated the output of grain according to the grid unit of cultivated land, the output of meat according to the grid unit of grassland, and the output of aquatic products according to the grid numbers of the water bodies. The specific methods were as follows. First, using conditional functions in ArcGIS, information on cultivated land, grassland, and water bodies were extracted from the land-use data of the study area. Second, for cultivated land and grassland, we allocated the outputs based on the ratio of the NDVI value of each grid to the total NDVI value of different land types, and ultimately allocate the outputs of grain, meat, and milk products to the cultivated land and grassland grid units, respectively [42]. For the allocation of aquatic products, we used an average allocation method based on the total number of water body pixels-that is, the ratio of the output of aquatic products to the total number of water grid cells-to evenly allocate the number of aquatic products. Finally, all land-use types were merged with the support of ArcGIS, resulting in a grid map of food supply services. The formula used for NDVI allocation is as follows: where G i is the output of grain, meat and milk products distributed by cell i, G sum is the total output of grain, meat, and milk in the DBM of western Anhui Province, NDVI i is the normalized vegetation index of cell i, and NDVI sm is the sum of the NDVI values of cropland or grassland in the study area.
(b) Soil retention (sr): In this study, we used the sediment delivery ratio model, which is a submodel in the InVEST software, to estimate the amount of potential soil erosion (RKLS) and actual soil erosion (USLE) of each cell unit based on the Universal Soil Loss Equation (USLE), and then calculate the difference between these values to obtain soil conservation (SR) data. The formula used is as follows: where Sr is the amount of soil retention (t·hm −2 ·a −1 ) and R is the rainfall erosivity (MJ·mm·hm −2 ·h −1 ·a −1 ). R was calculated based on the monthly average precipitation and annual average precipitation in the study area by using Wischmeier's formula [43]. K is the soil erodibility (t·hm 2 ·h·hm −2 ·MJ −1 ·mm −1 ), which was calculated by using the EPIC formula [44] on the basis of the relevant soil data, LS is the slope length gradient factor, which was extracted by the watershed slope length and slope factor calculation method [45] based on the DEM, P is the factor of water and soil conservation measures, which was determined by the slope index α [46], C is the vegetation coverage and crop management factor, which was calculated by the relationship formula between vegetation coverage and C value [47], and the value range is from 0 to 1. (c) Water yield (wy): The water yield submodel in the InVEST software was applied to estimate the water yield. The water yield module in the InVEST model is mainly based on the Budyko curve and water balance principle [44]. The water yield capacity is equal to the precipitation on each grid in the study area minus the actual evapotranspiration (including soil evapotranspiration and plant evapotranspiration). Compared with other hydrological models, the InVEST model has advantages in the spatial expression and visualization of calculation results [48,49]. The formula is as follows: where Y (x) is the annual water yield (mm) at cell x, AET (x) is the annual actual evapotranspiration (mm) at cell x, and P (x) is the annual average precipitation (mm) at cell x. (d) Carbon storage (cs): Carbon storage mainly includes four basic carbon pools: aboveground biomass, underground biomass, soil, and dead organic matter. The ES of carbon storage is to estimate the amount of carbon currently stored in the landscape on the basis of the average carbon density of carbon pools of different land-use types in the study area multiplied by the area of each land-use type [50,51]. In this study, the carbon storage module in the InVEST software was used to evaluate the amount of carbon currently in the study area. The formula is as follows: where C _total is the total carbon storage (t), C _above is the carbon storage (t) of the aboveground biomass, C _below is the carbon storage (t) of the belowground biomass, C _soil refers to the carbon storage in soil, and C _dead refers to the dead organic carbon storage.
(e) Biodiversity maintenance (bd): The InVEST software calculates habitat quality by combining landscape sensitivity and external threat intensity and evaluates biodiversity service function based on habitat quality [52]. The Habitat Quality module in InVEST was applied to evaluate the habitat quality index to reflect the function of providing biodiversity services (the InVEST model assumes that areas with good habitat quality have high biodiversity). The habitat quality index is a dimensionless comprehensive index to reflect the habitat suitability and habitat degradation degree of different land-use types in the study area. The calculation formula and model parameter table of the habitat quality index were derived from relevant references [44,[53][54][55]].

Spatial Expression and Measure of Trade-Offs and Synergies among ESs
(a) Expression of spatial heterogeneity for typical ESs: For better comparison and calculation in the next step, those five aforementioned ESs were standardized primarily to eliminate this dimension. In this paper, the Fuzzy membership method, which is mainly based on the overlay toolset of the spatial analysis toolbox in the ArcGIS 10.2 software, was used to normalize five typical ESs in the study area. After making them dimensionless, the Raster Calculator tool was used to overlay the normalized grid map of ESs to obtain the grid map of the total value of ESs in the study area. The amounts of change in ESs from 2005 to 2020 were obtained by calculating the differences between the ESs in 2005 and 2020. In addition, the Zonal Statistics function in ArcGIS was used to calculate the amounts of the five ESs in the different administrative units, and then the characteristics of spatial distribution and spatial heterogeneity of the total value, and the amounts of change in ESs were discussed by calculating the local autocorrelation, i.e., the Moran's I index, among ESs and drawing the map of local indicators of spatial association (LISA) based on GeoDA 10.0.
(b) Measure of trade-offs/synergies among ESs: The trade-offs and synergies among the five ESs were explored by using the spatial analysis tools in ArcGIS and the correlation analysis function in the R language software. First, random points were created by using the Create Random Points tool in ArcGIS to take raster values from the supply maps of various ESs to these points. Second, the attribute data of random points with the values of the ESs were imported into the Excel software. Finally, the data were imported into the R language, and the maps of the lower triangle and the maps of correlation heat with the coefficient matrix were produced based on the R language, which can measure the trade-offs/synergies among the five ESs.
(c) Spatial patterns of trade-offs/synergies between ESs pairs and identification of ES bundles: In this study, we aimed to provide a clearer visualization of the spatial distribution characteristics of the trade-offs/synergies between ESs pairs. First, the mean values of ESs at the township administrative scale were computed using the Zonal Statistics function in ArcGIS 10.2. Then, we used the bivariate global and local Moran's I indices to further explore the spatial expression of trade-offs/synergies between ESs pairs [49]. Second, to classify the homogeneity of ESs, ecosystem service bundles at the scale of township administrative units were identified by using the SOM method. The SOM method changes the network parameters and structure in a self-organizing and adaptive manner by automatically finding the intrinsic laws and essential properties in samples. It is an unsupervised learning neural network method that can allocate each township administrative unit into ecosystem service bundles based on the similarities of ecosystem service co-occurrence in space [56]. In this study, the identification of ES bundles, by using the SOM method, was executed in the R language.

The Research Framework
In the study, the correlation coefficient method, exploratory spatial data analysis (ESDA), bivariate spatial autocorrelation, SOM and InVEST models, R software, ArcGIS software, and GeoDA software were applied. First, we evaluated ESs from 2005 to 2020 and analyzed their spatiotemporal patterns. Second, ecosystem service contribution structures in different land-use types, the slope effect, and trade-offs/synergies correlations between ESs were analyzed. Last, the spatial expressions of trade-off/synergy relationships between ESs pairs, ecosystem service bundles, and spatial distribution and composition structures of ecosystem service bundles were identified and analyzed. The utilized research framework is shown in Figure 2.

Spatial Patterns and ESs Changes
The five ESs, i.e., food supply, soil retention, water yield, carbon storage, and biodiversity maintenance, were standardized after the assessment and then were classified into five levels in the ArcGIS software, forming the lowest, lower, medium, higher, and highest levels, respectively. The total ESs in the DBM of western Anhui Province showed a changing trend of "three increases and two decreases." Among the five ESs, food supply, soil retention, and water yield had increased, with annual growth rates of 2.83%, 6.50% and 2.99%, respectively from 2005 to 2020. Conversely, the other two ESs, carbon storage and

Spatial Patterns and ESs Changes
The five ESs, i.e., food supply, soil retention, water yield, carbon storage, and biodiversity maintenance, were standardized after the assessment and then were classified into five levels in the ArcGIS software, forming the lowest, lower, medium, higher, and highest levels, respectively. The total ESs in the DBM of western Anhui Province showed a changing trend of "three increases and two decreases." Among the five ESs, food supply, soil retention, and water yield had increased, with annual growth rates of 2.83%, 6.50% and 2.99%, respectively from 2005 to 2020. Conversely, the other two ESs, carbon storage and biodiversity maintenance, had decreased with annual decline rates of 0.03%, and 0.07%, respectively, showing weak downward trends from 2005 to 2020. Since the DBM was situated in the national rainfall center in 2020, water yield increased from 682.22 mm in 2005 to 987.67 mm in 2020, with a growth of 44.77%. Meanwhile, soil conservation changed fastest, rising from 1.74 × 10 8 t in 2005 up to 3.44 × 10 8 t in 2020, with a growth of 97.49%; food supply also showed an upward trend, from 664.07 × 10 4 t in 2005 up to 945.67 × 10 4 t in 2020, with a growth of 42.41%. Spatially, the low-level region of food supply was mainly distributed in the southern parts of the study area, such as Jinzhai County, Huoshan County, and the south of Shucheng County. The reason for this is that the southern mountainous areas are mainly forest land and slope farmland, and hence the food supply is low, whereas the high-level food supply region was mainly distributed in the northern parts of the study area, such as Huoqiu County and Shou County. It must be specially noted that Shou County is a major grainproducing county in China. This region has excellent irrigation water quality, moderately fertile soil, and a fine climate condition; hence, the region's food supply capacity is the highest in the study area. The statistical description of food supply also showed that the maximum, mean, and variance of food supply on the cell unit had increased, indicating that the measures of cropland consolidation in the DBM of western Anhui Province had played a certain role in promoting grain production, and that the measures of cropland concentration and large-scale management had optimized the spatial layout of food supply carriers and then improved the food supply. The areas with high levels of water yield were found to be distributed in the south of the study area, whereas the low-level areas were in the north of the study area; the reasons for this are mainly related to the heavy rainfall and relatively low evaporation [57] in the southern mountainous areas (Figures 3 and 4).
Land 2023, 12, x FOR PEER REVIEW 9 of 23 major grain-producing county in China. This region has excellent irrigation water quality, moderately fertile soil, and a fine climate condition; hence, the region's food supply capacity is the highest in the study area. The statistical description of food supply also showed that the maximum, mean, and variance of food supply on the cell unit had increased, indicating that the measures of cropland consolidation in the DBM of western Anhui Province had played a certain role in promoting grain production, and that the measures of cropland concentration and large-scale management had optimized the spatial layout of food supply carriers and then improved the food supply. The areas with high levels of water yield were found to be distributed in the south of the study area, whereas the low-level areas were in the north of the study area; the reasons for this are mainly related to the heavy rainfall and relatively low evaporation [57] in the southern mountainous areas (Figures 3 and 4). From the perspective of spatial patterns, the spatial patterns of ESs had obvious characteristics. The low-level region of soil retention was mainly located in the middle and north of the study area. The high-level region of soil retention was mainly distributed in the south of the study area. This area is the intersection of nature reserves, scenic spots, and forest parks. Thus, here, the regional vegetation coverage is high, and the soil retention is relatively high. The high-level region of carbon storage was mainly located in the southern mountainous area, which has a humid climate, high vegetation coverage, and high soil organic matter content, whereas the vegetation coverage and soil organic matter content in the northern area are relatively low due to the presence of cropland and builtup land. Therefore, the carbon storage in the northern parts was generally lower than that in the south of the study area. The high-level region of biodiversity maintenance was mainly distributed in the south of the study area, whereas the low-level region was mainly located in the scope of built-up land with a low density of vegetation coverage and strong interference of human activities.

Spatiotemporal Heterogeneity of ESs
Revealing the spatial clustering characteristics of the total ESs and the temporal change of ESs in the study area will help us understand the spatial heterogeneity of ESs and perform effective regulation for ESs in practice. In this study, ArcGIS was used to derive spatial statistics on the total ESs in 2005 and 2020 and the temporal change of ESs from 2005 to 2020 according to township units, then local spatial autocorrelation analysis was carried out. The Moran's I index and the map of local indicators of spatial association (LISA) were generated, as shown in Figure 5 by using the GeoDA software. The total of ESs in 2005 and 2020 and the change amounts of ESs from 2005 to 2020 in the DBM of western Anhui Province had reasonably significant spatial heterogeneity. The results showed that the Moran's I indices of the total ESs in 2005 and 2020 were 0.3995 and 0.4305 with Z-Score over 2.58, respectively, indicating that they passed the significance test at the 0.01 level and the spatial clustering effect was highly significant. The map of LISA shows that the scope of High-High clustering and Low-Low clustering area in 2005 and 2020 were also mainly located in the central north and the south of the study area, respectively. Overall, the spatial distribution of the total ESs had still changed in the past 15 years. The number of units in the High-High region increased from 27 in 2005 to From the perspective of spatial patterns, the spatial characteristics of ESs were obvious. The low-level region of soil retention was mainly located in the middle and north of the study area. The high-level region of soil retention was mainly distributed in the south of the study area. This area is the intersection of nature reserves, scenic spots, and forest parks. Thus, here, the regional vegetation coverage is high, and the soil retention is relatively high. The high-level region of carbon storage was mainly located in the southern mountainous area, which has a humid climate, high vegetation coverage, and high soil organic matter content, whereas the vegetation coverage and soil organic matter content in the northern area are relatively low due to the presence of cropland and built-up land. Therefore, the carbon storage in the northern parts was generally lower than that in the south of the study area. The high-level region of biodiversity maintenance was mainly distributed in the south of the study area, whereas the low-level region was mainly located in the scope of built-up land with a low density of vegetation coverage and strong interference of human activities.

Spatiotemporal Heterogeneity of ESs
Revealing the spatial clustering characteristics of the total ESs and the temporal change of ESs in the study area will help us understand the spatial heterogeneity of ESs and perform effective regulation for ESs in practice. In this study, ArcGIS was used to derive spatial statistics on the total ESs in 2005 and 2020 and the temporal change of ESs from 2005 to 2020 according to township units, then local spatial autocorrelation analysis was carried out. The Moran's I index and the map of local indicators of spatial association (LISA) were generated, as shown in Figure 5 by using the GeoDA software. The total of ESs in 2005 and 2020 and the change amounts of ESs from 2005 to 2020 in the DBM of western Anhui Province had reasonably significant spatial heterogeneity. and perform effective regulation for ESs in practice. In this study, ArcGIS was use derive spatial statistics on the total ESs in 2005 and 2020 and the temporal change of from 2005 to 2020 according to township units, then local spatial autocorrelation anal was carried out. The Moran's I index and the map of local indicators of spatial associa (LISA) were generated, as shown in Figure 5 by using the GeoDA software. The tota ESs in 2005 and 2020 and the change amounts of ESs from 2005 to 2020 in the DBM western Anhui Province had reasonably significant spatial heterogeneity. The results showed that the Moran's I indices of the total ESs in 2005 and 2020 w 0.3995 and 0.4305 with Z-Score over 2.58, respectively, indicating that they passed the nificance test at the 0.01 level and the spatial clustering effect was highly significant. map of LISA shows that the scope of High-High clustering and Low-Low clustering in 2005 and 2020 were also mainly located in the central north and the south of the st area, respectively. Overall, the spatial distribution of the total ESs had still changed in past 15 years. The number of units in the High-High region increased from 27 in 200 of these over 80% were located in Yu'an District, Jin'an District, Huoqiu County, and Shou County, and the spatial distribution here formed a basin structure of "low, middle and high outside". The local autocorrelation analysis showed that the Moran's I index of the ecosystem service change was 0.6548; meanwhile, the Z-score was greater than 2.58, indicating that the spatial correlation had reached a highly significant level. The LISA cluster map showed that the Low-Low area of the ecosystem service change was mainly distributed in Huoqiu County, Shou County, and Jin'an District in the middle and north of the study area, with a total of 38 township units, whereas the High-High region was mainly located in Jinzhai County and Shucheng County in the south of the study area, with a total of 36 township units ( Figure 5). Local spatial autocorrelation analysis accurately expressed the clustering characteristic of temporal change of ESs from 2005 to 2020 in the study area, and this can provide a decision-making basis for township ecological planning and ecosystem function regulation to some extent.

Impact of Land-Use Types on the Supply of ESs
Land use/cover change is an important reason for the spatial differences between surface ESs. Different land-use types affect ESs distinctly [58]. In this study, the Condition Analysis Function in the ArcGIS software was used to extract the cells of cropland, woodland, grassland, water bodies, and built-up land, and then the cells of the five ESs in the study area were extracted in accordance with the positions of cells in different land-use types.
Taking 2020 as an example, the results showed that the amounts of ESs from cropland and woodland accounted for 37.42% and 42.70% of the total ESs in the DBM of western Anhui Province, respectively. Therefore, cropland and woodland were the main contribu- tors to the total ESs in the study area, whereas water bodies accounted for only 3.01% of the total ESs in the study area because of their relatively small area proportion of the study area. In terms of the composition of each ecosystem service, the sources of ESs were found to be significantly different, as can be seen in Figure 6. (a) The total food supply mainly came from cropland, accounting for 73.41%. (b) The total of soil retention mainly came from woodland and grassland, accounting for 58.44% and 34.29%, respectively. (c) The total of water yields mainly came from cropland, woodland, and grassland, accounting for 53.75%, 19.22%, and 12.81%, respectively. (d) The total of carbon storage mainly came from cropland and woodland, accounting for 26.91% and 64.77%, respectively. (e) Similarly, over 80% of the total biodiversity maintenance came from cropland and woodland. These results show that the supply sources of ESs in the study area were mainly cropland and woodland, whether with regards to the composition of total ESs or the composition of single ESs, as shown in Figure 6.
Analysis Function in the ArcGIS software was used to extract the cells of cropland, woodland, grassland, water bodies, and built-up land, and then the cells of the five ESs in the study area were extracted in accordance with the positions of cells in different land-use types.
Taking 2020 as an example, the results showed that the amounts of ESs from cropland and woodland accounted for 37.42% and 42.70% of the total ESs in the DBM of western Anhui Province, respectively. Therefore, cropland and woodland were the main contributors to the total ESs in the study area, whereas water bodies accounted for only 3.01% of the total ESs in the study area because of their relatively small area proportion of the study area. In terms of the composition of each ecosystem service, the sources of ESs were found to be significantly different, as can be seen in Figure 6. (a) The total food supply mainly came from cropland, accounting for 73.41%. (b) The total of soil retention mainly came from woodland and grassland, accounting for 58.44% and 34.29%, respectively. (c) The total of water yields mainly came from cropland, woodland, and grassland, accounting for 53.75%, 19.22%, and 12.81%, respectively. (d) The total of carbon storage mainly came from cropland and woodland, accounting for 26.91% and 64.77%, respectively. (e) Similarly, over 80% of the total biodiversity maintenance came from cropland and woodland. These results show that the supply sources of ESs in the study area were mainly cropland and woodland, whether with regards to the composition of total ESs or the composition of single ESs, as shown in Figure 6. From the area proportion of various land-use types, it is clear that the proportion of cropland accounted for 48.89% of the total area in 2020, making it the largest land-use type in the study area. The woodland area accounted for approximately 28.16%, the grassland area accounted for roughly 10.38%, the total area of built-up land accounted for about 7.84%, whereas the total area of water bodies accounted for only 4.72%, making this the smallest land-use type in the study area. The area proportion of the five land-use types in the study area was significantly different. Therefore, further analysis of the composition of ESs per unit area can reflect the supply capacity of different land-use types to the total of ESs in the study area. The results showed that in terms of supply capacity per unit area, the main contribution to the total ESs in the study area in 2020 was woodland and grassland, accounting for 32.42% and 26.01%, respectively, which can be seen in Figure 7. From the area proportion of various land-use types, it is clear that the proportion of cropland accounted for 48.89% of the total area in 2020, making it the largest land-use type in the study area. The woodland area accounted for approximately 28.16%, the grassland area accounted for roughly 10.38%, the total area of built-up land accounted for about 7.84%, whereas the total area of water bodies accounted for only 4.72%, making this the smallest land-use type in the study area. The area proportion of the five land-use types in the study area was significantly different. Therefore, further analysis of the composition of ESs per unit area can reflect the supply capacity of different land-use types to the total of ESs in the study area. The results showed that in terms of supply capacity per unit area, the main contribution to the total ESs in the study area in 2020 was woodland and grassland, accounting for 32.42% and 26.01%, respectively, which can be seen in Figure 7.
Based on the composition of each ecosystem service per unit area, the following were evident: cropland and grassland contributed over 70% to food supply per unit area; soil retention per unit area mainly came from woodland and grassland, accounting for 37.04% and 58.96%, respectively; built-up land was the main source of water yield per unit area, accounting for 28.66%; carbon storage per unit area mainly came from woodland, accounting for 61.47%; and woodland, grassland, and water bodies had strong supply capacities of biodiversity maintenance per unit area, accounting for 34.21%, 28.79%, and 26.57%, respectively (Figure 7). The results revealed that in view of the supply composition of the total amount of ESs and the supply capacity per unit area comprehensively, the ESs in the study area mainly came from the woodland, which accounts for 42.70% of the total amount of ESs in the study area and 32.42% of the ESs per unit area. Therefore, the contribution of the land-use type of woodland to ESs in the DBM of western Anhui Province reached approximately 1/3. Based on the composition of each ecosystem service per unit area, the following were evident: cropland and grassland contributed over 70% to food supply per unit area; soil retention per unit area mainly came from woodland and grassland, accounting for 37.04% and 58.96%, respectively; built-up land was the main source of water yield per unit area, accounting for 28.66%; carbon storage per unit area mainly came from woodland, accounting for 61.47%; and woodland, grassland, and water bodies had strong supply capacities of biodiversity maintenance per unit area, accounting for 34.21%, 28.79%, and 26.57%, respectively ( Figure 7). The results revealed that in view of the supply composition of the total amount of ESs and the supply capacity per unit area comprehensively, the ESs in the study area mainly came from the woodland, which accounts for 42.70% of the total amount of ESs in the study area and 32.42% of the ESs per unit area. Therefore, the contribution of the land-use type of woodland to ESs in the DBM of western Anhui Province reached approximately 1/3.

Impact of Slope on Trade-Offs and Synergies among ESs
The terrain is an important factor affecting land use, ecosystem, and soil erosion. It controls the redistribution of surface water, heat, and other energies, and dominates the intensity and direction of surface runoff [39]. Revealing the influencing mechanism of slope on trade-offs and synergies among ecosystem services is helpful to establish multiscale ecosystem function regulation and ecological optimization decision-making schemes [36]. Therefore, the effects of slope on trade-offs and synergies among ESs, based on the correlation analysis functions in the R language, were discussed in the study. In this study, the study area was reclassified into five grades based on slope, namely, Grade I slope area (slop1 ≤ 2°), Grade II slope area (2° < slop2 ≤ 6°), Grade III slope area (6° < slop3 ≤ 15°), Grade IV slope area (15° < slop3 ≤ 25°), and Grade V Slope area (25° < slop5) according to the technical specification for land-use status investigation by using the ArcGIS software.
Taking 2020 as an example, random points were created within the scope of different grades; then, the attribute values of the five ESs were extracted at the corresponding random points, and then the attribute values were imported into the R language to generate the maps of the lower triangle, and finally, the correlations between ESs pairs were analyzed on the basis of the maps. In a map of the lower triangle, the red round cake indicates the trade-offs with negative correlations between these pairs of ESs, and the blue round cake indicates synergies with positive correlations. The darker the color and the larger the round cake, the greater the correlation coefficient and the more significant the trade-offs and synergies between ESs. The results show that different trade-offs and synergies were

Impact of Slope on Trade-Offs and Synergies among ESs
The terrain is an important factor affecting land use, ecosystem, and soil erosion. It controls the redistribution of surface water, heat, and other energies, and dominates the intensity and direction of surface runoff [39]. Revealing the influencing mechanism of slope on trade-offs and synergies among ecosystem services is helpful to establish multi-scale ecosystem function regulation and ecological optimization decision-making schemes [36]. Therefore, the effects of slope on trade-offs and synergies among ESs, based on the correlation analysis functions in the R language, were discussed in the study. In this study, the study area was reclassified into five grades based on slope, namely, Grade I slope area (slop1 ≤ 2 • ), Grade II slope area (2 • < slop2 ≤ 6 • ), Grade III slope area (6 • < slop3 ≤ 15 • ), Grade IV slope area (15 • < slop3 ≤ 25 • ), and Grade V Slope area (25 • < slop5) according to the technical specification for land-use status investigation by using the ArcGIS software.
Taking 2020 as an example, random points were created within the scope of different grades; then, the attribute values of the five ESs were extracted at the corresponding random points, and then the attribute values were imported into the R language to generate the maps of the lower triangle, and finally, the correlations between ESs pairs were analyzed on the basis of the maps. In a map of the lower triangle, the red round cake indicates the trade-offs with negative correlations between these pairs of ESs, and the blue round cake indicates synergies with positive correlations. The darker the color and the larger the round cake, the greater the correlation coefficient and the more significant the trade-offs and synergies between ESs. The results show that different trade-offs and synergies were found among ESs in the study area based on the different slope grades, and the numbers of trade-offs were more than those of synergies on the whole. The trade-offs and synergies among ESs had different sensitivities to the slope grade changes, i.e., trade-offs between wy and bd, wy and cs were always found under five slope grades. As we can see in Figure 7, the correlation coefficients of the trade-offs between wy and bd in five slope conditions were −0.78, −0.75, −0.62, −0.65 and −0.73, respectively, while the correlation coefficients of the trade-offs between wy and cs in five slope conditions were −0.65, −0.72, −0.83, −0.89 and −0.89, respectively, with all correlation coefficients reaching significant levels, and the trade-offs became increasingly significant with increase in slope grade. Similarly, weak synergies were found between fs and wy at the slope Grade I (slop ≤ 2 • ) and II (2 • < slop ≤ 6 • ), and the synergies became more significant at the slope Grades III (6 • < slop ≤ 15 • ), IV (15 • < slop ≤ 25 • ), and V (25 • < slop), which reveal a typical slope effect. Conversely, bd and cs always showed significant synergies at different grades of slope. As we can see in Figure 8, the correlation coefficients of the synergistic relationships between bd and cs in five slope conditions were 0.76, 0.80, 0.66, 0.69 and 0.80, respectively. The relationships between fs and sr, sr and wy, sr and bd, and sr and cs at Grades III (6 • < slop ≤ 15 • ), IV (15 • < slop ≤ 25 • ), and V (25 • < slop) were unclear because their correlation coefficients were small, despite the fact that there were some trade-offs and synergies with weak correlation below the slop Grade II.
The results unveiled that when the slope increased, the correlation coefficients between most ESs increased gradually: that is, the significance of trade-offs and synergies increased with slope increase, such as the pair of ESs between wy and cs, and between cs and bd. Conversely, there remained a few correlations between ESs that were less sensitive to slope change, such as these pairs: fs and sr, wy and sr, cs and sr, and cs and bd. Ultimately, the slope effects on trade-offs and synergies were significant between typical ecosystem service pairs in the study area (Figure 8). The diversity of ESs and the identification of trade-offs and synergies at different slope levels in the study area can provide a practical decision-making basis for ecological management. For example, on a slope greater than 25 degrees, we found that the trade-off relationship between carbon storage and water yield reached the maximum and the synergies between carbon storage and biodiversity maintenance also reached the maximum. Therefore, the planning of ES demand in areas with complex terrain must be more cautious than the planning in flat areas, and thereby, based on the corresponding trade-offs and synergies, a management and regulation strategy for maximizing ES supply in the study area can be proposed. However, the influencing mechanism of the slope effect on the trade-off/synergy changes among ESs warrants further analysis. the trade-offs between wy and cs in five slope conditions were −0.65, −0.72, −0.83, − and −0.89, respectively, with all correlation coefficients reaching significant levels, an trade-offs became increasingly significant with increase in slope grade. Similarly, w synergies were found between fs and wy at the slope Grade I (slop ≤ 2°) and II (2° < ≤ 6°), and the synergies became more significant at the slope Grades III (6° < slop ≤ IV (15° < slop ≤ 25°), and V (25° < slop), which reveal a typical slope effect. Conversel and cs always showed significant synergies at different grades of slope. As we can s Figure 8, the correlation coefficients of the synergistic relationships between bd and five slope conditions were 0.76, 0.80, 0.66, 0.69 and 0.80, respectively. The relations between fs and sr, sr and wy, sr and bd, and sr and cs at Grades III (6° < slop ≤ 15° (15° < slop ≤ 25°), and V (25° < slop) were unclear because their correlation coeffic were small, despite the fact that there were some trade-offs and synergies with weak relation below the slop Grade II.
The results unveiled that when the slope increased, the correlation coefficient tween most ESs increased gradually: that is, the significance of trade-offs and syne increased with slope increase, such as the pair of ESs between wy and cs, and betwe and bd. Conversely, there remained a few correlations between ESs that were less s tive to slope change, such as these pairs: fs and sr, wy and sr, cs and sr, and cs and Ultimately, the slope effects on trade-offs and synergies were significant between ty ecosystem service pairs in the study area (Figure 8). The diversity of ESs and the ide cation of trade-offs and synergies at different slope levels in the study area can prov practical decision-making basis for ecological management. For example, on a s greater than 25 degrees, we found that the trade-off relationship between carbon sto and water yield reached the maximum and the synergies between carbon storage an odiversity maintenance also reached the maximum. Therefore, the planning of ES dem in areas with complex terrain must be more cautious than the planning in flat areas thereby, based on the corresponding trade-offs and synergies, a management and reg tion strategy for maximizing ES supply in the study area can be proposed. However influencing mechanism of the slope effect on the trade-off/synergy changes among warrants further analysis.

Temporal Characteristics of Trade-Offs/Synergies among ESs
Research on the temporal characteristics of trade-offs/synergies is beneficial to understand the law of mutual conservation between ESs to some extent. Based on the R language, correlation heat maps were generated in this study. In a scatter matrix of correlation heat maps, the histogram and kernel density curve of five services are on the main diagonal. The histogram can reflect the distribution characteristics of horizontal axis data (ESs), and the kernel density curve reflects the concentration degree of horizontal axis data distribution. The part above the main diagonal is the correlation coefficient between ESs. The part below the main diagonal is the scatter diagram and smooth-fitting curve between ESs.

Temporal Characteristics of Trade-Offs/Synergies among ESs
Research on the temporal characteristics of trade-offs/synergies is beneficial to understand the law of mutual conservation between ESs to some extent. Based on the R language, correlation heat maps were generated in this study. In a scatter matrix of correlation heat maps, the histogram and kernel density curve of five services are on the main diagonal. The histogram can reflect the distribution characteristics of horizontal axis data (ESs), and the kernel density curve reflects the concentration degree of horizontal axis data distribution. The part above the main diagonal is the correlation coefficient between ESs. The part below the main diagonal is the scatter diagram and smooth-fitting curve between ESs.
The results uncovered that the correlation coefficients between these pairs of ESs in 2005 and 2020 mostly reached significance at the 0.01 level except for these pairs of ESs between fs and sr, and fs and wy, in 2020. Synergies with positive correlation were found between sr and cs, sr and bd, and cs and bd in 2005 and 2020. Among them, the positive coefficients between cs and bd reached above 0.90 in the two periods, and the correlation between these two was the strongest. Therefore, the relationship between cs and bd was identified as the dominant synergy in the study area. In comparison, these pairs of ESs, i.e., wy and cs, wy and bd, fs and cs, fs and bd, sr and wy, and fs and sr, showed the relationships of trade-offs with negative coefficients, of which the pair of ESs including wy and cs had the strong negative correlation, changing from −0.72 in 2005 to −0.50 in 2020; this was identified as the dominant trade-off relationship in the study area (Figure 9). In addition, in recent 15 years, the number of trade-off types and synergy types in the study area had remained unchanged. Among them, the type of synergies still accounted for 40%, and the type of trade-offs accounted for 60%, but the trade-off and synergy relationships between all pairs of ESs were weakening. Among them, the maximum value

Spatial Distributions of Trade-Offs/Synergies between ESs Pairs
The study of trade-off/synergy relationships based on quantitative methods o lation coefficients can easily and conveniently identify the types of trade-offs and gies in the interrelationships between typical ESs in the DBM of western Anhui Pr as well as identify the strengths and weaknesses of correlations between pairs of E hough the correlation coefficient matrix can be used to quickly identify trade-off/s relationships between two pairs of ESs, it is weak in expressing the spatial distrib of the trade-off/synergy effects between ESs specifically. Therefore, strengthening t tial expression of trade-off/synergy relationships between ESs is beneficial for the regulation of ecosystems.
This study analyzes the spatial patterns of trade-offs and synergies of ESs ba

Spatial Distributions of Trade-Offs/Synergies between ESs Pairs
The study of trade-off/synergy relationships based on quantitative methods of correlation coefficients can easily and conveniently identify the types of trade-offs and synergies in the interrelationships between typical ESs in the DBM of western Anhui Province, as well as identify the strengths and weaknesses of correlations between pairs of ESs. Although the correlation coefficient matrix can be used to quickly identify trade-off/synergy relationships between two pairs of ESs, it is weak in expressing the spatial distributions of the trade-off/synergy effects between ESs specifically. Therefore, strengthening the spatial expression of trade-off/synergy relationships between ESs is beneficial for the precise regulation of ecosystems.
This study analyzes the spatial patterns of trade-offs and synergies of ESs based on the bivariate local Moran's I. Figure 10 shows that the spatial patterns of trade-offs and synergies between different ecosystem service pairs differed significantly, with the High-High and Low-Low values of the bivariate local Moran's I expressing the High-value and Low-value areas of synergies between ESs, respectively, and the High-Low-and Low-High -value areas of Moran's I expressing trade-offs between ESs, respectively. The results showed that the spatial distribution pattern of trade-off/synergy relationships among ESs in the study area had not changed much in the past 15 years, and presented a more stable spatial pattern. This was strongly related to the progress of conservation and development in the DBM as well as land-use changes and the level of socio-economic development in the region. In the years 2005 and 2020, sr and cs, sr and bd, and cs and bd showed synergistic relationships in many places spatially, in which the High-value synergistic areas were distributed in the south of the study area, while the Low-value synergistic areas were mainly located in the north-central part of the study area, showing significant spatial differences.
Bivariate local spatial autocorrelation clearly reveals the spatial trade-off/synergy relationships of ecosystem services, but ignores potential patterns of spatial association of ESs and does not provide insight into the interdependence between multiple ESs [59]. A bundle of ESs is a group of ESs that recur in time and space, reflecting the interdependence of those ESs [30]. SOM is an unsupervised learning neural network method that measures the similarity between ESs of different township units based on the similarity principle, and classifies township units with high similarity into the same ES bundles. In this study, the SOM method was applied to identify ES bundles and their spatial distribution at the township scale, and the results of SOM self-organizational mapping analysis can be used to tailored management measures and solutions for different ecological environments. The SOM cluster analysis of five typical ESs in the DBM of western Anhui Prov finally identified six types of ecosystem service bundles. Based on their clustering cha teristics and combined with the realities of the study area, these six types of ecosys service bundles were named 'ecological forest conservation areas', 'urban ecologic fragile areas', 'plain flood prone areas', 'ecological regulating functional areas', 'impor water reservoir area', and 'irrigation and main grain production areas' in this study (  The SOM cluster analysis of five typical ESs in the DBM of western Anhui Province finally identified six types of ecosystem service bundles. Based on their clustering characteristics and combined with the realities of the study area, these six types of ecosystem service bundles were named 'ecological forest conservation areas', 'urban ecologically fragile areas', 'plain flood prone areas', 'ecological regulating functional areas', 'important water reservoir area', and 'irrigation and main grain production areas' in this study ( Figure 11   As shown in Figure 11, bundles of similar ESs are obviously clustered in space, while the composition structure of each ecosystem service bundle has obvious differences, and the rose diagram in this paper shows the relative size of each ecosystem service in each ecological subdistrict of the DBM of western Anhui Province. In Bundle 1, the composition structure of ESs is relatively balanced, and the supply of ESs is high, mainly including sr, bd, cs, and wy, and this type of service is mainly located in the southern forest area with high topography, so it is named 'ecological forest conservation areas'. In Bundle 2, the amount of ESs is scarce, and these types of service are mainly located in the area of urban construction land, so the bundle is named 'urban ecological fragile areas'. In Bundle 3, the ecosystem service type is mainly wy, and the area is mainly located in the north-central part of the relatively gentle topography, so it is named 'plain flood prone areas'. Bundle 4 As shown in Figure 11, bundles of similar ESs are obviously clustered in space, while the composition structure of each ecosystem service bundle has obvious differences, and the rose diagram in this paper shows the relative size of each ecosystem service in each ecological subdistrict of the DBM of western Anhui Province. In Bundle 1, the composition structure of ESs is relatively balanced, and the supply of ESs is high, mainly including sr, bd, cs, and wy, and this type of service is mainly located in the southern forest area with high topography, so it is named 'ecological forest conservation areas'. In Bundle 2, the amount of ESs is scarce, and these types of service are mainly located in the area of urban construction land, so the bundle is named 'urban ecological fragile areas'. In Bundle 3, the ecosystem service type is mainly wy, and the area is mainly located in the north-central part of the relatively gentle topography, so it is named 'plain flood prone areas'. Bundle 4 is mainly located in the mountainous area in the south of the study area, and the ecosystem service types are mainly cs, bd, and sr, so the bundle is named 'ecological regulation functional areas'. In Bundle 5, the ecosystem service types are mainly wy, bd, cs, and a small amount of sr. There are five important reservoirs in the DBM, so this bundle is named 'important water source storage areas'. In Bundle 6, the ecosystem service types are mainly food supply services fs and wy, and the distribution is mainly located in the Pishihang Irrigation area and the main food production area, so it is named 'irrigation and main grain producing areas'. The ecosystem service bundles in the DBM are better identified by using the SOM method, and thus different ecological function zones can be formed. In practice, the ecosystem service bundles can be used to provide some reference for ecological regulation and planning.

Discussion
This study discussed the spatial heterogeneity of ESs in the study area to a certain extent and analyzed the impact of regional land-use change and topography on the tradeoffs/synergies among ESs. Finally, we obtained the results of the total amount and temporal and spatial patterns, trade-offs/synergies relationship characteristics, and scale dynamic changes of five typical ESs in the study area. The study showed that there is a more obvious topographic gradient difference in the trade-off synergistic relationships of ecological services in the DBM of western Anhui Province. As the slope changed, there were significant differences in trade-offs/synergies between ESs. For example, the trade-off correlations of water yield and habitat quality, and water yield and carbon storage, under the five slope conditions were highly negative and showed highly significant trade-offs, with the highest values (−0.75 and −0.89, respectively). Biodiversity maintenance and carbon storage under the five slope conditions showed synergistic correlations of 0.76, 0.80, 0.66, 0.69 and 0.80, respectively. Due to the fact that the development of the forestry economy and ecological economy industries in the DBM is more active in the mountainous areas with complex topography, food supply and water yield in the study showed a significant synergistic relationship at higher slopes, with the maximum correlation reaching 0.83 when the slope was at Grade IV (15 < slop ≤ 25 • ). In addition, the number of pairs of ESs in a trade-off relationship gradually increased with an increase in slope. Among them, significant tradeoffs were found between food supply and carbon storage and food supply and biodiversity maintenance as the slope increased to Grade III or higher. Therefore, for the DBM with complex topography, ESs should be regulated according to local conditions to enhance the overall ecosystem benefits. Decision-makers should take into account the trade-off/synergy relationship between ESs for two reasons: on the one hand, to enhance the gains of the synergistic relationship between ESs, and on the other hand, to mitigate the interaction of trade-off relationships between ESs in order to maximize the supply of ESs. According to the results of this study, woodlands and grasslands in the DBM are the main providers of ecosystem service contribution, and therefore, land-use optimization and allocation should be reasonably carried out according to the topographic distribution characteristics, such as forestry and grass, as appropriate in middle-and high-altitude areas, to improve the carbon storage capacity and habitat quality of the study area.
However, some limitations remain in the study. (a) We used data on land use and DEM with a spatial resolution of 30 m in this study. Thereby, compared with the need for refinement in the study of ESs, the accuracy of relevant parameters should have been further improved in this study. Concurrently, the research on ecosystem service tradeoffs and synergies must still strengthen the fine empirical and simulation analysis with high resolution and large scale [60]. Due to the complexity of ecological processes, relevant studies, such as those on the simulation and optimization of ecological processes, ground experiments, micro monitoring, and large-scale ecological network monitoring, should be further strengthened [49]. Research on trade-offs/synergies among ESs with micro-to large-scale includes fine experimental research under the combination of coupling multi-dimensional factors (such as slope, vegetation coverage, rainfall, soil type, land-use type, and so on), and this is also an important theoretical basis for ecological problem identification, ecological restoration, and ecological compensation in the future [61,62]. (b) Scale effect is the core difficulty of ecosystem research. Currently, ESs are widely used in the research on watersheds, the Loess Plateau, mountainous areas, agricultural and pastoral areas, oceans, and so on. The research conclusions on trade-offs and synergies also showed some differences in different regions. Limited by the data acquisition of the corresponding time series in the study area, this study only selected the rapid urbanization stage as the study period and made a scale comparative analysis on a shorter time series scale. Nonetheless, revealing the detailed characteristics of its dynamic change trend by analyzing the correlation between long-time series ecological services is easy [63,64]. We will further strengthen multi-period and multi-objective prediction and scenario simulation research [65] and carry out the comparative analysis of key time nodes and periods on a large time scale [63], combined with analysis of the important periods of social and economic development. The corresponding analysis can better reveal the coupling mechanism between ecosystem and human factors such as population, social economy, urbanization, and so on. (c) In general, the results in this paper concerning grain, oil, and vegetables based on the arable land NDVI matching method are more satisfactory, but the method of spatialization of meat products through grassland NDVI will have certain shortcomings due to the fact that the current meat output is dominated by family captive farming and factory farming. Therefore, the results of matching meat output and grassland may still have some uncertainty. Despite the shortcomings of this method, the main food sources in this study area are grains and vegetables, while the production of meat and aquatic products accounts for a small proportion of the total food production, and so the spatialized results of meat and aquatic product production have a small impact on the overall study results in terms of production, but if the study area is the main production area or an important base for meat product supply, i.e., when meat production accounts for a large proportion of the food supply, it is necessary to optimize the method according to the scale of pastures and farming patterns in the study area.
Currently, research related to mountain ESs is becoming richer and research methods are being enhanced. Coupled ecosystem service evaluation can be applied to the construction of regional ecological security patterns [66]. The results of supply and demand analysis based on the spatial dimension of ESs can also be applied to the research on the optimization of territorial spatial planning and ecological governance [67]. The results of this study can provide technical and theoretical support for ecological function zoning, ecological restoration, and precise regulation in the Dabie Mountains and similar regions. The DBM is an old revolutionary base in China. Presently, the comprehensive revitalization of rural areas and the control of relative poverty are important goals and tasks for the region's future development. Relevant research will provide a basis for the formulation of "win-win" goals and strategies for the optimization of ecosystem function and the coordinated improvement of ecological quality and resident wellbeing in the DBM.

Conclusions
By using the InVEST model and spatial autocorrelation analysis, the spatial patterns and spatiotemporal heterogeneity of five typical ESs in the DBM of western Anhui Province from 2005 to 2020 were studied, and then the impacts of land-use type and slope on the ESs were analyzed. The spatial and temporal scales of trade-offs/synergies among ESs were also discussed. Some meaningful results and findings were obtained from this study: (a) From 2005 to 2020, the total value of ESs in the DBM showed a changing trend of "three increases and two decreases". Among the ESs, only food supply, soil retention, and water yield showed upward trends, with annual growth rates of 2.83%, 6.50% and 2.98%, respectively, whereas carbon storage and biodiversity maintenance showed downward trends, with annual decline rates of 0.03% and 0.07%, respectively. Additionally, in view of the supply composition of the total amount of ESs and the supply capacity per unit area comprehensively, the ESs mainly came from the woodland, which accounts for 42.70% of the total amount of ESs and 32.42% of the ESs per unit area. Therefore, the contribution of the woodland to ESs in the DBM reached approximately 1/3. In practice, the regulation of land-use types with different total contributions and contribution rates of ESs can be targeted and strengthened. Overall, the development trend of ES evolution in the DBM is good, while to obtain the sustainable improvement of the comprehensive capacity of ES, one would need to further strengthen the monitoring of ecosystem change in the DBM and the analysis of the dynamic mechanisms of influencing factors. (b) The Moran's I indices of the total of ESs in 2005 and 2020 were 0.3995 and 0.4305, respectively, and the Moran's I index of the ecosystem service changes in the two periods also reached 0.6518, indicating that both the total of, as well as changes in, the ecosystem services had significant spatial characteristics. The Low-Low cluster regions with decreasing change were mainly in the central and northern parts of the study area, whereas the High-High cluster regions with increasing change were mainly in the south of the study area. The impact of land-use type on ESs was obvious, i.e., cropland and woodland were the main contributors to the total of ESs in the DBM, but the supply capacity of woodland per unit land area was the strongest. The analysis of spatial agglomeration characteristics reveals the spatial distribution pattern and change characteristics of ESs in the study area, helping the regulation and precise policy application of ecosystem management in practice. This study only discussed the spatial patterns and clustering characteristics of various ESs by using the spatial autocorrelation analysis method, but the reason for the basin structure of "low in the middle and high outside" formed by the changes of ESs in recent 15 years has not been deeply discussed. The complex driving mechanism of the spatial patterns of ESs should be further analyzed with regards to social and economic development. (c) The significance of trade-offs and synergies increased with slope increase. Fundamentally, the slope effect on trade-offs/synergies was significant between typical ecosystem service pairs in the study area. The number of trade-offs relationship among ESs was more than that of synergy relationships in 2005 and 2020. Among them, the relationships between carbon storage and biodiversity maintenance and water yield and carbon storage were the two dominant relationships of synergies and trade-offs among ESs in the study area. Based on the results of the cluster analysis, the DBM can be divided into six types of ecosystem service bundles, and similar ecosystem service bundles show obvious clustering in space, while there are obvious differences in the composition structure of each ES bundle. The changing trend of ES bundles was relatively stable for 15 years. The identification of ecosystem service bundles is of great significance for the division of ecological functional zones and ecological regulation in the DBM.
The identification of dominant relationships in the trade-off/synergy relationships between ESs in the DBM is important for the judgment of the regional direction of ecosystem service regulation. In practice, the identification of dominant relationships in tradeoffs/synergies of ESs, combined with the slope effect of trade-offs/synergies, can not only reveal the most important relationship categories in trade-offs/synergies of ESs, but also provide information on the spatial topographic scales at which trade-offs/synergies interactions occur. Additionally, six ES bundles were recognized, and their spatial distributions and structures were revealed in the study; this can be helpful for the precise regulation of ESs and the enhancement of the total values of ESs. These are imperative to the precise regulation of ESs and the enhancement of the total benefits of the ecosystem in DBM.  Data Availability Statement: Data sharing is not applicable to this article as no new data were created or analyzed in this study.