Spatial covariance of ecosystem services and poverty in China

ABSTRACT Ecosystem services (ESs) are known to be particularly important to the rural poor globally and effective management of such services is argued to be a sustainable pathway out of poverty. However, there is as yet no clear evidence as to how important ESs are for poverty alleviation, partly as there are very few large-scale studies addressing this issue. Here, we examine patterns of spatial covariation of income poverty and provisioning services and biodiversity using county-level data across China (n = 1924). We conduct our analyses both at the national scale and at the subnational scale. At the national scale, poor counties have significantly lower levels of agricultural provisioning services and water availability, but significantly higher levels of forest-related provisioning services and biodiversity. This finding supports the hypothesis that in general, high levels of poverty co-occur with areas with high levels of non-agricultural ESs. However, in the forest-dominated counties in southern China, low poverty, high densities of forest-related provisioning services and high levels of natural forest cover co-occur. Our results highlight the scale and context dependency of patterns of co-occurrence of poverty and ESs, and the importance of large-scale analyses for understanding the relationships between poverty and ESs. EDITED BY Berta Martín-López


Introduction
It is now widely accepted that sustainably managing ecosystem services (ESs) is a major global societal challenge (e.g. TEEB 2010) and that a key requirement for meeting this challenge is the mapping and quantification of the spatial distributions of both the supply of and demand for ES (Burkhard et al. 2012). This large and rapidly growing literature (reviewed in Crossman et al. 2013) is largely focused on spatial covariation of ES with each other, both singly and as part of ES bundles (e.g. Raudsepp-Hearne et al. 2010a). Recent work has frequently focused on explicitly considering the supply and demand of ES (Burkhard et al. 2012) separately, most often with the goal of identifying the most important areas for ES supply and demand (reviewed in Martinez-Harms et al. 2015). However, despite widespread recognition that relationships between both the supply and demand of ES are scale dependent (e.g. Anderson et al. 2009;Scholes et al. 2013), most existing ES mapping studies are carried out at a single spatial and temporal scale (Martinez-Harms et al. 2015). In addition, the majority of ES maps are based on land use and land-cover data (Martinez-Harms et al. 2015). Such land-cover-based maps are useful proxies of potential supplies of ES but are not well suited in of themselves for linking such supplies to human beneficiaries. In addition, land-cover-focused mapping means that disaggregation of beneficiaries by social group is difficultthis is potentially a major issue as different socio-economic groups derive different benefits from different groups of ES (Daw et al. 2011).
Of all segments of society, it is the rural poor who are believed to be most directly dependent on ESparticularly those from forestsfor their livelihoods (e.g. Sunderlin et al. 2008;Persha et al. 2011). It has also been argued that good management of ES is an important mechanism for sustainable development and poverty alleviation (TEEB 2010). While there is some evidence that increasing agricultural yields can be achieved without adversely impacting non-agricultural ES (Pretty et al. 2006), partly through smallscale sustainable agriculture and agroforestry (Perfecto and Vandermeer 2010), the majority of case studies suggest the oppositeincreases in human well-being usually result in degradation of non-agricultural ES (Raudsepp-Hearne et al. 2010b). Moreover, while heavily forested areas tend to be associated with high levels of rural poverty globally (Sunderlin et al. 2008), there is as yet no consensus about the importance (or lack thereof) of ES for alleviating poverty (reviewed in Suich et al. 2015). The literature that does exist is largely based on case studies examining the relationships between poverty and environmental resources, though there have been some recent large-scale analyses of this relationship (e.g. Angelsen et al. 2014;Ferraro et al. 2015).
The main objective of this study is to help address this paucity of large-scale understanding of the relationships between income poverty and ES by examining the spatial covariation of poverty and key provisioning services across all rural counties (n = 1924) in China. Our approach is to first quantitatively identify the spatial distributions of ES bundles, as this enables us to identify the major socioecological subsystems that exist across rural China (Raudsepp-Hearne et al. 2010a). We then explore the relationship between poverty and ES both at the national scale and within each ES bundle. Our rationale for this approach is that patterns of ES co-occurrence and trade-offs are known to be scale and context dependent (e.g. Anderson et al. 2009) and, as such, our a priori prediction was that relationships between poverty and ES would vary both between the scale (national vs. subnational) and socioecological context (bundle) considered.
We focus on China here, as despite extraordinary recent economic growth that reduced rural poverty from 75% to 13% between 1980(De Janvry et al. 2005) China continues to have a large poor rural population. The remarkable reductions in rural poverty between 1980 and 2001 have been attributed to not only increases in agricultural production but also large increases in non-farm incomes (De Janvry et al. 2005;Imai and You 2014). However, the remaining rural poor population may be less well placed to move out of poverty via these mechanisms, as they disproportionately live in forested mountainous areas (Rodriguez et al. 2009;FAO 2012). This suggests that ES from forests may have an important role to play in these regions of China. Moreover, China has the largest system of payment for ecosystem service (PES) schemes in the world, meaning that ES research in China has major policy relevance. There is a large and rapidly increasing literature on ES in China, including national-scale analyses (e.g. Ouyang et al. 1999Ouyang et al. , 2016Chen and Zhang 2000;He et al. 2005) as well as case studies examining the impacts of PES on both livelihoods and environmental outcomes (e.g. Cao et al. 2010;Li et al. 2011;Liu et al. 2008;Yang et al. 2013aYang et al. , 2013b; see also Suich et al. 2015). However, to our knowledge, ours is the first analysis to examine the spatial covariation of poverty and ES at the scale of China.

