Land cover and precipitation controls over long-term trends in carbon gains in the grassland biome of South America

Carbongainsareakeyaspectofecosystemfunctioningsince theyrepresent theenergyavailable for upper trophic levels. Carbon gains (or primary production) are strongly correlated with other ecosystem attributes such as secondary production and they are also the support for the provision of many ecosystem services. Given the documented dependency of primary production on precipitation, we expect that altered precipitation regimes, such as those projected by climatemodels,will have a significant impact on carbon gains. Landuse and land cover changes are also expected tohave a significant impact on thedynamics of carbongains. Wegenerated a spectral database of the fraction of photosynthetically active radiation interceptedbyvegetation (fPAR), inorder tostudylong-termtrends(i.e.,decades)incarbongainsanditsspatialandtemporalrelationships with precipitation and land cover patterns inUruguay,which is part of theRio de la PlataGrasslands, one of the largest temperate grasslands biome of the world. We found that carbon gains of native forests and grassland afforestationexhibited thestrongestpositive spatial response toprecipitation,whereas cropsandrangelands the weakest. In addition, we found that the temporal response of carbon gains to precipitation was strong and positiveforall landuses.Althoughtherewerenotclear trends inprecipitation,wefoundstrongnegativetrendsin carbongains throughtime,particularly inrangelandsof the ‘‘NorthernCampos’’ofUruguay,where these trends represent a decrease between 10% and 25% of the annual aboveground net primary production. On the other hand, positive trends in carbongains through timewere associated to grasslandafforestation andnative forests. Therefore,duringtheperiodanalyzed, landcoverhadastrongerinfluenceontheobservedtrends incarbongains than precipitation. These patterns emerged as a consequence of the interaction among precipitation, temperature, edaphic factors and management. Present trends in the controlling factors of C gains would exacerbate the observed patterns with serious consequences for the provision of ecosystems services.


