Drivers of Productivity Trends in Cork Oak Woodlands over the Last 15 Years

Higher biodiversity leads to more productive ecosystems which, in turn, supports more biodiversity. Ongoing global changes affect ecosystem productivity and, therefore, are expected to affect productivity-biodiversity relationships. However, the magnitude of these relationships may be affected by baseline biodiversity and its lifeforms. Cork oak (Quercus suber) woodlands are a highly biodiverse Mediterranean ecosystem managed for cork extraction; as a result of this management cork oak woodlands may have both tree and shrub canopies, just tree and just shrub canopies, and just grasslands. Trees, shrubs, and grasses may respond differently to climatic variables and their combination may, therefore, affect measurements of productivity and the resulting productivity-biodiversity relationships. Here, we asked whether the relationship between productivity and climate is affected by the responses of trees, shrubs, and grasses in cork oak woodlands in Southern Portugal. To answer this question, we linked a 15-year time series of Enhanced Vegetation Index (EVI) derived from Landsat satellites to micrometeorological data to assess the relationship between trends in EVI and climate. Between 2000 and 2013 we observed an overall decrease in EVI. However, EVI increased over cork oaks and decreased over shrublands. EVI trends were strongly positively related to changes in relative humidity and negatively related to temperature. The intra-annual EVI cycle of grasslands and sparse cork oak woodland without understorey (savannah-like ecosystem) had higher variation than the other land-cover types. These results suggest that oaks and shrubs have different responses to changes in water availability, which can be either related to oak physiology, to oaks being either more resilient or having lagged responses to changes in climate, or to the fact that shrublands start senesce earlier than oaks. Our results also suggest that in the future EVI could improve because the rate of increase in minimum EVI is greater than the rate of decrease in maximum EVI, and that this is contingent on management of the shrub understorey as it affects the rate of decrease in maximum EVI. This will be the challenge for the persistence of cork oak woodlands, their associated biodiversity and social-ecological system.


Introduction
One of the key ecological paradigms is the relationship between biodiversity and productivity [1,2].Higher biodiversity leads to more productive ecosystems [3,4], and more productive environments support more biodiversity [5,6].Two mechanisms relate biodiversity to productivity.On one hand, resource acquisition can be done by obtaining the necessary resources along a productivity gradient, commonly referred to as niche complementarity [7].On the other hand, species may improve their understory shrub species.Therefore, we expect no changes in EVI in patches with dense and sparse cork oak forests without understorey (Figure 1c,e).

(a)
We expect a decrease in EVI in dense and sparse cork oak with understorey (Figure 1b,d) because the negative EVI of the shrub layer effect will dominate over the stable effect of the cork oak.This is because the life-cycle of shrubs operates in a seasonal time span while oaks operate in a yearly time span.(b) Sparse cork oak woodlands with shrub understorey (Figure 1d) are more resistant to changes in temperature and humidity than sparse cork oak woodlands without shrub understorey (Figure 1e) because in the absence of oaks, shrub leaves are relatively more resilient to changes in humidity and temperature than annual or perennial grasses.(c) Grasslands show substantially stronger intra-that inter-annual variability in EVI compared to sparse and dense cork oak woodlands.This is because the life-cycle of grasses operates in a shorter time span than that of shrubs and oaks.
Remote Sens. 2016, 8, 486 3 of 13 (a) We expect a decrease in EVI in dense and sparse cork oak with understorey (Figure 1b,d) because the negative EVI of the shrub layer effect will dominate over the stable effect of the cork oak.This is because the life-cycle of shrubs operates in a seasonal time span while oaks operate in a yearly time span.(b) Sparse cork oak woodlands with shrub understorey (Figure 1d) are more resistant to changes in temperature and humidity than sparse cork oak woodlands without shrub understorey (Figure 1e) because in the absence of oaks, shrub leaves are relatively more resilient to changes in humidity and temperature than annual or perennial grasses.(c) Grasslands show substantially stronger intra-that inter-annual variability in EVI compared to sparse and dense cork oak woodlands.This is because the life-cycle of grasses operates in a shorter time span than that of shrubs and oaks.
Topography is moderate, with gentle slopes and low altitude (159 to 238 m).Climate is Mediterranean with an Atlantic influence, mean annual precipitation of 500 mm, and maximum average daily temperatures ranging from 5 °C to 35 °C [46].The area holds a high richness of species with varying conservation status (Table 1; [47,48]).
Topography is moderate, with gentle slopes and low altitude (159 to 238 m).Climate is Mediterranean with an Atlantic influence, mean annual precipitation of 500 mm, and maximum average daily temperatures ranging from 5 ˝C to 35 ˝C [46].The area holds a high richness of species with varying conservation status (Table 1; [47,48]).We used a time series of EVI to describe the phenological evolution in our study area.EVI is particularly well suited for cork oak woodlands, because it is generally less affected by the saturation of dense canopy covers [49], which are characteristic of cork oak woodlands.Further, EVI has been particularly useful for ecosystem modeling [50], monitoring vegetation phenology [51], and assessing vegetation responses to droughts [52].We generated our time series from all available Landsat imagery between the years 2000 and 2013 for the two footprints that cover our study area (Path/Row 204/33, 204/34), resulting in a time series of 296 images.We pre-processed all images by converting digital numbers into surface reflectance values using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS; [53]), and (b) masking all clouds and cloud shadows using Fmask [54].We chose this approach because for calculating EVI we need surface reflectance and we chose to use a standard method across all the samples.We then calculated the EVI for each image (Equation (1)): where ρx is the surface reflectance for band ˆ(NIR = shortwave infrared band, R = red band, B = blue band), L is the canopy background adjustment, C1 and C2 are the coefficients of the aerosol resistance term, and G is a gain or scaling factor [31].

