Human Land-Use Practices Lead to Global Long-Term Increases in Photosynthetic Capacity

1 Department of Biology, University of Maryland, College Park, MD 20742, USA; E-Mail: bfagan@umd.edu 2 Biodiversity and Climate Research Center, Senckenberg Gesellschaft für Naturforschung, Senckenberganlage 25, D-60325 Frankfurt, Germany; E-Mail: katrin.boehning-gaese@senckenberg.de 3 Department of Biological Sciences, Goethe University, Max-von-Laue-Straße 9, D-60438 Frankfurt, Germany 4 Smithsonian Conservation Biology Institute, National Zoological Park, Front Royal, VA 22630, USA; E-Mail: leimgruberp@si.edu 5 Department of Ecological Modeling, Helmholtz Center for Environmental Research-UFZ, D-04318 Leipzig, Germany; E-Mail: gunnar.dressler@gmail.com 6 Code 610.9, NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA; E-Mails: compton.j.tucker@nasa.gov (C.J.T.); jorge.e.pinzon@nasa.gov(J.E.P.) 7 Department of Geographical Sciences, University of Maryland, College Park, MD 20771, USA; E-Mails: dubayah@umd.edu (R.O.D.); gchurtt@umd.edu (G.C.H.) 8 National Socio-Environmental Synthesis Center (SESYNC), 1 Park Place, Suite 300, Annapolis, MD 21401, USA


Introduction
Climate change research often focuses on long-term productivity trends in vegetation at regional [1][2][3] to continental and global scales [4,5].A key variable in these studies of vegetation productivity is the satellite-derived Normalized Difference Vegetation Index (NDVI), that is either used directly as a measure of photosynthetic capacity [6], or serves as an input variable for climate change models predicting terrestrial net primary production (NPP [5,7,8]).On a global scale, these studies often ignore that NDVI may not be driven by climate alone.Because the majority of the Earth's ice-free land surface has been altered by human action [9][10][11][12] and because these same areas account for ~90% of terrestrial net primary production [13], trends in NDVI may also be linked to human land-use practices, such as land conversion, irrigation, and nitrogen deposition.Although several regional and local analyses demonstrate significant impacts on NDVI trends from changes in human land-use practices [14][15][16], global-scale analyses relating trends of NDVI to these direct anthropogenic effects have received little attention [5].
Here, we provide one of the first global analyses of direct anthropogenic effects on long-term trends in NDVI from 1981 to 2010 (however see [17] for an analysis that also include land use but on a different scale).First, we quantified variation in NDVI trends across "anthromes" (anthropogenic biomes, after [13]).Anthromes were introduced by Ellis and Ramankutty [13] and classify the earths land surface based on a combination of data on population density, land-use, and land cover [13].They characterize land surface based on global patterns of sustained, direct human interaction with ecosystems.Comparisons of trends in NDVI among anthromes thus represent a useful probe of direct anthropogenic effects on long-term trends in NDVI.We also examined the effect of human-population density, as a surrogate for the degree of anthropogenic effects.Lastly we tested for specific mechanistic linkages that may underlie spatial variation in trends of NDVI among ecoregions and investigated effects of human land conversion (measured as percentage village, urban areas and croplands per ecoregion), nitrogen deposition, and irrigation.

