Greater loss and fragmentation of savannas than forests over the last three decades in Yunnan Province, China

Yunnan Province, southwest China, has a monsoonal climate suitable for a mix of fire-driven savannas and fire-averse forests as alternate stable states, and has vast areas with savanna physiognomy. Presently, savannas are only formally recognised in the dry valleys of the region, and a no-fire policy has been enforced nationwide since the 1980s. Misidentification of savannas as forests may have contributed to their low protection level and fire-suppression may be contributing to vegetation change towards forest states through woody encroachment. Here, we present an analysis of vegetation and land-use change in Yunnan for years 1986, 1996, 2006, and 2016 by classifying Landsat imagery using a hybrid of unsupervised and supervised classification. We assessed how much savanna area had changed over the 3 decades (area loss, fragmentation), and of this how much was due to direct human intervention versus vegetation transition. We also assessed how climate (mean annual temperature, aridity), landscape accessibility (slope, distance to roads), and fire had altered transition rates. Our classification yielded accuracy values of 77.89%, 82.16%, 94.93%, and 86.84% for our four maps, respectively. In 1986, savannas had the greatest area of any vegetation type in Yunnan at 40.30%, whereas forest cover was 30.78%. Savanna coverage declined across the decades mainly due to a drop in open parkland savannas, while forest cover remained stable. Savannas experienced greater fragmentation than forests. Savannas suffered direct loss of coverage to human uses and to woody encroachment. Savannas in more humid environments switched to denser vegetation at a higher rate. Fire slowed the rate of conversion away from savanna states and promoted conversion towards them. We identified remaining savannas in Yunnan that can be considered when drafting future protected areas. Our results can inform more inclusive policy-making that considers Yunnan’s forests and savannas as distinct vegetation types with different management needs.