Methods
We use a variety of data sets and analytical techniques in this analysis, which are outlined in a flowchart ( Figure 1).

Creation of data sets
We use a variety of data sources in this manuscript, which we describe and justify in two sectionspoverty and ESs. Our base geographic information system (GIS) layer for all analysis is the rural subset of Chinese counties (n = 1924), which we derived from the county-level administrative map of China. All GIS analyses were conducted in ArcGIS 10 (ESRI, Redlands, California) based on the Albers Equal Area projection.

Overview and justification of indicator of poverty
Povertythe inability of an individual or family to meet basic needsis a multidimensional issue that has both monetary and non-monetary aspects (Sen 1992;Schleicher et al. 2017). The monetary aspect of poverty is relatively straightforward to measure as it can be assessed based on a single measure (minimum level of income). By contrast, the non-monetary aspectwhich relates to issues such as human capabilitiesis considerably more complicated to measure as it is highly multidimensional (Schleicher et al. 2017). Our indicator of poverty in this analysis is whether or not a county is designated as one of 592 'key poverty-stricken' counties (hereafter 'poor county') by the Chinese government (State Council Leading Group Office of Poverty Alleviation and Development 2012). This measure is exclusively of income (monetary) poverty. In general, a county with a mean per capita annual income of less than 2300 RMB in 2011 (approximately $365 US as of 31/ 12/2011) was classified as 'poor' by the Chinese government, with some regional adjustment of this criteria in eastern and western parts of China (State Council Leading Group Office of Poverty Alleviation and Development 2012). This classification has been revised four timesin 1986, 1994, 2006 and 2012. We used the 2012 version of the list of the 'key poverty-stricken' counties for our analyses. Our rationale for focusing only on income poverty here is twofold. First, we use the binary indicator of 'poor county' due to its simplicity and policy relevance. Income poverty is widely used to study poverty in China (Imai and You 2014), and the list of 'key poverty-stricken counties' we use here is a key indicator used by the Chinese government to target government assistance (State Council Leading Group Office of Poverty Alleviation and Development 2012). Second, income poverty is the aspect of poverty most directly linked to the ES for which data are available at the scale of China (provisioning ES), as multiple studies have shown that wealth in rural areas is directly linked to natural resources derived from the environment (Cavendish 2000;Porro et al. 2015).