Results and Discussion
Globally, anthromes with intensive land use (dense settlements, villages, and croplands) all had significant (p < 0.05) and steep increases in annual mean NDVI over a 29-year period compared to other anthromes.Forests exhibited marginally significant increases (β = 0.0006, p = 0.06), whereas we found no significant trends for wildlands or rangelands (Figure 1A,C,E,G,I,K).We emphasize that wildlands in the anthropogenic biomes classification used here [13] simply lump all land surfaces that are free or almost free of any human footprint.They thus comprise extremely varied ecosystems across the globe where some may show positive and other may show negative trends.Our statement only refers to the average across all these diverse ecosystems.
We also calculated robust linear regression Theil-Sen (TS) slope estimators for each spatial location in the global NDVI data (using the GIMMS3g NDVI data).Across all anthromes, the medians of the per-pixel TS estimators agreed closely with the slope estimates of the annual mean (compare slope estimates in Figure 1(left) with medians in Figure 1(right)).Based on the average NDVI value in the respective anthromes, the TS estimators showed that NDVI of villages changed the most over the 29 year study period (5.9% increase), followed by dense settlements (4.3%) and croplands (3.7%).With a 2.5% increase, global net NDVI of forests increased less than half as much as villages.Global net NDVI increases in wildlands and rangelands were close to zero.We show that in years where global NDVI decreased, the global average NDVI of land types such as villages or croplands appear to have dampened minima compared to wildlands or forests (compare minima in Figure 1A,C,E,G,I,K).Variation induced by seasonal changes in NDVI, such as those associated with the interplay between human activities and monsoon dynamics [18], complicate interpretation of this pattern.
In our second analysis, we tested for the effect of human-population density, as a surrogate for the degree of anthropogenic effects, and found strong linkages between trends in NDVI and human population density.Here, we used ecoregions as spatial units [19].Ecoregions are distinct with respect to natural communities and resident species and thus ideal to test for spatial variation in anthropogenic effects on NDVI.In ecoregions with low human population densities, NDVI trends were small and close to zero, whereas in ecoregions with high population densities, NDVI increased significantly over time (Figure 2).For example, densely populated ecoregions with an average population density of ~500 people/km 2 showed an average increase of ~0.025 NDVI units over the 29 year study period (a 6% increase compared to the global mean).Except for Australia, average annual increases in NDVI were also greater in more densely populated ecoregions when data were analyzed at the continental scale.Europe, which has the greatest average population density across ecoregions, showed the strongest relationship between trends in NDVI and population density (Figure 3).
Results from multiple regression analysis of the importance of possible drivers of trends in NDVI.Analyses were performed using ecoregional means of per-pixel TS estimators of trends in NDVI (n = 770, coefficient determination: 0. Lastly, we tested for specific mechanistic linkages that may underlie spatial variation in trends of NDVI among ecoregions.In particular, we tested for human land conversion (measured as percentage village, urban areas, and croplands per ecoregion), nitrogen deposition, and irrigation and found that all three of these factors had significant effects.The percentage converted lands per ecoregion, the percent irrigated area per ecoregion, and the amount of nitrogen deposition all were positively related to increasing trends of NDVI (Table 1).Together these three factors explained 21% of the variation in trends in NDVI.

North America
Population density (km −2 ) Annual change NDVI 1 10 100 1000 −0.001 0.000 0.001 0.002 These results were robust to several methodological concerns and alternative hypotheses.For example, analyses excluding the moist tropical biome yielded similar results, suggesting the uncertainty of remote sensing estimates in the moist tropics [20] did not drive our findings (Table 2).Results were also robust to the hypothesis that humans have settled preferentially in areas with high NDVI (alternatively, that areas with high NDVI had greater absolute increases in NDVI) because standardizing the per-pixel TS estimators by the mean NDVI did not significantly change model outcomes (Table 2).Finally, results were robust to the spatial unit of analyses, because substituting the biologically defined ecoregions with square grid cells of length 500 km also did not substantially change model outcomes (Table 2).Ecoregions with the fastest increases in NDVI were predominately located in temperate Europe, the Mediterranean, the Sudano-Sahelian belt in Africa, western India, and northern China (Figure 4A).Ecoregions with the fastest decreases were predominately located in boreal forests of North America, as well as dry forests in South America (mostly Paraguay) and drylands of northeastern Australia.When comparing ecoregions with the fastest NDVI increases (the 5% of land surface with the fastest increases in NDVI, n = 73 ecoregions) against those showing the largest NDVI decreases (the 5% of land surface with the fastest decreases in NDVI, n = 38 ecoregions), we found that the ecoregions exhibiting the fastest NDVI increases had a greater average population density (median of 94 vs. 7 persons per km 2 , Figure 4), greater percentage of converted lands (88% vs. 9%, Figure 4B), greater percentage of irrigated lands (49% vs. 4%, Figure 4B), and greater deposition of nitrogen (573 vs. 209 mg N/m 2 /year, Figure 4B).
Ecoregions with the fastest changes in NDVI corresponded well with decreases and increases in NDVI shown in other global and regional studies [5,21].Decreasing trends in NDVI have been reported for ecoregions with the fastest decreases, like the Dry Chaco region in South America [14], northeastern Australia [22], and boreal Canada and Alaska [23,24].Likewise, increasing trends in NDVI have been reported for ecoregions with the fastest increases across Europe [25,26], northern China [27], central India [2], and the Sudano-Sahelian belt in Africa [28,29].
Although trends in NDVI in some of these ecoregions, especially those with negative trends in remote areas like boreal Canada and Alaska or northeastern Australia, have been predominately attributed to climatic effects [22][23][24], there is growing evidence for direct anthropogenic effects in almost all regions with significant human populations.For example, increasing trends in NDVI in India and northern China have been related to increasing rates of fertilization and irrigation [2,28].Likewise, while previous studies have chiefly related the Sahelian greening trend to climate, e.g., [7], new investigations point to land conversion to croplands as one of the main causes for changes in NDVI [29].In Europe, nitrogen deposition may be one of the main drivers of carbon sequestration in forests and-together with reforestation and climate changes-may well explain the large increases in NDVI observed across the continent [25,26,30].We explain 35% of the variability in trends in NDVI by population density alone (Figure 3D), supporting the hypothesis that anthropogenic drivers are contributing to NDVI increases in Europe.There are also examples for direct anthropogenic effects in ecoregions with rapidly decreasing NDVI.Results of a recent study on deforestation undertaken by the UN-FAO [31] show that forest loss in South America accounts for more than 50% of the global forest area loss in the past two decades.One region particularly affected is the Dry Chaco in Paraguay, a region with rapid land conversion from forest to soybean plantation [32], identified by our analyses as an extremum for decreasing NDVI.

NDVI Data and Pre-Processing
We used NDVI data from the Advanced Very High Resolution Radiometer (AVHRR) instrument developed into the GIMMS3g dataset (Global Inventory Modeling and Mapping Studies) [33,34].Previous versions of the GIMMS data set have been used extensively for testing temporal trend in NDVI [35], have been tested against ground controls points [36] and are thought to be the best long term data set for analyzing long term trends in vegetation dynamics [37].The current GIMMS3g dataset represents an improvement over previous AVHRR NDVI data sets because it was coprocessed with SeaWiFS and MODIS Terra data [20].The SeaWiFS and MODIS instruments have spectral bands specifically designed for vegetation monitoring, better atmospheric correction algorithms, reduced geometric distortions, and improved radiometric sensitivity [38].MODIS NDVI has been collected over the past 13 years from the Terra platform and the past 11 years from MODIS on the Aqua platform.For the respective MODIS overlap periods, the per-pixel trends in NDVI from GIMMS3g agree well with those from Terra and Aqua MODIS, except in the moist tropics, which are difficult for remote sensing due to humidity and cloud cover [20].MODIS Terra and Aqua NDVI data do not compare well for the tropics either.Using the longer-term GIMMS3g dataset, we conducted the analyses discussed below twice, once with all terrestrial pixels and once excluding the moist tropics.Since our analyses used yearly means (see below) and the studies cited above did not compare yearly means of NDVI between MODIS and GIMMS3g, we performed additional analyses comparing annual means between MODIS terra and Aqua, and the GIMMS3g data set.We sampled 34 1-degree areas across the globe and found a high linear correlation between the MODIS and GIMMS data sets.A Gaussian mixed model between annual means of MODIS Terra and GIMMS3g accounting for sampling area with a random effect showed very close relationships between both data sets (β = 0.74, t = 24.57,N = 298, coefficient det.= 0.99) as did the same comparison with MODIS Aqua (β = 0.75, t = 25.79,N = 247, coefficient det.= 0.99).If the data was standardized by the mean NDVI across years of each area to make areas more comparable, the relationship was still strong.MODIS Terra vs. GIMMS3g: (β = 0.68, t = 21.25,N = 298, coefficient det.= 0.65) and MODIS Aqua vs. GIMMS3g: (β = 0.69, t = 22.62, N = 247, coefficient det.= 0.65).
The GIMMS3g dataset of bi-monthly NDVI composites provides global coverage at a spatial resolution of 8 km × 8 km.For our analyses, we used the 672 composites for the 29 year period 1982-2010.We rescaled the data to the original floating-point NDVI scaling of -1 to 1, we identified and removed falsely included water pixels with a correction layer based on vector layers of continental coastline, lakes and reservoirs [39], and we applied a minimum threshold that set all pixels with NDVI values less than 0.05 to 0.05.NDVI values below 0.05 indicate non-vegetated areas like bare soil, water, snow, and ice.This 0.05 threshold is based on an analysis of stable pixels in the Sahara Desert that had an average NDVI of 0.094 and has been used previously in studies examining trends in NDVI [40].

Trend Analysis & Anthropogenic Effects
We estimated linear trends of annual means of NDVI for the 29 point time series from 1982 to 2010 for each spatial location in the data set (n = 1,975,140 pixels) using the Theil-Sen (TS) median estimator (function mblm in package mblm in R).The TS estimator represents a robust linear regression method that calculates the median slope among all pairs of observations.Commonly used for studies in vegetation trends based on NDVI time-series [20,41,42], the TS estimator is relatively insensitive to outliers and is accurate even for skewed and heteroscedastic data.
To relate these estimates to anthropogenic effects we performed three different analyses.First, we evaluated trends in NDVI for each of the six groups of anthropogenic biomes identified by Ellis and Ramankutti [13]; wildlands, rangelands, forests, croplands, villages, and dense settlements.We calculated histograms and medians for the TS estimators for each group of anthropogenic biomes and the global trend for each group based on the yearly average of NDVI for all pixels of a group using generalized least squares (GLS).We incorporated the temporal correlation structure with a 1st order autoregressive moving average.For all GLS analyses, we calculated the coefficient of determination based on the improvement of the log-likelihood from the fitted null model (intercept only) to the full fitted model (function r.squaredLR in R package MuMIn, [43]).
Our second and third analyses were carried out on the scale of ecoregions.We used the World Wildlife Fund's classification scheme, which distinguishes ecoregions [19] that are distinct with respect to natural communities and resident species.We first tested whether human-population density, as a surrogate for the degree of anthropogenic effects, was related to NDVI trends.We estimated human population density for each ecoregion as the mean of the 1990 and 2010 population estimates [44].We next tested for the effects of land conversion, nitrogen deposition, and irrigation on NDVI trends by performing multiple regressions that included these factors as predictors.As a proxy for converted lands, we used the percent area covered by dense settlements, villages, and croplands per ecoregion [13].
For nitrogen deposition we used the mean nitrogen deposition per ecoregion in 1993 ( [45], an average from the 1980s to the 2010s would have been preferable but was not available).For irrigation we calculated the percentage of irrigated area per ecoregion based on the UN-FAO irrigation map [46].
For both analyses we used the ecoregional means of the per-pixel TS estimators as the response variable.We used GLS accounting for spatial autocorrelation with an exponential spatial correlation function and accounting for differences in sizes of ecoregions using the inverse of ecoregion size as a weight in the variance function.For the multiple regression, we tested predictors for collinearity and used a cutoff correlation coefficient of 0.7 [47].
For robustness, we first repeated the analyses excluding all pixels located in the moist tropical broadleaf biome, because of unresolved and ongoing debates about the validity of remote sensing estimates in moist tropical regions [20,48] and the poor agreement between (1) GIMMS3g and MODIS Terra NDVI and (2) MODIS Terra and Aqua NDVI data in the moist tropics found by Fensholt and Proud [20].Second, to ensure that our results do not stem from humans choosing areas with greater NDVI or from areas of greater NDVI exhibiting greater absolute increases of NDVI, we repeated our analyses by standardizing the per-pixel TS estimators by the mean NDVI value of each pixel.Third, even though ecoregions represent ecologically distinct units and as such are appropriate geographic units for examining the effects of covariates on trends in NDVI, we also repeated our analyses substituting squared grid cells of 500 km by 500 km for ecoregions.Lastly, to examine whether our results hold true also within continents and not just on the global scale, we performed the analyses for population density separately for each continent.Because some continents have only relatively few ecoregions, we did not perform the by-continent analyses for the multiple regression analyses.

Conclusions
In conclusion, our study provides one of the first global synopses quantifying the effect of the Earth's human footprint on trends in photosynthetic capacity.Globally, more than 20% of the variability in NDVI trends can be associated with human land-use practices.Consequently, global climate change models and analyses of primary productivity should incorporate direct anthropogenic effects.Land use practices in densely populated areas, and other areas featuring major human impacts through land conversion, irrigation, and nitrogen deposition, may contribute more to global changes in photosynthetic capacity than previously anticipated.

Figure 2 .
Figure 2. Annual change of NDVI of ecoregions, based on ecoregional means of per-pixel Theil-Sen estimators) vs. human population density.Circle diameters are proportional to the size (area) of each ecoregion.Red line indicates trend based on generalized least squares model with exponential spatial autocorrelation structure and accounting for differences in sizes of ecoregions in the variance function (β: 0.00032, p < 0.0001, Δ AIC: 105, coef.determination: 0.13).For corresponding analyses by continent, see Figure 3.

