Ecosystem Services and Their Driving Forces in the Middle Reaches of the Yangtze River Urban Agglomerations, China

The impact of human activities on ecosystems can be measured by ecosystem services. The study of ecosystem services is an essential part of coupled human and natural systems. However, there is limited understanding about the driving forces of ecosystem services, especially from a spatial perspective. This study attempts to fill the gap by examining the driving forces of ecosystem services with an integrated spatial approach. The results indicate that more than US$430 billion of ecosystem services value (ESV) is produced annually in the Middle Reaches of the Yangtze River Urban Agglomerations (MRYRUA), with forestland providing the largest proportion of total ESV (≥75%) and hydrological regulation function accounting for the largest proportion of total ESV (≥15%). The average ESV in the surrounding areas is obviously higher than those in the metropolitan areas, in the plains areas, and along major traffic routes. Spatial dependence and spatial spillover effects were observed in the ecosystem services in the MRYRUA. Spatial regression results indicate that road density, proportion of developed land, and river density are negatively associated with ecosystem services, while distance to a socioeconomic center, proportion of forestland land, elevation, and precipitation are positively associated with ecosystem services. The findings in this study suggest that these driving factors and the spillover effect should be taken into consideration in ecosystem protection and land-use policymaking in urban agglomerations.


Introduction
Urban agglomerations in China have become the main spatial carrier of urbanization [1]. Human activities in urban agglomerations have increasingly intensified land-use/land-cover change (LULCC) and exacerbated the evolution of ecosystem services [2,3]. Clearly identifying the evolutionary mechanisms of ecosystem services in urban agglomerations carries great significance for policy implications in land-use planning and ecosystem management. However, there is limited understanding about the determinants of ecosystem services from a spatial perspective. Within this context, it is crucial to characterize the spatial pattern of ecosystem services and explore their driving forces in urban agglomerations.
More specifically, an accurate assessment of ecosystem services is the foundation for studying the evolutionary mechanisms of ecosystem services, yet the assessment of ecosystem services and the However, we live on an increasingly human-dominated planet, and the impact of human activities on the terrestrial ecosystem has dramatically exceeded those of the physical elements [24]. The human factors affecting the evolution of ecosystem services are both direct and indirect. The chief direct factors are land-use change [25], afforestation [26], deforestation [27], and land reclamation [28], among others. The primary indirect factors include urbanization [29,30], economic growth [16], population migration [31,32], traffic accessibility [33], and land-use policies [34].
Urbanization is a complex process, as are economic agglomeration, population agglomeration, and the increased demand for developed land and ecological land in urban areas [35]. Urbanization-led LULCC has profoundly transformed ecosystem services [3] and the proportion of developed land should be an essential consideration in the evolution of ecosystem services. High population density and intense economic activities in urban areas have resulted in large-scale modification of ecosystem services [36]. The construction of infrastructure, such as highways and railways, promotes economic growth and population migration [37], thereby affecting LULCC and the provisioning capacity of ecosystem services. Proximity factors, such as distance to city centers, are closely associated with other socioeconomic factors that should also be considered [38].