Overview and justification of ES and biodiversity data
We considered nine provisioning services, and two measures of biodiversity, as well as change in forest cover in this study. As discussed, we focus on provisioning ES due to their direct relevance to rural livelihoods (Yang et al. 2013a;Hamann et al. 2015) and hence income poverty, as well as issues of data availability. We included two measures of biodiversity for its role both in regulating the ecosystem processes that underpin all ES as well as its value as an ES in its own right (Mace et al. 2011). We had no a priori expectation that biodiversity would relate directly to income poverty; however, understanding the degree to which poverty co-varies with biodiversity is of policy interest as a positive association would suggest that biodiversity conservation measures are likely to work best if combined with poverty alleviation measures. Finally, we quantify changes in forest cover between 2001 and 2009 as a proxy for changes in the potential supply of ES from these forests, as all of our other proxies of ES represent a single snapshot in time; this metric gives us some indication of the degree to which incomes that depend on forest provisioning ES can be sustained.
The proxies for the nine ES we consider in this study are crop net income density (per km 2 ), density of large animals, densities of goats and sheep, densities of agricultural production from orchards and plantations ('woody agriculture'), density of non-timber forest product ('NTFPs') harvests, density of timber harvests, density of fuel wood harvests, density of bamboo harvests and natural fresh water availability. These proxies differ considerably in the degree to which they represent realized ES due to issues of data availability. Densities of large animals, goats and sheep, 'woody agriculture', and NTFPs, timber, fuel wood and bamboo harvests are partial measures of ecosystem benefits (Fisher et al. 2009) as these proxies do not distinguish between contributions of the ecosystems and further capital inputs (e.g. fertilizer). However, crop net incomewhich we take as our proxy for agricultural productivityis a good measure of the ecosystem benefits of crop production, as by considering the net rather than the gross value we exclude all capital inputs except labour. By contrast, natural water availability is an indicator of the potential supply of water by the environment. We also calculated the percentage of available natural water used by humans (percentage of available water used), as this gives us an indicator of the degree to which human demand for water is met by the naturally available supply (Burkhard et al. 2012).
Our first measure of biodiversity is extant woody plant richness per county, as this is the only taxa for which reliable data are available at sufficient resolution for China (Fang et al. 2011), Plants are arguably also the taxa that is most directly relevant to ES (Cardinale et al. 2012). Our second measure of biodiversity is the number of threatened species (endemic, endangered and nationally protected species) per county, which we obtained from a recent Chinese study (Ouyang et al. 2016). This latter measure is complementary to plant diversity, as it quantifies both the existence value of biodiversity placed on it by Chinese society (e.g. Anderson et al. 2009), as well as the counties where biodiversity is most vulnerable to anthropogenic change. We calculated densities (in km 2 ) of all ES to adjust for the large differences in the areas of counties (mean area 4177 km 2 ± sd 10,127 km 2 ). A summary of the types of data sets underpinning our biodiversity and ES indicators is given in Table 1; detailed methods are as follows: (1) Crop net income, (2) numbers of large live animals and (3) numbers of live goats and sheep. We obtained these statistics directly from the county-level data available in the Chinese rural statistic yearbook 2009 (National Bureau of Statistics 2009a) as this is the best available nationally consistent data set. All main food and oil crops as well as cotton were included in the estimate for crop net income; the net unit incomes for these main crops are also available in the same yearbook. Large animals include cattle, horses, yaks, mules and donkeys.
We obtained yield statistics for all 31 provinces in China for all these forest and orchard-derived indicators of ES from the China forestry statistics yearbook 2009 (National Bureau of Statistics 2009b) as this is the best freely available nationally consistent data set. We considered all tea, coffee, tea oil, orchard fruit and nut production as woody agriculture production; and all food, condiments, medicine and industrial material that are collected in the forest rather than grown to be NTFPs. Fuel wood, timber, small bamboo and large bamboo are all separate categories in the forest statistics data and were considered as such here.
As the above data are only available for provinces, we had to make a simplifying assumption to spatially disaggregate these data to the county level. This assumption is that a county's share of each of these services is directly proportional to amount of the type(s) of forest in which it grows in each county. We calculated the amount of four relevant types of forest in each countyplantation forest, plantation bamboo forest, natural forest and natural bamboo forestby digitizing the Atlas of forest resources in China (Xiao 2005), which is based on the latest National Forest Inventory of China (1999China ( -2003. Plantation forests include orchards as well as plantations planted to provide timber; so, we assumed that this type of forest cover provides all orchard-derived 'woody agriculture', all timber and 50% of fuel wood. In the absence of primary data, we assumed that natural forests provide all NTFPs and small bamboo (as bamboo forests are classified of consisting only of large bamboo), and 50% of fuelwood. We also assumed that 50% of large bamboo comes from natural bamboo forests and 50% from plantation bamboo forests. For example, if county X had 10% of the plantation forest in Fujian province, we assumed it would therefore also have 10% of Fujian's orchard-derived ESs, 10% of Fujian's timber production and 5% of Fujian's fuelwood production.
(9) Plant biodiversity and (10) threatened biodiversity. We used extant woody plant richness and threatened biodiversity richness as our two surrogates for biodiversity. The richness of woody plants by county was calculated based on the Atlas of woody plants in China (Fang et al. 2011), which is the most comprehensive database of species distributions in China. This atlas documents the county-level distributions of all 11,405 woody plants native to China, and is a compilation of all published country-level and provincial floras (more than 300 volumes in total), and a large number of additional published local floras and inventory reports. We defined species richness of woody plants as the total number of species in the county. Our measure of threatened biodiversity was the recent county-level assessment of the number endemic, endangered and nationally protected species per county (Ouyang et al. 2016).
(11) Natural water availability and percentage of available water used.
We defined natural blue water availability as the water provided by the ecosystem in each county (mm/ year, which is equivalent to litres/m 2 /year), and percentage of available water used by calculating the percentage of this naturally available water that is appropriated by humans (human water use/natural water availability). Counties are considered to be water scarce if the percentage of available water used by humans is between 20% and 40%, and severely water scarce if the ratio of human water use to availability exceeds 40% (Alcamo, & Henrichs, 2002;Müller Schmied et al. 2014), while values above 100% signify that humans are appropriating more water than is naturally available in the county (i.e. through imports or by relying on aquifers). We calculated both human water use and natural water availability at the county level by using the widely used global integrated water model WaterGAP3 (Alcamo et al. 2003;Verzano et al. 2012;Flörke et al. 2013). WaterGAP3 consists of two main components: (1) a water balance model to simulate the characteristic macro-scale behaviour of the terrestrial water cyclethis gives the estimate of natural water availability; and (2) a human water use model based on estimates of water withdrawals and consumptive water uses for agriculture, industry and domestic purposes. The detailed methodology for WaterGAP3 is available in the Supplementary materials.
(12) Calculation of the change in forest cover from 2001 to 2009.
In addition to calculating the current amount of different types of forest cover based on 1999-2003 National Forest Inventory data described earlier, we also used 500 × 500 m MODIS land-cover data (Justice et al. 1998) from the years 2001 to 2009 to calculate the change in total forest cover for each county. We first extracted any pixels of bad quality (i.e. where cloud cover etc. were a problemthese constituted less than 0.2% of all pixels) and then projected and transformed the data to Albers Equal Area projection from the Sinusoidal projection. We then reclassified Evergreen Needleleaf Forest, Evergreen Broadleaf Forest, Deciduous Needleleaf Forest, Deciduous Broadleaf Forest and Mixed Forest into a single forest landcover class and calculated the percentage forest cover in each county for both 2001 and 2009.

