Spectral-Based Monitoring of Climate Effects on the Inter-Annual Variability of Different Plant Functional Types in Mediterranean Cork Oak Woodlands

: Using remotely sensed data to estimate the biophysical properties of vegetation in woodlands is a challenging task due to their heterogeneous nature. The objective of this study was to assess the biophysical parameters of different vegetation types (cork oak trees, shrubs and herbaceous vegetation) in cork oak woodland through the analysis of temporal trends in spectral vegetation indices (VIs). A seven-year database (from 2011 until 2017) of in situ observations collected with a ﬁeld spectroradiometer with a monthly basis was used and four VIs were derived, considered as proxies for several biophysical properties of vegetation such as biomass (Normalized Difference Vegetation Index—NDVI); chlorophyll content (MERIS Terrestrial Chlorophyll Index-MTCI), tissue water content (Normalized Difference Water Index—NDWI) and the carotenoid/chlorophyll ratio (Photochemical Reﬂectance Index—PRI). During the analyzed period, some key meteorological data (precipitation, temperature, relative air humidity and global radiation) were collected for the study site, aggregated at three different time-lags (short period (30 d), medium period (90 d) and hydrological period (HIDR)), and their relationship with VIs was analyzed. The results showed different trends for each vegetation index and vegetation type. In NDVI and NDWI, herbaceous vegetation showed a highly marked seasonal trend, whereas for MTCI, it was the cork oak and Cistus salvifolius , and for PRI, it was Ulex airensis that showed the marked seasonal trend. Shrubs have large differences depending on the species: the shallow-rooted Cistus salvifolius showed a higher seasonal variability than the deep-rooted Ulex airensis . Our results revealed the importance of temperature and precipitation as the main climatic variables inﬂuencing VI variability in the four studied vegetation types. This study sets up the relationships between climate and vegetation indices for each vegetation type. Spectral vegetation indices are useful tools for assessing the impact of climate on vegetation, because using these makes it easier to monitor the amount of “greenness”, biomass and water stress of vegetation than assessing the photosynthetic efﬁciency. Proximal remote sensing measurements are fundamental for the correct use of remote sensing in monitoring complex agroforest ecosystems, largely used to inform policies to improve resilience to drought, particularly in the Mediterranean region.