INTRODUCTION
Global climate models predict an increased incidence of extreme weather events, such as droughts and floods and to a lesser extent, in the mean annual temperature and precipitation (Ra ¨isa ¨nen 2002, IPCC 2013).Given the documented regional relationship between precipitation and primary production (Sala et al. 1988, Lauenroth and Sala 1992, Paruelo et al. 1999), we expect that grasslands and shrublands, from semi-arid to sub-humid, be among the most affected biomes (Hsu et al. 2012).The effects of these projected changes on carbon gains are expected to be also lagged as vegetation, particularly in these systems, shows inertial dynamics associated to its dependence on water resources and soil characteristics (Oesterheld et al. 2001, Wiegand et al. 2004, Fabricante et al. 2009, Sala et al. 2012).Therefore monitoring of C gains is a first step to understand the response of vegetation to altered atmospheric composition, climate or land use/land cover.
Among the techniques used to estimate carbon gains and its dynamics, those derived from remotely sensed data have been shown particularly appropriate, given their spatial (from regional to global) and temporal (potentially from 1980) coverage.Moreover, the availability and quality of remote sensing data has increased dramatically over the past 20 years.Other related techniques, such as the empirical up-scaling of eddy covariance measurements (Baldocchi 2008, Jung et al. 2009) constitute a new and powerful alternative for estimating carbon, water and energy fluxes between the biosphere and the atmosphere.However, the absence of eddy covariance towers in some regions of the world (southern South America, our focal region, Middle East, India and central Africa) restricts its use, particularly in these regions (see Baldocchi 2008 and http://fluxnet.ornl.gov/maps-graphics).
Vegetation indices have been widely used to study vegetation dynamics in response to climate and human activities (e.g., Tucker et al. 2001, Baldi et al. 2008, Beck et al. 2011, Pettorelli 2013).However the virtual absence of continuous and congruent databases, in terms of radiometric, spatial and temporal resolution, limits its use for mid-to-long term studies.Some recent advances have been made in order to inter-calibrate vegetation indices and fully exploit their potential in analyses of mid and long-term (see for example Mao et al. 2012).The requirement of long, continuous and congruent records arises from the urgent need of reliable estimates of trends in carbon uptake at regional, continental and global scales, as we need to assess the role of countries, regions and biomes as net sources or sinks of carbon.
The Uruguayan territory is part of the Rio de la Plata grasslands (Fig. 1), one of the largest areas of temperate grasslands of the world (Soriano et al. 1992).The Rio de la Plata grasslands occupy more than 70 3 10 6 ha in southern South America, including the Campos in Uruguay and southern Brazil and the Pampas in Argentina.In Uruguay, natural grasslands, though modified by grazing, still cover most of the area of the country.These grasslands are co-dominated by C 3 and C 4 grasses, whereas woody species are by far less important than the herbaceous component (Altesor et al. 2006).The Rio de la Plata grasslands, and the Uruguayan Campos in particular, have experienced profound changes mainly associated to the expansion of agriculture and afforestation (Jobbagy et al. 2006, Paruelo et al. 2006, Vega et al. 2009).
Global trends in carbon gains have been described in several studies, using satellite platforms with diverse spatial resolutions and temporal coverage, but to our knowledge, none of them spanning more than 20 years.Using spectral information from Global Inventory Monitoring and Modeling Studies (GIMMS) and Pathfinder Advanced Very High Resolution Radiometer Land (PAL), Nemani et al. (2003) found a global increase in net primary production (NPP) for the period 1982-1999.According to these authors 25% of the Earths vegetated land surface exhibited significant increases in NPP.The Rı ´o de la Plata grasslands were one of the regions that exhibited significant increments in C gains.However Fensholt and Proud (2012) found that for a different period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) and depending on the dataset used, negative (normalized difference vegetation index [NDVI] from GIMMS-3g) or non significant (NDVI from MODIS) trends in C gains were prevailing in Uruguay.Beck et al. (2011) found contrasting results in C gains in Uruguay using four normalized differv www.esajournals.orgence vegetation index datasets derived from the Advanced Very High Resolution Radiometer (AVHRR) namely PAL, GIMMS, LTDR Version 3 (Land Long Term Data Record) and FASIR (Fourier-Adjustment, Solar zenith angle corrected, Interpolated Reconstructed).The degree of consistency among trends derived from these data sets was moderate to poor.For the Rio de la Plata grasslands, inconsistencies were greater in southern Uruguay and in the humid Pampa (Argentina).Increases in the fraction of photosynthetically active radiation intercepted by vegetation (fPAR) derived from PAL data sets were documented for South America by Paruelo et al. (2004).This study showed that in Uruguay the incidence of significant trends (either positive or negative), was minimal.In a continental comparison of AVHRR derived datasets (PAL, GIMMS and FASIR), Baldi et al. (2008) found similar results to those reported by Beck et al. (2011), i.e., continental increases in C gains but with moderate to poor levels of consistency among trends from different datasets and low occurrence of significant trends in our focal Although there is spectral information from sensors onboard different satellite platforms since 1980, currently we lack a reliable, continuous, congruent and unified record of vegetation indices (as NDVI) in terms of temporal, spatial and spectral resolution.A strong candidate to study trends in ecosystem functioning is the GIMMS dataset, which covers the period 1981-2007. However Baldi et al. (2008) and Alcaraz-Segura et al. (2010) showed that this dataset had serious problems to capture changes in carbon gains associated to known land cover changes.Therefore, combining and inter-calibrating datasets from different sensors and periods would be particularly relevant to achieve more accurate and consistent estimates of mid and long term C gains for focal regions (i.e., Uruguay) and worldwide.
Our first goal was to generate a continuous and congruent NDVI and fPAR based C gains estimates database for Uruguay from 1981 to 2011, by means of the splicing of spectral information from two satellite generated databases: Moderate Resolution Imaging Spectroradiometer (MODIS) on board NASA's Terra satellite (from 2000 to 2011) and Land Long Term Data Record (LTDR) on board NOAA's satellites (from 1981 to 1999).The NDVI is a spectral index related to fPAR (400-700 nm) (Sellers et al. 1992), the main determinant of NPP (Monteith 1972(Monteith , 1981)).
With this database at hand, our objectives were 1.To estimate the spatial relationship between mean annual carbon gains and mean annual precipitation for the main land uses covers in Uruguay (''spatial models'' based on multiyear averages on different locations, as proposed by Sala et al. 1988, 2012, Paruelo et al. 1999).2. To estimate trends in C gains from 1981 to 2011. 3. To estimate the temporal relationships between carbon gains and precipitation (''temporal models'' based on year-to-year change at single locations; Sala et al. 1988, 2012, Paruelo et al. 1999).4. To establish if temporal trends in carbon gains and it relationships with precipitation vary among different land covers.