Statistical analyses
We used cluster analysis to identify ES bundles of densities of co-occurring ES, following the approach first outlined by Raudsepp-Hearne et al. (2010a). As results of cluster analysis are known to be sensitive to both the clustering algorithm used and the number of clusters selected, we used the 'clValid' package in R 3.31 (R Development Core Team 2016) to select the most appropriate method and number of clusters by comparing multiple methods using multiple validation measures (Hamann et al. 2015). This showed that 'pam' is the most appropriate method, and that the most appropriate number of clusters for our data is two. 'Pam' implements partitioning around medoids, a method which is more robust to outliers than the older K-means clustering methodology (Kaufman and Rousseeuw 1990). We only considered provisioning services in the cluster analysis as our goal was to characterize bundles of ES that were directly relevant to poverty alleviation. These are crop net income density, large animal density, goat and sheep density, woody agriculture density, timber production density, fuelwood production density, NTFP density, large bamboo production density, small bamboo production density and natural water availability. We rescaled all variables to be between 0 and 1 to ensure equal weighting of all variables in the cluster analysis.
We then employed Wilcoxian non-parametric tests to identify the degree to which densities of ES differ between poor and non-poor counties. We ran these analyses both at the regional scale (within each of our two clusters of counties) and for China as a whole. We did not run all possible comparisons to minimize the possibilities of Type I errors and adjusted for multiple comparisons using the Holm adjustment (Holm 1979).
Finally, we used random forest analyses (Breiman 2001) to identify the extent to which the poverty status of a county (poor vs. non-poor) could be predicted based on the density of the provisioning ES we include in the PCA analysis. Again, we ran these analyses both at the regional scale (clusters of counties) and for China as a whole. Random forests are a powerful machine-learning statistical classifier that is well suited to large, complex data sets with large numbers of correlated predictor variables with nonlinear and complex interactions and make no assumptions about the distributions of predictor variables (Cutler et al. 2007). Random forests are becoming increasingly used in ecological studies due to their good predictive performance compared with traditional statistical modelling approaches and robustness (Prasad et al. 2006;Cutler et al. 2007) and have recently been used for ES analyses (e.g. Meacham et al. 2016) to which they are well suited, given the prevalence of the issues (correlations, interactions, non-linearities) in ES data that random forests are ideally suited to overcoming. All random forest analyses were carried out using the randomForest package in R 3.31 (R Development Core Team 2016).
More detailed methods for the random forest analysis are available in the Supplementary materials.