Figure 3 .
Figure 3. Annual change of NDVI of ecoregions (Africa (A), Asia (B), Australia (C), Europe (D), North America (E), South America (F), based on ecoregional means of per-pixel Theil-Sen estimators) vs. log10 of human population density.Circle diameters are proportional to the size (area) of each ecoregion.Red line indicates trend based on generalized least squares model with exponential spatial autocorrelation structure and accounting for differences in sizes of ecoregions in the variance function.β: coefficient of log10 of human population density, d: coefficient of determination, significance codes: *** p < 0.0001, ns: not significant.

Figure 4 .
Figure 4. (A) Trends in NDVI across the globe, 1981-2010.Ecoregional [19] extremes for NDVI increase (defined as the 5% of land surface with the fastest increases in NDVI, n = 73 ecoregions) are in red, whereas ecoregional extremes for NDVI decrease (defined as the 5% of land surface with the fastest decreases in NDVI, n = 38 ecoregions) are in blue; (B) Boxplots contrast the ecoregional extremes in A for increases (n = 73) and decreases (n = 38) in terms of NDVI trends, population density, percent converted lands [13], percent irrigated lands, and nitrogen deposition.

Table 2 .
Additional models relating trends in NDVI to (A) population density and (B) to nitrogen deposition, irrigation, and converted lands in order to ensure robustness of results against potential methodological concerns (see main text).