Introduction

Land degradation is a systemic global problem1,2,3,4 but the scale of the problem is disputed, with global estimates of degraded areas ranging from <10 to >60 million km2 5. Changes in vegetation in drylands are predominantly caused by two factors: (i) anthropogenic climate change (ACC), which includes both changes in water availability driven by trends in precipitation and increases in temperature6,7, as well as increased water use efficiency (carbon gain per unit of water lost) in response to rising atmospheric CO28; and (ii) land use (LU) practices, including grazing, cropping and deforestation2,9. Unsustainable LU is considered the primary negative driver of dryland degradation9,10,11. The impact of climate change (CC) on drylands is also generally thought to be negative, with some studies suggesting that anthropogenic forcing has already increased arid areas12,13,14.

Despite evidence for LU-induced degradation and the studies that find increased aridification over drylands, satellite estimates of vegetation greenness (a proxy for net primary productivity (NPP)) show a significant global increase since 198010. The key drivers of this global increase in apparent vegetation productivity are the vegetation’s response to rising CO28,15, increases in rainfall and temperature16,17 and LU10. Model simulations which prescribe LU, attribute almost all of the trend in satellite-derived greening to CO2 fertilization15, while satellite-derived models that do not account for CO2, explicitly find either climate or LU as the dominate factor10,17. Neither approach explicitly accounts for rapid ecosystem change (break points) in their proportioning of the relative contributions of each driver. This can lead them to miss or underestimate rapid changes driven by processes like extreme fires, deforestation, reforestation, changes in agricultural policy, etc.18,19,20,21. Disentangling the roles of climate (temperature and precipitation), CO2 and LU thus remains a key challenge22 and has been identified as a key knowledge gap by the United Nations Convention to Combat Desertification2 (UNCCD), the Intergovernmental Panel on Climate Change23 (IPCC), and the Intergovernmental Science–Policy Platform on Biodiversity and Ecosystem Services3.

Here we quantified the scale of global desertification, which both the UNCCD and IPCC define as degradation in arid, semi-arid, and dry sub-humid areas23. These definitions further define degradation as the long-term reduction or loss of biological productivity among other things. Here we identify areas undergoing long-term reductions in vegetation in dryland areas, hence the desertification according to relevant international conventions, using the satellite-based GIMMSv3.1g Normalized Difference Vegetation Index (NDVI) data. We calculated the overall vegetation change using a non-parametric trend analysis applied to peak growing season NDVI (NDVImax). We then attributed this change to CO2, climate variability (CV), CC, and LU using a modified version of the Time Series Segmented Residual Trends (TSS-RESTREND) method19,24. This approach quantifies the effect of interannual CV as well as long-term changes in climate and CO2 fertilization in addition to ecosystem break points caused by LU (see “Methods”). To quantify uncertainties, we used a 12-member ensemble made up of statistical model runs performed using a combination of observation-based gridded datasets (four precipitation and three temperature datasets). We show that 6% of dryland areas have undergone desertification since 1982 with a further 20% of dryland areas being at high risk of future desertification as a result of unsustainable LU practices and ACC.

Results and discussion

The extent and drivers of dryland vegetation change

Globally, of the 44.5 million km2 of drylands, 6% of these areas experienced desertification (i.e., significant negative change in NDVImax), 41% showed significant greening (i.e., significant positive change), and 53% had no significant change between 1982 and 2015 (Fig. 1a). The mean (±1 SD) of the area-weighted dryland vegetation change, as represented by the change in NDVImax was 0.031 ± 0.053. We estimated the scale of desertification to be 2.70 million km2, which is significantly below a previous estimate of ~10.5 million km2 over the same region, but over a different time window (1982 and 2003)1. A large part of this discrepancy can be attributed to climatic differences in the end dates of the studies (2003 vs. 2015), with increased rainfall over regions including the Sahel and India25,26. This large difference between our estimate and this existing dryland degradation estimate highlights that time-series vegetation trend analysis is sensitive to the start and end conditions17,18. For this reason, understanding what is driving the observed vegetation change is more important than the current directions of vegetation change for projecting future changes and vulnerabilities. It should also be noted that although the amount of land we estimate to have experienced desertification has decreased by ~70% between these two estimates, the number of people impacted has only decreased by ~25% (250 million to 189 million)27.