Results
Cluster analysis identifies two broad socioecological subsystems with distinctive bundles of ES (Figures 2  and 3). One cluster (Northern China [NC]; n = 1061) consists of the arid, largely unforested counties of NC, and is dominated by intensive agricultural landscapes, grasslands, steppes and deserts; however, it also encompasses the boreal forests found in northeastern China. The other cluster (Southern Forests [SFs]; n = 863) consists mostly of the subtropical and tropical forest-dominated counties of southeast China ( Figure 2); a few central and northern Chinese forested counties also fall within this cluster. The percentage of poor counties within these two clusters is 28% (NC), 32% (SF) and 29% for all rural Chinese counties (n = 1924).
There are large variations in both the distributions of ES and biodiversity at the national and bundle scale and within individual bundles (Table 1; Figure 3), and in how the distributions of ES and biodiversity vary between poor and non-poor counties (Table A1; Figure 3). The NC cluster has high levels of agricultural production when compared to the SF cluster, but very low water availability, lower forest cover and lower production of most forestrelated ES and much lower biodiversity (Table 1; Figure 3). At the national scale, poor counties in China have significantly lower (45%) crop income densities, goat and sheep densities and levels of human appropriation of naturally available water than non-poor counties, but significantly higher natural and plantation forest cover as well as increases in forest cover, significantly higher harvest densities of small bamboo and NTFPs, and significantly higher biodiversity (Table A1; Figure 3). In contrast to the national-scale pattern, densities of harvests of forestrelated provisioning servicesparticularly large bambooare generally much higher in the non-poor than poor SF counties, as is the amount of bamboo forest and the availability and human appropriation of naturally available water. However, as at the scale of China, biodiversity and increases in forest cover between 2001 and 2009 are much higher, and crop net incomes lower, in the poor than non-poor counties in the SFs. In the NC cluster, both crop income densities and natural water availability are much lower in the poor than in the non-poor counties (Table A1; Figure 3).
The random forest analyses show that high densities of harvests of key provisioning services are a good predictor of a county being classified as 'not poor', both at the scale of China as a whole, and within the SF and NC clusters (Table 2, Figure 4). Crop net income density is by far the most important single predictor of a county being classified as not poor at the scale of China and within both the NC and SF cluster. The density of large bamboo production is the second most important predictor of not being poor at the scale of China; this is due to its importance as a predictor in the SF cluster (Figure 4).
While classification error in the random forest analyses in terms of predicting if a county was 'not poor' is low (14% or less in all three analyses), we are unable to accurately predict if a county was classified as 'poor' (38% or greater classification error in all three analyses; Table 2).