Introduction
Savannas are ecosystems characterized by a ground layer dominated by grasses and an open canopy of trees that allows direct sun to penetrate to the understory (Scholes and Archer 1997). Competition between trees and grasses asymmetrically favours trees, therefore, savannas are dependent on disturbance to remove tree biomass and open up canopies, either by fire (Higgins et al 2000), herbivory (Sankaran et al 2013), or drought (Fensham et al 2009). In Asia, savannas are extensively found across the drier parts of the monsoonal belt of its tropical and subtropical regions (Ratnam et al 2016). Fossil and molecular evidence indicate that savannas and their associated species have been widely present in the region since at least the early Pleistocene (Shen et al 2009, Edwards et al 2010, Zhang et al 2011, Ratnam et al 2016, Chu et al 2021. Additional pieces of evidence include stochastic gradient boosting models that predict the existence of savannas in Asia using savanna climate envelopes from other continents (Ratnam et al 2016) and dynamic global vegetation models that simulated vegetation patterns using current climate ; the presence of high endemism and richness of C 4 grass species, which typify savannas (Ratnam et al 2016, Welker et al 2020, Chu et al 2021; and several in situ studies that have characterized the physiognomy, diversity, and fire-adapted traits of these communities in South and Southeast Asia (Khaing et  Our research focuses on the parkland (PRK) and woodland (WDL) savannas of Yunnan Province in southwest China (figure 1), an elevated region with a strongly monsoonal and relatively dry climate (Shi et al 2017, Shi andChen 2018) suitable for sustaining savannas and forests as alternative stable states (Bond et al 2005). In China, open-canopied PRK savannas were, until quite recently, only formally recognised in the dry valleys that incise Yunnan, and these have been subject to intense investigation (Jin andOu 2000, Zhu et al 2020). The vegetation outside the dry valleys of Yunnan was previously classified as 'forest' with small patches of treeless grasslands recognised at higher elevations (Hou 2001). Their classification followed the Food and Agricultural Organization of the United Nations (FAO) definition of forests, in which vegetation with >10% tree cover is considered forest (Bastin et al 2017), but it is now well understood that the tree cover at which savanna switches to forest is between 60% and 80% (Hirota et al 2011, Staver et al 2011, Griffith et al 2017, Khaing et al 2019. A new vegetation map of China concluded that much larger swathes of Yunnan are comprised of vegetation with a grass-dominated ground layer (Su et al 2020). Northern Yunnan is presently dominated by pine and oak WDL savannas with denser canopies and a ground layer dominated by C 4 grasses, a vegetation widely distributed across higher elevation tropical Asia (Jin 2002, Ratnam et al 2016. Diversity in the valleys of Yunnan is well-researched with over 2000 plant species recorded (Jin 2002) among which there are at least 188 grass species, of which >80% are C 4 (Osborne et al 2014) dependent on open habitats, and 79 species belong to the Andropogoneae clade, which dominate fire-prone savannas worldwide (Ripley et al 2015). We are unaware of any published community studies of the savanna vegetation on the tablelands intervening the valleys. Nevertheless, there is strong evidence that C 4 savannas have been widely present in Yunnan since at least the early Pleistocene, placing their existence prior to any possible impact by hominins. The antiquity of the valley savannas is supported by fossil evidence indicating substantial C 4 presence and a switch from forest to grassy biomes at about 3.4-2 Ma (Biasatti et al 2012, Yao et al 2012. A phylogeographic study of two common C 4 Andropogoneae grasses in Yunnan has demonstrated that savannas have been present extensively across Yunnan in the valleys and across the intervening tablelands since the early Pleistocene, but possibly as early as the late Miocene (Chu et al 2021). At present, land conservation in Yunnan protects the Hengduan Mountains Biodiversity Hotspot in the northwest and different forest types found throughout the province (Zhao et al 2019a). Only one nature reserve has been specifically gazetted to protect the dry valley PRK savannas of Yuanjiang, and some montane grasslands are protected in the Hengduan Mountains and in limestone karsts in eastern Yunnan. The vast landscapes of WDL savannas are not viewed as conservation priorities, and as a result, their rates of loss under human impacts have not been evaluated.
Numerous assessments of landscape change using remotely sensed imagery have demonstrated rapid rates of human-driven vegetation loss in Yunnan (Xu et al 2007, Liu et al 2014b, Ning et al 2018, Zhang et al 2019a, but these studies have not specifically assessed the spatiotemporal dynamics of PRK and WDL savanna coverage at the provincial scale. In addition, several national environmental policies have been put in place in China aiming for landscape protection, but with possible negative consequences for remaining intact savannas in Yunnan. Firstly, a national-level fire suppression policy has been effectively enforced since the late 1980s (Yi et al 2017), thereby limiting the occurrence of fires that are crucial for suppressing woody vegetation and maintaining savannas as open ecosystems in more humid environments (Bond 2019, Fogarty et al 2020. Second, a national-level ecological restoration project introduced during the 1990s called 'Grain To Green Program' (GTGP) halted the decrease of forests, which may have inadvertently led to significant afforestation of grassland areas in the southwest (Liu et al 2014b). In addition to direct land-use changes and landscape policies, increased atmospheric CO 2 may have been fertilizing the trees in Yunnan, in turn increasing their growth rates and causing increased woody encroachment (Buitenwerf et al 2012). Tree growth is suppressed by increased aridity and can be enhanced by increased temperature when water is not a limiting factor (Williams et al 2010, Panthi et al 2020, so it is probable that the changes in density in intact vegetation are happening faster in wetter and or warmer parts of Yunnan. We analysed vegetation and land-use change in Yunnan for the period 1986-2016 (30 years of change) using Landsat imagery, which is the period over which the fire suppression policy and GTGP became effective (late 1980s and 1990s, respectively). Our specific objectives were to understand: (a) how much of the savanna vegetation had been lost and fragmented relative to the other major vegetation type in Yunnan, namely forest, (b) where savanna changes happened spatially across Yunnan and how much remained intact, and (c) which environmental drivers, either natural or anthropogenic, were associated to the changes in savanna across Yunnan.

Methods
To address the above objectives, we had five goals in our analysis: first was to assess the full extent of savanna physiognomy in Yunnan in 1986. Second was to quantify the rates of attrition and fragmentation of savannas versus forests in the province. Third, we sought to determine whether the savanna vegetation had undergone woody encroachment by explicitly defining two savanna vegetation classes that differ in woody plant canopy cover (PRKs and WDLs), and then measure the rates of change in these classes as well as forests (higher woody canopy cover). Fourth, we tested whether climate variables, fire, and landscape accessibility were correlated with changes in savanna canopy cover and transitions to other vegetation (forest) and land cover types. Fifth, aware of the low level of conservation of savannas in Yunnan, we identified savannas that have remained mostly intact over the last 3 decades as potential sites for future protected areas.

Generating land cover maps and detecting inter-decadal changes
Land cover maps of Yunnan were generated for years 1986, 1996, 2006, and 2016, encompassing 3 decades of change. Since we wanted to classify Landsat imagery as early as 1986, when no adequate reference data to train a supervised classifier was available, we used a hybrid of unsupervised and supervised techniques combined with knowledge-based interpretation (figure 2; see description of Yunnan and detailed methods in supplementary (available online at stacks.iop.org/ERL/17/014003/mmedia)). We first defined nine land cover classes in the province using a combination of field observations, expert knowledge coupled with visual inspection of high-resolution imagery from Google Earth (www.google.com/earth/desktop, accessed on 5 July 2019) (Olofsson et al 2014), and existing land cover maps (Liu et al 2003, 2014b, Xu et al 2005, Diallo et al 2009, Zhao et al 2012, Lu et al 2015, Ning et al 2018, Zhang et al 2019a, Su et al 2020. We defined nine land cover classes, comprised of non-vegetated classes, namely water bodies (WATs), snowed regions (SNO), and built-up and bare rock areas (BURs); non-natural vegetation areas, namely croplands (CROs), tree plantations (TRPs), and bare ground (BAG); and natural vegetation areas, which are forests (FOR), which have closed tree canopies, and two savanna classes, which have open canopies: grassy, sparse-canopied PRKs and denser-canopied WDLs. See supplementary tables 1 and 2 for their full descriptions.
Image composites used in the classification were generated for each mapping year from Google Earth Engine (GEE; http://earthengine.google.com) using a median ee.Reducer function that composited tier 1 calibrated top-of-atmosphere reflectance images obtained by Landsat-5 Thematic Mapper and the Landsat-8 Operational Land Imager available in the GEE image collection (Loveland and Dwyer 2012). The composites were comprised of seven spectral bands (BLUE, GREEN, RED, NIR, SWIR1, SWIR2, TIR) and were resampled from the original 30 m spatial resolution to 100 m using nearest neighbour method. Five spectral indices, which include normalized difference vegetation index (NDVI), enhanced vegetation index, soil-adjusted total vegetation index, normalized difference tillage index, and land surface water index, were also calculated for each study year and used to improve classifications (see supplementary S1.3). Each image stack per mapping year contained 12 image layers, consisting of the seven spectral bands and five indices.
To delimit the classification algorithms only to vegetation classes, we masked the non-vegetated classes (WAT, SNO, BUR) from the image stacks. These classes were separated using their NDVI (Running et al 1995) and thermal band (Price 1981) values determined using a decision tree algorithm (De Alban et al 2018; see supplementary S1.4.1). To aid in our delineation of regions-of-interest (ROI) polygons, we first ran unsupervised classification on the masked image stacks to generate a computerautomated classified image with 100 clusters in GRASS GIS (GRASS Development Team 2020). We then assigned each cluster to the remaining vegetated land cover classes (both natural and non-natural) through visual interpretation conducted by crosschecking the pixels of each cluster with true-colour (RED-GREEN-BLUE) and false-colour (NIR-RED-BLUE) image composites, as well as high-resolution imagery from Google Earth (available for 2006 and 2016) (Olofsson et al 2014).
The resulting maps from the unsupervised classification were then used to delineate ROIs for supervised classification. We performed stratified random sampling over the pixels from the unsupervised map to generate 500 points used as guides in drawing the ROIs. Each ROI (at least 1 km 2 ) was drawn around the points, with their extents delineated based on the potential homogeneity of the target land cover class. When delineating ROI polygons for 1986 and 1996, we utilised additional mask layers generated from 2016 image statistics to assist in locating TRP and WDL classes due to challenges in discerning them during visual interpretation of the older Landsat image composites and lack of high-resolution imagery for those years (see supplementary S1.4.4). ROIs were then divided for training (70%) and testing (30%).
For supervised classification, we used random forest (RF) (Gislason et al 2006, Rodriguez-Galiano et al 2012 (see supplementary S1.4.5 for RF parameterisation settings). We performed pixel-based accuracy assessment on the RF products by generating confusion matrices for each year, and computed the overall accuracy, user's and producer's accuracies, for each land cover class. We computed error-adjusted area estimations and the confidence intervals for each land cover class to quantify estimation uncertainties (Olofsson et al 2013). Vegetation maps from RF and the non-vegetation maps masked out earlier were combined, and isolated pixels were smoothed out using a Majority Filter tool to produce final maps. The raster operations, supervised classification, and accuracy assessments were performed in R (R Core Team 2021), while smoothing was performed in ArcMap (ESRI 2011).
To detect inter-decadal land cover change dynamics, we utilised the final maps to produce crosstabulation matrices that summarised the land cover changes from 1986 to 2016 using the Semi-Automatic Classification Plugin in Quantum GIS (Congedo 2014). We used the older map as the 'reference map' and more recent map as the 'classification map' (i.e. to detect change between 2006 and 2016, we used the 2006 map as 'reference' and the 2016 map as 'classification'). To visualize the land cover changes, we created Sankey diagrams using the networkD3 package (Allaire et al 2017) in R.

Assessing the fragmentation and legal protection levels of savannas and forests
We calculated fragmentation statistics for natural vegetation classes PRK savanna, WDL savanna, and forest (FOR) for each mapping year. These included number of fragments, size of fragments (mean patch area, largest patch index), geometric complexity (mean patch shape ratio), physical connectedness (patch cohesion index), and edge density. These were calculated using the landscapemetrics package in R (Hesselbarth et al 2019).
We also determined the level of protection that each land cover types received in 2016 by calculating their areal extents within Yunnan protected areas provided by CAS-RESDC (www.resdc.cn/ data.aspx?DATAID=272, accessed on 16 May 2021). We also located persistent savanna areas by identifying unchanged pixels of PRK and WDL since 1986, and subsequently derived the number of decades a pixel had remained as PRK or WDL.

Determining drivers of change in savanna coverage
We conducted logistic regression analyses to assess whether retention or conversion of savannas to other land cover types were due to selected environmental parameters, including climate (mean annual temperature (MAT), aridity), topography (slope), presence of fire (hereafter referred to only as 'fire'), and accessibility (distance from the nearest road (hereon referred to only as 'distance')). MAT and aridity both impact plant growth rates (Williams et al 2010, Panthi et al 2020, and hence rates of vegetation change. Fire reduces tree biomass, which helps to maintain open ecosystems such as savannas (Higgins et al 2000, Lehmann et al 2009. Slope and distance both constrain human accessibility and thus direct human impacts (Liu et al 2014a, Alphan 2017. To assess drivers of woody encroachment, we analysed transitions of savanna to denser vegetation (PRK to WDL, PRK to FOR, and WDL to FOR). To assess drivers of fragmentation from land use conversion, we investigated transitions of PRK and WDL to/from farmlands and non-vegetation classes. In these analyses, we lumped together CRO, TRP, and BAG under one 'farmland' class, and BUR, WAT, and SNO under one 'non-vegetation' class. To quantify land cover change, we reclassified the transition maps into binary according to the transition being analysed. Pixels that retained their classification in the next time step were scored '0' , while those that transitioned to another were scored '1' .
As we were interested in state changes between two levels (i.e. 0 and 1) in each analysis, we tested predictors of land cover transitions using logistic regression in R. Environmental data of both '0' and '1' pixels were first extracted using the extract() function of the raster package in R. The pixels environmental data served as predictors, while their land cover transition scores represented the change variable in the regression models. We used MAT, aridity, slope, and fire in analysing PRK and/or WDL conversion to denser-canopied WDL and FOR to test whether these predictors were responsible for inducing or preventing canopy closure. We added distance in assessing PRK and WDL conversion to/from farmland to investigate whether accessibility influenced land cover change. In testing PRK and WDL conversion to/from non-vegetation, we used slope, distance, and fire. We transformed predictors as needed to approximate normality prior to regression. To assess the robustness of predictor coefficients, we ran the model for each type of transition 1000 times using randomly subsampled equal populations of 1's and 0's for each run, and constructed 95% quantile intervals for each coefficient; 95% quantile intervals which did not overlap zero were considered significant. The explanatory power for each significant predictor was computed from the median of all the partial R-squared values calculated using the partR2 package (Stoffel et al 2021) across all runs. We analysed changes across 3 decades (1986-2016) and decadal changes (1986-1996; 1996-2006; 2006-2016). Analyses were conducted at 1 km resolution to match the native resolution of available environmental data.
Sources of environmental data are listed in supplementary S2.

Net changes over 3 decades, 1986-2016
In 1986, the majority of Yunnan was covered by natural vegetation (figure 3). The combined savanna classes covered 40.30% of the landscape, with PRK savannas at 27.09% and WDL savannas at 13.21%, while FOR was at 30.78%. There were also substantial CROs at 17.73%, while TRPs had low cover at 3.33%. BAG covered 5.38%, whereas BUR was at 1.72%. WATs and SNO regions had the lowest coverage at 0.48% and 0.28%, respectively. By 1996, savannas had a net increase of 2.26%, as PRK cover increased by 4.29% and WDL declined by 2.03%. FOR cover increased by 1.69%. CROs dropped by 4.85%, but BAG increased 1.67%. In 2006, savannas declined slightly due to a 6.32% drop in PRK areas despite a 5.35% increase in WDLs. Decreases in FOR (0.92%), CRO (0.74%), and BAG (1.47%) coverage occurred while TRPs (3.91%) and BURs (0.27%) increased.

Losses and gains in savanna coverage
Extensive land cover changes were observed in Yunnan between 1986 and 2016 (figure 4; see supplementary S5 for the land cover change matrices), and a total of 316 800 km 2 (74.7% of Yunnan's total land area) underwent change. The gross gains and losses  within each land cover type indicated dynamic transitions throughout the study period. The observed decline in PRK savanna cover was largely due to woody cover increases to WDL savanna (5.60%) and FOR (3.66%), and direct land conversion to CRO (2.88%), BAG (4.04%), plantations (0.98%) and BUR (0.31%). Most PRK savanna that was persistent until 2016 was on the eastern and northwest sections of Yunnan (figures 5(a) and 6). Gains in PRK savanna were scattered across the upper section of province, with a concentration of converted WDL savanna and CRO to PRK savanna in the east ( figure 5(b)).
Total WDL savanna coverage was very dynamic. Persistent WDL savanna was ∼6% of 13.21% in 1986, but the remainder changed substantially mostly due to tree cover changes between PRK savanna and FOR states, and also through conversion to TRPs in the south and northeast (figure 5(c) and (d)).

Fragmentation and protected area coverage analyses of natural vegetation cover types
Fragmentation generally increased for all natural vegetation land cover types from 1986 to 2016, but was more severe for savannas (table 1). The number of patches increased during this period for PRK savannas (29.3%), WDL savannas (59.8%), and FOR (23.9%). Fragment sizes also generally decreased, most severely for PRK, with mean patch area decreasing 52.9% (from 20.25 ha to 9.53 ha) between 1986 and 2016. The mean patch areas of WDL (25.3%) and FOR (17.8%) decreased less. The largest patch index for PRK also decreased dramatically, from 6.84 to 0.33 between 1986 and 2016. In contrast, the largest patch index values of WDL (0.17 vs 0.15) and FOR (2.78 vs 2.97) had changed less during the same period. The edge density of WDLs (21.30 vs 29.77) and FORs (24.08 vs 25.57) increased, but decreased for PRKs (27.16 vs 25.56). Fragment connectivity of savannas decreased, with PRK (99.65%-96.23%) experiencing greater loss than WDL (93.28%-91.87%), whereas it remained stable for FOR (99.52%-99.43%). The geometric complexity of the fragments, as indicated by the mean patch shape ratio, was largely stable across time for savannas and forests.
The protected area coverage analysis indicated low legal protection percentages for all natural vegetation land cover types for 2016. Of the three, FOR had the most legal protection, with 7.60% (10 127 km 2 out of 133 160 km 2 total 2016 coverage) inside protected area boundaries. The savannas had lower protection, with PRK at 3.09% (2161/70 022 km 2 ), and WDL at 1.84% (1233/66 824 km 2 ). In addition, most of these nature reserves do not cover the persistent savanna fragments (figure 6). PRK protection is only 11.96% (2161 km 2 out of 18 064 km 2 total nature reserve area), while WDL is at 6.82% (1233/18 064 km 2 ). In contrast, FOR protection is at 56.06% (10 127/18 064 km 2 ), a clear indication that savannas are underrepresented in most of the protected areas.

Drivers of change in savannas
The logistic regression models indicated varying significant associations between the different environmental parameters and land cover changes between and across the 3 decades (table 2). In the PRK savanna to WDL savanna and FOR transitions, the general trend in between decades indicated that PRK savannas located in warmer and wetter (more positive aridity index) areas had a greater probability to transition to denser-canopied vegetation (WDL savanna or FOR). Although the overall analysis suggested that PRK in drier areas transitioned to WDL, aridity only weakly explained the variation (lowest partial R 2 ) in this case. For WDL savanna to FOR transitions, WDL in wetter regions converted to FOR. The effect of temperature differed in between decades; however, the overall transition suggested that WDL in colder areas became FOR. While the effect of slope varied per decade, the overall analyses suggested that PRK in gentler slopes (negative trend) converted to WDL, and WDL in steeper slopes (positive trend) converted to FOR. Fire had a negative effect on transitions to denser canopy cover. In the majority of these transitions, MAT was the predictor explaining the most variation.
Fire was positively related to conversion of WDL and FOR to more open PRK and WDL, but the rest of the predictors had variable effects in these transitions. The analyses suggested that WDLs in colder, wetter, and steeper regions converted to PRKs, and that FORs in drier and flatter places switched to PRK and WDL. MAT had a generally positive trend across the decades, suggesting that transitions to more open vegetation happened in warmer areas. Overall, FOR to PRK transitions happened in colder places, while FOR to WDL occurred in warmer spots. MAT and aridity were the most important predictors in explaining variation. PRK and WDL that converted to farmlands during 1986-2016 were generally in warmer, wetter areas close to road networks and had no fire occurrence. PRK that converted to farms were in flatter places, while WDL that did were in steeper slopes. On the reverse side, farmlands that became PRK or WDL were in colder, steeper regions far from roads and had experienced fire. Moreover, there is a clear trend of wetter farmlands becoming WDL, while the trend for PRK is unclear. In these transitions, MAT, aridity, distance, and slope explained most of the variation in the data.
PRK conversion to non-vegetation areas were also in gentler slopes and near roads. Transitions to nonvegetation classes from 1986 to 2006 happened in colder PRKs. However, the overall transition suggested that PRK in warmer and wetter areas converted to non-vegetation. The overall WDL to nonvegetation transition had similar results, although there were insufficient data to analyse inter-decadal transitions for 1986-2006 due to the small percentage of WDL areas that were converted to non-vegetation areas. Conversely, non-vegetation areas that became PRK and WDL were located in steeper areas and far from roads. In addition, non-vegetation areas in warmer areas became WDL overall. Distance and slope explained most of the variation, whereas fire yielded non-significant associations.

Accuracy
While the overall accuracy numbers we obtained were acceptable for our analyses and comparable to other large-scale land cover studies (Bruggeman et al 2016, Su et al 2020), we observed variations among decades. We attribute these to the different locations of the training ROIs used for each time step, despite ensuring the stability of their spectral signatures (see boxplots in supplementary S1.4.3). It was easier to discern between classes for 2006 and 2016, with high-resolution imagery available to guide ROI sampling, but more challenging for 1986 and 1996. Nevertheless, the UA and PA numbers we obtained for PRK and WDL savannas were well within reasonable range, and thus reliable for analysing trends in savanna cover. Focusing on UA/PA of target land cover classes, which were savannas in our case, is included in the 'good practice' recommendations for area change and accuracy assessment (Olofsson et al 2014). The confidence intervals from the area-adjusted accuracy computations indicate that misclassified areas only comprise a small percentage of the total area of the savanna classes. These misclassified areas were between PRKs and CRO/BAG (see confusion matrices in supplementary S4), which could be due to similarity in spectral signatures of these classes, as grasses are the main vegetation in both classes. During field observations, we noticed some farmlands occur interspersed in between PRKs, which may have contributed to the confusion given the 100 m resolution used in the study. A finer-scale classification, however, would be less useful as our aim is to depict trends in savanna dynamics of a large area (∼400 000 km 2 ); the environmental data are also only available in coarser resolution, thus, the land cover data would always be aggregated to accommodate the drivers of change analyses. We improved the separation of these classes by including spectral indices in the classification.

Dynamic land cover changes in Yunnan exposed the loss of savanna coverage
Our spatial analyses showed that savannas were much more extensive in Yunnan than previously mapped, accounting for 40.3% of Yunnan's land area in 1986, and that savannas were not only limited to the hot Table 2. ANOVA results for generalized linear regression analyses testing whether conversion of land cover types are driven by MAT, aridity, slope, presence of fire, or distance from roads. To achieve a normal distribution required by the model, aridity and slope were transformed using square root, while distance from roads was log10-transformed. For each land cover transition, columns under the time periods list each significant (P < 0.001) predictor's coefficient sign and partial R 2 , while ns means not significant, and nd means no data, thus not included in the analysis.
Overall 1986-1996 1996-2006 2006-2016   dry valleys of the province. This is direct land cover evidence that savannas indeed exist in Yunnan, as postulated by spatial models derived from climate spaces of savannas from other continents (Ratnam et al 2016) and by larger scale dynamic global vegetation model simulations of tropical east Asian vegetation . The more detailed canopybased classification scheme we used divided the previously mapped 'forests' based on the FAO definition, and thus gave more resolution to the vegetation types existing in Yunnan. PRKs occupied almost the entire northern and eastern regions of the province, then extended westward to Kunming Prefecture and the Dali-Lijiang area in the northwest, and southward to Lincang. WDLs were extensive in the northeast along the border with Guizhou province, in the Wenshan area on the southeast, and occurred as transitions between PRKs and FORs in other parts of Yunnan. Savannas declined by 8% in coverage over the 30 years due to the overall decrease in PRK areas, despite a small increase of WDLs. We observed two notable losses in PRK cover. First was the massive loss of PRKs converted to BAGs and CROs, which occurred concurrently with conversions of cropland to urban areas and plantations. This is corroborated by local-scale studies in the Hengduan mountains during 1990-2010 (Wang et al 2018(Wu et al 2015, and in the Dianchi Lake watershed where a 62.1% decline in grassland cover during 1974-2008 was observed (Zhao et al 2012), mostly through conversion to urban and agricultural land. These PRK to farmland conversions occurred in warmer, gently sloped areas accessible via road networks, whereas PRK to built-up area conversions took place in plateaued areas near roads. The most significant drivers explaining these transitions were in areas with relatively warmer temperature and closer proximity from roads.
Second, substantial areas of PRK converted to WDL and then to FOR, indicating woody plant encroachment across Yunnan Province. These transitions occurred throughout the study period, but most were apparent during 2006-2016, and in more humid environments that experienced no fire. This was observed in a local-scale land cover assessment of Greater Kunming, where 1200 km 2 of grassland were converted to forest and thickets from 2003 to 2010 (Lu et al 2015). These transitions may possibly be due to afforestation, fire suppression, atmospheric CO 2 enrichment, or a combination of all. In our analysis, PRKs in steeper regions converted to WDLs during 1996-2006, coinciding with previous observations that open grasslands converted to woodlands in montane areas of southwest China due to implementation of the GTGP during this period (Liu et al 2014b). Fire was negatively associated with transitions to denser canopy cover vegetation and positively associated with transitions to more open vegetation, suggesting fire could promote the persistence of open PRKs in Yunnan, as in other ecosystems (Fogarty et al 2020). Atmospheric CO 2 enrichment increases the growth efficiency of C 3 trees (Ehleringer et al 1997, Hoffmann et al 2000, Kgope et al 2010, and is suggested to be responsible for woody encroachment in savannas in Africa (Buitenwerf et al 2011). Model simulations of tropical and subtropical Asian vegetation have also suggested that elevated CO 2 is supporting woody plant biomass increases and vegetation transitions .

Fragmentation of savannas in Yunnan is greater than for other vegetation types
Although all three natural vegetation types experienced greater fragmentation over time, PRKs suffered the worst degradation, breaking into smaller patches with greater exposed edges and less connectivity. WDLs suffered a massive increase in patch number during 2006-2016, after being relatively stable from 1986 to 2006; however, their connectivity slightly increased throughout the decades. Fragment numbers also increased in FORs, but the largest forest fragments increased and patch cohesion barely changed throughout the decades. The decrease in connectivity of savanna fragments can negatively affect grassland species, which respond sensitively to fragmentation introduced by woody species (Fuhlendorf et al 2002, Cunningham andJohnson 2019), and can introduce disruptions in genetic flow between populations (Honnay et al 2007, Helm et al 2009, affecting the biodiversity supported by these habitats (Jin 2002, Simon et al 2009, Ratnam et al 2016. To our knowledge, these are the first fragmentation statistics of Yunnan's savannas. Existing fragmentation research in the region only focused on forests (Liang et al 2014, 2020, Zhang et al 2019a, Zhao et al 2019b. Our results suggest that fragmentation of savannas is severe and requires urgent attention.

Conservation recommendations
Our analysis suggests that Yunnan's forest cover was largely stable, and even experienced a slight increase over the 3 decades. Forests also possessed a higher level of protection among the two vegetation types. These, combined with the fact that forest fragmentation was also less compared to savannas, imply that the implementation of the GTGP and the protected area systems in place are effectively protecting forests, despite the low protection coverage relative to the total forest land area. The province might also continue to observe an uptick in stable forest cover, as China continues its existing programs to conserve and expand forests (Chen et al 2019). Policies affecting Yunnan include the sloping land conversion program (SLCP) to prevent cultivation on land with slopes steeper than 25 • (Xu et al 2005), and the GTGP established in 1999 to restore natural ecosystems through the return of former croplands back to forests or savannas (Chen et al 2019). A good policy model for savannas could be the SLCP, as our drivers analysis showed that steeper slopes did in fact contribute positively to PRK and WDL conversion from farmlands and non-vegetation areas. Wellinformed implementation of the GTGP on former farmlands rather than PRKs would also be helpful. A policy regarding the controlled use of fire would also be beneficial for savanna protection, as its presence was demonstrated to be favourable for retaining savannas. The use of fire for maintaining savannas is well-researched in other countries and demonstrated to be helpful for suppressing woody encroachment (Wilgen et al 2004, Andersen et al 2005, Schmidt et al 2018. Considering the low protection presently given to savannas, it is vital to establish protected areas that are savanna-inclusive. We located several areas in the central and eastern side of the province where savannas have persisted as candidates for potential protected areas (figure 6). While there are intact savannas that are currently unprotected, such as in central and southeast Yunnan, there are PRKs in the north that would benefit from the expansion of adjacent or nearby existing protected area boundaries.
Finally, a study that looked at how land cover change is tied up with ecosystem service values (ESVs) in northeast Yunnan stated that increasing grassland proportion with respect to forestland would increase ESVs in certain counties (Wang et al 2018), which underscores the contributions of savannas and forests for human well-being and biodiversity conservation. A Yunnan-wide assessment of ESVs is recommended, for understanding potential synergies and trade-offs among multiple ecosystem services in savannas and other land cover types can help improve conservation priorities (Eastburn et al 2017), as well as land and fire management strategies that may encourage the participation of communities as stewards of savannas (Sangha et al 2021).

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).