Land-Cover
We selected only land-cover types of interest to assess EVI trends over stable land-cover types.To do so, we used an existing land-cover map and selected only areas which did not show any land-cover change.The map was based on visual interpretation of a 1 m spatial resolution aerial photograph of the year 1995 [46], and distinguished nine land-cover types within our study area: dense cork oak woodlands with understory (DCoW), dense cork oak woodlands without understory (DCoWt), sparse cork oak woodlands with understory (SCoW), sparse cork oak woodlands without understory (SCoWt), olive yards and orchards (O), riparian vegetation (RV), grasslands (G), eucalyptus plantations (E), and human settlements (S; Figure 2).We chose to focus on five of these land-cover types to test our hypotheses: (i) dense and sparse cork oak without understorey (DCoWt, SCoWt) because these areas represent the response of cork oaks; (ii) dense cork oak with understorey (DCoW) because it represents the combined response of oaks and shrubs; (iii) sparse oak with understorey (SCoW) because it represents the response of shrubs (mostly Cistus spp.); and (iv) grasslands (G) that represent the response of the herbaceous layer.We used data from this station aggregated to daily averages.We chose to analyze only the effect of temperature (T), relative humidity (RH), and heat index (HI) on the EVI trend because these weather parameters are more tightly linked to productivity than the parameters directly measuring wind or metrics that depend on wind.Due to malfunctioning, there were data gaps between the end of March 2006 and May 2007, December 2008 and October 2009, March and July 2011, and April 2012 to February 2014.Unfortunately, it was not possible to complement these data gaps with other meteorological data from the surroundings because the closest station was in pine plantations.

Data Analysis
We based our analysis on a sub-region of ~1.8 km 2 in the immediate vicinity of the weather station so that we could test effects of local weather on local productivity measurements.Further, this close match between the scales of the measurements ensured that the potential relationships are interpretable at the local level, before scaling up to more general patterns for larger areas.This choice reduced the number of available data points, however, still within what is considered a sufficient sample size if significant differences are to be detected and we do not think that it hinders the relevance of the findings.These locations were randomly chosen across all land-cover types, and the random selection provided proportional representation of each land-cover.We then extracted the aggregated weather data for each day in our time series for and merged this information to the EVI values of our Landsat time series.This resulted in a time series containing, overall, 107 observations (i.e., time points) between 2000 and 2015, which served as the input for our analysis (Table 2).We used data from this station aggregated to daily averages.We chose to analyze only the effect of temperature (T), relative humidity (RH), and heat index (HI) on the EVI trend because these weather parameters are more tightly linked to productivity than the parameters directly measuring wind or metrics that depend on wind.Due to malfunctioning, there were data gaps between the end of March 2006 and May 2007, December 2008 and October 2009, March and July 2011, and April 2012 to February 2014.Unfortunately, it was not possible to complement these data gaps with other meteorological data from the surroundings because the closest station was in pine plantations.