Introduction
Oak woodlands are agro-silvopastoral ecosystems in which evergreen oak species (Quercus ilex L. ssp. rotundifolia and/or Quercus suber L.) coexist with shrubs and annual herbaceous vegetation [1]. They provide several ecosystem services, such as cork production, pasture, the prevention of severe wildfires, carbon sequestration and biodiversity conservation [2]. In Portugal, they represent 33% of the total forested area [3] and are heterogeneously distributed in the Western Mediterranean region [4]. This region is characterized by mild, wet winters and hot, dry summers, with a high susceptibility to aridity increase [5]; therefore, it is considered a hot spot for climate change [6]. Climatic projections for the Mediterranean region foresee longer summer droughts, with more frequent extreme events, such as heatwaves and an increase in annual mean temperature between 2 and 4.5 • C, with a substantial change in the rainfall distribution [7]. The coexistence of several plant functional types (PFT, i.e., trees, herbaceous vegetation and shrubs) in this ecosystem makes it difficult to understand how each one, individually, is affected by climate drivers. PFTs are expected to play different roles in ecosystems [8] (p. 709). PFTs are defined as non-phylogenetic sets of species that have similar effects on ecosystem functioning and show close similarities in their response to environmental and biotic controls [9]. Better knowledge of climate-PFT interactions is needed to fully understand climate-ecosystem interactions.
Several studies have been conducted to understand these dynamics in cork oak woodland: climate-cork growth relationships [10]; differences in the vegetation indices representing the fraction of photosynthetically active radiation (fPAR) and light use efficiency (LUE) among different PFTs [11]; the effect of a dry winter season on understory species [12]; how cork oak woodlands respond to low temperatures in winter [13]; the role of photochemistry in the resilience to severe summer stress of Cistus salvifolius shrubs growing in a cork oak woodland [14]; and a meta-analysis on herbaceous and woody plants fine-root responses to precipitation changes [15]. Many of these studies used traditional fieldwork to conduct their analyses, which usually involves intense and expensive soil and/or plant sampling on a closely spaced grid [16]. An alternative to ground field sampling is the use of indicators that are known to be proxies of vegetation properties, such as vegetation indices (VIs). These can be defined as spectral data algorithms used for quantitative and qualitative evaluations of vegetation cover, or as indicators of biophysical and biochemical properties of vegetation (vigor, growth dynamics, canopy cover, leaf composition, abiotic/biotic stress) [17]. These spectral data can be retrieved for long periods through satellite images in a more cost-effective way than with traditional fieldwork [18]. Regardless of this advantage, there are two main problems associated with satellite data retrieval: (1) pixel resolution, which varies among different satellites (for example, 30 m for Landsat and 500 m for MODIS); and (2) the revisit time (for example, 16 days for Landsat) [17]. These characteristics lead to a loss of temporal information. MODIS and Landsat were referred given their importance for long-term monitoring and retrospective analyses of vegetation properties. When using this type of platforms, one must find a compromise between spatial and temporal resolution of spectral variables to obtain data in useful timescales and with an appropriate resolution, which will obviously depend on the ecosystem being studied, its size and the purpose of the study. Thus, it is an important task to disentangle the spectral trends of different components of the ecosystem (PFT) with field data, such as hyperspectral data at canopy level, to upscale it to satellite platforms. Hyperspectral sensors allow the collection of contiguous spectral bands throughout the visible, near-infrared, mid-infrared and thermal infrared portions of the electromagnetic spectrum enabling the construction of an almost-continuous reflectance spectrum [19].
One of the most used VIs is the Normalized Difference Vegetation Index (NDVI), measuring the "greenness" of vegetation, which is positively correlated with canopy density [20]. It is widely used at different spatial and temporal scales, assessing largescale vegetation trends [21,22]. Ref. [23] (pp. [17][18][19] showed that long-term NDVI trends, based on 34 years of Landsat imagery, are correlated with cork oak woodland health and productivity. Ref. [24] (pp. 2384-2387) assessed a relationship between NDVI and climatic factors according to latitude. They showed that in the northern latitudes, temperature was the primary climatic factor of vegetation increase, whereas in southern latitudes, the decreases in greenness were strongly correlated with both the increase in temperature and decreases in precipitation. Another frequently used VI, the Normalized Difference Water Index (NDWI), was developed to estimate plant water contents [25]. Significant correlations between NDWI and soil moisture content, pasture moisture content, surface pasture temperature, biomass productivity and pasture quality degradation index were established, suggesting that NDWI can be used in site-specific management for monitoring pasture deterioration/restoration [26]. The MERIS Terrestrial Chlorophyll Index (MTCI) is related to the chlorophyll content of vegetation. This is reflected in the abrupt reflectance change (red edge) that occurs in healthy leaves between red and near-infrared (NIR) wavelengths [27]. The existence of a strong relationship between seasonal temperature variation and MTCI has also been demonstrated [28] (p. 4), from which it was inferred that the chlorophyll content was sensitive to temperature. Another VI, the Photochemical Reflectance Index (PRI), is used as an indicator of the de-epoxidation state of xanthophyll pigments, and hence is related to photosynthetic processes [29]. In the long term, the PRI is sensitive to changes in carotenoid pigments which are also involved in photosynthetic carbon uptake [30]. The PRI has been used to predict light use efficiency (LUE), the amount of carbon sequestered per unit of photosynthetically active radiation (PAR) absorbed, and can be usefully employed to assess the photosynthetic response of vegetation to climate drivers over time [31]. Vegetation may not have an immediate biophysical response, represented by the vegetation indices, to variations in climate, showing a delay that can vary depending on the PTF and on the biophysical/structural property represented by the vegetation index [11]. To the best of our knowledge, there are few studies which consider different time-lags of climatic variables as independent variables to assess their influence on the various types of vegetation in cork oak woodland. Therefore, in this study, we tested different time aggregation periods (30 days (30 d), 90 days (90 d)) and all the period elapsed since the beginning of the hydrological year (HIDR) in the relationship between the VI temporal series and climatic variables for each considered vegetation type. The main goal of this study was to assess the suitability of using spectral vegetation indices to monitor the impacts of climate drivers on the biophysical properties of three PFTs (cork oak, herbaceous vegetation and shrubs) coexisting in a cork oak woodland. The specific objectives are: To compare the seasonal and inter-annual variability of a set of VIs, selected from the literature, among different PFTs; To assess which are the most determinant climatic variables influencing the variability of the VIs; To assess the sensitivity of each PFT to the temporal aggregation period of those climatic variables, i.e., to infer time lag effects on the impacts of climate factors on vegetation status.

Study Area
The study was conducted in Herdade da Machoqueira do Grou estate (Figure 1), located in Coruche, central mainland Portugal (39 • 08 20.9 N, 9 • 19 57.7 W, 165 m altitude). The site is a Mediterranean cork oak woodland in which three different PFTs coexist: trees, shrubs and herbaceous vegetation [11]. The tree density is about 177 trees/ha with a Leaf Area Index (LAI) of about 1.5. The herbaceous layer is dominated by grasses (44.5%) and legumes (28.7%) [12]. The more abundant shrub species are the semi-deciduous cistus and the evergreen sclerophyll Ulex airensis. The area is characterized by a Mediterranean climate with hot and dry summers and a mildly cold wet season. The mean annual precipitation recorded at the meteorological station, located at Santarém (39 • 12 N, 8 • 44 W), for the period of 1981-2010 was 652 mm, with a mean daily temperature of 17 • C. The soil is a cambisol with 81% sand, 5% clay and 14% silt [11].

Spectral Data
Hyperspectral field measurements were collected monthly from April 2011 to December 2017. Four cork oak trees, four individuals of each shrub species (cistus and ulex) were sampled for monitoring, whereas for the herbaceous vegetation, eight quadrats (50 × 50 cm 2 ) were selected. The two shrub species were considered separately due to their different structure and physiology. Our data have a gap in shrub spectra collected between 2014 and 2016, due to clear-cutting that occurred in 2014. The individuals/samples of each PFT were randomly chosen, with the condition that they were not located in the shadow cast by trees at the time of measurements. All measurements were taken at nadir and 1 m above the target. To minimize atmospheric and bidirectional effects, the measurements were performed under clear sky conditions between 12:00 and 14:00 solar time. Each measurement of spectral reflectance was the average of 5 replicates and was taken with a portable FieldSpec3 spectroradiometer (ASD Inc., Boulder, CO, USA). Remote sensing of vegetation studies with a focus on biochemical and biophysical vegetation properties have been widely carried out in the visible and near-infrared region (VNIR) ranging from 0.3 to 1.0 µm, and short-wavelength IR (SWIR) ranging from 1.0 to 2.5 µm regions [32]. The VNIR was measured with a spectral resolution (full-width half-maximum) of 3 nm and a sampling interval of 1.4 nm, whereas the SWIR was measured with a spectral resolution of 10 nm and a sampling interval of 2.0 nm. A white reference of known reflectance (Spectralon panel, Labsphere, Inc., North Sutton, NH, USA) was used to normalize for variations in atmospheric conditions and to convert the measurements into absolute reflectance (ρ). These reflectance measurements were used to calculate the four VIs considered in this study (Table 1). Climatic variables were measured at the top of a 22 m tall tower located in the study site. Data were collected at 10 Hz with a laptop connected to a data logger (CR1000, Campbell Scientific, Inc., Logan, UT, USA). Precipitation (ARG100, Environmental Measurements Ltd., Gateshead, UK), global radiation (NR Lite2 Kipp&Zonen, Delft, The Netherlands), humidity and air temperature (CS215, Campbell Scientific, Inc., Logan, UT, USA) were collected continuously and averaged at 30 min time steps [28]. Then, the data were averaged (maximum temperature, TMAX; minimum temperature, TMIN; relative air humidity, RH; and global radiation, RAD) or summed (precipitation, PP) at different time-steps: 30 days (30 d), 90 days (90 d) and all the periods elapsed since the beginning of the hydrological year (HIDR-considered to start on 1 September, according to the rainfall pattern in Portugal) prior to the date of each spectral measurements. The aggregation period of climatic data, matched to spectral measurements collected in June or July, was larger than that matching the spectral measurements conducted in January or February.

Data Analysis
The main steps of the methodology are summarized in Figure 2. Temporal time series of spectral data were plotted by vegetation indices and vegetation type for the study period and then compared. To explore the seasonal (winter, spring, summer and autumn) effects in vegetation-type VIs, we performed an ANOVA type II assay for each of the vegetation indices using the R package "car" [35] with two factors (season and vegetation type), using Sidak's method for p-value adjustment to multiple comparisons. A value of p < 0.05 (two-tailed) was considered. The seasonal periods were divided, according to the rainfall pattern in Portugal [36] into winter (December, January and February), spring (March, April and May), summer (June, July and August) and autumn (September, October and November). To measure the strength of the linear association between the VIs (NDVI, NDWI, MTCI and PRI) and climatic variables (TMAX, TMIN, PP, RH and RAD), for each type of vegetation (cork oak, herbaceous vegetation, cistus and ulex) and for each time step (30 d, 90 d and HIDR), correlation analysis for the 7-year study period was performed, using Pearson's correlation coefficient (r). To detect potential lags or the cumulative effects of climate variables on VIs, we performed Generalized Linear Models (GLM) analyses using the R package "glmulti" [37]. The selection of the best models was based on the Akaike Information Criterion (AIC) [38], whereas multicollinearity was tested through variance inflation factor (VIF) analyses [39]. All climatic variables were used as input data of each model, but two main limitations were applied: (1) only the main effects were returned (the interactions of the variables were not explored); and (2) only the two most significant variables were chosen. These limitations are intended to reduce the entropy for the interpretation and clearly identify the more influential climatic variables for the models.  Figure 3 shows the time series of the meteorological variables considered in this study, namely, total annual precipitation (mm) ( Figure 3A), average of the maximum and minimum daily temperature ( • C) ( Figure 3B), average daily global radiation (W·m −2 ) ( Figure 3C) and relative humidity (%) ( Figure 3D) recorded at the study site between April 2011 and December 2017. In this period, there were four years (2012, 2015, 2016 and 2017) with total precipitation below the annual average for the study period (592 mm) ( Figure 3A).

Temporal Trends of Vegetation Indices
In this study, we assessed different temporal trends: (1) seasonal trends (to study the differences among seasons: winter, spring, summer and autumn); and (2) over the study period (to assess the trend in each VI for the different PTFs). Figure 4 shows the analysis of variance for each VI (NDVI, NDWI, MTCI and PRI) with two factors (season-winter, spring, summer and autumn-and vegetation type-cork oak, herbs, cistus and ulex). For NDVI ( Figure 4A), cork oak and ulex showed no significant differences between the seasons. For herbaceous and cistus, winter and spring differed from summer and autumn with higher values for the winter/spring season and lower values for the summer/autumn season. The same trend was also observed in the NDWI ( Figure 4B) for the herbaceous vegetation and the PRI ( Figure 4D) for cistus. In winter, the NDWI of the cork oak was statistically different from that of the summer, with higher values in summer. In autumn, the NDWI of cistus was statistically different from the other seasons, presenting lower values ( Figure 4B). In cork oak, summer and autumn MTCI values were significantly higher than those observed in winter and spring, whereas herbaceous vegetation showed no significant differences among the seasons except for winter, where MTCI presented lower values than for other seasons (this result correlates with the observed trends in the section above). In MTCI ( Figure 4C), winter values for cistus stand out from the other seasons, as do spring values for the ulex. The PRI ( Figure 4D) for cork oak and ulex was notable in winter. Summer values of PRI observed in cistus were statistically higher than those observed in winter and spring. Regarding the VI time series (Figure 5), our results showed that the NDVI cork oak had the lowest amplitude of values ranging from 0.71 to 0.87 (Table 2). In contrast, herbaceous plants were the PTF that showed the highest amplitude ranging from 0.16 to 0.63. Despite the low amplitude, the NDVI values in the cork oak varied from the lowest recorded in April, reaching the highest values immediately afterward in June ( Figure 5A). In the herbaceous species, the NDVI reached the highest values in the winter (December/January) and the lowest in the summer (July/August). Concerning shrubs, cistus exhibited the larger variation in NDVI values, being the lowest in the summer (August) and prolonged until October, when they begin to rise, reaching their maximum in February. Ulex presented a seasonal behavior similar to cistus but with a smaller variation in values ( Table 2). For this shrub species, the lowest values were registered in the summer (around August extending until October) and the highest in the winter (January/February) ( Figure 5A).  Similarly to the NDVI, the NDWI of the cork oak exhibited the lowest variation over the study period, ranging from −0.09 to 0.04. The lowest values were recorded in March/April and the highest around May/July. The temporal trends in NDWI for ulex also showed a small range of values, ranging from the lowest values in October to the highest values in March/April. A larger variation in NDWI was observed for herbaceous vegetation and cistus. In herbaceous vegetation, the lowest values were recorded in July, extending until October, whereas the highest values were observed in winter (January/February). In cistus, the highest values were observed in March and the lowest in autumn (October).
Concerning MTCI, cork oak and cistus showed a larger variability ( Table 2). The minimum values of MTCI of cork oak were generally recorded in May, whereas for cistus, the minimum was in July. The maximum value for cork was observed in June, similarly to that observed for the NDVI, whereas the MTCI of cistus reached the maximum values in February. The Ulex airensis registered the highest values in March after it decreased until reaching the minimum in October. Concerning herbaceous plants, MTCI presented an unexpected behavior, as it would be expected that MTCI presented a similar trend to NDVI, because they are indices related to chlorophyll content and the "greenness" of plants, respectively. What we observed in the results obtained for herb NDVI and MTCI values ( Figure 5A,C) was that the MTCI of herbaceous plants had higher values in the summer and autumn months than in the winter and spring months, a behavior that did not occur with the NDVI of herbaceous plants. In fact, the NDVI of herbaceous plants presented higher values in winter than in summer, which can be explained by the absence of herbaceous vegetation in the summer season. In order to understand this result and taking into account that, during the summer and fall, an absence of green herbaceous vegetation was observed, we plotted the regression between the NDVI and MTCI of herbs per season (Figure 6), and also evaluated the contribution (inference) of the soil per se (Figure 7) (we measured the reflectance of bare soil separately from the vegetation types over the study period) to the herbaceous NDVI and MTCI results. We designated these outputs as "soilNDVI" and "soil MTCI". It showed that in winter and spring ( Figure 6A,B), these two vegetation indices of herbs followed the same trend. On the other hand, when it comes to summer and autumn ( Figure 6C,D), NDVI and MTCI showed the opposite trend, i.e., with MTCI increasing when NDVI decreases.  The PRI ( Figure 5D) was the vegetation index that presented the smallest variation in values among the vegetation indices. This variation did not differ much for the various PFTs, but even so, ulex stands out as being that with the highest amplitude. The cork oak reached the minimum value in October and the maximum value in February. The The PRI ( Figure 5D) was the vegetation index that presented the smallest variation in values among the vegetation indices. This variation did not differ much for the various PFTs, but even so, ulex stands out as being that with the highest amplitude. The cork oak reached the minimum value in October and the maximum value in February. The herbaceous vegetation and cistus presented the minimum values in winter (December/January) and the maximum values in July, whereas the ulex recorded the maximum value in February, falling to the minimum value immediately after, in April.

The Influence of Climate
The correlation results for each vegetation type, between VIs and climatic variables aggregated at different periods (precipitation was summed and temperature, relative humidity and global radiation were averaged), are shown in Figure 8. Regarding cork oak, the NDVI and the NDWI presented very similar correlation results, despite the correlation values in the NDWI being higher than those of the NDVI. These two vegetation indices had a positive correlation with temperature (maximum and minimum) and radiation, whereas with precipitation and relative humidity, the correlation was practically null or negative. One exception should be noted: for accumulated precipitation for the whole hydrological year, there was only one period of aggregation verified with positive correlation (NDVI with r = 0.42 and NDWI with r = 0.48). Both presented higher correlation values with shorter periods of aggregation (30 days) for temperature and radiation. In the cork oak MTCI, the highest values were found for the aggregation period of 90 days for practically all climatic variables, except for rainfall, where the hydrological period had a higher correlation value than the other periods. As for the NDVI and NDWI, the cork oak MTCI had a positive correlation with temperature (maximum and minimum) and radiation, whereas with precipitation and relative humidity, the correlation was practically null or negative. The PRI of cork oak showed different trend values. The correlation values were negative with temperature, radiation and precipitation, whereas for relative humidity, they were positive. Regarding temperature (maximum and minimum), the cork oak PRI presented the highest correlation value over 30 days, whereas precipitation, relative humidity and radiation, showed the stronger correlation over the period of 90 days.
Herbaceous vegetation and cistus had very similar correlation values. We observed that for NDVI and NDWI, the correlation values were negative for temperature (maximum and minimum) and radiation and positive for precipitation and relative humidity. In turn, the PRIs of these two types of vegetation presented positive values for temperature (maximum and minimum) and radiation, and negative values for precipitation and relative humidity. The exception was the correlation value between the PRI of herbaceous plants and hydrological precipitation, where the value was positive (r = 0.45). In general, all vegetation indices obtained from herbaceous plants showed stronger correlations with the climate variables when aggregated at shorter periods (30 and 90 days), whereas the MTCI and PRI of cistus presented higher correlations when aggregated at a shorter period (30 days) and the NDVI and NDWI with the longer periods (90 days and hydrological year). This similar trend in the correlation values between herbaceous and cistus was not verified in the MTCI where the correlation values of the herbaceous with climate variables were closer to the correlations observed for the MTCI of cork oak. Thus, herbaceous plants showed positive correlation values with temperature and radiation, whereas those with precipitation and relative humidity were negative. The MTCI of cistus maintained the trend in the correlation values of the other vegetation indices.
Ulex was the vegetation type that had the weakest correlation values between NDVI or NDWI and meteorological variables. In NDVI, these were practically null, except for precipitation, where the highest value was with 90-day precipitation (r = 0.30). Ulex had negative correlation values with temperature (maximum and minimum) and radiation with NDWI, MTCI and PRI, with higher values for the aggregation period of 90 days in NDWI and MTCI and the aggregation period of 30 days for the PRI. Regarding precipitation and relative humidity, it had positive values, except for the PRI, whose correlation with precipitation was negative. With global radiation, ulex was strongly correlated with longer periods of aggregation of climatic variables in the NDWI and MTCI and shorter periods in the PRI. Temperature (maximum and minimum) was strongly correlated with the 90-day aggregation of climatic variables in the NDWI and MTCI, and with the 30-day aggregation in the PRI. In precipitation, this correlation was stronger with the longer aggregation periods of the climate variables in NDWI, MTCI and PRI, whereas in relative humidity, it was over 90 days. Table 3 shows the selected GLM models in which the dependent variables are the VIs (NDVI, NDWI, MTCI, PRI) and the independent variable are meteorological variables (TMAX, TMIN, PP, RH, RAD) at different time lags (30 d , 90 d, HIDR). Cork oak presented better results at the annual scale with NDWI (R 2 = 0.58, p-value = 0.0000), whereas NDVI was the worse (R 2 = 0.33, p-value = 0.0000). Minimum temperature and precipitation were the most important factors for estimating NDVI, MTCI and PRI of cork oak, although in the latter, the influence was negative. In NDWI, the best model selected radiation aggregated at 30 days (short period).
In herbaceous plants, the best model obtained was with the NDVI (R 2 = 0.68); in this case, the radiation was aggregated at 90 days, the climate variable selected in the model. Additionally, in the PRI, for the second-best herbaceous model (R 2 = 0.65), the 30-day radiation was selected in the model alongside relative humidity aggregated at the hydrological year. For NDWI and MTCI, the selected variables were temperature and precipitation, but at different aggregation time periods: longer for NDWI and shorter for MTCI.
In cistus, the best model was the NDWI with R 2 = 0.78. Maximum temperature and hydrological precipitation proved to be important for estimating the NDVI, NDWI and PRI, although for the latter the precipitation was selected for the shortest period (30 days). In cistus MTCI, the selected variables were the hydrological relative humidity and the 30-day radiation.
Ulex showed the weakest performance for all vegetation indices. The best model was obtained for the estimation of MTCI (R 2 = 0.38). The temperature and precipitation were the climatic variables selected for estimating the vegetation indices. In the ulex NDVI, the precipitation at 90 days was selected, and in the NDWI and PRI, the minimum temperature and precipitation (TMIN_90 and PP_HIDR; TMIN_30 and PP_30, respectively). In the MTCI of ulex, the maximum temperature at 90 days (with a negative influence) and the precipitation at 30 days (with a positive influence) were selected.

Discussion
Our results show the trends in several VIs for the seven years of study, as we explored the differences in the VIs among cork oak, herbaceous vegetation and two shrub species (cistus and ulex) for the different seasons. Due to forest management practices, the shrubs were cut in 2014. Therefore, there is a gap in the data for this PFT and the obtained results must be interpreted taking into consideration the use of a shorter data time series.

Cork Oak
Our VI/climate oak model results showed minimum temperature and precipitation to be the most important climatic factors in the modeling of cork oak NDVI (canopy greenness), MTCI (chlorophyll content of the cells) and PRI (light assimilation in leaves), but at different time scales. For NDVI and PRI, a short time-scale of minimum temperature was selected (30 days), revealing that vegetative vigor and the carotenoid/chlorophyll ratio of cork oak are more sensitive to minimum temperature differences (positively for NDVI and negatively for PRI) in short periods. The results showed that cork oak is sensitive to low temperatures, which is in accordance with previous studies [13] (p. 1089), and suggests that this sensitivity of cork oak to low temperatures may explain the change in the geographic distribution of the species, causing the replacement of Quercus suber by Quercus ilex, another evergreen oak more tolerant of cold. Similar results for the oak PRI relationship with minimum temperature in a Mediterranean oak forest (holm oak forest) were found by [40], in which the minimum values of maximum photochemical efficiency were reached when cold was combined with drought conditions. Additionally, [13] (p. 1088) observed that cork oak is highly sensitive to low temperatures through the reduction in Fv/Fm values as a result of photoinhibition caused by low winter temperatures. As for NDWI, the cork oak model only selected the climatic variable with a high correlation value (r = 0.76)short-term radiation (30 days). This result can reflect how well cork oak is adapted to the Mediterranean climate in which high irradiance values are commonly observed [41]. Although in the NDWI cork oak model the minimum temperature is selected, its influence assumes a longer period (hydrological period) and negative effect, meaning that the higher the minimum temperature in the hydrological year, the lower cell's water content of cork oak. In MTCI, as for NDVI, the minimum temperature positively affects this index but at a longer aggregation period (90 days), meaning that the chlorophyll content of the cells is more slowly affected by the minimum temperature than the vegetative vigor of the canopy.
Cork oak also showed low susceptibility to water stress, as precipitation aggregated at hydrological (longer period) year was selected in NDVI and MTCI. In fact, cork oak is a very resilient tree concerning water deficit, due to the characteristics of its root system which is profound and well developed [42] and fits into a drought avoidance strategy, conferring trees' high adaptability under scarce and variable water resources [43]. On the other hand, in the cork oak PRI model, the 30-day aggregation period of precipitation was selected (with a negative influence, i.e., higher 30-day precipitation, lower PRI values), indicating that cork oak photosynthetic efficiency is sensitive to water stress.

Herbaceous Vegetation
NDVI was the best model for herbs (R 2 = 0.68), in which the selected variable was that with high correlation values (r = −0.82) with radiation at 90 days, negatively affecting the "greenness" of annual herbs. Although NDVI and MTCI showed the same correlation trends (in terms of being positive or negative) for cork oak and shrubs, in herbaceous vegetation, that does not happen (Figure 8). We attribute that result to the dry season (summer which is prolonged into autumn).We observed that, in winter and spring, these two vegetation indices follow the same trend for herbs, but when it comes to summer, when the Leaf Area Index (LAI) of herbs is diminishing until it almost disappears when the rainy season starts, NDVI and MTCI exhibit the opposite trend, i.e., when NDVI decreases due to an absence of herbs, MTCI increases ( Figure 6). As MTCI, a proxy for the chlorophyll content of vegetation, we attribute this result to the soil reflectance in the dry season ( Figure 7). Accordingly, with this result, we consider that MTCI is not a good proxy to track herbaceous vegetation phenology in summer due to the bias of soil reflectance that affects MTCI in this period of the year. This soil background sensitivity of MTCI has been reported by other studies where, despite a strong correlation between chlorophyll content and MTCI, it was significantly weakened when the Leaf Area Index (LAI) was low [44]. For the NDWI model, precipitation at the 90-day aggregation period and hydrological minimum temperature were selected. These variables have correlations of r = 0.71 and r = −0.42, respectively, with NDWI. Thus, according to our results, the water contents of herbaceous cells react slowly in a negative way for variations at minimum temperature (higher minimum temperature; smaller NDWI values) and react positively to precipitation at 90 days. In the PRI model for herbaceous vegetation, hydrological relative humidity (r = −0.3) and radiation at 30 days (r = 0.79) are important for model carotenoid/chlorophyll ratio.

Shrubs
Overall, in all models, cistus obtained the best results with NDWI (R 2 = 0.78). It showed high sensitivity to high temperatures at 90 days (affecting negatively) and hydrological precipitation (affecting positively) mainly in NDVI and NDWI, although in PRI, short-term precipitation was selected. MTCI models did not select any of these climatic variables. Instead, hydrological relative humidity and short-term radiation were chosen. Similar results were achieved by other studies where the results showed how adapted cistus is to the Mediterranean precipitation seasonality [12]. In ulex, the best model was for MTCI (r = 0.38), despite it being very weak. Temperature and precipitation were important to model this vegetation type. The observed PRI peaks in ulex corresponded to the bloom of the yellow flowers, which had a distinctive spectral signature. The blooming stage lasts from mid-January until March when flowering is at its highest [45].

Cork Oak
The evergreen nature of cork oak was reflected in the NDVI trend, in which very low amplitudes of values were observed ( Figure 5A). This "greenness" maintenance throughout the year is reflected in the consistency of the NDVI values. A small variation in NDVI, related to the leaf renewal of cork oak, occurs in late spring/beginning of summer, which has also been observed in other studies [46]. Our results revealed no statistical differences among seasons ( Figure 4A), nor significant difference in the inter-annual NDVI amplitude ( Figure 5A). This seasonal stability of NDVI was used in previous studies to derive long-term trends in the productivity of cork oak woodlands in Portugal [23]. Regarding NDWI, the observed low amplitude of values in cork oak ( Figure 5B) might be because trees minimize the loss of water in leaves with strong stomatal limitations when facing drought [43,46]. In seasonal analysis ( Figure 4B), cork oak NDWI showed differences between winter and summer. MTCI is the vegetation index that showed the highest amplitude of value variation ( Figure 5C), showing the same pattern of NDVI but with a high amplitude of values (high variability). This might be due to the higher sensitivity of MTCI to chlorophyll content compared with NDVI [47]. Seasonally, our results show that summer stands out from other seasons ( Figure 4C). In PRI, cork oak showed a minimum value around April, which quickly started to increase, coinciding with the renewal of the cork oak canopy ( Figure 5D). This result in PRI is corroborated by other studies in which it was observed in late spring/early summer, and the highest values of LUE in accordance with the maximum LAI value of cork [42].

Herbaceous Vegetation
Among all vegetation types, the herbaceous vegetation showed the greatest annual variability in NDVI ( Figure 5A). The two observed well-defined statistically seasonal groups (winter/spring and summer/autumn) for this type of vegetation in NDVI and NDWI are due to the absence of green herbaceous vegetation during the summer. Indeed, the annual growth cycle of the herbaceous species begins with the first rain of the year, the beginning of the hydrological period, and lasts until late spring and early summer when plants set seed and senesce [48]. MTCI, due to soil reflectance in summer and part of autumn, is not a good proxy for herbaceous vegetation trends in the dry season. In PRI, herbs stand out with higher values in summer and the lowest in winter.

Shrubs
Cistus and ulex are two species very well adapted to the Mediterranean climate [12]. These species deal with severe water stress, high temperatures and high irradiance during summer [49]. In our results, cistus had a higher annual and seasonal variability than ulex across all Vis. This difference between these two species may be because, although they belong to the same PFT, they have different adaptation strategies to face climatic variability. Cistus, a semi-deciduous shrub, loses a considerable portion of leaf area to limit water losses by transpiration during the summer months. However, cistus maintains a fraction of leaf area in this period to quickly expand with autumn rains and start to photosynthesize before the other species [12]. On the other hand, ulex adjusts its leaf morphology to avoid drought stress, reducing the transpiring surface area through the modification of leaves into spines, hence reducing the Leaf Area Index (LAI) [12].
Cistus MTCI showed a high amplitude of values for the annual trend, in which the highest values were reached in early winter (December) and the lowest in summer (august). At the same time, cistus PRI showed the highest values in summer and the lowest in winter. These results are corroborated by other studies in which the total chlorophyll content of cistus leaves fell gradually over the summer months and increased in November [14]. As a result of the declining chlorophyll content during summer, the total carotenoid/total chlorophyll ratio increased through the summer, which led to the highest values of PRI in cistus that we obtained in summer.

Conclusions
Our study aimed to analyze how different vegetation types (trees, shrubs and herbaceous vegetation) that coexist in cork oak woodlands are influenced by climate. The methodology adopted relies on vegetation indices that are known to be proxies for several biophysical properties of vegetation. In our study, indices were calculated from spectral data on a monthly basis from canopy-scale field measurements with a spectroradiometer device. The chosen climatic variables (temperature, precipitation, relative humidity and radiation) were aggregated in three time steps (30 days, 90 days and hydrological period) before the spectral measurements to assess the existence of different time lags. Then, we analyzed the vegetation index trends along with the 7-year study and their differences across seasons (winter, spring, summer and autumn) for each vegetation type. Finally, we propose 4 models for each type of vegetation, where each one has the VI as a dependent variable and the climatic variable (NDVI, NDWI, MTCI and PRI) as an independent variable. According to the amplitude of values of VIs (Table 2), the results of this study showed that MTCI is the best VI to assess the phenology of cork oak, cistus and ulex, whereas for herbaceous vegetation, it is NDVI. Although MTCI had the highest range of values in herbaceous vegetation, we found that MTCI is not suitable to track herbaceous vegetation phenology in summer due to the influence of soil reflectance, which is very high in this season (Figure 7). NDVI showed the evergreen nature of cork oak through the non-statistical differences among seasons ( Figure 4A). In herbs, NDVI and NDWI results highlight the statistical differences between winter/spring and summer/autumn ( Figure 4A,B). In cistus, similar differences can be seen in NDVI and PRI ( Figure 4A,D). In ulex, as for cork oak, NDVI does not show a marked annual phenology trend ( Figure 4A). Instead, PRI seems to be a good proxy for the flowering of ulex compared with other seasons ( Figure 4D).
Temperature and precipitation were the variables found to be more important to explain the observed VI variability in the different vegetation types of cork oak woodlands. In some VIs, radiation assumes an important role in modeling, such as in the NDVI and PRI of herbs, NDWI of oak and MTCI of cistus.
In Mediterranean environments, it is crucial to understand the role played by climate, specifically by precipitation and temperature, given the projected scenarios of climate change. Our study contributes to evaluating the ability of spectral vegetation indices to represent vegetation dynamics in Mediterranean ecosystems, promoting the use of remote sensing in ecology studies.
The inclusion of spectral indexes into models can greatly improve the actual estimates of gross primary productivity (GPP) at local and global scales. The linkages between field, airborne and spaceborne measurements are needed to map fine-scale patterns across large areas, over large periods. Field data are key for constraining how to generalize earth observation measurements and models used in vegetation monitoring, especially in heterogeneous forests, where different types of vegetation play different roles and respond to climate factors at different time lags. Our results give an insight particularly useful for assessing the sensitivity of vegetation to climate change and to defining adaptation strategies to these changes.