Title Examining the linkage between shrub encroachment and recent greening in water-limited southern Africa

. There have been reports of widespread increases in satellite-derived vegetation indices in drylands across the globe. At the same time, there are numerous field-based observations of increases in the density of small woody plant species, a process known as shrub encroachment. We compare ground-based observations of shrub encroachment with information from 14 years of Tropical Rainfall Measuring Mission (TRMM) precipitation data, Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) and MODIS land cover data for southern Africa, focusing on areas that receive less than 1000 mm of precipitation per year. Cumulative precipitation and maximum NDVI were computed for each year and linear regression analysis was used to correct for interannual variability in precipitation. Accounting for precipitation effects, we found an underlying, overall greening in southern Africa, with NDVI increasing þ 3.34 % on average over the 2000–2013 study period. Over 15 % of the region is undergoing statistically significant change, strongly biased towards greening. The strongest greening was in a coherent band stretching from northern Namibia to Lesotho. Ground-based reports of shrub encroachment tend to overlap with strong greening; we suggest that these processes are linked. Reports that intersect our period of record showed overwhelmingly positive trends, indicating that our method is suitable for detecting shrub encroachment. Shrublands and savannas accounted for 78 % of the study region. Savannas did not change appreciably ( (cid:2) 0.21 % total change) from 2000–2013 while shrublands were greening faster than any other land cover type over the same period ( þ 8.93 % ). Fast rates of greening in this biome may be indicative of shrub expansion. Large-scale shrub encroachment will have important consequences for dryland degradation and global carbon dynamics.