Spectral data and processing
Image splicing.-Toobtain a continuous series of carbon gains estimates (i.e., NDVI and fPAR) we spliced spectral information from two satellite products with different characteristics.The LTDR product (Long Term Data Record) is generated from data collected by the AVHRR sensor on board of different NOAA satellites (Pedelty et al. 2007).This product combines a high temporal resolution (daily) with a moderate spatial resolution (0.058 or 2500 ha at the equator).LTDR NDVI is derived from channels 1 (visible, 580-680 nm) and 2 (near infrared, 725-1100 nm) of the AVHRR sensor.This product also provides associated quality information (''quality flags'' in the LTDR terminology) that allows estimating the usefulness of the vegetation indices.The NDVI MOD13Q1 product is derived from data of bands 4 (red, 545-565 nm) and 5 (infrared, 1230-1250 nm) of the Terra sensor on board of the NASA earth observing system (EOS-NASA) MODIS (''Moderate Resolution Imaging Spectroradiometer'').This sensor combines a high spatial resolution (250 m) with a moderate temporal resolution (fortnightly), and as LTDR platform, provides additional quality information.
In order to obtain time series of vegetation indices with monthly temporal resolution we applied the algorithm presented in Fig. 2. In the case of LTDR, given that this product provides a daily image of global coverage, we first clipped the image to the study area, and then applied the 15-day maximum value composites (MVC).This technique selects as representative of the fortnight, the maximum daily NDVI value of the period.Then the monthly value was obtained as a weighted average of the fortnightly values.The MODIS product provides an image every 16 days, for a grid of cells that encompass the whole earth surface.From this grid, we select the cell (''tile,'' in the MODIS terminology) that covers Uruguay.In order to obtain monthly values of NDVI from MODIS, we applied a weighted average on the fortnight NDVI values as done with the LTDR product.Once we obtained the monthly time series of NDVI from each satellite product, we proceeded to the degradation of the MODIS monthly values to the LTDR spatial resolution.This was done by averaging the 400 MODIS pixels encompassed by each LTDR pixel.The NDVI series obtained are congruent in terms of temporal (monthly) and spatial (2500 ha) resolutions.
Given the differences in spatial, temporal and spectral resolution, we expect differences of scale between monthly NDVI values from MODIS and LTDR (i.e., differences in the mean and range of variation between NDVI from spatially degraded MODIS and LTDR for a given spatial location).In order to establish the existence of these inconsis-tencies between the NDVI series, we compared the mean and standard deviation of LTDR and MODIS derived vegetation indices corresponding to each LTDR pixel (2500 has).As each LTDR pixel consists in approximately 400 MODIS pixels, a monthly value of NDVI derived from MODIS corresponds to the mean of the 400 values that fall inside the area corresponding to an LTDR pixel.With this information, we calculated and compared the mean and standard deviation of the NDVI from MODIS and LTDR.In order to compare the means of the series we v www.esajournals.orgcalculated the difference of means as l NDVI MODIS À l NDVI LTDR .For comparing scales of variability, we calculated the quotient between the standard deviations of the series as r NDVI MODIS =r NDVI LTDR .
In the case of differences between means, the minimum and maximum values were 0.501 (l NDVI MODIS .l NDVI LTDR ) and 0.313 (l NDVI MODIS , l NDVI LTDR ), respectively.The differences in the means as percent of the mean for indices from both satellite platforms for each pixel were expressed as 100 3 ((l NDVI MODIS À l NDVI LTDR )/ l NDVI MODIS-LTDR ).We found that 97.2% of the pixels (6796) showed percent differences in the range of (À5%, þ5%) of the general mean, whereas 97.7% (6833) of the pixels exhibited differences in the range (À10%, þ10%).For the case of the quotient between standard deviations, the extreme values were 0.002 (which implies that r NDVI LTDR is 500 times bigger than r NDVI MODIS ) and 1.66 (which implies that r NDVI MODIS is more than a time and a half bigger than r NDVI LTDR ).The quotients were categorized in eight intervals, r NDVI MODIS / r NDVI LTDR 0.25 (r NDVI LTDR 4r NDVI MODIS ) to r NDVI MODIS /r NDVI LTDR ! 4 (4r NDVI LTDR r NDVI MODIS ).There was 86.2% of the pixels (6029) in the classes 0:67 r NDVI MODIS =r NDVI LTDR 1:5, indicating that most of the series exhibited a variability around a half or double of the other platform series variability.Although these checks showed that the congruence of the NDVI series from LTDR and MODIS platforms seems to be fairly good, we re-scaled the LTDR NDVI to make its variance and mean identical to those from MODIS NDVI, given that technical improvements and quality assessments given by this platform, make it more reliable (Fensholt et al. 2009).
We re-scaled the LTDR derived NDVI series with the standard deviation of the MODIS NDVI from all the NDVI values and pixels encompassing Uruguay (excluding pixels from urban areas, coastlines and water bodies), r NDVI MODIS ALL (0.0743), according to where NDVI M LTDR is the value of the NDVI from LTDR re-scaled according to the standard devi-ation of the NDVI from MODIS platform and z LTDR is the standardised value of NDVI derived from LTDR platform (considering the mean and standard deviation of LTDR NDVI values from all months and pixels from Uruguay, excluding those pixels from urban areas, coastlines and wa ter bodies, l NDVI LTDR ALL ¼ 0.5807 and r NDVI LTDR ALL ¼ 0.0780, respectively).With the values of NDVI M LTDR and NDVI MODIS (both series in MODIS scale of variation) we proceeded to calibrate the relationship between NDVI and fPAR with the parameters estimated by Grigera et al. (2007;see below).The rescaling process was done in this way (with the mean and standard deviation from all NDVI monthly values from all pixels), and not pixel to pixel (with the mean and standard deviation from all NDVI monthly values from each pixel), in order to eliminate the differences in indices that could have arisen as an artefact of differences in sensors and, at the same time, not to lose the capacity to detect the effect on NDVI of possible changes in land use of pixels between 1999 and 2000 (the last year of LTDR and the first year of MODIS).Other rescaling options were tested.For example, scaling on a per pixel basis, or considering mean and standard deviation values from those pixels without trends in any of the NDVI series.This alternative gave results almost identical to those presented below.