Study Area and Data
The MRYRUA in this study consists of three provinces: Hunan, Hubei, and Jiangxi ( Figure 1). The MRYRUA is located in the transitional zone between China's second and third steps. The types of landforms in the study area are complex and diverse, with significant spatial differences. The overall terrain is high in the west and low in the north. There are many types of landforms, including plains, hills, and mountains in the MRYRUA. The MRYRUA is surrounded by the Wu Mountains and Xuefeng Mountains in the west, Nanling Mountains in the south, and Wuyi Mountains in the east. The area is rich in farmland, forestland, water resources, and biodiversity, all of which play an important role in grain-producing and ecosystem services provisioning in China. With the Wuhan megalopolis, Changsha-Zhuzou-Xiangtan urban agglomerations, and Poyang Lake city group as its development centers, the MRYRUA has become an emerging national-level urban agglomeration in China. Rapid development in the MRYRUA greatly promoted the land-use transition and the deterioration of ecosystem services. During the period of 1995-2015, deforestation, afforestation, land reclamation, and returning sloping farmland to forestland occurred frequently in the MRYRUA, causing a dramatic effect on regional ecosystem services. In the context of the important strategic position of the MRYRUA and the environmental problems, it is necessary to study the supply capacity of the ecosystem services and their driving mechanism.
The 30-m resolution LULCC dataset was downloaded from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn) [39]. Landsat TM/ETM+ remote sensing imagery was the primary source of data for the LULCC dataset. Using the human-computer interactive interpretation method, Liu et al. constructed the national LULCC dataset in China at 5-year intervals [39]. According to specific research needs, the land uses are divided into six first-level types-farmland, forestland, grassland, water area, construction land, and unused land. To ensure the quality of the LULCC datasets, nationwide field surveys were conducted; the overall accuracy was over 90% [39][40][41]. The 30-m resolution LULCC dataset was downloaded from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn) [39]. Landsat TM/ETM+ remote sensing imagery was the primary source of data for the LULCC dataset. Using the human-computer interactive interpretation method, Liu et al. constructed the national LULCC dataset in China at 5-year intervals [39]. According to specific research needs, the land uses are divided into six first-level types-farmland, forestland, grassland, water area, construction land, and unused land. To ensure the quality of the LULCC datasets, nationwide field surveys were conducted; the overall accuracy was over 90% [39][40][41].