Fig. 1: Primary drivers of vegetation changes between 1982 and 2015.
figure 1

a The observed dryland vegetation change from 1982 to 2015 measured using the change in p Normalized Difference Vegetation Index (ΔNDVImax). Non-dryland and hyper-arid regions are masked in dark gray and areas where the change is insignificant (αFDR = 0.10) or smaller than the error in the sensors (±0.001) are masked in white. b Map showing the largest absolute attributed driver (CO2, land use, climate change, and climate variability). Non-dryland regions are masked in dark gray and the magnitude of other factors is not considered. c The mean and magnitude (mean absolute value) of the different drivers of change in NDVImax between 1982 and 2015. The error bars show the SD of the area-weight grid cells. d Frequency distribution function for the change in NDVImax (1982–2015) attributed to land use, which shows that there are opposing changes at local scales which cancel out when averaged globally.

Full size image

Figure 1b shows that globally, CO2 fertilization was the largest absolute attributed driver of dryland vegetation change in 44.1% of areas, followed by LU (28.2%), CV (14.6%), and then CC (13.1%). However, when averaged globally (Fig. 1c), the per-pixel contribution of CO2 (0.021 ± 0.011) was much larger than the contribution from CV (0.006 ± 0.020), CC (−0.002 ± 0.023), or LU (0.005 ± 0.032). The relative contribution (67.8% CO2, −5.6% climate, 15.5% LU) fall within the range of global estimates calculated using a model based factor analysis (70.1 ± 29.4% CO2 fertilization, 8.1 ± 20.6% climate, 3.7 ± 14.7% LU)15, despite a difference in the study domains. It should be noted here that model based factor analysis did not quantify the role of CV15, which we find accounts for 19.4% of the observed global dryland greening between 1982 and 2015.

If only the global mean effect size is considered, climate and LU seem to have a very small impact compared with CO2 fertilization, which seemingly contradicts well-documented evidence of LU and climate impacts9,11,16,17,28, and the spatial patterns shown in Fig. 1b. For example, a recent satellite-based study, which did not consider the role of CO2, attributed 60% of observed global land changes to LU activities and the remaining 40% to other factors including climate10. For comparison, we repeated our ensemble analysis, removing the role of CO2 fertilization (Supplementary Fig. 1), and found that 60.4% of global dryland vegetation change would have been attributed to LU and 39.6% to CC and variability combined, which is consistent despite a difference in the study domains and the attribution methods (see Supplementary Text 1). This result underlines the need to explicitly consider the positive role of CO2 as a driving mechanism for change in dryland ecosystems8. When CO2 is included in the analysis, the mean effect of LU and climate are small, because there are roughly equal areas of positive and negative change that largely cancel out when averaged globally (Fig. 1c, d), a result which holds even if we assume a different level of vegetation response to elevated CO2 (Supplementary Fig. 2). This means it is also important to consider the magnitude of the different drivers (Fig. 1c).

The impact of anthropogenic climate change

Combining the change due to CO2 fertilization and CC provides a quantification of the role of ACC in recent desertification (Fig. 2a). Globally, we found that ACC had a positive (greening effect) over the study period (NDVImax: 0.019 ± 0.027). Although broadly positive, ACC also had a desertifying effect across 12.55% (5.43 million km2) of drylands areas. Hotspots where ACC had a desertifying effect include parts of the western United States, eastern Brazil, Iraq, Syria, Jordan, Kazakhstan, Uzbekistan, Mongolia, and Australia. Crucially, the negative effects of ACC are disproportionately felt by poorer nations with 85% of the 213.4 million people impacted living in developing or newly industrialized countries. However, a negative ACC forcing does not guarantee an area experienced desertification (Fig. 3a). Only 13.8% (0.75 million km2) of areas with a negative ACC forcing, experienced significant desertification (αFDR = 0.10) and in only 2.27% (0.015 million km2) of the areas experiencing desertification, did we find that climate was the sole negative driver.

Fig. 2: The contribution of anthropogenic climate change to vegetation change from 1982 to 2015.
figure 2

