China’s Urban and Rural Development Signiﬁcantly Affects the Pattern of Human Appropriation of Net Primary Production

: Increasing human activities have greatly inﬂuenced the ecosystem and the use of ecological resources, and the unbalanced urban–rural development in China (urban and rural areas being two major bases of human activities) has always been accompanied by heterogeneous ecological effects. Human appropriation of net primary production (HANPP) is an integrated indicator quantifying the human domination of productivity and harvest in the biosphere. Identifying the unbalanced constraints of urban and rural development on HANPP has become necessary for improving human– land relationships. This study analyzed the spatial distribution and regional differentiations of the HANPP in China in 2015 and investigated how HANPP and its components responded to unbalanced regional urban–rural development. The results show that the total amount of HANPP was 2.68 PgC and gradually decreased from the southeast to the northwest of China in 2015, representing 60.33% of the NPP pot . In addition, HANPP luc , harvest through cropland, livestock grazing, and forestry contributed 60.70%, 29.86%, 8.53%, and 0.91%, respectively, to the total HANPP, with HANPP luc playing the dominant role in 21 provinces. There was a signiﬁcant differentiation ( p < 0.05) in the spatial distribution of HANPP (gC/m 2 ), HANPP harv (gC/m 2 ), and HANPP luc (gC/m 2 ), especially between the Huanyong Hu Line and the western–eastern part of China, fundamentally resulting from uneven regional development. In addition, biomass production–consumption decoupling existed in most regions in China, 17 provinces were identiﬁed as consumption type, and a universal positive correlation ( p < 0.05) was identiﬁed between the production–consumption ratio of occupied biomass and HANPP harv (%HANPP). Different drive mechanisms were found between urban–rural development and HANPP, and each HANPP index was more likely to be affected by urban economy (UE), rural population (RP), and rural agricultural technology (RA) in China. The higher regional average nighttime light intensity, the proportion of the built-up area, and the urban road area corresponded with a large HANPP luc value. Conversely, HANPP would decrease as the proportion of urban green spaces increased. Furthermore, HANPP (%NPP pot ) and HANPP (gC/m 2 ) mostly depended on the rural development index, while HANPP luc and HANPP harv were mainly controlled by urban and rural development, respectively. Our ﬁndings help understand, ﬁrst, how unbalanced regional development inﬂuences human-induced biomass occupation, the comprehensive urban ecological construction, and rural ecological restoration and, second, that the overall planning of urban–rural integration development must be strengthened to face greater ecological pressures in the future.