Data Analysis
We based our analysis on a sub-region of ~1.8 km 2 in the immediate vicinity of the weather station so that we could test effects of local weather on local productivity measurements.Further, this close match between the scales of the measurements ensured that the potential relationships are interpretable at the local level, before scaling up to more general patterns for larger areas.This choice reduced the number of available data points, however, still within what is considered a sufficient sample size if significant differences are to be detected and we do not think that it hinders the relevance of the findings.These locations were randomly chosen across all land-cover types, and the random selection provided proportional representation of each land-cover.We then extracted the aggregated weather data for each day in our time series for and merged this information to the EVI values of our Landsat time series.This resulted in a time series containing, overall, 107 observations (i.e., time points) between 2000 and 2015, which served as the input for our analysis (Table 2).To measure intra-annual variability in productivity, we calculated the coefficient of variation (CoV) for the EVI time series and for the EVI for each land-cover type individually.CoV measures the dispersion in the values across a period of time, therefore, allowing for the comparison of the relative variability of EVI for each land-cover type [55].High CoV values indicate a strong intra-annual variability in EVI, while low CoV values indicate a more homogeneous EVI and, therefore, a lower intra-annual variability.
To test whether there is inter-annual variability within the EVI time series, we used Kendall's tau (τ; [56]), a commonly used test to detect trends in time series (for examples see [57,58]).We assessed whether there were significant trends in overall EVI, and by land-cover type.addition we tested whether EVI trends were varying in the same direction or in different directions by correlating minimum and maximum EVI using Spearman's correlation coefficient.
Finally, to test whether temperature (T), relative humidity (RH), and heat index (HI) were related to EVI trends, we correlated the two trends for the dates for which we had both weather and EVI data.We repeated the analysis for the average EVI combined for all land-cover types, and for average EVI per land-cover type.To do so we correlated maximum and minimum EVI with the three weather parameters using Spearman's correlation coefficient.

Intra-Annual Variability
EVI in the study area followed a yearly phenological cycle for all land-cover types.EVI was highest in spring (February to May; EVI Feb = 0.355 ˘0.05, EVI Mar = 0.397 ˘0.05, EVI Apr = 0.415 ˘0.04, EVI May = 0.375 ˘0.06) and lowest during summer (July, August, and September; EVI Jul = 0.259 ˘0.03, EVI Aug = 0.238 ˘0.04, EVI Sep = 0.229 ˘0.04).Overall, the onset of productivity occurs from October to November, at the time of the first rains, and the downfall from April to May, as temperatures increase.Cork oak woodlands with understorey have the smoothest EVI yearly cycle, i.e., have the lowest CoV, whereas grasslands have the least smooth cycle with the highest CoV (Figure 3).Sparse cork oak woodlands without understorey follow closely the CoV of grasslands, likely because of their savannah-like morphology.