Dependent Variables
The theoretical framework for the measurement of ecosystem services value (ESV) proposed by Costanza et al. (1997) has greatly promoted the study of the assessment of ESV globally [8]. Xie et al. (2003Xie et al. ( , 2008 revised the ecosystem services classification and equivalent table based on the expertise of ≥700 ecologists in China using the method of Costanza et al. (1997) [42,43]. Specifically, ecosystem services were reclassified into four main types and nine subtypes. The equivalent table of ESV was modified with the concept of equivalent value per unit area. The economic value of grain production per unit area of farmland was assumed to be 1. The equivalent value per unit area of other ecosystem services was identified based on their relative importance to grain production of farmland. The ESV equivalent factor per unit area was defined as equal to 1/7 of the average economic value of grain production per unit area of farmland. Chen et al. (2019a) calculated the ESV equivalent value ($US344.927/(hm 2 a)) based on the grain yield data and the grain price data in the MRYRUA [44] (Eq. (1)). Ecosystem services have been localized based on biomass, but biomass was not completely positively correlated, particularly in the water areas and wetlands [42]. However, the water area and wetlands in the MRYRUA are vast, contributing more than 10% of the total ESV [44] Chen et al.

Dependent Variables
The theoretical framework for the measurement of ecosystem services value (ESV) proposed by Costanza et al. (1997) has greatly promoted the study of the assessment of ESV globally [8]. Xie et al. (2003Xie et al. ( , 2008 revised the ecosystem services classification and equivalent table based on the expertise of ≥700 ecologists in China using the method of Costanza et al. (1997) [42,43]. Specifically, ecosystem services were reclassified into four main types and nine subtypes. The equivalent table of ESV was modified with the concept of equivalent value per unit area. The economic value of grain production per unit area of farmland was assumed to be 1. The equivalent value per unit area of other ecosystem services was identified based on their relative importance to grain production of farmland. The ESV equivalent factor per unit area was defined as equal to 1/7 of the average economic value of grain production per unit area of farmland. Chen et al. (2019a) calculated the ESV equivalent value ($US344.927/(hm 2 a)) based on the grain yield data and the grain price data in the MRYRUA [44] (Equation (1)). Ecosystem services have been localized based on biomass, but biomass was not completely positively correlated, particularly in the water areas and wetlands [42]. However, the water area and wetlands in the MRYRUA are vast, contributing more than 10% of the total ESV [44] where EV is the economic value of grain production per unit area of farmland (ESV equivalent value; dollars/(hm 2 · a)); h r is the sown area of the rth grain crop (hm 2 ); p r is the average price of the rth grain crop ($/ton); q r is the per unit area yield of the rth grain crop (ton/hm 2 ); M is the sown area of all the grain crops (hm 2 ); and ESV corrected indicates the corrected ESV. VCI indicates the average vegetation condition index computed from the normalized difference vegetation index of the farmland; VCI indicates the annual average VCI; VCI k is the average VCI of the kth county; EV ij is the jth category of ESV equivalent value for the ith land use type; and LUC i represents the area of the ith land-use type.  Table 1). Several important trends and transitions between different land-use types are also worth noting during 1995-2015. The most significant transitions happen between farmland and forestland: 5906.48 km 2 of forestland converted to farmland, and 5459.57 km 2 of farmland converted to forestland. For farmland and construction land, 5001.00 km 2 of farmland were converted to construction land, while only 781.48 km 2 of construction land were converted to farmland ( Table 2). A significant imbalance can be found in the amount of cultivated land transferred in and transferred out during 1995-2005 and 2005-2015. The possible reason is that rapid urbanization took over a large amount of farmland, resulting in a rapid decline in the proportion of farmland. Additionally, since the implementation of the Farmland Balance Policy, the same amount of farmland should be reclaimed to make up for the farmland occupied by construction land, resulting in a large amount of forestland being converted to farmland. At the same time, due to the policy of returning farmland to forestland and the implementation of a series of ecological engineering projects, a large amount of farmland was converted into forestland during the study period.

ESV in the MRYRUA from 1995 to 2015
The average ESVs in the MRYRUA in 1995, 2005, and 2015 were US$7724.03/hm 2 · a, US$7817.47/hm 2 · a, and US$7877.02/hm 2 · a, respectively, documenting a gradually increasing trend during the period studied. Forestland provided more than 75% of the ESV in the MRYRUA, while farmland and water areas provided more than 10% of the total ESV ( Table 3) In terms of spatial distribution, the spatial pattern of the average ESV in the MRYRUA is relatively stable, and the low-value areas of the average ESV are distributed primarily in the Jianghan Plain, the Poyang Lake Plain, and the Dongting Lake Plain, especially in the key cities and surrounding areas ( Figure 2). The high-value areas of the average ESV are mainly south of the Dabie Mountains and Wushan in the western area of Hubei, the Xuefeng Mountains in the central and western part of Hunan Province and Nanling in the southern part of Hunan Province, the Wuyi Mountains in the eastern region of Jiangxi Province, and the Luoxiao Mountains between Jiangxi Province and Hunan Province. The supply capacity of ESV provided by ecosystems varied significantly due to the significant differences in natural conditions and socioeconomic development levels in mountainous and plain areas.
the central and western part of Hunan Province and Nanling in the southern part of Hunan Province, the Wuyi Mountains in the eastern region of Jiangxi Province, and the Luoxiao Mountains between Jiangxi Province and Hunan Province. The supply capacity of ESV provided by ecosystems varied significantly due to the significant differences in natural conditions and socioeconomic development levels in mountainous and plain areas.

Independent Variables
Based on the driving factors found in previous studies that impact ecosystem services, for this study we selected population density, road density, and proximity to represent human factors, and selected elevation, slope, temperature, precipitation, and river density to represent physical factors [16][17][18]46]. The proportion of forestland and the proportion of developed land were chosen to represent the land-use pressure [28]. Because a large number of variables may lead to a multicollinearity problem, a variance inflation factor (VIF) test was performed on the independent variables. After the VIF test, a set of factors was selected as independent variables (Table 4). The spatial interaction and spatial diffusion assume that the attributes at one location are interdependent with those at other locations [47]. Global spatial autocorrelation and local spatial autocorrelation are commonly used to measure their spatial relationships. The regional spatial dependence of the average ESV was measured with the global Moran's I index and local Moran's I Index (LISA). Specifically, global spatial autocorrelation was used to verify the overall aggregation of the specific phenomena or attribute values in the spatial distribution of the entire region, while the local Moran's I was used to measure the degree of spatial difference between a certain area and its surrounding areas. The scatter plot of Moran's I provided a clear spatial clustering pattern of average ESV. The first and third quadrants in the scatter plot were, respectively, the high-high and low-low types; and the second and fourth quadrants were, respectively, the low-high and high-low types. The LISA map helps in understanding whether any of the spatial patterns are significant. Based on the statistical hypothesis testing, the Z score ≥ 1.96 or Z score ≤ -1.96 showed that the spatial autocorrelation was statistically significant at the 5% level. GeoDa095i software (University of Chicago, Chicago, IL, USA) was employed in this study to test the spatial autocorrelation of the county-level average ESV in the MRYRUA in the years 1995, 2005, and 2015 using the queen's contiguity weight method [48]. The equations are as follows: where AESV i , AESV j are the spatial observations in position i and j, respectively; AESV is the average value of ecosystem services; and W ij is the spatial weight connection matrix (i, j = 1,2,3 ... n). E(I) and Var(I) are, respectively, the expected value and the variance of Moran's I. The global Moran's I value generally ranged from -1 to 1. Moran's I values > 0 indicate positive spatial autocorrelation, Moran's I values < 0 indicate negative spatial autocorrelation, and Moran's I approaching 0 indicates a random distribution pattern.

Spatial Regression Analysis
The spatial effects of the independent variables have often been neglected in previous research about ecosystem services, and it is commonly assumed that the effects are spatially independent and identically distributed. Using an ordinary least squares (OLS) model (Equation (4)) may lead to an incomplete explanation of the regression results because it ignores the potential spatial effects [49,50]. In this study, we use three spatial regression models, including the spatial lag model (SLM) (Equation (5)), spatial error model (SEM) (Equation (6)), and spatial error model with lag dependence (SEMLD) (Equation (7)) [37,50]. The SLM assumes spatial autocorrelation occurs in the dependent variable, emphasizes the neighborhood effect, and considers the phenomenon of spatial diffusion (the spillover effect) in the dependent variable across geographic units [51]. The SEM focuses on the neglected and unobserved spatial interdependencies among variables. The SEMLD is a spatial autoregressive model augmented by adding the spatially lagged dependent variable to the SEM. All the models were conducted in GeoDa095i at the county level. We identified the performance of the four models based on the log-likelihood value, the Akaike information criterion (AIC), and the Schwarz criterion (SC). The equations are given in Equations (6) through (9): AESV t = X t β + ε, ε = λW 2 ε + ξ AESV t = X t β + ρW 1 AESV t + ε, ε = λW 2 ε + ξ where AESV t is the matrix of dependent variable in time t; X t is an n × k independent variables matrix in year t; n is the number of study units; k is the number of independent variables; β is a vector of coefficients of X t indicating the influence level of the independent variables on the dependent variables; ρ is a spatial lag parameter; ε is a random error term vector; λ is a spatial error parameter; and W 1 and W 2 are spatial weight matrices for the lag term and the error term, respectively.

Spatial Autocorrelation
The exploratory spatial data analysis method was used to further study the spatial distribution of aggregation and abnormalities of average ESV and the relationship of average ESV with the adjacent areas in the MRYRUA.  (I)) did not change significantly, and p was significant at 0.01%, indicating that the average ESV in the MRYRUA remained stable in magnitude and in spatial distribution pattern. The counties of high-high type were concentrated primarily in the mountainous areas (i.e., Wu Mountains in western Hubei, Xuefeng Mountains in western Hunan). The low-low type was distributed chiefly in the central areas of the Wuhan metropolis, Changsha-Zhuzhou-Xiangtan Metropolis, and Poyang Lake city group, which are surrounding counties of major cities. The distribution of the high-low type was apparent in the counties surrounding Nanchang, there were few low-high types distributed discretely during the study period, and the overall LISA clustering pattern of average ESV did not change significantly during the study period (Figure 3).

Spatial Regression
The Moran's I analysis indicates that spatial regression models for spatial estimation and verification should be established. After OLS estimation, a Moran's I index residual test was performed on the residuals; the Moran's I coefficients in 1995, 2005, and 2015 were 0.400, 0.426, and 0.398, respectively. The p-value was 0.001, indicating strong spatial autocorrelation among the residuals. The spatial dependence diagnostics by the OLS model further showed that statistically significant spatial lag and spatial error terms existed in the model residuals ( Table 5). The OLS estimation in the classical linear regression model may have errors in model design with no consideration of spatial autocorrelation. Based on the measures of model fit to the data-AIC, SC, and log-likelihood (Tables 5 and 6)-the performance of the spatial regression models was better fitted to the data than that of the OLS model. Among the three spatial regression models, the SEMLD had the best fit to the data, with the lowest AIC and SC values and highest log-likelihood.
The regression results of the SEMLD are shown in Table 6. It was unanticipated that population density was insignificant in all the models from 1995 to 2015. The possible reason is that other determinants, such as the proportion of developed land, have advantages over population density,

Spatial Regression
The Moran's I analysis indicates that spatial regression models for spatial estimation and verification should be established. After OLS estimation, a Moran's I index residual test was performed on the residuals; the Moran's I coefficients in 1995, 2005, and 2015 were 0.400, 0.426, and 0.398, respectively. The p-value was 0.001, indicating strong spatial autocorrelation among the residuals. The spatial dependence diagnostics by the OLS model further showed that statistically significant spatial lag and spatial error terms existed in the model residuals ( Table 5). The OLS estimation in the classical linear regression model may have errors in model design with no consideration of spatial autocorrelation. Based on the measures of model fit to the data-AIC, SC, and log-likelihood (Tables 5  and 6)-the performance of the spatial regression models was better fitted to the data than that of the OLS model. Among the three spatial regression models, the SEMLD had the best fit to the data, with the lowest AIC and SC values and highest log-likelihood.
The regression results of the SEMLD are shown in Table 6. It was unanticipated that population density was insignificant in all the models from 1995 to 2015. The possible reason is that other determinants, such as the proportion of developed land, have advantages over population density, and the impacts of population density were not obvious during these periods. Similar results can be found in Hu et al. (2015) [28]. There was a negative spatial association between road density and ecosystem services, but not in all models. Different levels of road density had different impacts on ecosystem services. Similar results can be found in other studies-i.e., that an increase in the density of roads is associated with a decrease in ecosystem services by leading to landscape fragmentation and strengthening socioeconomic activities [52,53].
Distance to a socioeconomic center proved to be an important spatial determinant for ecosystem services, and there was a positive relationship between ecosystem services and distance to a socioeconomic center. Elevation exhibited a significant positive association with ecosystem services in the years studied, indicating that counties at higher elevation tend to have higher ecosystem services. River density and ecosystem services were negatively associated in the OLS and SLM models. Numerous studies have provided important evidence that settlements tend to be formed closer to water bodies or rivers [54,55], which inevitably damages the ecosystem. In an additional Pearson's correlation analysis of river density and construction land, a significant negative correlation was The negative association between the proportion of developed land and ecosystem services indicates that counties with a higher proportion of developed land tended to have lower ecosystem services. An increase in the proportion of developed land means an increase in the consumption of ecosystem services related to human activities, and the literature has shown that an increase in the proportion of developed land leads to ecosystem deterioration [28,56,57]. A statistically significant positive association can be observed between the percentage of forestland and ecosystem services. Using SEMLD for interpretation, 1% of the percentage of forestland increase corresponds to a 0.410%, 0.340%, and 0.319% increase of average ESV, in 1995, 2005, and 2015, respectively. Spatial lag terms in the SEMLD were statistically significant in 1995, 2005, and 2015, indicating that the spatial spillover effect of ecosystem services was common in the MRYRUA. The ecosystem services in one county were not only closely related to the individual elements (physical and social factors) but also to other neighborhood factors (ecosystem services in adjacent counties). Each 1% increase in average ESV in a neighboring county correlated to 0.276%, 0.134%, and 0.121% increases in each county in 1995, 2005, and 2015, respectively. Moreover, the spatial error terms in the SEMLD were statistically significant in 1995, 2005, and 2015, documenting that ecosystem services were impacted not only by the driving forces studied in this study but were also influenced by other, omitted variables.

A Summary of the Findings
National strategies like the Central China Rising Strategy, the Yangtze River Economic Belt Strategy, and the Triangle of Central China Strategy have greatly promoted urbanization and industrialization in the MRYRUA. As a consequence, a considerable amount of farmland has been converted into construction land. In response, a series of land-use policies, such as the Basic Farmland Protection Regulation, the Farmland Balance Policy, and the increasing versus decreasing balance policy, have been implemented in an attempt to control urban expansion and prevent the decline of farmland [58,59]. However, because of the huge economic benefits gap between farmland and construction land, the targets for construction land expansion control and farmland protection were not met [60]. Additionally, land-use activities, such as deforestation and cultivation, as well as forestry engineering in the shelterbelt program in the upper and middle reaches of the Yangtze River basin, resulted in frequent conversion between farmland and forestland.
Spatial determinants identified in this study further showed negative associations among ecosystem services and road density, the proportion of developed land, and river density, while the distance to a socioeconomic center, proportion of forestland land, elevation, and precipitation were positively associated with ecosystem services. Previous studies have evidenced that developed land tends to expand along road networks [38]. The road network affects the flow of elements and resources within and among regions, and the circulation and accumulation of labor, capital, technology factors, and resources affect population redistribution and land development, directly or indirectly. The geographic characteristic (such as distance to a socioeconomic center) has an evident impact on land-use change [61]. The availability and accessibility of multiple resources in urban centers exacerbate the agglomeration of socioeconomic activities. Ecosystem services around the central counties of cities and surrounding counties are becoming increasingly important. Previous studies provide empirical evidence that human settlements are more likely distributed near water bodies because of invaluable natural amenities [55]. River density is therefore negatively associated with ecosystem services. The positive association between the proportion of forestland and ecosystem services shows that an increase in forestland promotes an increase in ecosystem services. Reforestation is an effective means for improving ecosystem services [28]. The ecosystem services in the MRYRUA display significant spatial autocorrelation and spillover effects, indicating that the ecological conservation and land-use policies in a single area cannot completely solve the problem of regional ecological issues and land-use problems. The prevention and control of ecological issues needs effective coordination across counties, especially in urban agglomerations.

Policy Implications
In practice, there always exist conflicts among various plans (e.g., economic and social development planning, environmental protection planning, and land-use planning). For example, land-use planning tends to ignore ecological impacts. In response, multiple planning integration efforts have been promoted in China to balance the differences among various spatial plans and avoid unfavorable outcomes. For example, the three designation zones (agricultural, ecological, and urban) and three lines (urban development boundary, ecological red line, and permanent basic farmland red line) are to alleviate the issues of the decline of farmland and the deterioration of ecosystem services due to the expansion of construction land [62]. Therefore, it is necessary to strategically optimize national spatial planning and ecosystem services management. Empirical research shows that the evolution of regional ecosystem services is impacted not only by physical and human factors but also by the charateristics of the neighboring areas. It is noteworthy that ecosystem services in the MRYRUA exhibit significant spatial autocorrelation and spatial spillover features. In other words, regional ecosystem services are affected by ecological regulation and industrial structure as well as by neighboring counties and more distant counties; the ecological conservation regulation policy for a single area cannot completely solve the problem of regional ecological issues [10]. The effective control of ecological issues requires intercounty joint efforts.

Limitations and Future Research Directions
The revised benefit transfer method and an integrated spatial regression were used to measure the spatial features and the influence mechanism of ecosystem services. However, this study focused only on the supply capacity of the ecosystem services; their demand capacity was an essential aspect of ecosystem services that also should be considered [9]. The study examined only the impact mechanism of total ESV; future research could pay attention to the driving forces of regulating services, supplying services, supporting services, and cultural services. Furthermore, an integrated cross-sectional spatial regression approach was used to identify the spatial determinants of ecosystem services; dynamic spatial panel data models that blend the inter-individual differences and intra-individual dynamics could be considered in future research [63].

Conclusions
In this study, we examined the spatiotemporal pattern of the ESV in the MRYRUA from 1995 to 2015 using a revised benefit transfer method. Then we identified the spatial autocorrelation of ESV with the global Moran's I index and local Moran's I index. We ultimately analyzed the spatial determinants of average ESV using an integrated spatial regression approach and found a gradually increasing trend in ESV. Forestland provided the greatest proportion of ESV (more than 75%) in the MRYRUA, and farmland and water areas provided more than 10% of the total ESV. The hydrological regulation function accounted for the largest proportion of the total ESV (more than 15%) among all the ecosystem functions. We observed evident spatial dependence and spatial spillover effects of average ESV in the MRYRUA from 1995 to 2015. The average ESV in an area was impacted not only by its factors but also by the characteristics of its neighboring areas due to the spatial spillover effects. The spatial regression results showed that road density, river density, and proportion of developed land were negatively associated with the average ESV, while elevation and precipitation were positively associated with the average ESV. Distance to a socioeconomic center was also found to be an important spatial determinant for average ESV. The results in this study provided important implications for the cross-regional collaborative management of ecosystem services and sustainable land use.