a The contribution of anthropogenic climate change (Climate Change + O2) component to the change in vegetation between 1982 and 2015 (ΔNDVImax). Non-dryland and hyper-arid regions are masked in dark gray and areas where the change is insignificant (αFDR  = 0.10) or smaller than the error in the sensors (±0.001) are masked in white. b The mean area-weighted pixel effect size of the drivers of observed vegetation broken down by observed vegetation change direction. Positive indicates greening and negative is desertification. Error bars show the SD. c Frequency distribution function for the change in NDVImax (1982–2015) attributed to by anthropogenic climate change.

Full size image
Fig. 3: Desertification risk and regional drivers.
figure 3

a Regions experiencing, or at risk of experiencing, desertification. Areas with a significant negative change in vegetation (αFDR = 0.10) are classified “Desertification,” areas where the anthropogenic climate change (ACC: CO2  + Climate Change) and land-use (LU) components are both negative but the change in vegetation is not significant are classified as “LU and ACC.” Areas where the anthropogenic climate change has had a negative effect are classified as “ACC” and areas where land use had a negative impact classified as “LU.” b map of the regional subdivisions used. c The mean per-pixel ΔNDVImax and its drivers broken down by regions shown in b. The error bars show the SD.

Full size image

Drivers of desertification

In the 2.70 million km2 of drylands that experienced desertification, a negative LU component was the primary driver in 79.9% and a contributing factor across 99.0% of areas (Fig. 2b). Even though the average impact of CC (NDVImax: −0.004 ± 0.030) and CV (−0.002 ± 0.024) are much smaller than LU (NDVImax: −0.040 ± 0.034), climate remains an important driver of desertification. Ecosystems that are experiencing reduced water availability or drought conditions are much more vulnerable to degradation from LU and vice versa, with the negative effects compounding7. For example, over parts of Central Asia we observed negative changes in both CC and LU (Fig. 4), which is consistent with the strong evidence of long-term degradation driven by unsustainable LU practices resulting in the well-documented Aral Sea disaster26. Similarly, the negative impacts of decreased rainfall over the semi-arid Caatinga forest of Brazil has amplified the effects of widespread deforestation and grazing intensification. The decrease in rainfall in South America results from both CC (mean precipitation anomaly 1982–2015 ≈ −0.2) as well as a negative phase of CV (mean precipitation anomaly ≤ 0 and mean temperature anomaly ≥ 0) from 2009 to 2014 (Supplementary Figs. 3 and 4). For further discussion of these regions and comparison with regional studies, see Supplementary Text 2.1. In addition to the drylands experiencing desertification, there are additional 12.0 million km2 and 507 million people, living in areas where the desertifying effect of LU has been offset by a positive ACC signal (Fig. 3). These regions, along with the 7.2% of areas with a negative CC, but no significant vegetation change, are at the highest risk of future desertification.

Fig. 4: The drivers of global vegetation change.
figure 4

The changes in NDVImax between 1982 and 2015 (ΔNDVImax) attributed to (a) CO2 fertilization, (b) climate variability, (c) climate change, and (d) land use. Non-dryland and hyper-arid regions are masked in dark gray. Areas where the change did not meet the multi-run ensemble significance criteria detailed in the “ Methods,” or are smaller than the error in the sensors (±0.001) are masked in white.

Full size image

Drivers of dryland greening

We also observed widespread global greening, with 18.0 million km2 of drylands having a significant positive vegetation change (Fig. 1a). CO2 was the largest driver of this change (Figs. 1b and 2b), in line with previous findings of a dominant CO2 fertilization effect on global vegetation greening8,15. Where our results differ from previous findings is in highlighting the importance of LU and CV. Unlike model based approaches which prescribe LU15 and have large variability in the simulated response of vegetation to climate29, our approach empirically determines the impact of climate and LU on a per-pixel basis using observations. In the regions that experienced greening, we find that CO2 was the largest attributed driver in ~40% of areas followed closely by LU (~38%), CV (~13%), and CC (~8%). The importance of CV and LU is especially apparent when considering regional drivers, with one or both playing a large role in the observed greening in the Sahel, India, China and Australia (Figs. 3c and 4, and Supplementary Discussion 2.2).

We used an ensemble approach to minimize the uncertainty caused by observational datasets and structural change detection to account for the ecosystem break points driven by processes like deforestation observed in ~20% of areas (Supplementary Fig. 5). Our estimate of the total attributable change (CO2 + LU + CC + CV) varied from the observed vegetation change by only 3% and when mapped spatially, the observed vegetation change and the total attributable change show very consistent patterns of greening and browning (Supplementary Fig. 6). Furthermore, our greening attribution results are consistent with regional studies done in the Sahel, India, China, and Australia (Supplementary Discussion 2.2).