Discussion
At the scale of China as a whole, our results are consistent with the global pattern of low levels of income poverty co-occurring in areas with low levels of non-agricultural ESs (Raudsepp-Hearne et al. 2010a). More specifically, our analyses show that high levels of provisioning ESparticularly those related to agricultural productiongenerally coincide with low levels of rural poverty in China. Our results support economic analysis of household-level time series data from China which shows the importance of farming as strategy for staying out of poverty for the rural population (Imai and You 2014). However, our results also show that having low levels of provisioning ES is in itself not a good predictor of a county being poor. Again, this is not a surprising finding, as income rural poverty is a multifaceted issue that is driven by many factors that we were unable to quantify in this study. These include not only access to markets and remittances from migrants to urban areas (Donaldson 2011), demographics and education (Imai and You 2014), income from sources other than natural resources (De Janvry et al. 2005) but also biophysical factors such as rainfall and elevation (Olivia et al. 2011). However, this broad pattern of co-occurrence of poverty and high forest cover that exists both in China and elsewhere (Sunderlin et al. 2008) can mask positive associations between high harvests of forest-based ES and low poverty, such as we find within the SF cluster in China. This latter finding raises two important questions: (1) Are these forest-  (NCP and NC, respectively) and Southern Forest (SFP and SF) counties. Boxes correspond to the interquartile range of values, while whiskers incorporate extend to the most extreme data point that is no more than 1.5 times the interquartile range away from the edge of the box. Non-overlapping 'notches' indicate 'strong evidence' (Chambers et al. 1983) for significant differences between clusters. The dotted line gives the median value for each service for all rural counties in China.
related ES a sustainable route for poverty alleviation in this region and (2) more generally, does the cooccurrence of high levels of crop net income, natural forest cover, harvests of forest ES and stored carbon in non-poor SF counties support the idea that an agroecological matrix (Perfecto and Vandermeer 2010) can be a sustainable path out of poverty?
Answering both these questions requires an understanding of the relationships between realized and potential ES in this region. The amount of natural forest cover and extant biodiversity of each county are both measures of potential ES, while the harvests of forest ES are measures of realized ES. At present, in non-agricultural Chinese counties, harvests of forest ES and amounts of bamboo forest are both highest in non-poor SFs, suggesting that the high realized ES are largely the result of there being a high potential supply (Burkhard et al. 2012) of forest ES in these counties. A particularly interesting result is the strength of the association between levels of large bamboo production and whether or not a county is classified as poor or not poor within the SF cluster. China has by far the largest bamboo forests in the world, and while the potential economic benefits of bamboo (i.e. for building materials, furniture production, biochar etc.) are well understood from case studies (Song et al. 2011), our study is the first to highlight the relationship between bamboo forests and poverty levels at very large spatial scales. Our results indicate an urgent need for further work into the linkages between bamboo forests and poverty both in China and globally, and for better integration of bamboo production into forest policy, from which it has and continues to be largely marginalized (Buckingham et al. 2011).
However, we do not have information on whether current levels of harvests of provisioning ES are sustainable, or the intensity of management of the forests; both have implications for the non-provisioning ES bamboo (Song et al. 2011) and non-bamboo forests (Carnus et al. 2006) can provide. More work is also needed to understand whether the increase in forest cover in the SFs (Table 1; Figure 3) is mainly due to an increase in forest plantations rather than natural forests (Ying et al. 2010) or bamboo forests, as remotely sensed data cannot reliably indicate forest type; this problem is compounded by the age (1999)(2000)(2001)(2002)(2003) of the forest inventory data on which our estimates of the different forest types are based.
An additional important issue for understanding the role of forests in alleviating poverty is to understand the degree of dependence of the people within these regions on provisioning ES (Yang et al. 2013a). We therefore run post-hoc in analyses where we calculate the per capita rural GDP from 'primary industries' (forestry, fishing and farming) for all counties to enable us to quantify the degree to which rural livelihoods in each country are directly dependent on provisioning ES (Supplementary Methods). The per capita rural GDP from 'primary industries' in the non-poor SF counties is significantly higher than in The classification error is calculated based on the ratio of mismatches to matches in classification. The relative importance of predictor variables is obtained by randomly permutating the values of each variable within each classification tree for the out of bag observations and comparing the difference in the misclassification rates between the randomly permutated variables and the actual valuesthis is the mean decrease in predictive accuracy. The full names of the ecosystem services considered are as follows: Crops: crop net income (RMB/km 2 ); LAnimal: large animals (individuals/km 2 ); GoatSheep: goats and sheep (individuals/km 2 ); Timber: timber production (m 3 /km 2 ); NTFP: NTFP production (t/km 2 ); Fuelw: fuelwood production (m 3 /km 2 ); WoodyAg: fruit, nuts, tea and coffee (t/km 2 ); LBamboo: large bamboo production (stems/km 2 ); SBamboo: small bamboo production (t/km 2 ); WaterAvail: natural water availability (mm/year).
the poor ones, but the percentage of the total per capita rural GDP that comes from 'primary industries' is significantly higher in the poor than nonpoor SF counties ( Figure A2; Supplementary Methods). These results support work in the Wolong Nature Reserve that not only shows that it is the poor who are generally most directly dependent on ES (Yang et al. 2013a) but also suggests that the relative wealth of the non-poor SF counties is a result of greater levels of both natural and non-natural capital in these counties than in poor forested Chinese counties. More work is needed to understand the relative contribution of the provisioning ES we consider here and other non-consumptive ES (such as tourism) to rural incomes in China.
Another important limitation of our results is that they are based on the current spatial co-occurrence of ES harvests and poverty, as not all of the data sets that underpin our analyses (e.g. the atlas of forest resources and the WaterGAP3 model outputs) are available at multiple time points. As such, our results cannot show causality or allow the marginal effects of changes in stocks of ES on changes in poverty levels to be ascertained. Such temporal analyses are required to identify whether it is indeed forest-related provisioning ES that are leading to poverty alleviation, or whether the higher incomes in non-poor counties stem from non-consumptive ES like tourism (e.g. Ferraro and Hanauer 2014) or other non-ESderived incomes sources.
Achieving sustainable rural livelihoods in NC is likely to prove even more difficult than in the SFs, due to the aridity of this region. NC has by far the highest densities of crop net income (the median crop income density in the SF cluster is less than half of that in the NC cluster; Table 1), but less than a third of the median natural water availability of the SFs. Indeed, the median level of appropriation of natural water by humans in NC is 47%, making this region severely water scarce (Jiang 2009;Müller Schmied et al. 2014) as a whole. This water scarcity limits the potential increase of sustainable practices such as rice-fish co-culture (Xie et al. 2011) and also means that this region is likely to become increasingly reliant on water imports from other parts of China; indeed, some counties in NC are already using far more water than is naturally available within their boundaries. Overall, our results add to the growing evidence that maintaining current levels of agricultural production using current practices is not sustainable over the longer term (e.g. Guo et al. 2010;Dearing et al. 2012) within much of NC.
More generally, our results show the importance of spatial context and scale of analysis in quantifying the relationships between ES and poverty, thereby building on the growing literature on the effects of scale on ES (e.g. Anderson et al. 2009;Willemen et al. 2010;Raudsepp-Hearne et al. 2010a;Haase et al. 2012). Our findingstogether with other recent large-scale work (Angelson et al. 2014;Ferraro et al. 2015)therefore highlight the importance of largescale studies for understanding the linkages between poverty and ES. Our work also shows that there is a great deal of spatial variability within rural 'greenloop' systems (regions in which the populations are highly directly dependent on ES) (Cumming et al. 2014;Hamann et al. 2015) in terms of which ESs are most important for rural livelihoods. However, further work is required to ascertain whether the linkages between ES and poverty we observe at the county scale for China hold at finer resolutions. For example, there is evidence that the 'poor county' designation does not adequately represent withincounty variation in poverty levels in China (World Bank 2001). Similarly, within-county ES trade-offs between forest and non-forest-related ES are also likely, as even within-village differences in locations, can affect the level of income that individuals receive from forests (Angelson et al. 2014).
Finally, it is important to note that we only consider one aspect of poverty (income poverty) in our analysis and focus on the subset of ES (provisioning ES) most relevant to this. Recent work in China has shown that the degree of overlap of multidimensional measures of poverty (including factors such as education and access to sanitation) and income poverty (as measured by the same indicator we consider here) is only 31% ). In addition, a recent review on the relationship between human well-being (of which monetary income is only one aspect) shows that most studies globally only focus on one aspect of human well-being and provisioning services (Cruz-Garcia et al. 2017). Indeed, some studies have shown that non-monetary aspects of human well-being can be intrinsically linked to cultural ES and are 'constituents' rather than 'determinants' of human wellbeing (Schleicher et al. 2017). As such, much more work is clearly required on examining the linkages between multidimensional measures of well-being and poverty and a wider suite of ES in China and elsewhere.