INTRODUCTION
Previous studies have reported systematic increases in remotely sensed vegetation indices, a phenomenon commonly referred to as ''greening,'' occurring in the past three decades in temperate and subtropical water-limited environments across the globe (Bai et al. 2008, Helldén and Tottrup 2008, De Jong et al. 2011, Donohue et al. 2013). According to some estimates, approximately 70% of drylands are undergoing greening (Bai et al. 2008, substantially more than exhibit decreases in vegetation indices (or ''browning''). Positive trends in greenness have been documented over Africa using different datasets and trend detection procedures (De Jong et al. 2011, Gaughan et al. 2012 1 September 2015 v Volume 6(9) v Article 156 2013) and increases in vegetation have been broadly linked to increases in atmospheric CO 2 (Donohue et al. 2013), increased water availability (Hickler et al. 2005) and changes in plant community composition (Bai et al. 2008). Such trends have been useful in tracking changes in regional and global productivity (Zhu and Southworth 2013) but the causes and phenomenological specifics of large-scale greening remain ambiguous; it is unclear whether these changes are the result of precipitation trends or changes in plant community composition induced by other drivers. It is also unclear whether such changes represent an enhancement of understory herbaceous vegetation, an increase in woody cover, or both.
At the same time, there have been field-based reports of increases in the density of small (generally ,2 m) woody plants over southern Africa (Wigley et al. 2010, Buitenwerf et al. 2012, Rohde and Hoffman 2012, O'Connor et al. 2014), a phenomenon sometimes referred to as ''bush encroachment'' (Meik et al. 2002), ''bush thickening'' (Joubert et al. 2008), or ''shrub encroachment'' (Maestre et al. 2009). This usually is accompanied by decreased herbaceous biomass (Eldridge et al. 2011). More variable ecosystem changes that have been linked to shrub encroachment are the loss of mature trees (Herrmann and Tappan 2013), decreased plant (Ratajczak et al. 2012) and animal biodiversity (Meik et al. 2002, Blaum et al. 2007, Kaphengst and Ward 2008, increase in bare soil (D'Odorico et al. 2012), and a loss of ecosystem services (De Klerk 2004). Thus the classic model of shrub encroachment is commonly associated with land degradation. More recent findings challenge this paradigm, suggesting that the effects of woody plant encroachment are more nuanced and in some cases can even reverse degradation (Maestre et al. 2009). Shrub encroachment enhances soil carbon sequestration (Eldridge et al. 2011), and woody plant expansion on a large scale may have important consequences as a global dryland carbon sink (Poulter et al. 2014). Despite divergent notions regarding the fate of encroached systems, it is clear that increased woody cover results in substantial shifts in ecosystem structure and function (Huxman et al. 2005, Scott et al. 2006, Eldridge et al. 2011. Just as browning has been explicitly linked to land degradation (Evans andGeerken 2004, Bai et al. 2008), greening trends have been interpreted as an improvement in ecosystem status (Helldén and Tottrup 2008). Thus widespread greening seemingly runs counter to the more prevailing paradigm of desertification in semiarid ecosystems. But this greening could result from large-scale woody plant encroachment, which is generally interpreted as a mechanism of land degradation (Ward 2005, Eldridge et al. 2011, D'Odorico et al. 2012. In this case increased vegetation indices would not necessarily be indicative of dryland rehabilitation (see Herrmann and Tappan 2013). Despite reports of both greening (Helldén and Tottrup 2008) and shrub encroachment in southern Africa (Bester 1999, Wigley et al. 2010), evidence of their linkage at the regional scale is still missing. Here we investigate the relationship existing between them, and show how greening coincides with reports of shrub encroachment in the region. Due to the substantial ecological (Eldridge et al. 2011) andeconomic (De Klerk 2004) consequences of shrub encroachment, it is necessary to characterize the extent of ongoing structural changes in woody vegetation. The specific goals of this study are to (1) examine the extent of recent greening and browning in southern Africa using a satellite time series from the high quality Moderate Resolution Spectroradiometer (MODIS), (2) identify and remove the influence of climatic variability on these trends, and (3) use land cover data and previous reports of bush encroachment to interpret these trends and critically review the evidence that increased greenness is a manifestation of increased woody cover.

Study area
This study focuses on continental Africa south of 108 S to examine the relationship between greening and shrub encroachment. This region contains large precipitation gradients and multiple biomes (grassland, shrubland and savannas) allowing us to characterize the prevalence of greening over several distinct vegetation types. As the focus of this study is on water-limited regions, we limit our study to regions where mean annual precipitation is less than 1000 mm. v www.esajournals.org 2 September 2015 v Volume 6(9) v Article 156 This is somewhat more inclusive than what is typically considered to be semiarid (with precipitation based characterizations varying from 500 to 700 mm as an upper limit; Nicholson 2011), but it includes regions that may be seasonally water-limited.

Precipitation data
Monthly precipitation estimates in the form of the TRMM 3B43 product were used in this study. These data give average monthly precipitation rates that are blended estimates of the TRMM satellite and gauge data in a 0.25 3 0.258 gridded product with global coverage from 408 N to 408 S. Validation of this dataset indicates excellent accuracy in semi-arid regions of Africa (Nicholson et al. 2003).

Vegetation data
At large scales, greenness is most commonly characterized by the Normalized Difference Vegetation Index (NDVI), a metric that reflects above-ground green biomass based on the differential reflectance of red and near infra-red spectral bands. Almost all studies that have identified greening trends in drylands have done so using data derived from the Advanced Very High Resolution Radiometer (AVHRR) such as the Global Inventory Modeling and Mapping Studies (GIMMS; Tucker et al. 2005) dataset (Helldén and Tottrup 2008, De Jong et al. 2011, Donohue et al. 2013). Although NDVI time series from MODIS products have a shorter period of record compared with those from AVHRR products, MODIS-based products have advantages with respect to spectral fidelity and sensor continuity. Featuring onboard calibration, MODIS measures over narrower spectral bands compared with AVHRR (Huete et al. 2002).
The MODIS Terra Earth observing satellite has been collecting data on a continuous basis since February 2000. While this lifespan is shorter than AVHRR, which extends back to 1981, systematic drift in the AVHRR sensors and the patchwork of multiple instruments make long-term trend analysis potentially hazardous (Kaufmann et al. 2000), despite best efforts to apply appropriate corrections. Long-term changes in NDVI are small relative to seasonal variability in NDVI, and different products can potentially lead to different answers with respect to detecting these subtle changes. For example, Fensholt and Proud (2012) reported trends in GIMMS and MODIS NDVI using 10 years of overlapping data, and showed nontrivial differences in the geographic distribution of trend direction and magnitude in southern Africa. Thus, using the available MODIS data to assess changes in vegetative indices, as we do here, is a current research priority. We enhance our ability to identify trends in this shorter record by removing effects caused by interannual variability in climatic forcing.
MODIS NDVI data (MOD13A3 Level 5 product, 1-km spatial resolution, 1-month temporal resolution) were acquired over the study region from the first available MODIS images in February 2000 through June 2013 (see http:// e4ftl01.cr.usgs.gov/MOLT/MOD13A3.005/). This product uses atmospherically-corrected bidirectional reflectances and are masked for water, clouds, and heavy aerosols. These data were coarsened to the TRMM 0.258 resolution by averaging all high and moderate quality (quality flags 0 and 1) MODIS pixels in each TRMM cell at each time step. This resampling procedure allows us to consider only pixels that are of sufficient quality by removing data that may be contaminated by cloud cover. Resampled pixels that contained only water, urban, and bare ground land cover types were excluded from the analysis. This procedure and all subsequent statistical analyses were performed in MATLAB (version 2012a; MathWorks 2012).

Finding yearly precipitation and vegetation metrics
The phenology of the study region is generally water-limited (Nemani et al. 2003), driven by a summer (November-April) wet season. As such, building a yearly vegetation metric based on calendar year, which bisects the wet season, is not appropriate (De Jong et al. 2011). Instead, for each pixel we identified the beginning of the vegetative year (VY) to be the month with the lowest average monthly precipitation over the whole period of record. Yearly maximum NDVI (NDVImax) over each VY was then calculated. Yearly precipitation metrics were created by summing monthly precipitation values from the beginning of each VY through the month of NDVImax, determined on a per-year basis. We do not account for proposed inter-annual memory effects (Richard et al. 2008). This phenomev www.esajournals.org 3 September 2015 v Volume 6(9) v Article 156 non is manifested during the greening-up phase of the early wet season, so we do not expect memory effects to significantly influence observed NDVImax.
The first year of record (1999)(2000) did not comprise a full vegetative year. We included this year in the analysis because the yearly maximum in NDVI is expected to occur after the beginning of available MODIS data in February 2000. This is especially likely in the most water-limited southern areas of the study area, where there is a strong uni-modal yearly precipitation and vegetation cycle.

Trend analysis
Trends in observed yearly NDVI (''raw trends'') were computed using a nonparametric Thiel-Sen estimate of linear slope with a Mann-Kendall significance test. The Sen slope estimates linear trend using the median slope of all pairwise combinations present in the time-series and provides an outlier-resistant, non-parametric estimate of linear trend (Theil 1950, Sen 1968. A trivial explanation of increased greenness would be a simultaneous increase in precipitation. As our goal is to investigate underlying trends in vegetation (not caused simply by precipitation, the main climatic constraint) we must remove the effect that yearly precipitation has on observed NDVImax. We correct observed time-series for precipitation using a residual trend analysis procedure proposed by Evans and Geerken (2004) and widely used to correct for climatic constraint on vegetation (Evans and Geerken 2004, Herrmann et al. 2005, Wessels et al. 2007, Li et al. 2012, Andela et al. 2013, Fensholt et al. 2013. In this method linear regression is used to model yearly maximum NDVI as a function of yearly cumulative precipitation (see Data and methods: Finding yearly precipitation and vegetation metrics) to derive a pixel-wise sensitivity to precipitation: where V obs,i is the observed NDVImax for each year i, P i is yearly cumulative precipitation starting from the beginning or the rainy season, s is the sensitivity of the pixel to precipitation, given by the linear least squares slope, b 1 is a fitted intercept and e i are the residuals. Yearly estimated NDVImax (V est,i ) is calculated from this relationship using P i , s and b 1 : Finally, the NDVI residuals are found by subtracting V est,i from V obs,i . Temporal trends in these residuals (''corrected trends'') are also determined using the Theil-Sen method. As Andela et al. (2013) point out, this procedure assumes that any temporal trends detected from these anomalies are independent of long-term, concurrent trends in precipitation. We use discrepancies in the raw and corrected trends to locate regions where climate fluctuations have played a central role in determining greening or browning trends.

Land cover data
Land cover data are used to aid in interpretation of regional trends. We use the MODIS land cover product (MCD12C1, yearly values from 2001 to 2011, 500-m resolution). The MODIS algorithm uses ensemble decision trees to categorize land cover into one of 17 International Geosphere-Biosphere Program classes (Friedl et al. 2010). This dataset has estimated classification accuracy of 74.8%. However, much of the misclassification of the land cover classes that dominate our study region is due to confusion between ecologically similar classes (e.g., savanna versus woody savanna; Friedl et al. 2010). Land cover was coarsened to 0.258 resolution by selecting the most abundant land-cover class in each TRMM pixel over the whole period of record, provided that this land cover type accounted for more than half of the pixels. In some pixels the most abundant land cover class amounted to less than half of the pixels. These mixed pixels are still considered in this analysis wherever ''all'' pixels are referenced. For this reason, some land cover types that commonly occur as fine-scale, heterogeneous patches among other land cover types, such as cropland mosaics or ecotones, are underrepresented in the coarsened dataset.

Observed (raw) trends
Trends in observed NDVImax over the whole study region between 2000 and 2013 were equivocal (48.4% of pixels greening, as opposed v www.esajournals.org 4 September 2015 v Volume 6(9) v Article 156 to 51.2% browning). The average pixel-wise change using uncorrected trends was a slight decrease (À0.12%) over 14 years of record. However, trend directions and magnitudes showed coherence at finer spatial scales (Fig.  1A). The strongest browning was found in Angola and in a band stretching from northern Botswana to Mozambique. Greening was present to the largest extent in central to northern Namibia, with weaker trends in the Western Cape Province of South Africa and north of Lesotho. The areas with significant greening trends ( p , 0.1) were generally found around the foci mentioned above (Fig. 1B). Pixels with positive NDVI trends accounted for 41.5% of significantly changing pixels, and 4.9% of pixels overall, while pixels with negative NDVI trends accounted for 58.5% of significantly changing pixels, and 6.9% of pixels overall.

Corrected trends
Fig. 2 depicts the trends corrected for precipitation in NDVImax over southern Africa. There is an overall increase in greenness (þ3.34%) over the study period and 61% of pixels show positive trends. The corrected trends reflect some of the regional variability observed in the uncorrected trends. Browning is still seen on the Angolan coast and more spatially extended, albeit weaker negative trends span central Botswana to coastal Mozambique ( Fig. 2A). The Cape and central regions of South Africa, southern Botswana and most of Namibia show extended swaths of strong (.5 3 10 À3 NDVI per year) greening. The Mann-Kendall test for significance identifies clusters within these regions as undergoing significant change ( p , 0.1). Overall, statistically significant greening pixels outnumber significantly browning pixels by a factor of 2.8 ( Table 1).
The magnitude of corrected trends differs from uncorrected trends as well (Fig. 3). Greening tends to be stronger than browning, which we illustrate by examining the relative number of greening and browning pixels in either tail of the trend distributions, representative of ''strong'' trends. For example, the ratio of greening to browning pixels with moderately strong trends exceeding 65 3 10 À3 NDVI per year is 3.5 to 1. When we further increase this threshold to 67.5 3 10 À3 NDVI per year (corresponding to very strong trends and only present in 2.4% of all pixels), the ratio of greening to browning ratio jumps to 5.0 to 1. Nine percent of the study region in southern Africa is greening at least 5.0 3 10 À3 NDVI per year over 14 years of record, corresponding to a total change of 0.07 NDVI units. Given that the study area has an average NDVI max of 0.59, it is clear that greening of this scope and magnitude reflects substantial changes in vegetation and requires explanation.
As a whole, correcting for the influence of precipitation enhanced greening and attenuated browning. The correction procedure reduced the number of significantly browning pixels from 370 to 219 and increased the number of significantly greening pixels from 263 to 606. This is illustrated in a map of differences between corrected and raw trends (Fig. 4). The largest difference between corrected and observed trends occurs in coastal Angola, and is indicative of a larger pattern where some of the browning trends that were observed were significantly reduced or removed once precipitation was accounted for. Although positive corrected trends are generally stronger than corresponding uncorrected trends, there are some regions where the corrected trends are weaker. One such example is northern Namibia, indicating that here a component of

Land cover and greening
Open shrubland and savanna were the most abundant land cover in the study region, together accounting for 78.5% of the study area (Table 1). Cropland and grassland made up 1.56% and 5.09% of the region, respectively. Different land cover classes were generally found in spatially coherent bands that reflect variability (i.e., spatial patterns) in mean precipitation (Fig.  5). Greening tended to occur in more arid areas and land cover classes (croplands and open shrublands) while browning was dominant in wetter savannas and woody savannas (Table 1, Fig. 6). These two land covers experienced 0.62% and 1.62% decreases in greenness from 2000-2013. Statistically significant changes in greenness occurred in 10.6% of woody savanna pixels (85% percent of which were in the form of significant browning). Shrublands showed an 8.93% increase in greenness from 2000 to 2013. Trends were statistically significant in 18.9% of shrubland pixels, 95.1% of which were in the form of significant greening.

Literature-derived evidence of shrub-encroached areas
We compiled 24 different studies that report woody shrub encroachment over our study region (Table 2), amounting to 28 distinct sites. Eighteen of these sites report shrub encroachment within the time period of our analysis. Reports of encroachment tend to intersect with areas that exhibit greening, especially in regions of significant greening. Two of the studies (Moll andTrinder-Smith 1992, Milton-Dean 2010) in the Cape Province of South Africa document invasive species. Other than these two cases, all other studies report expansion of endemic woody species, generally stretching in a band from central Namibia to Lesotho (Fig. 2B, pink markers). For each site we extracted the rainfallcorrected trends of the pixel whose centroid was closest to the coordinates presented in each study. Despite a mismatch in the scale of the 0.25 pixels at which trends are calculated with ground based reports (commonly hectare-scale plots), this provides a gross validation of our technique over the sites that intersect with our Year À1 wide bins to assess pixel frequency. Light grey and black represent uncorrected and corrected trends, respectively. Solid and dashed lines represent all pixels and only pixels with a significant ( p , 0.10) greening or browning trend, respectively. v www.esajournals.org 7 September 2015 v Volume 6(9) v Article 156 study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Of the 18 distinct reports of encroachment that intersect with our study period, 17 sites (94.4%) exhibited a positive trend (Fig. 7). The average trend at these sites was þ2.945 3 10 À3 NDVI per year. Sites that report woody plant encroachment before year 2000 had an average trend of þ2.289 3 10 À3 NDVI per year but only seven out of ten sites (70%) showed positive trends.

Patterns of change in African drylands
Greening is focused in a NW-SE band stretching from Namibia to South Africa (Fig. 2). This is contrary to previous studies utilizing AVHRRderived data, which identify central Botswana as a focus of greening and little change in surrounding regions (De Jong et al. 2011). Instead, we identify weak-to-moderate browning over Botswana, although this country has very few pixels (either greening or browning) with significant slopes (Fig. 2). The magnitude of the browning in coastal Angola is likely associated with decreased productivity due to severe drought that occurred towards the end of the record, most notably in 2012-13 (Rembold et al. 2013). Removing the influence of precipitation significantly attenuates browning in this region, as is expected to occur with trends driven by water availability. The magnitude of precipitation-corrected trends was, in some areas, substantially different than raw trends. This has two important implications. Firstly, and in line with other studies (Zhu and Southworth 2013), it indicates that precipitation plays a primary role in determining dryland productivity, but it cannot explain all trends. Secondly, identifying signals of structural change in drylands without specifically accounting for the interannual variability in precipitation may simply reflect concurrent changes in water availability. This may be the cause of some of the discrepancies in the location of greening compared with previous research that does not apply corrected trends, but this is not to say that other factors such as different periods of record (AVHRR: 1982-present;MODIS: 2000-present), differences in the sensors, or spectral post-processing do not also play an important role. Although we identify regional greening and browning in different locations than previously reported, in a more general sense we confirm the large-scale dryland greening using the high quality MODIS dataset, once interannual variability in precipitation is taken into account.  The linkage between greening and shrub encroachment Our explanation of the large-scale greening is based mainly on two lines of evidence. The first is the correspondence between the locations of greening trends, especially significant ( p , 0.1) greening, and compiled reports of shrub encroachment (Fig. 2B). South Africa's Cape region These studies did not explicitly observe temporal changes in woody cover, but nonetheless characterized the landscape as having undergone woody plant encroachment (either through comparison with a reference or through expert assessment). Fig. 7. Histogram of rainfall-corrected trends at sites observed to have undergone encroachment. Locations of sites reporting woody plant expansion before our study period (pre-2000) are shown as white circles. Locations of sites reporting woody plant expansion that intersect our study period are shown as black circles.
v www.esajournals.org 9 September 2015 v Volume 6(9) v Article 156 notwithstanding, the general pattern of greening is in an NW-SE band, approximately 1800 km in length. The ground-based evidence of shrub encroachment largely mirrors this spatial pattern. While this evidence is essentially correlational, it is supported by the observed greening trends at 17 of 18 reference encroached sites. This seems to indicate that the residual trend analysis presented here detects woody plant encroachment. Studies that have used remotely sensed proxies of woody cover complement our findings. For example, Andela et al. (2013) used a time-series of satellite derived vegetation optical depth as an indicator of woody plant biomass. Their results support our findings; they observed increases in vegetation optical depth over the regions of southern Africa where we observe greening and interpreted it as woody plant expansion. Mitchard and Flintrop (2013), on the other hand, did not find evidence of increased woody cover in the study region using yearly minimum AVHRR NDVI. We believe this discrepancy may stem from differences in methodology and the fact that they did not consider climatic influence on vegetation in southern Africa, which had experienced severe droughts over the period of their analysis . This simultaneously highlights the need to understand the inherent differences between available data sets (such as MODIS and AVHRR) and the usefulness of residual trend analysis in removing the influence of interannual climatic variability.
The relationship between land cover and greening provide another form of evidence that the strongest greening is a manifestation of shrub encroachment. The location of greening is heavily skewed towards shrublands, regions defined as having woody vegetation less than 2 meters tall with 10-60% canopy cover in the IGBP classification scheme that we use in this study. Of the strongest 5% of trends, 73% are shrublands, despite this land cover accounting for 36% of the study region. Regions with this land cover type do not support extensive understory herbaceous biomass. Greening here may be inferred as a signal of increased woody plant cover. The changes occurring in shrublands are unparalleled in other natural land cover types (i.e., excluding croplands and crop mosaics) in extent and magnitude. Although grasslands seem to be greening overall, the proportion of grassland pixels exhibiting increased greening (62%) and the average trend (8.61 3 10 À4 NDVI per year) do not differ substantially from the averages across the whole region (see Table 1).
Due to other environmental constraints, the progression of shrub encroachment is expected to asymptote eventually (Roques et al. 2001). Thus we might not expect to see greening trends coincide with the location of studies that were conducted before our satellite record. However, these sites are generally greening too (Fig. 7). This may suggest that shrub thickening is still occurring at some sites where encroachment was documented decades earlier. Follow-up studies are needed to discern whether these sites have indeed reached a steady state in woody cover, or if the thickening of woody plants is more persistent than previously surmised.

Potential drivers of shrub encroachment
Explanations of shrub encroachment vary widely in the literature and authors have assigned different importance to drivers based on the particular system they are studying. Nonlinear interactions of drivers and encroached systems can sometimes make it difficult to attribute primacy to any one factor (D'Odorico et al. 2012). The most commonly invoked drivers are overgrazing by livestock, changes in fire frequency, and carbon dioxide fertilization.
The role of fire in disrupting the dominance of woody plants is well established. While surface fires rarely outright kill woody plants, frequent fire may prevent juvenile woody plants from recruiting into larger size classes (Dublin et al. 1990, Anderies et al. 2002, D'Odorico et al. 2006, Hanan et al. 2008. This effect has been confirmed by long-term experimental fire exclusion experiments, which demonstrates the proliferation of small woody plants when fire is limited (Du Toit et al. 2003).
Herbivory has disparate impacts on the rate and pattern of shrub encroachment. There has been widespread loss of ''mega-herbivores'' from southern Africa along with increased populations of grazing herbivores. Browsers such as elephants suppress tree growth through damage to roots and branches (Dublin et al. 1990). On the other hand, it is well established that grazing can induce a compositional switch to shrubs and reduce positive consequences of encroachment v www.esajournals.org 10 September 2015 v Volume 6(9) v Article 156 (Eldridge et al. 2013). Grazing has been implicated at many of the sites in our study region (Moleele andPerkins 1998, Kgosikoma et al. 2012). In addition, high grazing rates may interact with fire by lowering fuel loads, intensifying the rates of woody plant expansion (Roques et al. 2001).
Changes in land use other than ranching can also contribute to the proliferation of woody plants. Shackleton et al. (2013) describe the Wild Coast region of South Africa where there has been a systematic trend of farm abandonment and forest revegetation. Here, recolonization of natural species has led to substantial increases in woody species richness and overall woody cover.
The carbon dioxide fertilization hypothesis has also been used to explain large-scale shrub encroachment. This phenomenon, in which elevated atmospheric carbon dioxide enhances photosynthetic efficiency, has been demonstrated to benefit C3 type woody plants more than C4 grasses in both controlled environments (Kgope et al. 2010) and in natural systems (Morgan et al. 2007). There are reports of both remotely sensed greening and ground measurements of woody cover where the only obvious driver is CO 2 (Buitenwerf et al. 2012, Andela et al. 2013, Donohue et al. 2013. Whereas grazing and fire are influenced by site-specific management and land use, carbon dioxide acts as a global driver. The nature of observed greening, in which coherent trends extend spatially through disparate communities, provinces and countries, may indicate a driver that acts on the same scale. For example, where there are sharp boundaries in average livestock stocking rates (such as the western Botswana-South Africa border; Andela et al. 2013) there is no accompanying difference in NDVI trends. Therefore, we suggest that carbon dioxide fertilization may play a primary role in promoting shrub encroachment and subsequent greening, following Andela et al. (2013) among others (Helldén andTottrup 2008, Donohue et al. 2013). Greening is widespread, but it is not universal. In some regions CO 2 fertilization may be offset by other factors, such as overgrazing (Wigley et al. 2010). Predicting regions where CO 2 fertilization may or may not be manifested, the strength at which it will occur, and how it will interact with changing climate and land use, will be a research challenge. Furthermore, Bond (2012) suggests that increasing CO 2 may act synergistically with other drivers to tip the balance in favor of woody species where single factors, such as herbivory, were previously not strong enough to promote an expansion of woody species.
One aspect of CO 2 fertilization is the enhancement of plant water use efficiency by allowing more efficient carbon exchange with the environment and limiting stomatal water loss (Bond 2012, Donohue et al. 2013. Given increased atmospheric CO 2 , we may expect to see more pronounced fertilization in more water-limited areas. Our results may support this conclusion. When aggregated across all land cover types, we find that average precipitation-corrected greening is higher in more arid environments (Fig. 8). This may indicate that the benefits of CO 2 fertilization are more pronounced in arid regions, while other factors (e.g., fire, land use) dominate in more mesic systems.

Other possible explanations of observed greening
Greening may be manifested in systems undergoing increases in grass cover, woody cover, or both. An alternative explanation to shrub encroachment could be increased herbaceous biomass. We do not believe that this alone could account for the observed satellite-detected trends. Firstly, in arid systems that do not support an herbaceous understory, greening is necessarily an indicator of woody plant expansion. Where grasses are present, there is limited ground-based evidence of woody dieback and conversion of shrub cover to grassland, whereas reports of increased woody cover are more numerous (Table 2). Masubelele et al. (2014) report increased grass cover using long-term data from eight sites in southern Africa, but in the most recent sampling period ) most relevant to our discussion, half of these sites show increases in the ratio of shrub cover to grass cover. Secondly, carbon dioxide fertilization, a possible large-scale driver of structural change, is expected to favor increases in woody plants rather than herbaceous cover (as discussed in the previous section). Finally, the grass component of NDVI is more variable than that of woody cover on an annual basis (Archibald and Scholes 2007), largely in response to available moisture (Scanlon et al. 2002). Gradual v www.esajournals.org 11 September 2015 v Volume 6(9) v Article 156 increases in mean herbaceous cover due to factors other than increased precipitation could contribute to the greening signal. Increases in fire frequency could also account for both short and long term increases in herbaceous biomass, instigating a feedback towards more fuel, fire, and grasses. However, fire suppression is a common management technique and fire frequency has declined in southern Africa over the past half-century (Mouillot and Field 2005). Increased grass cover may be spurred on by land rehabilitation following a shift to more sustainable management regimes (Masubelele et al. 2014). A common example of this is decreased cattle stocking rates (Angassa andOba 2010, Eldridge et al. 2013). Further research is needed to characterize the nature and extent of rehabilitation following land abandonment and shifts in management techniques. The extent to which these may contribute to a greening signal is not clear.
We cannot absolutely rule out climate as a driver of widespread greening. While our precipitation de-trending accounted for rainfall on a yearly scale, we do not consider how sustained drought and decadal oscillations indirectly influence greening. Sustained departures from longterm climate have been shown to influence vegetation composition (Hickler et al. 2005). For example, sustained drought in fire prone savannas may reduce fuel loads and subsequently fire frequency. Longer fire return intervals can act as a demographic window to allow higher shrub recruitment in drought-afflicted areas. On the other hand, multiple years of above average rainfall may induce a compositional switch to the recruitment of less xerophytic species (Kraaij and Ward 2006). For example, wetter-than-normal conditions during the 1970s are associated with rapid woody plant encroachment (O'Connor et al. 2014). In this way the cumulative effect of extended climate anomalies is not fully accounted for. Finally, total annual precipitation is a good proxy variable of the general precipitation regime, but other hydrologic variables such as the average rainstorm frequency and depth can be more precise determinants of maximum woody plant cover (Good and Caylor 2011). Fig. 8. Distribution of corrected trends for data grouped into 100 mm precipitation bins. Boxes represent the first and third quartile of the data. Vertical lines show the 5th and 95th percentiles of trends in each bin. The median of each bin is shown as a horizontal line inside each box. Bins with medians significantly higher or lower than zero are labeled with stars above and below the bins, respectively. Significance at 0.05 and 0.01 levels are shown with one and two stars, respectively. Statistical significance is determined using a one-sample sign test (Gibbons and Chakraborti 2011) with resulting p values adjusted for multiple comparisons using the sequential Bonferroni technique (Holm 1979).
v www.esajournals.org 12 September 2015 v Volume 6(9) v Article 156 Accounting for specific climatologies, whether storm depth or seasonality, is a natural extension of the residual trend procedure used here and may provide valuable information about the effects of forecasted climate change may have on ecosystem structure and function.

CONCLUSIONS
This study examines the extent of greening in southern Africa and its coincidence with reports of shrub encroachment. This research complements previous studies that use the popular AVHRR GIMMS dataset by calculating trends in MODIS and accounting for year-to-year variability in precipitation with TRMM data. These corrected trends reveal general greening over southern Africa, in agreement with previous studies. We find that reports of shrub encroachment show a strong association with observed greening. Open shrublands contain the most and strongest magnitude trends, while savannas exhibit both greening and browning. Due to the extensive increases in woody plant cover in this region, we suggest that ongoing increase in atmospheric CO 2 concentrations may play a primary role in promoting woody plant encroachment.