Introduction
The increase in human activity (urbanization, deforestation, cultivation, hunting, etc.) has had a tremendous impact on ecosystems and ecological resources. Excessive interference has led to a series of contradictions in human-land relationships and has become the core of geography research [1]. Intensified human-land conflict and the excessive transformation and occupation of ecological resources have posed severe challenges to sustainable global development. In view of this situation, there is an urgent need to study ecological pressure and the factors that lead to it. However, quantifying the magnitude of the substantial amount of human interference and its consequences for the ecosystem is a challenge [2]. Human appropriation of net primary production (HANPP) is a credible aggregated indicator that quantifies the magnitude of human-activity-related alterations in the ecological energy flows from the perspective of the carbon cycle [3]. HANPP is always represented as the ratio of the ecological sources used by human activities to the potential ones, and can measure regional ecological pressure and sustainability. Net primary productivity (NPP) refers to the amount of photosynthetic active radiation fixed by green vegetation as the carbon storage per unit area [4], which determines the ecological resource capacity that humans can occupy [5].
Urban and rural areas are considered as two major bases of human activities, and there is a big difference between their construction and economic levels [6][7][8]. As the most typical way by which the ecosystem is altered by humans, urbanization can not only directly consume multiple ecological resources through land conversions, but also drive the exchange of the regional population, economy, and materials [9,10]. Meanwhile, as suppliers of materials, labor, and various services to cities [11], expansive rural areas bear the pressure of ecological resource consumption resulting from urban population growth. With the continuous concentration of population in cities, there is a local and telecoupling phenomenon between the consumption and production of ecological resources in urban and rural areas [12]. The unbalanced development between urban and rural areas directly affects the supply-demand relationship of ecological resources and regional functions, especially the food security and ecological security in the areas where ecological resources are supplied are under increasing pressure. Whether human activities in some areas have exceeded the threshold of the ecosystems and what are the influencing factors are the topical questions that need to be answered urgently. Therefore, exploring the differentiation of HANPP under the influence of urban and rural development is helpful to identify the areas of human activity overload, and to clarify the contributions of urban and rural areas to the sustainable development of an ecosystem. In recent years, with the implementation of strategies for rural revitalization and urban-rural integration, China, as a large developing country, has entered a critical period of transformation in urban-rural relations. However, the long-standing "dual system" has resulted in the division of urban and rural areas and unbalanced development [13] accompanied by heterogeneous ecological effects of human activities. Therefore, determining the internal driving mechanisms of the unbalanced urbanrural development on HANPP is helpful to judge the constraints affecting the ecological sustainable development, which is a crucial topic to alleviate the ecological contradictions between urban and rural areas and to improve human-land relations globally.
HANPP is an effective method of assessing human-land relationships and research on HANPP has received increasing attention internationally; however, related research in China is still in its infancy. Relevant studies mostly focus on the improvement of the simulation model [14], comparison with other similar indicators [15], and the evaluation and spatial-temporal dynamic analysis of HANPP on global [16][17][18], national [5,[19][20][21], or regional scales [2,[22][23][24]. Embodied HANPP, an indicator closely related to the import and export of ecological resources, is an emerging research field [25]. Different results can be found based on various global-scale models and periods, i.e., HANPP was estimated to be 20.32% in 1995 [26], 23.8% in 2000 [16], and 13-25% from 1910 to 2005 [18]. Existing national-scale studies mostly focus on European countries and always show a relatively high level of HANPP (over 50%), with a tendency to decline over a long time series [27][28][29]. However, case studies in Asian countries show that the ecological pressure has significantly increased [5,30], representing the imbalanced effect of human activities on HANPP in densely inhabited regions. Moreover, as a biophysical indicator for measuring regional sustainable development, the relationship between HANPP and ecological footprint [31,32] as well as biodiversity [33,34], has always been focused on. Studies on the driving mechanism of HANPP are mostly limited to the single factor of social economy, among which the decoupling of HANPP and population or GDP has attracted more attention [35]. Population density has been found to have a significant influence on the spatial distribution of HANPP [36]. However, less attention has been paid to the comprehensive studies of how HANPP responds to multidimensional indicators for urban and rural development, such as land conversion, population growth, economy development, technology progress, and ecological construction, and studies on the variation in HANPP during urban-rural relationship transformation are also particularly lacking in recent years.
The past decades have been characterized by rapid economic development and intensive human activities in China, and the duality of rural-urban development is a common problem in most regions, inducing an obvious spatial heterogeneity of ecological resources consumed by human activities. If the relationship between urban-rural development and HANPP was quantified, it would fill research gaps and help answer the essential question of how to reduce the ecological pressure of human activities and promote harmonious human-land relationships through coordinated urban-rural development, and thus be referable for further studies in other regions. Therefore, with the ultimate goal of identifying the influence mechanism of urban-rural development on HANPP and investigating the factors driving its heterogeneity, three detailed objectives were defined in this study: (1) to evaluate HANPP in China in 2015 and its spatial imbalance characteristics; (2) to depict the relationships between HANPP and China's urban-rural development and investigate HANPP responses to the degree of regional urban-rural coordination; and (3) to identify the dominant driving factors of different dimensions of urban and rural construction on HANPP and its components.

Data Sources
In this paper, three categories of datasets were chosen from 31 regions in China (22 provinces, 5 autonomous regions, and 4 municipalities). This excluded unobtainable statistical datasets from Hong Kong, Macau, and Taiwan. The dataset included remote sensing images, statistical yearbook data, and vector datasets. To estimate the terrestrial potential NPP, yearly land cover data MCD12Q1 with the resolution of 500 m were obtained from NASA Earthdata (https://lpdaac.usgs.gov/products/mcd12q1v061/ (accessed on 12 May 2022)) and MOD17A3HGF data were used to characterize the actual NPP (https: //lpdaac.usgs.gov/products/mod17a3hgfv061/ (accessed on 23 May 2022)). The Monthly Normalized Difference Vegetation Index (NDVI) with a resolution of 1 km was obtained from the MOD13A3 dataset (https://lpdaac.usgs.gov/products/mod13a3v061/ (accessed on 6 June 2022)) and was mainly used to spatialize the harvested NPP of each province by virtue of its linear relationship with regional crop and livestock production. Among the statistical data, the annual average temperature and the annual total precipitation data of 830 meteorological stations in China used in this study were collected from the China Meteorological Data Service Centre (http://data.cma.cn/ (accessed on 15 June 2022)) and used to drive the Thornthwaite Memorial model to evaluate the potential NPP in China. For the missing data of several individual stations, this study selected the threshold of 80% and used Lagrange's theorem to interpolate and fill in the gaps [37]. If the number of missing months exceeded 20%, the data of the station in that year would be taken as a null value. The dataset of 31 provinces in China on biomass harvested or destroyed during harvest for crop yield, wood harvest, and livestock, available in the China statistical yearbooks in the National Bureau of Statistics (http://www.stats.gov.cn/tjsj/ (accessed on 17 June 2022)), was quantified. Indicators representing urban and rural population growth, economic development, land expansion, ecological construction, and agricultural technology, which comprehensively reflect the level of urban-rural development in China, were also obtained from these annual statistical yearbooks.

Estimation Model of HANPP
On the basis of different research backgrounds and spatial scales, several definitions of HANPP have been used before [16,26,38,39], which caused a slight obstacle when the HANPP estimation results were being compared. In this study, we chose the definition proposed by Haberl et al. (2007) [16], which has been widely recognized in recent years. According to the definition of HANPP [16], two major parts were considered when human activities appropriated the ecosystem NPP, the harvested biomass and the HANPP caused by land use and land conversion. Specifically, agricultural harvest, forestry product, and livestock grazing are defined as principal combinations of the harvested part (HANPP harv ) [40]. The NPP of the potential natural vegetation without human disturbance is called NPP pot and the one of the currently prevailing vegetation is called NPP act , and the difference between them (HANPP luc ) was defined as the change in biomass caused by human-induced land use. Some amount of biomass still remains in the ecosystem after harvest (i.e., straw and stump), which is called NPP eco . Therefore, HANPP could be calculated using Formula (1) and transformed into more detailed forms using Formula (2). Figure S1 presents a conceptual model of HANPP.
Accurate estimation of NPP pot is the key to quantifying HANPP luc . In terms of the selection of estimation models, many studies have used the Lund-Potsdam-Jena dynamic global vegetation model (LPJ-DGVM) [16,41]. The Integrated biosphere simulator (IBIS) and LPJ-GUESS (general ecosystem simulator) models can also simulate the vegetation dynamics under different climate scenarios [42]. However, these models are mainly applied on the global scale, and the resolution is generally 0.5 • × 0.5 • , which is not suitable for the accurate characterization of small-scale HANPP when used alone. In addition, some studies have established the transformation relationship between NPP act and NPP pot in stable-vegetation areas based on climatic elements to achieve an extended estimation of non-vegetated areas. Among them, O'Neill et al. (2007) divided all pixels of land use data into ones disturbed and undisturbed by human activities and assumed that the NPP pot of undisturbed pixels is equal to NPP act , whereas the NPP pot of all disturbed pixels could be replaced by the average NPP act of the surrounding undisturbed pixels [2]. Models such as Miami, Thornthwaite Memorial, and Chikugo have been widely used for estimating NPP pot by establishing the relationship between meteorological indicators and potential vegetation [5,23,43]. To integrate the advantages of each model, improve the accuracy of NPP pot estimation, and represent the optimal NPP pot in their research, Erb et al. (2016) used the average value of the Miami model, the LPJ-DGVM model, and a vegetation-approach model [44]. Therefore, in this study, NPP pot was derived by combining the Thornthwaite Memorial model, the LPJ model, and the vegetation-approach model, among which the result of the LPJ model came from the potential NPP simulation with a resolution of 0.5 • × 0.5 • under the S2 scenario provided by the TRENDY program. In addition, according to the research of Erb et al. (2016) on global vegetation biomass [44], the conversion rules of NPP act and NPP pot of various land use types were extracted (Table S3) so as to drive the vegetation-approach model. The spatial patterns of NPP pot simulated by the three models are highly similar, and the results of each model and the average value of NPP pot in China are shown in Figure S2.
HANPP harv is the biomass harvested or destroyed during harvest (harvest of crops and crop residues, harvest through livestock grazing, harvest through forestry, etc.), which is derived from the statistics yearbooks from the Chinese National Bureau of Statistics database, and HANPP harv in this study was converted into carbon (gC) by assuming an average carbon content of dry matter (DM) biomass of 50% [18]. HANPP harv on cropland consists of crop harvest and the used or unused residues [45]. The annual yields of 15 kinds of crops widely cultivated in China were collected, and their harvested fresh weights were converted into dry matter equivalents with the moisture content (%). The moisture content of the specific crop was obtained from the International Network of Food Data Systems provided by FAO (http://www.fao.org/infoods/infoods/tables-and-databases/ (accessed on 20 March 2022)) and the literature [3,5]. Furthermore, the crop residues were calculated from crop yield data with the harvest factors and recovery rates (the former is the ratio of crop residual to primary crop harvest, and the latter is the ratio of used crop residues to available residues), which were obtained from the literature [18]. The grazed biomass was calculated as the difference between feed demand and supply [45]. In this study, 5 grass-feeding livestock species (asses, horses, mules, cows, and sheep) were transformed into standard sheep units and the feed demand was calculated by the forage consumption (1.8 kg of dry matter coarse forage per sheep unit per day). In addition, the feed supply was derived by the conversion coefficients from croplands to straw from 6 crops species (i.e., rice, wheat, maize, millet, sorghum, and soybeans) and the 25% straw utilization factor [5]. The forestry harvest includes the removal of timber and felling losses. The total timber and bamboo production was collected from the Chinese National Bureau of Statistics and the National Forestry Statistical Yearbook; 49.5% of wood density (0.495 t/m 3 ) was used to convert fresh weight to dry matter [46], and the felling losses were estimated using recovery rates defined as the ratio of wood removal to felling [18]. In addition, considering the method of Haberl et al. (2007) [16], HANPP harv on built-up land was assumed to be 50% of NPP act . The detailed calculation methods of HANPP harv are provided in Supplementary File S1.

Quantification of Urban and Rural Development
In this study, the focus was how HANPP responds to regional development in China, so multi-indicators were selected to quantify the unbalanced urban-rural growth. Urban areas expand spatially when people crowd into these areas as the economy improves. According to this, 10 urbanization indicators from three perspectives (economy, population, and spatial expansion) were chosen to detect their relationship with HANPP. The corresponding rural development could also be measured using those 10 indicators but from the perspectives of economy, population, and agricultural technology. The indicators from 31 Chinese provinces, presented in Table 1, which were collected from China's statistical yearbooks of 2015, provide a comprehensive reflection of China's urban-rural development. The level of coordinated urban-rural development in China was also assessed using these indicators. After standardizing with the data range, the weight of each indicator in Table 1 was calculated using the entropy method [47]. Then, the urban and rural development indexes (expressed as f (u) and f (r), respectively) in 31 provinces in China could be quantified. A larger value of f (u) or f (r) means a higher level of rural-urban development, with which the comprehensive urban-rural development index (T) was also calculated. In addition, referring to the coupled coordination model in physics, the coupling degree (C ∈ [0, 1]) expressed by f (u), f (r), and T was used to analyze whether the whole urban-rural system is ordered. Since the C value only reflects the interaction within each system, which cannot indicate the coordinated order [48], the coordinated development degree (D ∈ [0, 1]) in each province was evaluated by the coordination degree model, which was divided into 9 levels in this study [49] ( Table 2). The detailed formulas of all the indicators are presented in Supplementary File S2.

Statistical Methods
Spatial analysis in ArcMap 10.3 and multi-statistical methods were applied to the evaluation of spatial heterogeneity and driving factors of HANPP in China. In this study, the independent-samples t test was used to explore the differentiation of HANPP and its components and investigate whether HANPP differs significantly among different zones in China (i.e., western, eastern, northern, and southern area; coastal land; and inland area). Furthermore, stepwise regression analysis helped verify what kind of constraints the urban-rural development in China exerted on HANPP and then we extracted the factors dominating the differentiation of HANPP based on the standard partial regression coefficient in the model.

Spatial Distribution of HANPP in China
The total amount of human appropriation of biomass was 2.68 PgC, representing 60.33% of the NPP pot in China in 2015. Figure 1a shows that HANPP (gC/m 2 ) gradually decreased from the southeast to the northwest of China, showing an obvious difference on both sides of the Hu Huanyong Line (the mutation line of China's population density). Hunan, Hubei, Henan, Jiangxi, Guangxi, Guangdong, and Jilin Provinces had high HANPP values (greater than 1000 gC/m 2 ). The negative value of HANPP (gC/m 2 ) came from the contribution of HANPP luc , mainly concentrated in central Xinjiang, western Yunnan, and southern Tibet, reflecting that the actual carbon sequestration of vegetation was higher than the potential sequestration under the same climatic conditions. Among these, the negative HANPP values in Xinjiang Province mainly appeared in the cultivated land. In addition, two high-value strips were identified in the east-west direction of HANPP (gC) in each province (Figure 1b Figure 1c shows obvious positive and negative differences in HANPP luc (gC/m 2 ) in China in 2015, indicating that land use change played a bidirectional role on ecosystem NPP. An HANPP luc (gC/m 2 ) value of less than 0 means that the production capacity of vegetation after human intervention could be higher than the potential capacity under the same climatic conditions. The negative values of HANPP luc (gC/m 2 ) in Xinjiang, Ningxia, and Gansu Provinces were mainly the result of cultivated land, indicating that in regions with a relatively insufficient natural background, cultivation based on artificial management could break through the climate limit of biomass accumulation to a certain extent, which has been confirmed by an existing study [18]. On the contrary, positive values of HANPP luc (gC/m 2 ) were dominant across the whole country, indicating that human activity has destroyed the original green spaces to a certain extent, and HANPP luc (gC/m 2 ) values in big cities (i.e., Beijing, Tianjin, Shanghai, Shenzhen, and Guangzhou Provinces) were significantly higher than those in surrounding areas. Total HANPP luc (gC) values of over 0.1 PgC are mainly concentrated in Inner Mongolia (0.148 PgC), Tibet (0.130 PgC), Heilongjiang Province (0.118 PgC), and Sichuan Province (0.112 PgC) (Figure 1d). Although China has continuously strengthened the protection of grasslands and vigorously promoted the policy of "returning grazing land to grassland" in recent years, most grassland was still overloaded and overgrazed [50]. Inner Mongolia and Tibet are the key pastoral areas of China, and the high values of HANPP luc (gC) in Inner Mongolia and Tibet were mainly related to grassland degradation caused by grazing.
The North China Plain and northeast provinces in China had relatively high HANPP harv (gC/m 2 ) values ( Figure 1e), with an average value of more than 600 gC/m 2 , and the biomass harvest in these regions was mainly provided by cultivated land. Among them, Shandong and Henan Provinces had the most extensive distribution of cultivated land with stable production. In addition, there was significant spatial heterogeneity of HANPP harv (gC/m 2 ) among the provinces in southern China, where the harvest of crops and the amount of crop residues were much higher than those in the surrounding woodlands. The limited ecological land in northern Xinjiang Province showed relatively high HANPP harv (gC/m 2 ), mainly resulting from grassland and cultivated land. Figure 1f shows that the high-value areas of HANPP harv (gC) in each province in 2015 expanded from the North China Plain to the entire northeast of China. Among them, Henan (0.097 PgC), Heilongjiang (0.089 PgC), and Shandong (0.077 PgC) were the top three provinces of HANPP harv (gC), while Beijing, Tianjin, and Shanghai Province were at the bottom since they did not undertake the production function and their area is relatively small. HANPP (%NPP pot ) can reflect the ecological pressure caused by the use of resources resulting from human activities, and the low-value areas of HANPP (%NPP pot ) were mainly distributed in the southeast coastal areas and southwest China (Figure 1g). Specifically, the mountainous area accounted for more than 80% of the total area in Yunnan Province, which limited the destruction of natural vegetation and land reclamation by humans. In addition, the HANPP harv in southeastern coastal provinces mainly resulted from the lowlevel harvest of forestry, which also led to a relatively low HANPP (%NPP pot ). It is worth noting that Tibet is a major pasturing area in China, but the NPP act in Tibet was only 137.38 gC/m 2 , while the harvest through livestock grazing was already close to the NPP pot , so the HANPP (%NPP pot ) in most areas of Tibet was generally high in 2015. Estimating the HANPP (%NPP pot ) of each province helps identify the regions with oversaturated ecological occupation in China, which provides scientific support for formulating urban development strategies and ecological construction measures in corresponding areas. It can be seen from Figure 1h that HANPP was high (%NPP pot ) in the grain-producing provinces in the North China Plain, while HANPP was low in the southwestern provinces with complex terrain. With the major grain-producing provinces benefiting from improved agricultural technology and efficient manual management, the high value of HANPP (%NPP pot ) in these provinces was actually a sign of a significant increase in land production efficiency.

Contributions of HANPP Components
The total amounts of HANPP luc and HANPP harv were found to be 1.63 PgC and 1.05 PgC, respectively, in China in 2015. Specifically, the harvest through cropland, livestock grazing, and forestry contributed 29.86%, 8.53%, and 0.91% to the total HANPP, respectively, while the land-use-induced biomass appropriation contributed 60.70%. However, there were great differences in the contributions of HANPP components among provinces ( Figure 2). HANPP luc played the dominant role among 21 provinces (over 50% of HANPP), especially in the cities with rapid urban expansion (i.e., Shanghai and Beijing, with 74.82% and 74.78% of HANPP luc (%HANPP), respectively) and some coastal cities (i.e., Fujian, Zhejiang, and Guangdong, with 86.07%, 85.99%, and 79.07% of HANPP luc (%HANPP), respectively). On the contrary, the harvest on cropland (HANPP harv (crop)) contributed most to HANPP in Shandong (62.33%), Henan (57.30%), Jilin (51.66%), Jiangsu (50.58%), and Hebei (50.46%), the traditional grain-producing provinces, and livestock grazing in the traditional pastoral regions (i.e., Gansu, Xinjiang, Ningxia, and Qinghai) resulted in 20% of HANPP. It is worth noting that different urban functions and development planning resulted in various ecological occupation structures, and expansive impervious surface replaced almost all ecological land in Shanghai, the economic center in China, where HANPP luc contributed the highest to NPP pot in China (50.55%). Compared with the high HANPP (%NPP pot ) due to harvested biomass through cropland, livestock grazing, and forestry, the substantial compression of ecological land due to excessive urban development needs more attention.  The amount of HANPPharv can reflect the production capacity of an ecosystem, and the regional population directly determines the consumption demand of the local biomass. Without considering the local import and export, we built a production-consumption ratio, which was calculated as the ratio of HANPPharv to the multiplication of national biomass consumption per capita and the regional population. The production-consumption ratio was divided into three types: an area was self-sufficient if the ratio was close to one, a ratio over one indicated adequate biomass production, and a ratio under one indicated inadequate biomass production. The latter two values lead to ecological resources being circulated among regions. On the basis of HANPPharv in 2015, a significant imbalance could be found between the supply of and demand for biomass in China (Figure 3a). In all, 14 provinces in the north and southwest China were defined as the production type, and ratios over the 95th percentile were found in the sparsely populated Tibet (6.36) and Inner Mongolia (4.07). Conversely, the remaining 17 provinces in the middle and southeastern coastal areas in China needed external biomass supplies. Specifically, the production-consumption ratios of rapidly growing cities such as Shanghai (0.06) and Beijing The amount of HANPP harv can reflect the production capacity of an ecosystem, and the regional population directly determines the consumption demand of the local biomass. Without considering the local import and export, we built a production-consumption ratio, which was calculated as the ratio of HANPP harv to the multiplication of national biomass consumption per capita and the regional population. The production-consumption ratio was divided into three types: an area was self-sufficient if the ratio was close to one, a ratio over one indicated adequate biomass production, and a ratio under one indicated inadequate biomass production. The latter two values lead to ecological resources being circulated among regions. On the basis of HANPP harv in 2015, a significant imbalance could be found between the supply of and demand for biomass in China (Figure 3a). In all, 14 provinces in the north and southwest China were defined as the production type, and ratios over the 95th percentile were found in the sparsely populated Tibet (6.36) and Inner Mongolia (4.07). Conversely, the remaining 17 provinces in the middle and southeastern coastal areas in China needed external biomass supplies. Specifically, the productionconsumption ratios of rapidly growing cities such as Shanghai (0.06) and Beijing (0.09) were under the fifth percentiles. Furthermore, it should be clarified that a telecoupling of ecological resources always occurs in remote spaces, and the gradient differences of production-consumption ratios in Figure 3a represent interregional resource flows; more attention should be paid to this in future studies. (0.09) were under the fifth percentiles. Furthermore, it should be clarified that a telecoupling of ecological resources always occurs in remote spaces, and the gradient differences of production-consumption ratios in Figure 3a represent interregional resource flows; more attention should be paid to this in future studies. Five dominant production types were identified among 31 provinces with production-consumption ratios of HANPPharv through agriculture, forestry, and grazing ( Figure  3b), and a significant spatial aggregation could be found among all these types. A great number of ecological resources were needed among consumption regions with ratios under 1, which included the metropolis with high population densities (i.e., Shanghai, Beijing, and Tianjin) and several provinces in the middle of China. On the contrary, only two production regions (i.e., Inner Mongolia and Jilin) were found in the northeast of China where the production-consumption ratios of the components of HANPPharv were over 1. Furthermore, huge gaps existed among the three ratios in the remaining provinces and the highest ratio was selected as the dominant production type. Specifically, five agriculture-dominant provinces, situated mostly in the circum-Bohai Sea region, and five provinces in south China with better climate conditions, provided more biomass through forestry. HANPPharv through grazing was found as the dominant production type among eight traditional husbandry provinces in the west of China. Figure 3c shows that a universal positive correlation could be found between the production-consumption ratio and the HANPPharv (% of HANPP) through agriculture, forestry, and grazing in 31 provinces in China. The aggregation characteristics of sample points declared whether the HANPPharv was consumed locally, showing that the interregional flows and the exchange of agricultural and livestock products were more complicated than the forestry products. Five dominant production types were identified among 31 provinces with productionconsumption ratios of HANPP harv through agriculture, forestry, and grazing (Figure 3b), and a significant spatial aggregation could be found among all these types. A great number of ecological resources were needed among consumption regions with ratios under 1, which included the metropolis with high population densities (i.e., Shanghai, Beijing, and Tianjin) and several provinces in the middle of China. On the contrary, only two production regions (i.e., Inner Mongolia and Jilin) were found in the northeast of China where the productionconsumption ratios of the components of HANPP harv were over 1. Furthermore, huge gaps existed among the three ratios in the remaining provinces and the highest ratio was selected as the dominant production type. Specifically, five agriculture-dominant provinces, situated mostly in the circum-Bohai Sea region, and five provinces in south China with better climate conditions, provided more biomass through forestry. HANPP harv through grazing was found as the dominant production type among eight traditional husbandry provinces in the west of China. Figure 3c shows that a universal positive correlation could be found between the production-consumption ratio and the HANPP harv (% of HANPP) through agriculture, forestry, and grazing in 31 provinces in China. The aggregation characteristics of sample points declared whether the HANPP harv was consumed locally, showing that the interregional flows and the exchange of agricultural and livestock products were more complicated than the forestry products.

Regional Differentiations in HANPP
An obvious spatial differentiation of HANPP was found in China. However, the focus should be to identify whether there is a significant difference among regionalization zones, and nine HANPP indexes were chosen for further discussion (Table 3). Here, the independent-samples T test was used to detect the differentiation of nine HANPP indexes across the four zones (Zone 1 to Zone 4; Figure S3). The result showed that HANPP (gC/m 2 ), HANPP harv (gC/m 2 ), and HANPP luc (gC/m 2 ) were significantly different on both sides of the Huanyong Hu line (Zone 1) but their total amount and contribution to NPP pot did not pass the test. The same regional differentiation of HANPP indexes could be found between the western and eastern parts of China (Zone 2) as Zone 1, and the ecological pressure caused by human occupation (HANPP (%NPP pot )) was also significantly different in Zone 2. HANPP harv and HANPP luc have different sensitivities to meteorological elements, so the corresponding indexes showed various characteristics between the north and south parts of China (Zone 3), which was divided on the basis of regional temperature and precipitation. There is a difference between economic development and population agglomeration in China's coastal and inland regions, so the HANPP harv indexes related to land productivity did not show significant differentiations in Zone 4. Furthermore, from the perspective of the nine HANPP indexes, it is worth noting that HANPP (gC/m 2 ) showed significant differences among all the regionalization zones ( Table 3). The complex characteristics of HANPP indexes in China indicated that they were not only restricted by population density and terrestrial natural conditions but that a complex mechanism of regional comprehensive development (i.e., economy, urbanization, and agrotechnique) also had an impact on them.  Figure S3. *: p < 0.05, **: p < 0.01, and ***: p < 0.001.

Credibility of HANPP
In existing studies, HANPP (%NPP pot ) has attracted more attention than HANPP (gC/m 2 ) or HANPP (gC) due to its ability to represent regional ecological sustainability and the ecological pressure caused by human occupation. However, the lack of a uniform standard in the selection of estimation models of NPP pot and NPP act hinders the comparison of HANPP results among different regions. For instance, in related studies, a variety of models have been used to estimate NPP pot , e.g., LPJ-DGVM [16,41], the Guangsheng Zhou model [5], the new calibrated Miami model [45], and multiple-iterations modification based on the disturbed pixels in NPP act [2]. Rather than comparing the similarities of HANPP (%NPP pot ) between regions, it is better to focus on whether its spatial distribution and long-term variation characteristics are similar.
In China, HANPP (%NPP pot ) was found to be 60.33% in 2015, which is higher than it was in 2010 (57.80%), estimated by Chen et al. (2015) [5], as well as the global HANPP (%NPP pot ) values in 2000 (23.80%), 2005 (25.00%), and 2010 (20.70%) [16][17][18]. Huang et al. (2020) found that HANPP (%NPP pot ) was 72.4% in the Yangtze River Delta in 2015 [51], and we found a slightly lower value (68.92%) in the same region in this study. From the perspective of the spatial distribution characteristics of HANPP, the grid-scaled HANPP (%NPP pot ) in this study is highly similar to the results of Haberl et al. (2007) and Kastner et al. (2021) [16,17]. In general, the simulation of HANPP among provinces and cities in China is still in its infancy and the deviations are mainly due to different methods and datasets used during NPP pot and HANPP harv estimation within limited studies.

Relationships between HANPP and China's Urban and Rural Development
Significant spatial differentiation of HANPP was found under the impact of unbalanced human activities and the ecological environment in China, especially the duality of rural-urban development. In this study, 20 indictors reflecting regional economy, population, urban expansion, and rural agricultural technology were chosen, as shown in Table 1, and the correlation coefficients and nine HANPP indexes (Table 3) of 31 provinces were summarized (Figure 4). From the perspective of urban economy, significant positive correlations were identified between the added value of the secondary and tertiary industries (UE1) and HANPP (gC/m 2 ), HANPP harv (gC/m 2 ), and HANPP luc (gC/m 2 ). The added value of transportation, warehousing, and postal services (UE3) had a significant influence on most HANPP indexes (Figure 4a), which guarantees the circulation of regional ecological resources. However, UE2 and UE4 showed a negative correlation with the total amount of HANPP, HANPP harv , and HANPP luc , which is because a larger proportion of UE1 and a higher consumption expenditure of urban residents always appears in highly urbanized areas, which usually have a relatively weak biomass production (low HANPP harv ); the HANPP luc would also not increase significantly when urban construction is saturated, so the HANPP indexes would be suppressed instead. High urban population (UP1) resulted in a great amount of biomass occupation per unit area, whereas the negative correlations between urban population proportion (UP2) and HANPP (gC), HANPP harv (gC), and HANPP luc (gC) reflected the polarization characteristics of the urban population agglomeration and ecological resource production. Among the urban spatial expansion indicators, the higher regional average nighttime light intensity (US1), the proportion of the built-up area (US2), and the proportion of the urban road area (US4) always corresponded with large HANPP luc (% NPP pot ) and HANPP luc (gC/m 2 ) values, which also directly led to the significant positive correlations (p < 0.01) and relatively high correlation coefficients with HANPP (% NPP pot ) and HANPP (gC/m 2 ). Conversely, urban ecological construction gained valuable space to counter the degradation of natural vegetation caused by urban expansion, and most of the HANPP indexes, especially the HANPP luc indexes, decreased as the proportion of urban green spaces increased (US3).  HANPP indexes responded to rural development in different ways (Figure 4b). Stronger positive correlations existed between the added value of the primary industry (RE1) and all HANPP harv indexes than UE1, and the total amount of HANPP, HANPP harv , and HANPP luc in each province was more sensitive to the proportion of the added value of the primary industry (RE2) than other HANPP indexes. Higher consumption expenditure of rural residents (RE3) appeared not only in highly urbanized areas but also in the agriculture, forestry, and animal husbandry developed areas, which showed relatively lower correlations with the total amount of HANPP, HANPP harv , and HANPP luc than UE4. The higher the Engel coefficient of rural residents (RE4), the lower the HANPP (% NPP pot ) and all HANPP harv indexes. This shows that the main article consumed by residents was food, the regional economy was relatively lagging in terms of development, and the influence of human activities on ecology was relatively weak. Rural population (RP1) had significant and higher positive correlations with all the HANPP harv indexes than UP1 but no statistical correlations with HANPP luc indexes. Moreover, the biomass occupied through land use per unit area and its percentage in NPP pot would be reduced if the rural population was the main component (RP2) in a region. All the rural agricultural technology indicators showed similar impacts on the nine HANPP indexes, especially in terms of the significant positive effects (p < 0.01) on and high correlation coefficients among HANPP harv (gC), HANPP harv (% NPP pot ), and HANPP harv (gC/m 2 ). Adequate effective irrigated area (RA3) and total sown area of crops (RA4) are necessary for adequate food supply, the total power of the agricultural machinery (RA1) reflects the regional agricultural machinery ownership, and the harvest biomass of grain and forestry would benefit from these elements as well as the amount of agricultural fertilizer (RA2) applied.
An interaction effect always exists between regional economy improvement and biomass occupation, and industrial structure would alter the HANPP proration. Overdevelopment of the secondary and tertiary industries would restrain harvest biomass, but adequate ecological supply depends on the primary industry, so balanced industrial development should be the focus among regions. Harvested biomass (HANPP harv ) is mostly determined by the rural population rather than the urban population, while the urban or rural population proportion has hardly any influence on it, indicating that urban people mainly play the role of the consumer of ecological resources. Urban sprawl, a result of land surface modification by humans, usually leads to the irrecoverable loss of NPP. Cultivated land offers additional products along with the improvement of agricultural technology, although it replaces part of the natural vegetation. Moreover, worldwide, the increasing populations and economies are inevitably accompanied by rapid urbanization, and the corresponding requirements of ecological resources must be sustainable. Thus, through policies and agriculture techniques, it is necessary to improve agriculture production and the utilization efficiency of NPP per area, as well as properly control urban expansion.

HANPP Responses to Coordinated Regional Urban-Rural Development
As two major systems indispensable to regional development, urban areas and rural areas play an important role in promoting sustainable use of ecological resources. Therefore, coordinated urban-rural development and the impact on the HANPP in 31 provinces in China was quantified. The quantification of the urban development index (f (u)) indicated that highly urbanized areas (Figure 5a) were mostly concentrated in the coastal provinces but high values of the rural development index (f (r)) mostly appeared in southern China (Figure 5b). Specifically, Shanghai (f (u) = 0.90), Beijing (f (u) = 0.43), Henan (f (r) = 0.73), and Shandong (f (r) = 0.68) showed the peak values over the 95th percentiles of f (u) and f (r). In contrast, the valley values under the fifth percentiles were in Qinghai (f (u) = 0.04), Tibet (f (u) = 0.02), Tianjin (f (r) = 0.09), and Beijing (f (r) = 0.08). Furthermore, the degree of coordinated urban-rural development (D) among provinces was evaluated by the comprehensive urbanrural development index (T) and the coupling degree (C) (Table S4) and showed an obvious spatial differentiation between coastal and inland areas (Figure 5c). In this study, the D value was divided into nine levels ( Table 2) to map coordinated urban-rural development types in China (Figure 5d). Only a few coastal provinces had achieved primary coordination (Jiangsu and Guangdong) and barely coordination (Shandong and Zhejiang). Inversely, most regions showed a universal imbalance and urban-rural development gap. Notably, there was severe imbalance in Heilongjiang Province, Tibet, Gansu Province, and Xinjiang Province. In short, the ideal urban-rural coordination did not exist in China in 2015.
5c). In this study, the D value was divided into nine levels ( Table 2) to map coordinated urban-rural development types in China (Figure 5d). Only a few coastal provinces had achieved primary coordination (Jiangsu and Guangdong) and barely coordination (Shandong and Zhejiang). Inversely, most regions showed a universal imbalance and urbanrural development gap. Notably, there was severe imbalance in Heilongjiang Province, Tibet, Gansu Province, and Xinjiang Province. In short, the ideal urban-rural coordination did not exist in China in 2015.  (a,b) and coordinated urban-rural development degree (c) and types (d) in China. HANPP responses to coordinated regional urban-rural development (e). (PC, primary coordination; BC, barely coordination; NI, near imbalance; II, inchoate imbalance; MI, moderate imbalance; SI, severe imbalance; **: p < 0.05, and ***: p < 0.01).
The relationship between HANPP indexes and coordinated urban-rural development ( Figure 5e) helped create verifications and supplements for Figure 4. It is worth noting that f(u) and f(r) showed opposite driving effects on HANPP (gC), HANPPharv (gC), and HANPPluc (gC) and the regional ecological pressure characterized by HANPP (%NPPpot) and HANPP (gC/m 2 ) mostly depended on f(r). The highly urbanized regions with high f(u) usually appeared in provinces with limited territorial areas (i.e., Shanghai, Beijing, and Tianjin), in which the total amount of occupied biomass was at a relatively low level. On the contrary, provinces with high levels of f(r) had relatively large territorial areas in China and developed production capacity. As a component of HANPP, HANPPluc was mainly controlled by the urban development index (f(u)). Conversely, f(r) made significant positive impacts on HANPPharv. All the occupied biomass indexes per unit area showed the most sensitive response to coordinated urban-rural development degree (D) and good urban-rural coordination obviously improved HANPP (gC/m 2 ) and HANPPharv (gC/m 2 ) but relatively weaker urban-rural coordination influenced HANPPluc (gC/m 2 ) variation. From the perspective of six coordination types among provinces in China, the sample points of HANPPharv (gC/m 2 ) tended to be discrete, along with a higher degree of coordination (Figure 6), and the samples with moderate imbalance (MI) urban-rural coordination showed the highest degree of dispersion of HANPP (gC/m 2 ) and HANPPluc (gC/m 2 ). The relationship between HANPP indexes and coordinated urban-rural development ( Figure 5e) helped create verifications and supplements for Figure 4. It is worth noting that f (u) and f (r) showed opposite driving effects on HANPP (gC), HANPP harv (gC), and HANPP luc (gC) and the regional ecological pressure characterized by HANPP (%NPP pot ) and HANPP (gC/m 2 ) mostly depended on f (r). The highly urbanized regions with high f (u) usually appeared in provinces with limited territorial areas (i.e., Shanghai, Beijing, and Tianjin), in which the total amount of occupied biomass was at a relatively low level. On the contrary, provinces with high levels of f (r) had relatively large territorial areas in China and developed production capacity. As a component of HANPP, HANPP luc was mainly controlled by the urban development index (f (u)). Conversely, f (r) made significant positive impacts on HANPP harv . All the occupied biomass indexes per unit area showed the most sensitive response to coordinated urban-rural development degree (D) and good urban-rural coordination obviously improved HANPP (gC/m 2 ) and HANPP harv (gC/m 2 ) but relatively weaker urban-rural coordination influenced HANPP luc (gC/m 2 ) variation. From the perspective of six coordination types among provinces in China, the sample points of HANPP harv (gC/m 2 ) tended to be discrete, along with a higher degree of coordination (Figure 6), and the samples with moderate imbalance (MI) urban-rural coordination showed the highest degree of dispersion of HANPP (gC/m 2 ) and HANPP luc (gC/m 2 ).

Dominant Driving Factors of Each HANPP Index
Bidirectional effects were found among urban-rural development indictors of biomass appropriation, and the dominant driving factors of each HANPP index were detected by stepwise regression to eliminate the indirect correlation. Table 4 shows the fitting models and correlation coefficients from three perspectives (urban, rural, and urban-rural development), in which stronger decisive effects were displayed by the driving factor with larger standardized coefficients, and a higher R referred to a better fitting model. Among all the driving factors in the 27 fitting equations in Table 4, the rural development indicators (34 times) appeared more frequently than the urban development indicators (21 times). At the same time, in different dimensions of urban-rural development, each HANPP index was more likely to be affected by urban economy (UE) than rural economy (RE), more likely to be affected by rural population (RP) than urban population (UP), and more likely to be affected by rural agricultural technology (RA) than urban spatial expansion (US). From the analog effect of each equation, driven by 20 urban and rural development indicators, the analog effect of HANPPharv indexes was the best, followed by that of HANPP indexes, and HANPPluc indexes were the weakest; the selected urban and rural development indicators could explain their variation characteristics only to a certain extent. From the perspectives of urban, rural, and comprehensive urban-rural development,

Dominant Driving Factors of Each HANPP Index
Bidirectional effects were found among urban-rural development indictors of biomass appropriation, and the dominant driving factors of each HANPP index were detected by stepwise regression to eliminate the indirect correlation. Table 4 shows the fitting models and correlation coefficients from three perspectives (urban, rural, and urban-rural development), in which stronger decisive effects were displayed by the driving factor with larger standardized coefficients, and a higher R referred to a better fitting model. Among all the driving factors in the 27 fitting equations in Table 4, the rural development indicators (34 times) appeared more frequently than the urban development indicators (21 times). At the same time, in different dimensions of urban-rural development, each HANPP index was more likely to be affected by urban economy (UE) than rural economy (RE), more likely to be affected by rural population (RP) than urban population (UP), and more likely to be affected by rural agricultural technology (RA) than urban spatial expansion (US). From the analog effect of each equation, driven by 20 urban and rural development indicators, the analog effect of HANPP harv indexes was the best, followed by that of HANPP indexes, and HANPP luc indexes were the weakest; the selected urban and rural development indicators could explain their variation characteristics only to a certain extent. From the perspectives of urban, rural, and comprehensive urban-rural development, the analog effects of all the nine HANPP indexes were the best when only considering rural development indicators, followed by the regression models of comprehensive urbanrural development indicators, and the ability to explain HANPP indexes was relatively limited when only considering urban development indicators. From the perspective of the explanatory ability of each independent variable, the impact of urban economy (UE) on occupied biomass was much stronger than that of urban population (UP) and urban spatial expansion (US). On the contrary, rural population (RP) and rural agricultural technology (RA) had more significant explanatory effects on each HANPP index than rural economy (RE). Among all the independent variables, UP2, US1, US2, US3, and RA3 did not appear in any of the fitting models. Comprehensively considering the dominant driving factors and standardized coefficients of each HANPP index, it was found that UE3 had the strongest explanation for the differentiation characteristics of HANPP (gC/m 2 ) and that the restriction of RP and RA indicators on HANPP (%NPP pot ) was much stronger than that of urban development indicators. The total amount of HANPP (gC) was affected by the opposite effects of RA4 and US4, and RA4 showed a stronger control of it. In addition, RP and RA indicators showed absolute driving advantages in HANPP harv (%NPP pot ) and HANPP harv (gC) and RA indicators had a stronger positive driving effect on them. From the perspective of comprehensive urban-rural development, urban development indicators had no significant driving effect on the above two indexes. The influence of each independent variable on HANPP harv (gC/m 2 ) was relatively complicated, in consideration of RA, UE, and US indicators, and a better-fitting model could obtain HANPP harv (gC/m 2 ); the standardized coefficient showed that the urban development indicators played a more important role. Table 4 shows that the fitting equations between the independent variables and HANPP luc were relatively simple. US4, RE2, and UE1 were the dominant driving factors for HANPP luc (gC), HANPP luc (%NPP pot ), and HANPP luc (gC/m 2 ), respectively, and none of the population factors were correlated. It should be noted that the ability of selected independent variables to interpret all HANPP luc indexes was significantly weaker than their ability to interpret HANPP and HANPP harv , indicating that there are driving factors besides those selected in this study. Further analysis of the positive or negative driving differences of each independent variable showed that UE3, RA1, RA2, and RA4 had a significant positive influence in each fitting model, which could be verified and compared with the results in Figure 4. In summary, HANPP (%NPP pot ), HANPP harv (%NPP pot ), and HANPP luc (%NPP pot ) were more dependent on rural development, while urban development was the decisive factor for all HANPP indexes per unit area, and the total amount of occupied biomass of each province was more affected by rural development than urban.

Conclusions
In this study, we estimated the HANPP in China in 2015 based on multi-source datasets and quantified the contribution characteristics of HANPP components as well as the biomass production-consumption types in 31 provinces. The regional differentiations of HANPP indexes among typical regionalization zones were also identified. In addition, an evaluation framework of urban-rural development was built to investigate how HANPP indexes respond to the regional unbalanced development and we identified the dominant driving factors of HANPP indexes from the perspectives of economy, population, urban expansion, and agricultural technology. The results showed that the total amount of HANPP was 2.68 PgC and gradually decreased from the southeast to the northwest of China, which represented 60.33% of the NPP pot . The grain-producing provinces had high HANPP (%NPP pot ) values, while the southwestern provinces with complex terrain had low HANPP (%NPP pot ) values. HANPP (gC/m 2 ), HANPP harv (gC/m 2 ), and HANPP luc (gC/m 2 ) were significantly different (p < 0.05) between the Huanyong Hu Line and the western-eastern part of China, while HANPP harv did not show significant differentiations in the coastal land and inland areas.
From the perspective of HANPP components, in China, the total amount and the average value of HANPP luc and HANPP harv were found to be 1.63 PgC and 216.97 gC/m 2 and 1.05 PgC and 126.53 gC/m 2 , respectively, and the harvest through cropland, livestock grazing, and forestry contributed 29.86%, 8.53%, and 0.91% to the total HANPP, respectively. Land use change had a bidirectional influence on HANPP luc (gC/m 2 ); the negative values were mostly concentrated in a small number of cultivated land patches in northwest China, which indicated that the biomass production under intensive artificial management could be higher than the potential production. Moreover, the North China Plain and northeast provinces had relatively high HANPP harv (gC/m 2 ) values, mainly provided by cultivated land. Different urban functions resulted in various ecological occupation structures. HANPP luc played the dominant role among 21 provinces (over 50% of HANPP), especially in the cities with rapid urban expansion and some coastal ones. With the production-consumption ratio of biomass based on HANPP, 17 provinces in the middle and southeastern coastal areas in China were defined as the consumption type. Specifically, the ratios in Shanghai (0.06) and Beijing (0.09) were under the fifth percentiles, and a great number of ecological resources were needed from surrounding regions. In addition, a universal positive correlation (p < 0.05) was found between the production-consumption ratio and the HANPP harv (% of HANPP) through agriculture, forestry, and grazing in 31 provinces in China.
In different dimensions of urban-rural development, each HANPP index was more likely to be affected by urban economy (UE) than rural economy (RE); the impact of rural population (RP) was much stronger than that of urban population (UP); and rural agricultural technology (RA) had more significant explanatory effects on each HANPP index than urban spatial expansion (US). Specifically, in the dimension of economic development, the added value of transportation, warehousing, and postal services (UE3) had a significant influence on most HANPP indexes. From the demographic perspective, rural population (RP1) had significant and higher positive correlations with all the HANPP harv indexes than UP1 but no statistical correlation with HANPP luc . Meanwhile, agricultural technology indicators showed significant positive effects (p < 0.01) on and high correlations with HANPP harv . The higher regional average nighttime light intensity (US1), the proportion of the built-up area (US2), and the proportion of urban road area (US4) always corresponded with large HANPP luc values. Conversely, most HANPP indexes decreased as the proportion of urban green spaces increased (US3). Furthermore, the regional ecological pressure characterized by HANPP (%NPP pot ) and HANPP (gC/m 2 ) mostly depended on the rural development index (f (r)) and HANPP luc and HANPP harv were mainly controlled by f (u) and f (r), respectively. A higher degree of urban-rural coordination (D) obviously improved HANPP (gC/m 2 ) and HANPP harv (gC/m 2 ) but had a relatively weaker effect on the variation in HANPP luc (gC/m 2 ).
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/land12051062/s1: Figure S1: Conceptual model of human appropriation of net primary production (HANPP); Figure S2: Spatial distribution of NPP pot in China; Figure S3: Typical regionalization of China to detect the HANPP differentiation; Table S1: Moisture contents (%), harvest factors and recovery rates of typical crops; Table S2: Conversion coefficients from crop harvest to straw; Table S3: Regulation for assessing NPP pot by vegetation-approach model; Table S4: Urban-rural coordinated development indexes; Supplementary File S1: Estimation of HANPP harv ; Supplementary File S2: Evaluation of urban-rural coordinated development.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.