Policy implications
These findings have a number of clear policy implications for ES management in China. First, the spatial distributions of the bundles of ES we have identified and their bundle-specific associations with poverty can provide guidance on priority settings for more detailed place-based studies, as well as policy planning and evaluation at the scale of China as a whole, thereby providing a complement to the recent China-wide assessment of changes in ES (Ouyang et al. 2016).
Second, our results have policy implications for the extensive program of PES schemes that exist in China, the majority of which focus on maintaining or increasing forest cover to reduce soil erosion and flooding (Liu et al. 2007(Liu et al. , 2008. Our results suggest that harvests of ES from natural forests may be important for poverty alleviation, so schemes designed to conserve natural forests need to carefully consider their impacts on livelihoods. If such schemes allow sustainable harvests of forest ES, they may be able to offer major policy win-winsreducing rural poverty and soil erosion in the uplands while also providing flood prevention services for the lowlands, maintaining biodiversity and increasing carbon stocks (Liu et al. 2008;Chen et al. 2009;Zhou et al. 2009). However, restrictions on the use of these forests in areas under these schemes mean that such multiple benefits may be difficult to achieve. The Natural Forest Conservation Programme, which aims to conserve natural forests and to encourage the planting of new forests, including bamboo in southern China, has been shown to have fewer positive socio-economic impacts than the Grain to Green afforestation schemes by depriving communities of timber revenues (Liu et al. 2008), and to have the greatest adverse effect on the poorest people (Cao et al. 2010). Achieving poverty reduction in rural China therefore needs to include measures that go beyond PES schemes such as improving medical and educational facilities in poor rural counties (Yang et al. 2013b). In addition, increasing decentralization of the administration and governance of forests may be a way to sustainably increase the poverty alleviation benefits of China's forests, as community forests have been shown to be significantly more likely to provide both local livelihoods and biodiversity benefits when local communities can participate in rule-making in forest governance (Persha et al. 2011).

Conclusions
China has been extraordinarily successful in reducing rural poverty between 1980 and the early 2000s, largely through increases in agricultural production and nonfarm incomes (De Janvry et al. 2005;Imai and You 2014). However, lifting the remaining rural poor out of poverty has proved more difficult as they frequently live within remote areas (Imai and You 2014). Our national-scale analysis confirms that at the scale of China, provisioning ES related to agriculture is a key driver of rural incomes, supporting the hypothesis that in general, high levels of poverty co-occur with areas with high levels of non-agricultural ESs. However, in one of the ES bundleswhich consists of forest-dominated counties in southern Chinalow poverty, high densities of forest-related provisioning services and high levels of natural forest cover co-occur. Further work is urgently required to understand the degree to which offtake of these forest-related ES (particularly bamboo) is sustainable, how important it is for livelihoods at the household level and, crucially, the degree to which it can help lift the remaining rural poor in China out of poverty. However, our inability to predict which counties are poor based on provisioning ES also highlights that much more work is required to understand how other aspects of poverty than income are related to ES, including on the likely important, but difficult to quantify role that cultural ES play in maintaining human well-being.