ARTICLE https://doi.org/10.1038/s41467-019-08602-6 OPEN Ecosystem structural changes controlled by altered rainfall climatology in tropical savannas Wenmin Zhang1,2, Martin Brandt 1, Josep Penuelas 3,4, Françoise Guichard5, Xiaoye Tong1, Feng Tian 1,6 & Rasmus Fensholt 1 1234567890():,; Tropical savannas comprise mixed woodland grassland ecosystems in which trees and grasses compete for water resources thereby maintaining the spatial structuring of this ecosystem. A global change in rainfall climatology may impact the structure of tropical savanna ecosystems by favouring woody plants, relative to herbaceous vegetation. Here we analysed satellite data and observed a relatively higher increase in woody vegetation (5%) as compared to the increase in annual maximum leaf area index (LAImax, an indicator of the total green vegetation production) (3%) in arid and semi-arid savannas over recent decades. We further observed a declining sensitivity of LAImax to annual rainfall over 56% of the tropical savannas, spatially overlapping with areas of increased woody cover and altered rainfall climatology. This suggests a climate-induced shift in the coexistence of woody and her- baceous vegetation in savanna ecosystems, possibly caused by altered hydrological condi- tions with significance for land cover and associated biophysical effects such as surface albedo and evapotranspiration. 1 Department of Geosciences and Natural Resource Management, University of Copenhagen, 1350 Copenhagen, Denmark. 2 International Institute for Earth System Sciences, Nanjing University, 210023 Nanjing, China. 3 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, 08193 Catalonia, Spain. 4 CREAF, Cerdanyola del Vallès, 08193 Catalonia, Spain. 5 Centre National de Recherches Météorologiques (CNRM), UMR CNRS 3589 and Météo-France, 42 Avenue Gaspard Coriolis, 31057 Toulouse, France. 6 Department of Physical Geography and Ecosystem Science, Lund University, 22362 Lund, Sweden. Correspondence and requests for materials should be addressed to W.Z. (email: wenminzhg@gmail.com) NATURE COMMUNICATIONS | (2019)10:671 | https://doi.org/10.1038/s41467-019-08602-6 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08602-6 T ropical savanna ecosystems are characterized by the area coverage). We documented a global increase in LAImax over coexistence of woody and herbaceous vegetation1,2, the tropical savannas from 1992 to 2012 (Fig. 1a). Annual minimum distribution and abundance of which are controlled by vegetation optical depth (VODmin), which has been shown to be herbivores, fire and rainfall3,4. Changing patterns in one of these closely related to woody cover in both drylands and humid controlling factors can lead to an imbalance of this coexistence, regions21, was used as a proxy for woody vegetation. VODmin ultimately causing changes in the structure and functioning of senses the water content of all woody vegetation components, savanna ecosystems. It has been well documented that increasing including the stem and branches, and the impact of herbaceous pressure by herbivores (in particular livestock) favours woody plants on the VODmin signal is weak21. We found that VODmin plants at the detriment of herbaceous vegetation, causing a phe- estimated woody cover increase was significantly higher nomenon termed woody encroachment, which is widespread (expressed by z-scores) than the increase in LAImax (p < 0.05), across global savannas5 and in particular in Sub-Saharan Africa6. especially for areas with a mean annual rainfall less than 1000 Moreover, a declining fire frequency has been observed7, also mm yr−1 (Fig. 1a and Supplementary Figs. 2 and 3) (spatial known to sustain the growth of woody vegetation in savanna distribution of global mean annual rainfall is shown in Supple- ecosystems5. Recently, increased atmospheric CO2 concentrations mentary Fig. 4). The increase in LAImax was evenly distributed were suggested to increase the growth rate of trees relative to across rainfall regimes above 500 mm yr−1, whereas woody cover grasses5,8–10. The impact from changes in rainfall on the vege- increases were strongest for arid/semi-arid regions (Fig. 1a, b). An tation production (comprising both herbaceous and woody foli- increase in annual rainfall was observed for all humidity zones, age mass) is well documented from analyses covering global with relatively higher changes in the drier environments (Fig. 1b). drylands11. Likewise, global syntheses of the impact on woody A decline in woody cover was observed for high-rainfall zones vegetation (based on benchmark tree cover maps) evidence the (>1800 mm yr−1), presumably caused by deforestation21, but also close relationship between rainfall climatology and the density of higher rainfall in moist tropical drylands can lead to reduced woody plants12. At local scale, it has been found that changing productivity23. We further analysed changes in the sensitivity of rainfall climatology (more specifically the increase in heavy LAImax to annual rainfall (βLAI-Rainfall), expressed by the slope of a rainfall events) in drylands has led to increased growth of woody linear regression between LAImax and annual rainfall over a 15- plants2. This is because woody plant vegetation does not require year moving-window period. We found that a majority (56%) of an even temporal distribution of rainfall over the wet season to the global savanna areas experienced a reduced sensitivity over the same degree as herbaceous vegetation and is able to use water the period 1982–2015, and 68% of these areas (38% of global from the deeper soil layers2. Consequently, altered rainfall cli- savannas) coincided with areas of increased woody cover. More matology supposedly can give woody plants an advantage over specifically, the higher the change in annual rainfall and increase herbaceous vegetation with shallow rooted systems13. This in woody cover, the stronger the observed decrease in sensitivity advantage might also be associated with a decrease in the sus- between LAImax and annual rainfall (Fig. 1c). In particular, the ceptibility of juvenile woody vegetation to the controls induced by semi-arid region showed a predominance of areas of positive fires and herbivores3. trends in woody cover and a reduced sensitivity between LAImax A global increase in rainfall extremes has been observed in and annual rainfall, whereas for arid and humid regions the recent decades4,14, which may be expected to affect the structure distribution of areas with reduced and increased sensitivity of and functioning of terrestrial ecosystems2,15,16, including water LAImax to annual rainfall was more balanced (being most pro- balance15, vegetation productivity17 and ecosystem composi- nounced for humid regions) (Fig. 1d–f). tion16. Yet, our current knowledge on the impact of rainfall extremes on the competition between woody and herbaceous plants is based on only a few plot-based experimental Rainfall climatology driving changes in woody cover. Taken studies2,18,19, as it is challenging to assess woody vegetation together, our results suggest that a shift in ecosystem structure trends from global long-term optical satellite imagery systems. has occurred under changed rainfall conditions in tropical Recently, a few studies exploited the complementary use of savanna areas during recent decades. The relative weights of the optical and passive microwave satellite data to study temporal underlying factors controlling the observed increase in woody changes in woody cover20,21. However, no global scale study has cover (Fig. 2a) were determined by a bootstrapping technique yet documented temporal long-term changes in both woody based on the LMG (Lindeman, Merenda and Gold) method24. vegetation and green vegetation production and examined the We found that factors related to inter-annual dynamics in underlying controls, including an analysis of changing rainfall rainfall climatology (including annual rainfall, heavy rainfall climatology. frequency and rainy days) were primary factors controlling Here, we use satellite data of rainfall, LAImax and woody cover changes in woody cover with relative weights of 47% for all for global tropical savannas (see Methods) and find a higher tropical savannas, 78% in arid, 67% in semi-arid and 36% in increase in woody vegetation as compared with the increase in humid regions (Fig. 2a). Factors related to changes in air LAImax. We further observe a widespread declining sensitivity of temperature and solar radiation, human management (human LAImax to annual rainfall for tropical savannas. These results population density and fire changes) as well as static biotic- and suggest that altered hydrological conditions have caused a shift in abiotic-site conditions (mean annual rainfall, rainfall varia- the coexistence of woody and herbaceous vegetation in savanna bility, soil organic carbon, elevation and sand fraction) together ecosystems. accounted for the remaining relative weights (53% for all tro- pical savannas, 22% in arid, 33% in semi-arid and 64% in humid savannas). Change in rainfall intensity was excluded Results from the analysis due to its close relationship with change in Changes in LAImax and woody cover. Here, we used the growing annual rainfall (Supplementary Table 1). For comparison, lower season maximum leaf area index (LAImax) derived from optical contributions of all explanatory factors to changes in LAImax satellite data to study changes in savanna photosynthetic primary were observed (Supplementary Table 2). The contribution of production. LAImax is an indicator of the total green produc- each explanatory variable to changes in woody cover is shown tion22, including the foliage mass from both herbaceous and in Supplementary Fig. 5. Although grazing has shown to impact woody vegetation (see Supplementary Fig. 1 for tropical savanna on the dynamics of woody cover in Southern Africa25, the lack 2 NATURE COMMUNICATIONS | (2019)10:671 | https://doi.org/10.1038/s41467-019-08602-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08602-6 ARTICLE a 0.12 b 0.3 c Woody Woody 0.02 βLAI-Rainfall LAImax 40 20 Area (× 1000 km2) LAImax Rainfall Area (× 10,000 km2) Area 0.2 Rainfall Area Trend (z-score) Trend (z-score) Trend (z-score) 0.08 30 0.01 10 0.1 20 0.04 0.00 0 10 0.0 0.00 0 –0.1 –10 –0.01 Arid Semi-arid Humid –0.02 0 500 1000 1500 2000 –0.2 –0.25 0.00 0.25 0.50 Arid Semi-arid Humid Mean annual rainfall (mm yr–1) Regions Trend in woody cover (% yr–1) d 0.08 e 0.08 f 0.08 Density Density Density 70 50 Trend in βLAI–Rainfall Trend in βLAI–Rainfall Trend in βLAI–Rainfall 50 0.04 50 0.04 30 0.04 30 30 (z-score) (z-score) (z-score) 10 10 10 0.00 0.00 0.00 –0.04 –0.04 –0.04 Arid Semi-arid Humid –0.08 –0.08 –0.08 –0.4 –0.2 0.0 0.2 0.4 0.6 –0.4 –0.2 0.0 0.2 0.4 0.6 –0.4 –0.2 0.0 0.2 0.4 0.6 Trend in woody cover (% yr–1) Trend in woody cover (% yr–1) Trend in woody cover (% yr–1) Fig. 1 Changes in the relationship between vegetation and annual rainfall. a Changes in normalized woody cover and LAImax along a mean annual rainfall gradient (50 mm yr−1 steps) for 1992–2012. Arid, semi-arid and humid regions were defined as < 300, 300–700 and > 700 mm yr−1, respectively. b Comparison of changes in normalized woody cover, LAImax and annual rainfall for arid (n = 3853), semi-arid (n = 8159) and humid regions (n = 15,824). c Changes in the sensitivity of LAImax to annual rainfall (expressed by the linear regression slope over 15 years: βLAI-Rainfall) and normalized annual rainfall for 1982–2015 along with changes in woody cover (10% yr−1 steps). Normalizations were calculated by the z-score: Xi = (xi–x)/σ(x) (where x and σ(x) are the mean and standard deviation of xi). d–f Spatial agreement of areas with changes in woody cover and areas with changes in βLAI-Rainfall for arid, semi-arid and humid regions. Lines in b from top to bottom represent the maximum, third quartile, median, first quartile and minimum values a 100 b c Arid Semi-arid Humid 100 Relative weight (%) 80 0.6 Rescale trend 80 Sensitivity 60 60 40 0.2 40 20 20 0 0 –0.2 All Arid Semi-arid Humid 0 500 1000 1500 2000 All Arid Semi-arid Humid Regions Mean annual rainfall (mm yr–1) Regions Woody cover Rainfall climatology T&S Rainfall Heavy rainfall frequency Rainfall Heavy rainfall frequency Management Static Rainy days Rainfall intensity Rainy days Fig. 2 Drivers of changes in woody cover in tropical savannas. a Relative weights of variables describing rainfall climatology as well as environmental variables explaining woody cover changes (1992–2012) from the LMG method based on a least squares regression model. Rainfall climatology: changes in annual rainfall, heavy rainfall frequency and rainy days; T&S: changes in mean annual air temperature and solar radiation (represented by cloud cover); management: changes in human population density and fire (burned area fraction); static includes variables: mean annual rainfall, coefficient of inter- annual rainfall variability, soil organic carbon, elevation and sand fraction. b Changes in rainfall climatology and woody cover were averaged along the mean annual rainfall gradient (50 mm yr−1 steps); all variables were rescaled to 0–100. c Sensitivity of changes in woody cover towards changes in rainfall variables (expressed as slope from a linear regression; βwoody-Rainfall) for each rainfall region. Error bars indicate the 5–95 % confidence intervals. All rainfall climatology variables were normalized of reliable global data on grazing intensities impeded this factor in woody cover. We found a similar pattern between changes in to be included in the model. woody cover and changes in rainfall variables along the entire We further examined the relationship of changes in different rainfall gradient, with the strongest changes happening in arid rainfall variables (representing rainfall climatology) and changes and semi-arid regions (Fig. 2b). The sensitivity of changes in NATURE COMMUNICATIONS | (2019)10:671 | https://doi.org/10.1038/s41467-019-08602-6 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08602-6 woody cover to changes in rainfall climatology was assessed by all continents (Fig. 3a) with large clusters of reduced sensitivity in the slope of a linear regression between changes in woody cover the Sudano–Sahelian zone, in southern Africa and northern and changes in rainfall variables (βwoody-Rainfall) (annual rainfall, Australia. Scattered areas showing an increased sensitivity of heavy rainfall frequency and rainy days) with results varying LAImax to annual rainfall were also found across all continents across rainfall regimes (Fig. 2c). For arid savannas, woody cover (for example woodlands in southern Africa, humid zones of changes were most sensitive to changes in annual rainfall and western Africa, South America and western Australia). As rainy days. In semi-arid savannas, the sensitivity of changes in increased heavy rainfall frequency has been documented to woody cover to annual rainfall, heavy rainfall frequency and rainy favour woody vegetation2, supported by a biophysical explanation days was of almost equal size, while in the humid region, woody of the phenomenon5, we analysed the impact of heavy rainfall cover was most sensitive to number of rainy days. The sensitivity frequency on changes in woody cover (changes in all variables of each explanatory variable to changes in woody cover is shown describing rainfall climatology are shown in the supplementary in Supplementary Table 3. These findings were further supported Fig. 6). The spatial correspondence of changes in woody cover by a partial regression analysis that produced comparable results and heavy rainfall frequency is shown in Fig. 3b, illustrating that (Supplementary Table 4). 24% of global savanna regions experienced an increase in both woody cover and heavy rainfall frequency (results for annual rainfall, rainy days and rainfall intensity are shown in the Sup- Spatial distribution of ecosystem structure changes. The plementary Fig. 7). Hot spot areas of change matched the patterns reduced sensitivity of LAImax to annual rainfall was evident over of Fig. 3a and coincided with areas where woody encroachment a 40°N 20°N 0° 20°S Trend in βLAI-rainfall 40°S –0.06 –0.03 0 0.03 0.06 b Not significant 40°N 20°N Trend in woody cover 60 Percent of counts (%) 0° + 40 – 20°S – + 20 Trend in heavy rainfall frequency 0 40°S c 40°N 20°N Trend in woody cover 60 Percent of counts (%) 0° + 40 – 20°S – + 20 Trend in LAImax 0 40°S 100°W 50°W 0° 50°E 100°E 150°E Fig. 3 Spatial distribution of ecosystem structure changes in tropical savannas. a Significant (p < 0.05) linear trends in βLAI-Rainfall (1982–2015). b Spatial agreement of trends in woody cover (1992–2012) and heavy rainfall frequency (1982–2015). Purple areas show increased trends in both heavy rainfall frequency and woody cover, light red areas show an increase in heavy rainfall frequency, but a decrease in woody cover, blue areas show an increase in woody cover but a decrease in heavy rainfall frequency, and grey areas show a decrease in both woody cover and heavy rainfall frequency. The bar plot shows the percent of count for each group with the same colour scheme as the map (total number of pixels = 17,717). c Spatial agreement of trends in woody cover (1992–2012) and LAImax (1982–2015). Purple areas show increased trends in both woody cover and LAImax, light red areas show a decrease in woody cover and increase in LAImax, blue areas an increase in woody cover and decrease in LAImax and grey areas show a decrease in both variables 4 NATURE COMMUNICATIONS | (2019)10:671 | https://doi.org/10.1038/s41467-019-08602-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08602-6 ARTICLE a 0.6 b 0.06 c 150 Density Density Woody Trend in woody cover (% yr–1) 120 LAImax Trend in βLAI-rainfall (z-score) 20 0.4 90 Area 15 0.03 60 Area (× 1000 km2) 10 Trend (z-score) 5 30 0.10 100 0.2 0.00 0.0 0.05 50 –0.03 –0.2 –0.4 –0.06 0.00 0 –0.3 –0.2 –0.1 0.0 0.1 0.2 –0.3 –0.2 –0.1 0.0 0.1 0.2 –0.3 –0.2 –0.1 0.0 0.1 0.2 Trend in heavy rainfall frequency Trend in heavy rainfall frequency Trend in heavy rainfall frequency (events yr–1) (events yr–1) (events yr–1) Fig. 4 Relationship between vegetation and heavy rainfall frequency in tropical savannas. a Spatial agreement of areas of changed woody cover and areas of changed heavy rainfall frequency. b The same for areas of changes in βLAI-Rainfall and areas of changed heavy rainfall frequency. c Changes in normalized woody cover and LAImax for 1992–2012 along a gradient of changes in heavy rainfall frequency (0.05 events yr−1 steps) has been documented by field studies and satellite data5,7,18,19. higher as compared with increases in LAImax, and this was, at Areas with an increasing woody cover but decreasing heavy global scale, largely controlled by rainfall climatology. Increases in rainfall frequency were located in African woodlands north and woody cover were mainly found in areas with a decreasing sen- south of the equatorial rainforest and may be related to non- sitivity between LAImax and annual rainfall, suggesting that symmetric responses to wet years or changes in human woody vegetation benefits more from the increased annual management26. rainfall as compared with LAImax. A possible explanation is found Most areas with an increased woody cover also experienced an in the altered rainfall climatology (distribution of annual rainfall). increase in LAImax (54%) (Fig. 3c). However, areas of increased Here, we found a strong relationship between increases in woody LAImax and a decreased woody cover were found in the Brazilian cover and increases in the frequency of heavy rainfall events, cerrado and Nigeria, which spatially agree with areas of large- which has also been demonstrated at experimental sites2,19. Since scale forest to farmland conversions27,28. For areas of both herbaceous vegetation is known to be vulnerable to an uneven increased woody cover and LAImax, particularly in Southern seasonal distribution of rainfall34, altered rainfall climatology is Africa, Northern Australia and Central African Republic large likely to affect the coexistence of woody and herbaceous vegeta- clusters of pixels showing higher trends in woody growth as tion, e.g., favouring the expansion of the ranges of woody compared with LAImax is observed (Supplementary Fig. 3) plants2,13. This effect would imply that large-scale woody consistent with the substantial woody encroachment reported in encroachment reported in savanna ecosystems2 is at least partly refs. 5,6. driven by changes in the rainfall regime. This further implies a The majority of areas where heavy rainfall frequency increased pronounced change in the functioning of global savanna eco- also showed an increased woody cover (Fig. 4a), and areas of systems within the last 30 years, which would have a profound increased heavy rainfall frequency also aligned with areas of a impact on the biodiversity, carbon-storage capacity and provi- decreased sensitivity of LAImax towards annual rainfall (Fig. 4b). sioning of these ecosystems. The magnitude of increases in heavy rainfall frequency further Moreover, the results may contribute explaining the phenom- appeared to be associated with the magnitude of increased woody enon of the greening of global drylands11,29. A strong influence of cover, which was not the case for changes in LAImax (Fig. 4c). El Niño/La Nina events on vegetation has been reported in several regions of the world, notably via the impact on rainfall13,16. From Discussion a wider climatic perspective, the increasing frequency of La Niña/ Our study shows global increases in LAImax and annual rainfall, El Niño events with global warming suggested by some recent in particular in arid and semi-arid regions. This increase implies studies35,36 would lead to an increase in the frequency of heavy that the total production of green vegetation is increasing in rainfall episodes and droughts37, which could further favour recent decades, being in line with studies by Zhu29 and Fen- woody vegetation. Experimental studies have demonstrated that sholt11. However, at the same time, the sensitivity of LAImax to increased CO2 positively affected vegetation productivity, parti- annual rainfall declined over the period analysed, supporting that cularly the growth of woody plants9,38. However, the contribution a shift in the relationship between vegetation productivity and of CO2 fertilization on the global greening and large-scale woody rainfall is taking place30,31. As LAImax comprises both herbaceous encroachment is difficult to quantify and global evidence based and woody foliage production, the decreased sensitivity between on observational data has not yet been established5,9,10,39,40. Our LAImax and annual rainfall suggests a change in the ecosystems results do not exclude the possibility of a “CO2 effect”, but sug- structure, and/or changes in the rainfall climatology. It has been gest, from a global perspective, woody cover increases in tropical shown that increases in rainfall benefit herbaceous plants at the savannas to be influenced by altered rainfall climatology, expense of woody plant growth18, whereas woody plant pro- including a higher annual rainfall, higher frequency of extreme liferation often occurs in periods following lower rainfall, when rainfall events and number of rainy days induced by climate the competition with herbaceous plants is low32. However, the change. decreased sensitivity between LAImax and annual rainfall does not point towards a strongly increased herbaceous mass, since her- Methods baceous vegetation (and also woody foliage production) in Climatic data. We used monthly climatic data (rainfall, temperature and cloud cover) with a spatial resolution of 0.5° (1982–2015) from the University of East savannas reacts mainly to the rainfall conditions of a given year33. Anglia’s Climate Research Unit CRU TS 441 and CHIRPS satellite-based daily At the same time, we find that increases in woody vegetation (as rainfall estimates at 0.05° spatial resolution42, both resampled to 0.25° (bilinear and measured by the water content of the entire woody plant) were nearest neighbour, respectively). CRU rainfall was used to calculate the relationship NATURE COMMUNICATIONS | (2019)10:671 | https://doi.org/10.1038/s41467-019-08602-6 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08602-6 between LAI and annual rainfall. Heavy rainfall frequency were extracted from including arid (< 300 mm yr−1), semi-arid (300–700 mm yr−1) and humid regions CHIRPS using the number of rainy days (≥ 1 mm day−1) for 1982—2015 above the (> 700 mm yr−1) divided according to mean annual rainfall during 1992–2012. 90th percentile of observed daily rainfall, the number of rainy days and rainfall Each variable xi was normalized as Xi = (xi –x)/σ(x) (where x and σ(x) are the mean intensity were calculated for each year. and standard deviation of xi) before modelling the data and all analysis were conducted in R. Vegetation data. We used the latest version of GIMMS3g LAI v1 that was gen- erated by training the GIMMS3g v1 NDVI (Normalized Difference Vegetation Reporting summary. Further information on experimental design is available in Index) with MODIS LAI using an artificial neutral network43. This dataset was the Nature Research Reporting Summary linked to this article. provided with a bi-monthly temporal resolution in a 1/12° spatial resolution from 1982 to 2015. LAImax is an indicator of the total green vegetation production22, Data availability including the foliage mass from both herbaceous and woody vegetation. The growing GIMMS3g LAI data are available at http://cliveg.bu.edu/modismisr/lai3g-fpar3g.html. season maximum leaf area index (LAImax) was used to study changes in savanna The VOD data is provided by Y. Liu, Nanjing University of Information Science & photosynthetic primary production. The LAImax was resampled to a 0.25° spatial Technology, China. CRU TS 4 datasets (rainfall, temperature and cloud cover) are resolution using nearest neighbour interpolation to match the VOD pixel size. available from the University of East Anglia (http://www.cru.uea.ac.uk/). CHIRPS Contrary to LAI, vegetation optical depth (VOD) is sensitive to the water datasets can be downloaded from ftp://ftp.chg.ucsb.edu/pub/org/chg/ products /CHIRPS- content of both the green and non-green woody parts of vegetation44 and has been 2.0. The GFED4s datasets are available at http://www.globalfiredata.org/data.html. The shown to better reflect the woody vegetation signal20. The data were retrieved from Gridded Population of the World data can be downloaded from http://sedac.ciesin. satellite passive-microwave observations and several merged sensors at a spatial columbia.edu/data/collection/gpw-v4. The soil maps are available at http://worldgrids. resolution of 0.25˚ from 1992 to 201245. We followed the approach of Brandt org/doku.php/wiki:layers. The SRTM 90 is available at http://www.cgiar-csi.org/data/ et al.21 who found a strong linear correlation between annual minimum VOD srtm-90m-digital-elevation-database-v4-1. values and woody cover over Africa (r² = 0.81, slope = 64.46, regression forced through 0 to avoid negative values) and transformed the annual minimum VOD values to the unit woody cover (%) by multiplication with 64.46. By using the Received: 1 October 2018 Accepted: 15 January 2019 seasonal minimum, we also minimized the influence of annual herbaceous vegetation and avoided saturation effects over dense forests21. Environmental data. Additional datasets, including dynamic environmental/ demography variables (changes in air temperature, solar radiation, human popu- References 1. Sankaran, M. et al. Determinants of woody cover in African savannas. Nature lation density and burned area fraction) and static environmental variables (mean annual rainfall, coefficient of inter-annual rainfall variability, soil organic carbon, 438, 846–849 (2005). elevation and sand fraction), were used to study the relationship between changed 2. Kulmatiski, A. & Beard, K. H. Woody plant encroachment facilitated by herbaceous/woody vegetation and rainfall climatology changes. Trends in air increased precipitation intensity. Nat. Clim. Change 3, 833–837 (2013). temperature and solar radiation (represented by cloud cover) for 1992–2012 were 3. Staver, A. C., Archibald, S. & Levin, S. The global extent and determinants calculated using the CRU TS4 datasets41. The differences in human population of savanna and forest as alternative biome states. Science 334, 230–232 density between 2015 and 2000 were based on the Gridded Population of the (2011). World (GPWv4)46 dataset. Mean annual burned fraction for 1997–2012 from the 4. Sankaran, M., Ratnam, J. & Hanan, N. Woody cover in African savannas: Global Fire Emissions Database (GFED4s)47 was used to analyse the impact of the role of resources, fire and herbivory. Glob. Ecol. Biogeogr. 17, 236–245 changes in fire activity. Soil organic carbon (SOC) and sand fraction were extracted (2008). from the Harmonized World Soil Database (HWSD v1.2)48. The elevation 5. Stevens, N., Lehmann, C. E. R., Murphy, B. P. & Durigan, G. Savanna woody extracted from the Shuttle Radar Topography Mission data (SRTM 90 m)49. All encroachment is widespread across three continents. Glob. Change Biol. 23, datasets collected were resampled to a 0.25° spatial resolution using nearest 235–244 (2017). neighbour interpolation (except GPWv4 and GFED4s that were provided in 0.25° 6. Venter, Z. S., Cramer, M. D. & Hawkins, H.-J. Drivers of woody plant spatial resolution). encroachment over Africa. Nat. Commun. 9, 2272 (2018). 7. Andela, N. et al. A human-driven decline in global burned area. Science 356, 1356–1362 (2017). Additional data. We defined the areal extent of tropical savannas based on the 8. Niinemets, Ü., Flexas, J. & Peñuelas, J. Evergreens favored by higher ecoregion classification by Olson et al.50. To minimize the human footprint on the responsiveness to increased CO2. Trends Ecol. Evol. 26, 136–142 (2011). analyses conducted, global land cover (GLC-SHARE) data51 were used to mask out 9. Buitenwerf, R., Bond, W. J., Stevens, N. & Trollope, W. S. W. Increased tree cropland areas with a cropped percentage above 70%. ESA CCI land cover (https:// www.esa-landcover-cci.org/) from 2010 was used to mask irrigated areas. Addi- densities in South African savannas:>50 years of data suggests CO2 as a driver. tionally, ESA CCI land cover was used to limit our analysis to savanna vegetation Glob. Change Biol. 18, 675–684 (2012). by masking out water areas. 10. Devine, A. P., McDonald, R. A., Quaife, T. & Maclean, I. M. D. Determinants of woody encroachment and cover in African savannas. Oecologia 183, 939–951 (2017). Analysis. We performed a linear regression between herbaceous vegetation 11. Fensholt, R. et al. Greenness in semi-arid areas across the globe 1981–2007— (LAImax) and annual rainfall over a 15-year moving window to derive the beta an earth observing satellite based analysis of trends and drivers. Remote Sens. coefficient (βLAI-Rainfall), which is frequently used as a measure of the sensitivity of Environ. 121, 144–158 (2012). vegetation to variation in rainfall. Variables were normalized at each grid point by 12. Good, S. P. & Caylor, K. K. Climatological determinants of woody cover in the inter-annual standard deviation in each study period prior to linear regression Africa. Proc. Natl. Acad. Sci. 108, 4902–4907 (2011). analysis. A non-parametric Theil–Sen trend analysis was then applied to detect 13. Holmgren, M., Hirota, M., van Nes, E. H. & Scheffer, M. Effects of interannual trends in βLAI-Rainfall and the Mann–Kendall (MK) test was applied to evaluate βLAI- climate variability on tropical tree cover. Nat. Clim. Change 3, 755–758 (2013). Rainfall trends at the 95% significance level (p < 0.05) taking temporal autocorrela- 14. Donat, M. G. et al. Updated analyses of temperature and precipitation extreme tion into account. Trends in woody cover and rainfall climatology variables were indices since the beginning of the twentieth century: The HadEX2 dataset. J. also calculated using this method. Here, the MK significance test based on linear Geophys. Res.: Atmospheres 118, 2098–2118 (2013). trends was avoided as a linear response of changes in woody cover to changes in 15. Knapp, A. K. et al. Consequences of more extreme precipitation regimes for rainfall climatology variables was not to be expected. Moreover, due to the high terrestrial ecosystems. Bioscience 58, 811–821 (2008). inter-annual variability of climatology, significant trends could not be expected52. 16. Holmgren, M. et al. Extreme climatic events shape arid and semiarid Welch’s unequal variances t test was used to test for significant difference between ecosystems. Front Ecol. Environ. 4, 87–95 (2006). temporal changes in LAImax and woody cover. 17. Heisler-White, J. L., Blair, J. M., Kelly, E. F., Harmoney, K. & Knapp, A. K. The relative weights of the controlling factors on changes in woody cover were estimated using the LMG method (a bootstrapping technique)24, which determines Contingent productivity responses to more extreme rainfall regimes across a the explaining power of each variables in the model as a share of 100%. A least grassland biome. Glob. Change Biol. 15, 2894–2904 (2009). squares regression model and a partial regression model were used to relate 18. February, E. C., Higgins, S. I., Bond, W. J. & Swemmer, L. Influence of changes in woody cover to changed rainfall climatology variables and competition and rainfall manipulation on the growth responses of savanna environmental variables. Changed woody cover was used as the response variable. trees and grasses. Ecology 94, 1155–1164 (2013). The explanatory variables included: (1) change in annual rainfall; (2) change in 19. Gherardi, L. A. & Sala, O. E. Enhanced precipitation variability decreases heavy rainfall frequency; (3) change in rainy days; (4) change in annual air grass-and increases shrub-productivity. Proc. Natl. Acad. Sci. 112, temperature; (5) change in solar radiation; (6) change in human population 12735–12740 (2015). density; (7) change in burned area fraction; (8) mean annual rainfall; (9) coefficient 20. Tian, F., Brandt, M., Liu, Y. Y., Rasmussen, K. & Fensholt, R. Mapping gains of inter-annual rainfall variability; (10) soil organic carbon; (11) elevation; (12) and losses in woody vegetation across global tropical drylands. Glob. Change sand fraction. The model was run for tropical savannas (All) and its sub-regions Biol. 23, 1748–1760 (2017). 6 NATURE COMMUNICATIONS | (2019)10:671 | https://doi.org/10.1038/s41467-019-08602-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08602-6 ARTICLE 21. Brandt, M. et al. Human population growth offsets climate-driven increase in 47. Giglio, L., Randerson, J. T. & Van Der Werf, G. R. Analysis of daily, woody vegetation in sub-Saharan Africa. Nat. Ecol. & Evol. 1, 0081 (2017). monthly, and annual burned area using the fourth-generation global fire 22. Myneni, R. B. et al. Global products of vegetation leaf area and fraction emissions database (GFED4). J. Geophys. Res.: Biogeosciences 118, 317–328 absorbed PAR from year one of MODIS data. Remote Sens. Environ. 83, (2013). 214–231 (2002). 48. Wieder, W., Boehnert, J., Bonan, G. & Langseth, M. Regridded harmonized 23. Guan, K. et al. Seasonal coupling of canopy structure and function in African world soil database v1.2. Oak Ridge, Tennessee, USA. TN: ORNL Distributed tropical forests and its environmental controls. Ecosphere 4, 1–21 (2013). Active Archive Center. Retrieved from http://daac.ornl.gov (2014). 24. Grömping, U. Estimators of relative importance in linear regression based on 49. Jarvis, A., Reuter, H., Nelson, A. & Guevara, E. Hole-filled SRTM for the globe variance decomposition. Am. Stat. 61, 139–147 (2007). Version 4. Web publication/site, Retrieved from http://srtm.csi.cgiar.org 25. Skowno, A. L. et al. Woodland expansion in South African grassy biomes (2008). based on satellite observations (1990–2013): general patterns and potential 50. Olson, D. M. et al. Terrestrial Ecoregions of the World: A New Map of Life on drivers. Glob. Change Biol. 23, 2358–2369 (2017). Earth. Bioscience 51, 933–938 (2001). 26. Brandt, M. et al. Reduction of tree cover in West African woodlands and 51. Latham, J., Cumani, R., Rosati, I. & Bloise, M. Global Land Cover SHARE promotion in semi-arid farmlands. Nat. Geosci. 11, 328–333 (2018). (GLC-SHARE) Database Beta-Release Version 1.0-2014. Technical Report: 27. Audu, E. B. Fuel wood consumption and desertification in Nigeria. Int. J. Sci. Food and Agric. Organ, Rome, Italy (2014). Technol. 3, 1–5 (2013). 52. Fischer, E. M., Beyerle, U. & Knutti, R. Robust spatially aggregated projections 28. Lapola, D. M. et al. Pervasive transition of the Brazilian land-use system. Nat. of climate extremes. Nat. Clim. Change 3, 1033–1038 (2013). Clim. Change 4, 27–35 (2014). 29. Zhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Change 6, 791–795 (2016). Acknowledgements 30. Knapp, A. K., Ciais, P. & Smith, M. D. Reconciling inconsistencies in This study is jointly supported by the China Scholarship Council (CSC, 201506190076) precipitation–productivity relationships: implications for climate change. New and the Danish Council for Independent Research (DFF) project: Greening of drylands Phytol. 214, 41–47 (2017). (DFF-6111-00258): Towards understanding ecosystem functioning changes, drivers and 31. Estiarte, M. et al. Few multiyear precipitation–reduction experiments find a impacts on livelihoods. J.P. acknowledges the financial support from the European shift in the productivity–precipitation relationship. Glob. Change Biol. 22, Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P. F.G. 2570–2581 (2016). acknowledges funding from the NERC/DFID Future Climate For Africa programme 32. Riginos, C. Grass competition suppresses savanna tree growth across multiple under the AMMA-2050 project (Grant NE/ M019950/1). M.B. was supported by the demographic stages. Ecology 90, 335–340 (2009). AXA postdoctoral research grant. Also, we would like to thank Y.Y. Liu for providing the 33. Brandt, M. et al. Ground-and satellite-based evidence of the biophysical VOD data. mechanisms behind the greening Sahel. Glob. Change Biol. 21, 1610–1620 (2015). Author contributions 34. Zhang, W., Brandt, M., Tong, X., Tian, Q. & Fensholt, R. Impacts of the W.Z., M.B., R.F. and J.P. designed the study. W.Z. conducted the analyses with support seasonal distribution of rainfall on vegetation productivity across the Sahel. by M.B., R.F., J.P., F.G., X.T. and F.T. W.Z., M.B. and R.F. drafted the paper with Biogeosciences 15, 319–330 (2018). contributions by J.P. All authors discussed the results and commented on the paper. 35. Cai, W. et al. Increased frequency of extreme La Niña events under greenhouse warming. Nat. Clim. Change 5, 132–137 (2015). 36. Cai, W. et al. Increasing frequency of extreme El Niño events due to Additional information greenhouse warming. Nat. Clim. Change 5, 1–6 (2014). Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- 37. Mason, S. J. & Goddard, L. Probabilistic precipitation anomalies associated 019-08602-6. with ENSO. Bull. Am. Meteorol. Soc. 82, 619–638 (2001). 38. Yang, Y., Donohue, R. J., Mcvicar, T. R., Roderick, M. L. & Beck, H. E. Long- Competing interests: The authors declare no competing interests. term CO2 fertilization increases vegetation productivity and has little effect on hydrological partitioning in tropical rainforests. J. Geophys. Res. Biogeosci. 121, Reprints and permission information is available online at http://npg.nature.com/ 2125–2140 (2016). reprintsandpermissions/ 39. Bond, W. J. & Midgley, G. F. A proposed CO2-controlled mechanism of woody plant invasion in grasslands and savannas. Glob. Change Biol. 6, Journal peer review information: Nature Communications thanks Kelly Caylor and the 865–869 (2000). anonymous reviewer for their contribution to the peer review of this work. Peer reviewer 40. Higgins, S. I. & Scheiter, S. Atmospheric CO2 forces abrupt vegetation shifts reports are available. locally, but not globally. Nature 488, 209–212 (2012). 41. Mitchell, T. D. & Jones, P. D. An improved method of constructing a database Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in of monthly climate observations and associated high-resolution grids. Int. J. published maps and institutional affiliations. Climatol. 25, 693–712 (2005). 42. Funk, C. et al. The climate hazards infrared precipitation with stations—a new environmental record for monitoring extremes. Sci. Data 2, 150066 Open Access This article is licensed under a Creative Commons (2015). Attribution 4.0 International License, which permits use, sharing, 43. Zhu, Z. et al. Global data sets of vegetation leaf area index (LAI)3g and fraction adaptation, distribution and reproduction in any medium or format, as long as you give of photosynthetically active radiation (FPAR)3g derived from global inventory appropriate credit to the original author(s) and the source, provide a link to the Creative modeling and mapping studies (GIMMS) normalized difference vegetation index Commons license, and indicate if changes were made. The images or other third party (NDVI3G) for the period 1981 to 2. Remote Sens. 5, 927–948 (2013). material in this article are included in the article’s Creative Commons license, unless 44. Tian, F. et al. Remote sensing of vegetation dynamics in drylands: Evaluating indicated otherwise in a credit line to the material. If material is not included in the vegetation optical depth (VOD) using AVHRR NDVI and in situ green article’s Creative Commons license and your intended use is not permitted by statutory biomass data over West African Sahel. Remote Sens. Environ. 177, 265–276 regulation or exceeds the permitted use, you will need to obtain permission directly from (2016). 45. Liu, Y. Y. et al. Recent reversal in loss of global terrestrial biomass. Nat. Clim. the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Change 5, 1–5 (2015). licenses/by/4.0/. 46. Center for International Earth Science Information Network - CIESIN - Columbia University. Gridded Population of the World, Version 4 (GPWv4): Population Density, Revision 10. (2017). © The Author(s) 2019 NATURE COMMUNICATIONS | (2019)10:671 | https://doi.org/10.1038/s41467-019-08602-6 | www.nature.com/naturecommunications 7