In summary, understanding the causes of dryland degradation is an important and necessary step in targeting mitigating action that can reduce the impact of CC and prevent widespread desertification. Our change detection and attribution approach highlights the importance of accounting for the role of CO2 and accurately quantifying the impact of LU when considering potential drivers of change in dryland ecosystems. Our results show that, despite widespread vegetation greening, 6% of areas that have undergone desertification mostly over western Asia and South America. This desertification directly effects 190 million people. In addition, we showed that unsustainable LU practices or ACC has placed 20% of drylands at high risk of desertification. This impacts 580 million people with the risk experienced disproportionately by low socioeconomic countries. Overall, our results highlight the importance of understanding what is driving the vegetation change for projecting impacts, because without this understanding there is a high risk that mitigation strategies will fail to prevent desertification.

Methods

Quantifying desertification

There is no universally agreed upon definition of desertification5,23,30. Here we use the UNCCD definition of land degradation, which is a reduction or loss of the biological or economic productivity resulting from various factors, including climatic variations and human activities, with desertification being any land degradation in dryland ecosystems2. Drylands are defined by the UNCCD to be areas with an Aridity index < 0.05 or >0.652. Historically, this has been measured using a linear trend applied to a satellite-derived vegetation proxy1. In this study, we used the growing season maximum NNDVI (NDVImax) as a proxy of vegetation growth. For most regions, peak growing season NDVI was determined using the maximum value in a calendar year. However, pixels where peak the occurred in December, the January and February NDVI values of the subsequent year are considered part of the previous year’s growing season.

NDVImax has been found to have a highly significant correlation with NPP in a large range of different dryland ecosystem31,32,33. The Desertification chapter of the 2019 IPCC report on CC and LU both, the UNCCD definition of desertification, and NDVImax as a proxy of vegetation growth23. It should be noted that the UNCCD definition of desertification and the trend in vegetation data used to measure it will not identify processes such as shrub encroachment, over intensification of agriculture, or the invasions by non-native species, which have been linked to degradation but can cause increases in proxies such as NDVI23,30,34.

To produce a comparable estimate of desertification, we used a per-pixel non-parametric trend method (Theil–Sen slope estimator and Spearman’s ρ significance test), applied to the satellite-derived GIMMSv3.1g NDVImax dataset35. The global dryland vegetation change (Obs) was calculated as the difference between the expected values (E) at the start (1982) and the end (2015) of the time series \(({\mathrm{{Obs}}} = E_{2015} - E_{1982})\). It should be noted that this is identical to multiplying the annual trend by the length of the time series, in this case 34 years. We report all variables using the difference between expected values at the start and end of the time series (ΔNDVImax) rather than an as an annual trend to account for the ecosystem break points that are detected in some locations in the middle of the time series.

This study used version 3.1 of the 1/12° Global Inventory for Mapping and Modeling Studies (GIMMS)35 NDVI dataset, which spans 1982–2015. Although the shorter temporal but higher spatial resolution datasets from newer sensors such Moderate Resolution Imaging Spectroradiometer (MODIS) do offer advantages, the shorter temporal record poses a serious issue in dryland ecosystems. The natural variability in dryland ecosystems is greatly impacted by decadal climate modes the most significant of which is El Niño Southern Oscillation (ENSO)36,37. The intensity of ENSO events varies significantly and since 1980 there have been three extreme El Niño (1982, 1997, and 2015) events that significantly impacted dryland regions around the world38. Even with its almost 20-year record, MODIS has only captured one of these events (2015) compared with the three present in GIMMS record. It is for this reason that GIMMS remains the most widely used dataset for vegetation trend detection and attribution studies15.

Accounting for the CO2 fertilization effect

To attribute the change in NDVI to the CO2 effect on plant productivity between 1982 and 2015, we used a theoretical relationship that links the increase in photosynthesis to increasing CO239 (Eq. 1).