Inter-Annual Variability
During our study period (2000-2013) we found a small but significant downward trend in maximum EVI and a significant upward trend for minimum EVI.The increase in minimum EVI was slightly higher than the decrease in maximum EVI, with little correlation between maximum and minimum EVI (R 2 = 0.203; Spearman's rho = 0.211; Table 3).
Contrary to the expectation, cork oak woodlands without understorey showed a significant increase in minimum EVI and a trend towards a decreasing maximum EVI (Table 3).The average minimum EVI over the entire time series and for all land-cover types was 0.171, whereas, for DCoWt, it was 0.22.Further, and in line with the expectations, sparse cork oak woodland with understorey shrubs show a decreasing trend in maximum EVI.The average maximum EVI over the entire time series and for all land-cover types was 0.464, whereas, for SCoWt, it was 0.394.Despite the high intra-annual variability in grasslands EVI we found no significant trend in EVI over the entire time series, except for a potential increase in minimum EVI (average min EVI = 0.20, average max EVI = 0.434; Table 3).

Inter-Annual Variability
During our study period (2000-2013) we found a small but significant downward trend in maximum EVI and a significant upward trend for minimum EVI.The increase in minimum EVI was slightly higher than the decrease in maximum EVI, with little correlation between maximum and minimum EVI (R 2 = 0.203; Spearman's rho = 0.211; Table 3).
Contrary to the expectation, cork oak woodlands without understorey showed a significant increase in minimum EVI and a trend towards a decreasing maximum EVI (Table 3).The average minimum EVI over the entire time series and for all land-cover types was 0.171, whereas, for DCoWt, it was 0.22.Further, and in line with the expectations, sparse cork oak woodland with understorey shrubs show a decreasing trend in maximum EVI.The average maximum EVI over the entire time series and for all land-cover types was 0.464, whereas, for SCoWt, it was 0.394.Despite the high intra-     Since the relation between weather and EVI could be the opposite in growing and non-growing seasons we repeated the analysis just for the growing season (December to May).During the growing season, we found a strong negative correlation between EVI and temperature and heat index, and a less strong correlation with relative humidity (Table 5).We found the strongest correlation for sparse cork oak woodlands without understorey and grasslands, and the weakest correlation for sparse cork oak woodlands with understorey.

Discussion
Our study aimed at assessing whether trees, shrubs, and grasses in cork oak (Quercus suber) woodlands in Southern Portugal showed differences in the relationship between productivity and climate.We found that the selected metric of productivity-EVI-in areas with dense and sparse cork oak forest without understorey responded differently than in areas with understorey and of grasslands.
Cork oak woodlands without understorey maintained higher EVI and even a significant increase in minimum EVI.Further, this land-cover type showed a much stronger relationship with relative humidity then when there is a layer of understorey.During the growing season, however, EVI over this land cover type was mostly negatively correlated with temperature.Therefore, temperature and relative humidity seem to control EVI in the growing season and outside the growing season, respectively.Since the growing season in our study area occurs during the colder months and after the first rains, this explains the strong negative correlation with temperature.In this period there is no water limitation; thus, EVI increases, which mostly occurs during the hot and dry periods when relative humidity becomes the stronger control over EVI.A possible explanation is that the modification of the leaves in evergreen oaks allows them to maintain their photosynthetic rates despite water stress [59].Evergreen cork oaks have modified leaves called sclerophyllous leaves, which have a waxy layer on the upper side of the leaf and stomata in the lower side of the leaf.The stomata in these leaves are surrounded by hairs (trichoms), which capture and maintain moisture in the lower side of the leaf, therefore, allowing the stomata to stay open for longer periods of time, enhancing photosynthesis, resulting in the relatively constant EVI throughout the time series and a strong correlation with relative humidity.
We also found a significant downward trend in maximum EVI is sparse oak woodlands with understorey and a lower correlation with relative humidity.In these areas the dominant species in the shrub layer are Cistus spp.Contrary to cork oaks, Cistus spp.do not have sclerophyllous leaves and its growth is strongly related to precipitation and in dry periods these species tend to start senescence earlier than expected [60].This could justify why we find a significant decrease in EVI when patches are dominated by shrubs.This trend can be further exacerbated by the fact that these shrubs have a life expectancy of 10 years and that they may be replaced every 5-8 years, being two time periods within our time series.
We also found that EVI did not change in dense cork oak woodlands with shrub understorey and that the presence of understorey reduced the strength of the relationship between EVI and relative humidity and between EVI and temperature during the growing season.This could be because of the opposing direction of the trends in the oak and the shrub layers, thus, suggesting that the tree canopy and the shrub canopy layers have different physiological responses to RH, T, and HI.Despite Cistus species having an ability to restore soil and aid in oak regeneration [60,61], our results suggest that this is limited by the ability of the shrubs to persist.This could be because the dominant C. salvifolius produces allelopathic substances therefore maintaining a mono-specific shrub layer that is vulnerable to changes in precipitation/relative humidity.
These results suggest that the overall upward trend in minimum EVI is driven by the canopy of cork oaks, the downward trend in maximum EVI by the understorey canopy of shrubs, and that grasslands do not influence EVI trends over the time period we performed our analysis.Our results also suggest that in the future EVI could improve because the rate of increase in minimum EVI is greater than the rate of decrease in maximum EVI, and that this contingent on management of the shrub understory, as it affects the rate of decrease in maximum EVI.Grasslands showed a much stronger intra-annual variation in EVI and a stronger correlation between those and relative humidity, possibly because grasslands respond to changes in precipitation in a shorter cycle than shrubs and oaks.These results are corroborated by the trends in the sparse oak woodlands without understorey, which are very close to those of grasslands.Since maximum and minimum EVI trends are not correlated, i.e., the two trends are occurring independently of each other, it suggests that different processes may be operating simultaneously.We suggest that oaks are better adapted to Mediterranean climates and, therefore, are more resilient to changes, or that oaks take longer to respond to the ongoing changes than the understorey and their response will be delayed.
A few uncertainties are worth mentioning.First, our study area is comparatively local and, as such, we are unable to draw conclusions on the relationship across a larger scale.However, since our analysis relied on standardized methods (i.e., Landsat EVI data and climate station data) there is only little reason to believe that the relationship will not hold outside of our study area.Second, calculating EVI requires reflectance values, which we only achieved through a standardized method for atmospheric correction (i.e., LEDAPS).While some uncertainties may result from using LEDAPS, our EVI time series analyses suggest that using Landsat EVI time series data elsewhere in the world will be reliable.Last, a few uncertainties may remain from the input land-cover map, but we consider these uncertainties marginal since the land cover types used were hand-drawn and based on detailed regional knowledge.
Nevertheless, being able to link the remote sensing measurements to these processes is fundamental for ecosystem modeling and management.Our analysis shows one way in which we can take advantage of the existing archive of remote sensing data and analytical tools to understand ecosystem processes and further expand our understanding of the relation between productivity-biodiversity. While our results are focused on a local area, they are indicative that we can better link remote sensing metrics with local environmental conditions to understand the fundamental processes that may drive larger scale patterns.

Conclusions
We found different temporal trends in productivity for the tree canopy and the shrub canopy, and these trends in productivity were positively associated with relative humidity and negatively with temperature.Independently of whether oaks are better adapted to Mediterranean climates and, therefore, more resilient to changes, or oaks are taking longer to respond to the ongoing changes than the understorey, these results have implications for the persistence of cork oak woodlands and their associated biodiversity.While relative humidity seems to be stressing productivity, its effect is stronger in deciduous land-cover types than in cork oak dominated land-cover types.Cork oak dominated land-cover pixels were stable, in fact even improving productivity over this time period.However, when cork oak was replaced by shrubs or herbs or its density decreased, productivity also decreased.Therefore, if future land use preserves cork oaks, they may be sufficiently resilient to withstand changes in climate, or at least have yet to show major effects.However, if changes in land use favor replacement of cork oaks by other land-cover types or if cork oaks show lagged responses, then the effect of climatic changes may be stronger than expected, challenging the persistence of its associated biodiversity and social-ecological system.

Figure 1 .
Figure 1.(a) Study area location and distribution of oak woodlands in Portugal; (b) Dense cork oak woodland with shrub understory; (c) Dense cork oak woodland without shrub understory; (d) Sparse cork oak woodland with shrub understory; (e) Sparse cork oak woodland without shrub understorey.

Figure 1 .
Figure 1.(a) Study area location and distribution of oak woodlands in Portugal; (b) Dense cork oak woodland with shrub understory; (c) Dense cork oak woodland without shrub understory; (d) Sparse cork oak woodland with shrub understory; (e) Sparse cork oak woodland without shrub understorey.

Figure 2 .
Figure 2. Photo-interpreted land-cover types and respective extent in the study area.

Figure 2 .
Figure 2. Photo-interpreted land-cover types and respective extent in the study area.

Figure 3 .
Figure 3. Intra-annual variability in EVI per land-cover types.Land-cover types: dense cork oak woodlands with understory (DCoW), dense cork oak woodlands without understory (DCoWt), sparse cork oak woodlands with understory (SCoW), sparse cork oak woodlands without understory (SCoWt), and grasslands (G).Box plots represent the mean (bold line within the box), variance (size of the box), and the range of values (whiskers).Open circles and asterisks represent outliers, out and far out respectively.

Figure 3 .Table 3 .
Figure 3. Intra-annual variability in EVI per land-cover types.Land-cover types: dense cork oak woodlands with understory (DCoW), dense cork oak woodlands without understory (DCoWt), sparse cork oak woodlands with understory (SCoW), sparse cork oak woodlands without understory (SCoWt), and grasslands (G).Box plots represent the mean (bold line within the box), variance (size of the box), and the range of values (whiskers).Open circles and asterisks represent outliers, out and far out respectively.

Table 1 .
Species richness and their International Union for Conservation of Nature (IUCN) conservation status in the study area (CR-Critically Endangered, EN-Endangered, VU-Vulnerable, NT-Near Threatened).

Table 2 .
Number of Landsat pixels over each land-cover type in the study area.