Calibration of the relationship fPAR-NDVI
To estimate fPAR based on NDVI, we used an empirical approximation that assumes a nonlinear relationship between MODIS-NDVI and fPAR (Los et al. 2000, Piñeiro et al. 2006, Grigera et al. 2007).The non-linear relation between NDVI and fPAR accounts for the widely described saturation of NDVI at high Leaf Area Index (LAI) .3, and implies a linear relation between the simple ratio index (SR and fPAR.We parameterized the relation between NDVI and fPAR with data from Rio de la Plata grasslands (Grigera et al. 2007) assigning no absorption (fPAR ¼ 0) to NDVI values corresponding to pixels that had no green vegetation (bare soil) and maximum fPAR (fPAR ¼ 0.95) to NDVI values corresponding to pixels with high amount of green biomass (sown pastures with LAI .3 and high yielding wheat crops during anthesis, Grigera et al. 2007)

Land cover characterization for Uruguay
We generated a land cover characterization for Uruguay using MODIS high temporal resolution images and field observations, following the methodology proposed by Volante and Paruelo (in press) and Baeza et al. (2014).Six land cover classes were defined: rangelands (forage classes, including natural grasslands and perennial sown pastures), tree plantations and forests (and Pinus and Eucalyptus plantations and native forests), summer crops (maize, soybean, sunflower and sorghum), winter-summer crops (wheat-soybean double crop), water bodies and urban areas (Fig. 1).Based on this characterization, we estimated for each of the 6993 LTDR pixels that encompass the whole country the proportion of the pixel occupied by each land cover classes.
MODIS products were used to generate the land cover characterization.The Enhanced Vegetation Index (EVI) from MODIS product MOD13Q1, was used as a surrogate of the fPAR that is absorbed by green tissues (Sellers et al. 1992) and, hence, a key determinant of primary production (Monteith, 1981, Running et al. 2000).The spatial and temporal resolutions were 250 m and 16 days respectively.EVI temporal dynamics for a given pixel over a growing season was called ''phenological signature.''Each ''phenological signature'' described the seasonal dynamics of radiation interception by each land cover type.In order to do this, it was necessary to record the occurrence of land covers in the field (4309 ground truth points, taken during 2011 and 2012; Appendix A).Water bodies and urban areas were identified using MODIS product MCD12Q1.With the collected data along the whole territory, we built a geodatabase as a ''spectral library of phenological signatures'' including temporal and spatial variability within each land cover type and between different land cover associated with management practices, climate, topography, etc.Then, knowing the ''phenological signature'' of each land cover observed in field allowed us to determine what cover may be associated with the corresponding ''spectral signature'' and therefore identify the coverage for the whole territory.For that, we used See5 software to define the classification rules expressed as decision trees.The map generated showed the known spatial patterns of main land cover types in Uruguay, where tree plantation and forest land cover classes are present in all country except in the ''Cuesta Basa ´ltica'' region, crop areas dominate the west and east regions and the rest of the country is mainly occupied by rangelands (Fig. 1).The map has a global precision of almost 90%.Excluding water bodies and urban areas, precision associated to each land cover was higher than 85% (rangelands 90%, forests and tree plantations 97%, summer crops 86% and winter-summer crops 87%; Appendix A).Precision values above 0.8 are indicative of good to very good classification performance (Congalton 1991).

Climate data
We obtained precipitation records for the whole Uruguay from the TRMM mission product 3b43 (http://trmm.gsfc.nasa.gov/3b43.html).This product provides mean monthly precipitation rates in mm/h at a spatial resolution of 0.25 3 0.258, and from 1998 to 2011 with a monthly temporal resolution.The precipitation rates were converted to total monthly precipitation.In order to evaluate TRMM precipitation estimates we obtained monthly total ground precipitation measures from a network of 21 meteorological stations from Uruguay that cover the period from 1948 to 2011 (Direccio ´n Nacional de Meteorologı ´a and Instituto Nacional de Investigacio ´n Agropecuaria, INIA).The precipitation records from these meteorological stations were then correlated with the corresponding TRMM monthly precipitation estimates.All stations except one, exhibited correlations higher than 0.75, and all of them where highly significant (Appendix B).Mean correlation (standard error) was r ¼ 0.847 (SE r ¼ 0.034).Most of the stations (18-21) did not show significant differences in the mean of TRMM precipitation estimates and ground precipitation series.

Statistical analyses
Spatial relationship between fPAR and precipitation.-Thisrelationship was explored regressing the mean annual integral of fPAR against mean annual precipitation for each pixel with more than 75% of its surface covered by one of the three dominant land uses/covers: rangelands, crops and afforestation and native forests.In order to consider spatial autocorrelation, after fitting simple linear regression, we estimated the empirical residual semivariogram and fitted different theoretical variogram models to it (Beale et al. 2010).The best residual variogram model was selected by means of the weighted sum of squared residuals and then we included it in the error term of a new linear regression fitted by means of generalized least squares (Dormann et al. 2007).Finally for each land cover, the models with and without spatial autocorrelation were compared by means of Akaike's information criterion (AIC; Burnham and Anderson 2002) and the best models (those with the lowest AIC) were selected.
Temporal trends in precipitation (1998-2011) and fPAR .-Trends in TRMM precipitation estimates and fPAR were studied considering monthly and yearly temporal resolutions, and by both parametric and non-parametric statistical methods.For the monthly time series, in order to consider seasonal variability we analysed trends by means of linear regression with additive seasonal indicator variables (Cowpertwait and Metcalfe 2009) and then by means of the seasonal Kendall-Theil trend estimator (Hipel andMcLeod 1994, Helsel andHirsch 2002).For the yearly time series we estimated temporal trends by means of simple linear regression and Kendall-Theil regression.These approaches are complementary and were used together to add robustness to the analysis.Unlike linear regression, Kendall-Theil regression makes few assumptions about the underlying relationship between the response variable and the predictor, and is insensitive to the presence of outliers (Helsel and Hirsch 2002).The agreement between these two techniques, for each temporal resolution and for all pixels that encompass Uruguay, was assessed.
Temporal trends in precipitation use efficiency .-Weestimated trends in precipitation use efficiency (PUE; Vero ´n et al. 2005) for a subset of 12 ground meteorological stations with precipitation records from 1948 to 2011 and fPAR estimates from 1981 to 2011.PUE was estimated as the quotient between annually integrated fPAR (I-fPAR) and mean annual precipitation (MAP).In order to avoid precipitation carry-over effects on I-fPAR, we also estimated trends in 3 and 5 year integrated PUE.We applied simple linear regression (Draper and Smith 1981) and segmented linear regression (Muggeo 2008) in order to explore the existence of breakpoints in the dynamics of PUE.Best fit models were selected by means of Akaike's information criteria (AIC; Burnham and Anderson 2002) Patterns associated to land cover.-Inorder to establish if trends in C gains, and the correlations with precipitation were related to the main land uses in Uruguay (rangelands, afforestation/native forest and crops), we first selected the pixels with more than 75% of the area occupied by one of these land covers and at the same time presented significant trends or correlations with precipitation.We then compared the values of trends and correlations by means of one way analyses of variance with land use as factor.Given that ANOVAs were highly unbalanced (for example, there were 2075 pixels occupied mainly by natural grasslands that exhibited significant monthly trends, 136 by crops and 130 by tree plantations/native forests), we sampled randomly with replacement from each land cover, a number of pixels equal to the number of pixels in the least represented class (130 in the previous example), and then performed the analysis of variance.This process was repeated 5000 times, and then the percentage of significant results was obtained.

Spatial relationship between fPAR and precipitation
The strongest spatial relationship between mean annual integrate of fPAR (I-fPAR) and mean annual precipitation (MAP, the ''spatial v www.esajournals.orgv www.esajournals.orgmodel'') was found in afforestation/native forests, regardless of whether we consider or not spatial autocorrelation among pixels (Fig. 3 and Table 1).According to our classification, tree plantations and native forests were dominant in 142 pixels (more than 3500 km 2 , 2% of the countries area).The models without considering spatial autocorrelation were almost identical for rangelands (3139 pixels or more than 78000 km 2 , almost a half of the countries area) and crop classes (250 pixels or more than 6200 km 2 , or 3.5% of the countries area) but crops presented a higher slope when we included spatial autocorrelation.For all land cover classes considered, the statistical models including spatial autocorrelation were significant (Table 1), implying that some degree of spatial dependence was present in the relation between mean I-fPAR and MAP.

Temporal trends in precipitation
Only negative significant trends were detected in TRMM precipitation estimates, which affects an area of around 57.700 km 2 , mainly in Northern Uruguay (Fig. 4).These trends range from decreases in 2 mm/month to around 6 mm/ month.The percentage of agreement between results from seasonal Kendall-Theil regression and linear regression with seasonal indicators was modest (71.9%) indicating that in more than 4450 pixels both techniques detected the same temporal pattern.We did not find any significant trend in the annual accumulated TRMM precipitation estimates time series.This apparent inconsistence could be the result of very different sample size in the monthly resolution time series respect to the annual resolution time series (168 data values in the monthly time series vs. 14 data

Temporal trends in fPAR
Almost 60% of the area of the country exhibited significant trends in fPAR during the period 1981-2011 at a monthly temporal resolution (almost 4000 LTDR pixels or 100000 km 2 , Fig. 5a, c).Most of these significant trends were negative (2673 LTDR pixels or 66800 km 2 ) and almost perfectly coincident with the Cuesta Basa ´ltica geomorphologic unit, which is covered mainly by natural grasslands and devoted to cattle and sheep grazing.Mean observed trend (standard error) was À0.00215 (SE ¼ 1.53 3 10 À5 ).On the other hand, positive trends encompass an area of more than 30000 km 2 (1219 LTDR pixels) and were associated to tree plantations and forests in the North of the country and crops in the South.Mean observed trend (standard error) was 0.00272 (SE ¼ 1.72 3 10 À6 ).The percentage of agreement between trends estimated by means of seasonal Kendall-Theil regression and by linear regression with seasonal indicator variables was 92.8%, indicating that in more than 5780 pixels both techniques detected the same temporal pattern.
Trends in annual integral of fPAR were less abundant, and represented 25% of the countries area (.44000 km 2 or 1775 LTDR pixels; Fig. 5b, d).
As with the monthly fPAR trends, the spatial distribution of annual integrated C gains resembles closely the main land covers in Uruguay, i.e., rangelands exhibiting negative trends (n ¼ 1191, mean trend ¼ À0.00284, SE ¼ 6.2 3 10 À7 ) and

Temporal correlation between fPAR and precipitation
Significant positive correlations between monthly fPAR and monthly TRMM precipitation estimates were present in most of the Uruguayan territory (more than 140.000 km 2 , or 5704 LTDR pixels; Fig. 6a, c).Mean monthly correlation was r ¼ 0.3062 (SE r ¼ 0.0011) and ranged from 0.1563 to 0.5406.Positive significant correlations be-tween I-fPAR and accumulated annual TRMM precipitation estimates were higher than correlations at a monthly resolution, but were less extended in space (3303 LTDR pixels, or 82500 km 2 , almost a half of area of the country, Fig. 6b, d).Mean annual correlation was r ¼ 0.6978 (SE r ¼ 0.0014) and ranged from 0.5760 to 0.9442.At both temporal resolutions higher correlations were more frequent in northern Uruguay.

Temporal trends in precipitation use efficiency
According to AIC criterion, in 10 of the 12 sites analyzed, the most parsimonious models were linear (6 sites with negative trends and 4 sites with positive trends; Tables 2 and 3).Only one of these 10 sites exhibited a significant (a ¼ 0.05, positive) Fig. 6.Spatial (maps) and frequency (graphs) distribution of significant correlation beetwen monthly fPAR values and monthly precipitation estimates from TRMM (a, c) and significant correlations beetwen annual integral of fPAR and annual precipitation estimates from TRMM (b, d).
v www.esajournals.orglinear trend in PUE (''Estanzuela'' site in the SW of the country; Fig. 7).When we repeated the analyses for 3 and 5 year resolution, none of the sites exhibited significant trends (Appendix C).

Patterns associated to land cover
When analyzing fPAR trends and correlations with precipitation for the main land covers of Uruguay, some clear patterns emerged (Fig. 8).First of all, positive trends in C gains, as estimated by fPAR, were more common in afforestation/native forests, whereas negative trends were more frequent in rangelands and to a lesser extent in crops (Table 4).This pattern seems to be independent of the temporal resolution considered.Analyses of variance performed on trends from both temporal resolutions were highly significant (F 3,2338 ¼ 767.1, p , 0.001 for monthly resolution and F 3,1099 ¼ 441.6, p , 0.001 for annual resolution).Resampled ANOVAs reinforced these results: all replicates (5000) yielded significant results for both temporal resolutions.
When we examined the effect of land cover on correlations between precipitation and fPAR we found a slight trend to higher monthly correlations in grasslands and to higher annual correlations in grasslands and afforestation/native forests (Fig. 9).Results from ANOVA confirm these trends, showing that differences although weak, were significant (F 3,3265 ¼ 14938, p , 0.001 for monthly resolution and F 3,2034 ¼ 54271, p , 0.001 for annual resolution).As in the fPAR trends, all resampled ANOVAs yielded significant results.
Table 2. Results of the regression models of precipitation use efficiency against time for ten sites where the most parsimonious model was linear.b 1 represents the slope of the simple linear regression model, p, the p value associated and AIC, Akaike's information criterion.Latitude (Lat), longitude (long), mean of precipitation use efficiency during the period (mean PUE) and precipitation marginal response (PMR, estimated as the slope of (annually integrated fPAR) 3 1000 vs. mean annual precipitation) of sites are shown.Dominant land cover is also shown: R stands for rangeland, C stands for crops and TP-F for tree plantations and forests.In sites where there was not a clear dominance, the order of the two or three dominant land covers is shown.Significance codes are (.) for p , 0.05, (*) for p , 0.01 and (**) for p , 0.001.Table 3. Results of the regression models of precipitation use efficiency against time for two sites where the most parsimonious model was segmented linear.b 1fs represents the slope of the relationship between PUE and time in the first segment or period, whereas b 1ss represents the slope in the second period, BY represents the breakpoint year, p the represents p value, and AIC, Akaike's information criterion.Latitude (Lat), longitude (long), mean of precipitation use efficiency during the period (mean PUE) and precipitation marginal response (PMR, estimated as the slope of (annually integrated fPAR) 3 1000 vs. mean annual precipitation) of sites are shown.Dominant land cover is also shown: R stands for rangeland, C stands for crops and TP-F for tree plantations and forests.In sites where there was not a clear dominance, the order of the two or three dominant land covers is shown.Significance codes are (.) for p , 0.05, (*) for p , 0.01 and (**) for p , 0.001.

DISCUSSION
Our results describes for the first time, and making use of long term data (30 years), the relationship between C gains (using fPAR) and climate (precipitation) for the main land cover/ uses in a significant portion of the Rio de la Plata grasslands (Uruguay).The most remarkable finding is the prevalence of negative trends in C gains.When all trends (i.e., significant and non significant) are considered, around 65% of the area of the country exhibited decreases in monthly and annual integrated C gains.When only significant trends are considered, these percentages fall to 42.8% and 19.1% for monthly and annual trends, respectively.The decrease in fPAR (our estimator of C gains) was not uniformly distributed across the country, instead it was strongly associated to rangelands and to lesser extent, to crops.These land covers also exhibited a similar and weak spatial response to precipitation (Fig. 3, Table 1).On the other hand, afforestation and native forests exhibited the strongest, although weak, spatial fPAR response to precipitation among the land covers considered.This land cover was also associated to increases in fPAR during the period.The weak spatial relationships found could be a consequence of the relatively narrow mean annual precipitation gradient characterizing Uruguay (from 1000 mm to 1800 mm of accumulated annual precipitation), as compared with other studies estimating spatial relations between C gains and precipitation (from 200 mm to 3000 mm of mean annual precipitation; Sala et al. 1988, Huxman et al. 2004, Sala et al. 2012).
Positive temporal trends in fPAR associated to tree plantations and forests reveal a strong human control on ecosystem functioning.Tree plantations expansion was associated to government plans that promoted forestry at the beginning of 1990 decade which resulted in an accelerated increase in accumulated afforested area, mainly in Northern Uruguay (Torres Rojo  v www.esajournals.organd Fossati 2004).Our results also suggest that tree plantations (mostly monocultures of eucalyptus and pinus) have a greater spatial response to precipitation than rangelands and crops (Fig. 3, Table 1).Tree plantations showed a high increase in NPP in response to spatial variability in precipitation while in contrast, highly diverse natural grasslands or sown pastures did not increase their NPP, suggesting other limitations for NPP in these land covers.It also has been shown that tree plantations exhibited higher annual averages and lower seasonal ranges in light interception and delayed growing seasons (Vassallo et al. 2012).
Negative trends in C gains were observed mainly in the ''Cuesta Basa ´ltica'' region which is occupied by natural grasslands devoted mainly to extensive livestock and North Western Sedimentary basin (''Cuenca Sedimentaria del Noroeste'') in which crops and natural and implanted grasslands are dominant.Mean observed trend in fPAR (À0.00215 fPAR units/ month) represents important reductions in the productivity of these grasslands.Using average values of incoming photosynthetically active radiation (PAR) for the region and radiation use efficiency (RUE, or the ratio between primary production and the radiation intercepted by terrestrial vegetation)) from published calibrations for Rio de la Plata Grasslands (Piñeiro et al. 2006), this change in fPAR represents an aboveground net primary production (ANPP) reduction of 0.2 grDM .m À2. month À1 , or more than 700 Kg .ha À1 for the whole period analyzed .This ANPP value represents roughly between 10 and 25% of the annual ANPP of grasslands in the ''Northern Campos'' of Uruguay according to the values reported by Paruelo et al. (2010) and Baeza et al. (2010).It has been shown that RUE derived from eddy covariance C flux measurements is strongly associated to water stress and temperature for a series of sites along the major biomes of the world (Garbulsky et al. 2010).Specifically, RUE in temperate forests and savannas was negatively related to water stress (as assessed by the ratio of potential evapotranspiration/actual evapotranspiration), positively related to temperature in temperate forests and negatively related to temperature in grasslands and savannas.This could represent even greater reductions in ANPP, given the observed trends in monthly precipitation and temperature anomalies.
Several possible explanations arise for the prevalence of negative trends in C gains estimations.First, an increment in the livestock stocking rate during the period could have resulted in the degradation of rangelands, resulting from a reductions in plant cover and in the density of green tissues.Hilker et al. (2014) showed similar decreases in vegetation greenness of mongolian rangelands for the period 2002-2012, that were consequence of the dramatic increase in stocking rates observed during that period.However v www.esajournals.orgstock records for part of the period for the ''Cuesta Basaltica'' and North ''Western Sedimentary Basin'' geomorphological units, showed negative trends between 1990 and 2010 (Oyhantcabal 2014).These negative trends in animal stocks were accompanied by negative trends in the sheep/cow ratios.This could have resulted in an increment in standing senescent plant material and less intercepted radiation (Dibella et al. 2004).The observed reductions in livestock numbers were also accompanied by a reduction in the grazed area due to the transformation of rangelands into croplands (Paruelo et al. 2006) that may result in an increase in the stocking rate on the remaining grasslands.A proper identification of causes and the role of overgrazing on C gains trends would need additional data at finer spatial resolutions.
The second explanation could be associated to the interaction between temperature, soil types and vegetation response to water stress.Positive mean annual temperature anomalies (respect to the period 1981-2010) have been dominant in Uruguay during at least the last twenty years (see http://www.ncdc.noaa.gov/cag/mapping/global).This temperature rise could imply an increment in the atmospheric water demand, effect that could be reinforced if precipitation is reduced (Wu et al. 2011, Vicente-Serrano et al. 2013).Shallow soils in the ''Cuesta Basaltica'' geomorphological unit, could have experienced higher levels of water stress.On the other hand, vegetation on deeper soils, with a higher water holding capacity, as those in Southern Uruguay in general and in the ''North Western Sedimentary Basin'' in particular, may result less affected by water stress.This could result in lowered fPAR and C fixation in the Cuesta Basa ´ltica rangelands as opposed to increased or constant fPAR and C gains in the North Western Sedimentary basin.Projections of global climatological models predict continued increases in temperatures for different plausible scenarios of CO 2 and other greenhouse gasses emissions (IPCC 2013), in particular for our focal region.This could represent a deepening of present observed trends in fPAR and C gains.
Trends in fPAR were not accompanied by trends in precipitation use efficiency.Only one site in Southern Uruguay (site ''Estanzuela'') exhibited a significant increase in PUE for the period 1981-2011 (Fig. 7, Table 2).This result, imply that decreasing C fixation would not be the result of an altered response of vegetation to precipitation.
Some authors (Le Houerou 1984, Prince et al. 1998, 2007, Hein and de Ridder 2006) have hypothesized that degradation process would result in changes in the response of vegetation to precipitation, i.e., a decrease in precipitation use efficiency and precipitation marginal response (PMR, the slope of the regression between annual integrated fPAR and annual precipitation; see Vero ´n and Paruelo 2010).We did not find any pattern in PMR or PUE associated to land cover for the period 1981-2011 (Table 2, Fig. 7; Appendix C).The absence of trends in PUE could be associated to concomitant decreases in fPAR and precipitation, as evidenced by the high correlations between I-fPAR and annual precipitations, for the period 1998-2011, mainly in the ''Cuesta Basa ´ltica'' region (Fig. 6b, d).As discussed in previous paragraphs, this fact suggests observed trends in carbon gains are the result of the interaction between annual precipitation, soils and temperature.Other factors, in particular the interannual variability in the magnitude and distribution of precipitation events, soil nitrogen availability, edaphic characteristics and resource heterogeneity (Reynolds et al. 2004), could be plausible alternative explanations.
Our results highlight the importance of having long time series for detecting current changes in ecosystems.In particular, the splicing of spectral databases allowed us to resolve some inconsistencies in the reported trends that are consequence of the differences in temporal, spatial and spectral resolutions and the partial periods embraced by data from different sources (Nemani et al. 2003, Paruelo et al. 2004, Baldi et al. 2008, Zhao and Running 2010, Beck et al. 2011).
The main patterns detected are: (1) Forest/ afforestation where the principal land cover with positive trends in C gains and with the stronger spatial response to precipitation; (2) Rangelands and croplands presented negative trends in C gains; and (3) All land covers presented high correlations between C gains and precipitation for all the spatial or temporal resolutions considered resolution.These results must be placed in the context of future global change, for example the trends in precipitation (those detected by us in the monthly TRMM series and those predicted by global circulation models), temperature anomalies (see http://www.ncdc.noaa.gov/cag/mapping/global) and land use/cover changes (i.e., an increment in afforested and cultivated area), would deepen, the differences among land uses in terms of trends and spatial response to precipitation.Moreover, these findings point us to the fact that the same precipitation pattern (no trends, or negative trends in monthly precipitation) during the last three decades resulted in different trends in C gains depending mainly on soils depth (shallow in the North of the country vs. deep in the South) and on how humans use and transform the landscape.This could have serious consequences, given that C gains are one of the fundamental supports for the provision of many ecosystem services.

Fig. 1 .
Fig. 1.Land cover classification for Uruguay (left panel,) and geomorphological units (right panel).The maps are shown in WGS84 projection.

Fig. 2 .
Fig. 2. Schematic representation of the algorithm developed to splice NDVI from LTDR and MODIS products.

Fig. 3 .
Fig. 3. Relationship between mean annual integral of fPAR and mean annual accumulated precipitation for the main land uses in Uruguay.Broken lines represent the fit of the model without spatial autocorrelation in residuals, whereas solid lines represent the fit of the model with spatially autocorrelated errors.The maps in the right panel show the number and spatial distribution of pixels associated with these land uses.Limits of geomorphological units are shown in the maps.

Fig. 5 .
Fig. 5. Spatial (maps) and frequency (graphs) distribution of significant trend estimators from Kendall-Theil regression of monthly fPAR vs. time (change in fPAR/month, a, c) and annual integral of fPAR vs. time (change in annual integral of fPAR/year, b, d).Broken vertical lines represent zero.

Fig. 7 .
Fig. 7. Trends in precipitation use efficiency ((I-fPAR/annual precipitation) 3 1000) for the period 1980-2011.Red broken lines represent best fit models (selected by means of AIC).Vertical broken lines represent the period mean whereas horizontal broken lines represent the mean PUE.All models except those for site S9 (Estanzuela, marked with an asterisk) were nonsignificant at a ¼ 0.05.

Fig. 8 .
Fig. 8. Distribution of trend estimators from Kendall-Theil regression for monthly fPAR time series (a) and fPAR annual integral time series (b) as a function of land use: C (crops), F (forests) and R (rangelands).The horizontal broken line represents no change situation, i.e., trend estimator equal to 0.

Fig. 9 .
Fig. 9. Distribution of correlation between monthly fPAR time series and monthly precipitation estimates time series (a) and fPAR annual integral time series and annual accumulated precipitation time series (b) as function of land use: C (crops), F (forests) and R (rangelands).The horizontal broken line represents mean correlation across land uses, (r ¼ 0.3178 for monthly correlation and r ¼ 0.7033 for annual correlation).

Table 1 .
Results of the regressions beetwen fPAR anual integral and mean anual accumulated precipitation for the three main land uses in Uruguay.b 1 represents the slope of the relationship (change in fPAR per unit change in mean annual accumulated fPAR), AIC is Akaike's information criterion, R 2 the coefficient of determination and n the number of pixels.For the models with spatial structure, VM represents the best fit variogram model, and p, the probability associated with the comparison between the models with and without spatial structure.Significant slopes are indicated with an asterisk.
v www.esajournals.orgvalues in the annual time series).

Table 4 .
Number of pixels (and percentage respect total number of pixels in land use) from main land uses exhibiting positive (þ) or negative (À).Signif.represents the total number (percentage) of significant trends in each land use, whereas total represents total number of pixels with .75% of area covered by the corresponding land use.