$${\mathrm{{GPP}}}_{({\mathrm{{rel}}})} \approx \left[ {\frac{{\left( {c_a - {\Gamma}^ \ast } \right)\left( {c_{a0} + 2{\Gamma}^ \ast } \right)}}{{\left( {c_a + 2{\Gamma}^ \ast } \right)\left( {c_{a0} - {\Gamma}^ \ast } \right)}}} \right]$$
(1)

where GPP(rel) is the relative CO2 assimilation rate (%), ca is the atmospheric CO2 concentration (µmol mol−1), and Γ* is the CO2 compensation point in the absence of dark respiration (µmol mol−1). We set ca0 to the CO2 concentration in 198040 (~339 µmol mol−1) and Γ* = 40 (µmol mol−1).

Franks et al.39 argued that the longer term response of plants to increasing CO2 follows the ribulose 1,5-bisphosphate regeneration-limited rate (see also McMurtrie et al.41). Accordingly, this relationship implies a conservative response to CO2 (as plants may actually follow the Rubisco-limited rate when calculated on a intercellular CO2 concentration (Ci) basis during the period of 1982–201542) and ignores any indirect effects (i.e., increased water availability due to stomatal closure, which may extend the growth period in drylands, or interactions with seasonal rainfall43). This approach has previously been advocated as a plausible assumption to correctly estimate gross primary productivity (GPP) using satellite light-use efficiency models44. We then assume that there has been no change in ratio of GPP to autotrophic respiration (Ra) during this period (1982–2015) and, as a result, the relative change in GPP equates to the relative change in NPP based on Eq. 1. During the period 1982–2015, global air temperatures have risen, which may have led to an increase in Ra45. However, increasing temperature has also increased carbon uptake, and both GPP and Ra have been shown to acclimate to the prevailing temperatures46,47, meaning that it does not necessarily follow that the GPP : Ra ratio has changed.

We apply the nonlinear CO2 relationship (Eq. 1) to the raw NDVI data (NDVIobs) to produce a scaled NDVI estimate (NDVIadj) that excludes the CO2 fertilization effect using Eq. 2 to relate the relative change in NPP to a relative change in NDVI.

$$\frac{{\mathrm{{NPP}}}_{{\mathrm{{obs}}}}}{{\mathrm{{NPP}}}_{\mathrm{{base}}}} \approx \frac{{\mathrm{{NDVI}}}_{\mathrm{{obs}}}}{{\mathrm{{NDVI}}}_{\mathrm{{adj}}}}$$
(2)

where NPPobs is the NPP at the observed atmospheric CO2 concentration (ca), NPPbase is the NPP given the same climate conditions but an atmospheric CO2 concentration of ca0, NDVIobs is measured NDVI value, and NDVIadj NPP given the same climate conditions but an atmospheric CO2 concentration of ca0. Equation 2 was used to calculate a NDVIadj value for every in the full NDVIobs time series with the atmospheric CO2 concentrations taken from the IPCC historical forcing data40. This approach assumes that NPP and NDVI are linearly related: .. where b ≈ 0 and m varies spatially. The linear relationship between NPP and NDVI has been observed with both field estimates of NPP31,32,48 and estimates of NPP derived from remote-sensing platforms1,49,50. However, this assumption of linearity breaks down in densely vegetated regions where NDVI saturates and in biomes with very low above-ground biomass, where the spectral characteristics of the bare soil influences NDVI values51,52. As we have excluded hyper-arid and non-water-limited ecosystems, which is a standard practice for studies on desertification (for more information, see ref. 23), we expect this assumption of linearity to be robust for our analysis (areas with Aridity index < 0.05 or >0.65 are masked from analysis).

The change in vegetation (ΔNDVImax) attributed to the CO2 fertilization effect was calculated by first taking the difference between peak growing NDVI with and without CO2 fertilization effect (NDVIobs − NDVIadj). Similar to the calculation of the observed ΔNDVImax, the non-parametric Theil–Sen slope estimator and Spearman’s ρ test for significance was then applied to these values. A time-series plot of the mean global NDVIobs and NDVIadj is shown in Supplementary Fig. 7 to highlight temporal nature of this attribution.

Determining the impact of climate and land use

After the NDVI was scaled to remove the CO2 effect using the relationship from Franks et al.39, a statistical approach was used to attribute the change in NDVIadj to climate (both CV and CC) and LU. We used the recently developed TSS-RESTREND method, which allows vegetation changes due to LU to be separated from those driven by CC and variability19,24,53.

Dryland ecosystems have large natural interannual CV. To separate the impact of climate and LU, previous studies have generally fitted a statistical relationship between climate and vegetation, then used the trends in that relationship or its residuals to quantify LU impacts54,55. However, LU can have both a gradual impact on vegetation through processes such as grazing, which are captured by these methodologies54,56, and abrupt impact through processes such as deforestation, which cause these methods to break down18.

TSS-RESTREND differs from existing dryland trend attribution methods in that it is able to capture both the long-term trends and the step changes in NDVI that occur in regions where ecosystems have experienced significant structural changes19. To do this TSS-RESTREND incorporates a phenological change detection method57 to identify structural changes in the ecosystem, which manifest as break points in the NDVI time series. We used TSS-RESTREND v2.15, which has been updated to use both precipitation and temperature24, to calculate the Vegetation Climate Relationship (VCR) using the per-pixel optimal precipitation and temperature accumulation periods19,24. TSS-RESTREND was applied to the NDVIadj, with the LU driven component calculated using an ordinary least squared regression between the residuals of the VCR and time, accounting for any detected structural changes to the ecosystem. A similar approach was presented by the IPCC23, although that approach does not separate CC from CV and also assumes that all dryland plants follow a C3 photosynthetic pathway.

Separating climate change and climate variability

The TSS-RESTREND method separates the effects of LU from the combined effects of climate (variability and change) using VCR. To separate the effect of CC and CV, the observed climatology (the per-pixel accumulated precipitation and temperature data) was calculated for the period 1962 to 2015. A 20-year leading edge smoothing window was then applied to this observed climatology to remove the interannual CV. The long-term trend caused by CC was determined using the Theil–Sen slope estimator58 applied to the smoothed data and the results were used to detrend the observed climatology.

Using the per-pixel VCR, the total climate driven NDVI (NDVICL) and the NDVI due to CV (NDVIcv) were calculated using the observed climatology and detrended climatology, respectively. The difference between NDVICL and NDVICV is the change in NDVImax attributable to CC (NDVICC). The non-parametric Theil–Sen slope estimator and Spearman’s ρ test for significance was then applied to NDVICV and NDVICC to get the change attributable to CV and CC, respectively. The influence of Other Factors (OF), which could not be modeled, was calculated using \({\mathrm{{OF}}} = {\mathrm{{Obs}}} - ({\rm{CO}}_{2} + {\rm{LU}} + {\rm{CV}} + {\rm{CC}})\).

When discussing regional drivers of vegetation change in the main text and in Supplementary Discussion 2, we report the mean climate anomaly rather than the accumulated precipitation and temperature values to allow comparison of different accumulation and offset periods of different pixels. For the observed accumulated precipitation and temperature, the anomaly was calculated on a per-pixel basis using:

$$z_n = \frac{{x_n - \mu _{\mathrm{{obs.}}}}}{{\sigma _{\mathrm{{obs.}}}}}$$
(3)

where z = anomaly, n is the year, x = observed value, μ = mean of the per-pixel accumulated precipitation or temperature, σ = SD of the per-pixel accumulated precipitation or temperature. The temperature and precipitation anomaly attributed to CC was calculated using:

$$z_n = \frac{{\beta \times (n - n_0)}}{{\sigma _{\mathrm{{obs.}}}}}$$
(4)

where β is the per-pixel trend in accumulated precipitation or temperature, n0 is the first year of the analysis (1982). The temperature and precipitation anomaly attributed to CV was calculated using:

$$z_n = \frac{{x_{{\mathrm{{(adj)}}}n} - \mu _{\mathrm{{obs.}}}}}{{\sigma _{\mathrm{{obs.}}}}}$$
(5)

where x(adj) is the detrended precipitation or temperature. Regional breakdowns of the time series of observed, CC- and CV-driven precipitation and temperature anomaly are shown in Supplementary Figs. 3 and 4, respectively.

Calculating the impact of Anthropogenic Climate Change

In this study, we use the term ACC to refer to both changes in water availability driven by trends in precipitation and increases in temperature14,23, as well as increased water use efficiency (carbon gain per unit of water lost) in response to rising atmospheric CO28,15. The impact of ACC) is calculated using ACC = CO2 + CC. Although we acknowledge that this considering CO2 fertilization and CC together is not common in the literature, these two things are aspects of the same anthropogenic cause; hence, we have reason to discuss them together. It should be noted that although the methodology used in this study is able to capture the effects of other anthropogenic greenhouse gases like methane in the CC component, we cannot separately quantify any direct impact that changes in the amount of these gasses will have on vegetation.

Accounting for dataset uncertainties and different photosynthetic pathways (C3 vs. C4)

Burrell et al.53 showed that using an ensemble of TSS-RESTREND runs using different climate datasets improves the accuracy and minimizes the impact of errors and biases in the individual datasets. Climate data are relatively poorly sampled in dryland regions, which can amplify the documented discrepancies between different datasets59,60 We used two 12-member ensembles with matched runs made using TSS-RESTREND analysis performed using every combination of four precipitation and three temperature datasets (see Table 1 for details). The first ensemble assumes all plants are C3 and respond to elevated CO2, and the second ensemble assumes all plants are C4 with no eCO2 response (see Supplementary Figs. 1 and 2).

Table 1 Table of gridded datasets.
Full size table

A third 12-member ensemble, which accounts for the relative fraction of C3 and C4 plants, were calculated by taking the weighted mean of matched runs in ensemble one and two where the weights were the per-pixel fractions of C3 and C4 plants. Estimates of the relative fraction of C4 present in each pixel were derived from the matching 0.5° pixel in the North American Carbon Program (NACP) Global C3 and C4 SYNergetic land cover MAP (SYNMAP)61. The p-values for each ensemble member were combined using the same weights and the Stouffer’s Z-score method. All the results presented in the main paper are from this C4 adjusted ensemble.

We make the assumption that in dryland ecosystems, precipitation controls the amount of foliage cover. As a result, increases in water use efficiency with increasing CO2 will lead to increases in foliage cover in plants that use the C3 photosynthetic pathway8 while plants that use the C4 photosynthetic pathway are not expected to show this response62. Our assumption that C3 plant respond strongly to increased atmospheric CO2, and that C4 plants do not respond, which is consistent with theory and short-term CO2 enrichment experiments62. However, recent surprising results from a long-term CO2 manipulation experiment in a dryland ecosystem have shown much higher responses in C4 plants62. For this reason, the results of both the C3 and C4 12-member ensembles are included in the Supplementary Material for those interested parties (see Supplementary Figs. 1, 2, 8, and 9).

Determining statistical significance

For both Obs trend and CO2-driven change in NDVI, the Spearman’s rank correlation co-efficient test was used to measure statistical significance for each pixel63. In order to determine field significance and account for the multiple testing problem the Benjamini–Hochberg procedure was then applied to these p-values to control the false discovery rate (FDR) (αFDR = 0.10)64. As for the uncertainty in the approach we used to measure the CO2 fertilization effect, the C3 and C4 12-member ensembles included in the Supplementary Material provide an estimate of the upper and lower bounds of CO2 responses (see Supplementary Figs. 1 and 2).

For LU, CV, and CC there are 12 members in each ensemble. The p-values of these members were combined using the Fisher’s combined probability test and then the Benjamini–Hochberg procedure was applied to these p-values to determine statistical significance (αFDR = 0.10). In addition, we also applied the IPCC protocol for determining ensemble significance and agreement65. For a pixel to be significant under this protocol, >50% of ensemble members must find a significant change (αFDR = 0.10) and, of those significant runs, >80% must agree on the direction of change. If a pixel fails either the overall significance test or the IPCC protocol, the estimate of that component is masked for that pixel. Similar criteria have also been applied to ensemble breakpoint detection (>50% of runs must find a significant breakpoint, 80% of which must be in a three-year window) the results of which are included in supplementary material.

All the Climate datasets where interpolated from their native resolution to the 1/12° grid of the GIMMSv3.1 g datasets using the First-Order Conservative Remapping66 in CDO67. Additional datasets where used to aid the interpretation and discussion our results. All dataset where converted to the GIMMS grid to allow for per-pixel comparison. To separate CC from CV, a 20-year leading edge moving window was used where the value for a given year is the mean of the previous 20 years. For this reason, it was necessary to have climate data that goes back to 1960, which not all datasets do. Those climate datasets with insufficient temporal record were extended using TERRACLIMATE precipitation and temperature, and the Delta Bias Correction method68.