Nutrients and water availability constrain the seasonality of vegetation activity in a Mediterranean ecosystem

Anthropogenic nitrogen (N) deposition and resulting differences in ecosystem N and phosphorus (P) ratios are expected to impact photosynthetic capacity, that is, maximum gross primary productivity (GPPmax). However, the interplay between N and P availability with other critical resources on seasonal dynamics of ecosystem productivity remains largely unknown. In a Mediterranean tree–grass ecosystem, we established three landscape‐level (24 ha) nutrient addition treatments: N addition (NT), N and P addition (NPT), and a control site (CT). We analyzed the response of ecosystem to altered nutrient stoichiometry using eddy covariance fluxes measurements, satellite observations, and digital repeat photography. A set of metrics, including phenological transition dates (PTDs; timing of green‐up and dry‐down), slopes during green‐up and dry‐down period, and seasonal amplitude, were extracted from time series of GPPmax and used to represent the seasonality of vegetation activity. The seasonal amplitude of GPPmax was higher for NT and NPT than CT, which was attributed to changes in structure and physiology induced by fertilization. PTDs were mainly driven by rainfall and exhibited no significant differences among treatments during the green‐up period. Yet, both fertilized sites senesced earlier during the dry‐down period (17–19 days), which was more pronounced in the NT due to larger evapotranspiration and water usage. Fertilization also resulted in a faster increase in GPPmax during the green‐up period and a sharper decline in GPPmax during the dry‐down period, with less prominent decline response in NPT. Overall, we demonstrated seasonality of vegetation activity was altered after fertilization and the importance of nutrient–water interaction in such water‐limited ecosystems. With the projected warming‐drying trend, the positive effects of N fertilization induced by N deposition on GPPmax may be counteracted by an earlier and faster dry‐down in particular in areas where the N:P ratio increases, with potential impact on the carbon cycle of water‐limited ecosystems.


| INTRODUC TI ON
Gross primary productivity (GPP), the most important component of the terrestrial carbon cycle (Beer et al., 2010;Ciais et al., 2014), varies greatly over space and time. A better understanding of its temporal variability is crucial to accurately predict carbon fluxes under global changes (Xia et al., 2015). At the seasonal time scale, variability of GPP is controlled mainly by the dynamics of vegetation structure (Richardson et al., 2013), such as leaf area index (LAI), which is monitored using vegetation indexes (VIs) from remote (Zhang et al., 2003) or near-surface sensing measurements (Richardson, Braswell, Hollinger, Jenkins, & Ollinger, 2009). Therein, land-surface phenology (the timing of recurrent biological events such as the start of the green-up or dry-down of canopy), slopes and amplitudes of the seasonal cycle are often extracted from time series of VIs (Gerard, George, Hayman, Chavana-Bryant, & Weedon, 2020;Liu et al., 2016) used as metrics to describe the seasonal dynamics of canopy structural development (Soudani et al., 2012). Similarly, maximum gross primary productivity (GPP max ), representing the carbon sequestration potential of vegetation activity (i.e., photosynthetic capacity) and dominating the interannual variability of global net ecosystem exchange (NEE; Fu et al., 2019), in conjunction with extracted metrics from its time series, can also be used to describe seasonal variation of GPP.
Attention has been paid to canopy development and its phenology in terms of its response to both long-term climate (Cleland, Chuine, Menzel, Mooney, & Schwartz, 2007;Menzel et al., 2006;Wang et al., 2015) and short-term weather (de Jong, Verbesselt, Schaepman, & de Bruin, 2012). Likewise, seasonality of vegetation activity is also strongly affected by climate and meteorological factors (Richardson et al., 2010;Wingate et al., 2015). However, plants simultaneously balance uptake of multiple resources (i.e., carbon, nutrients, and water), and alter allocation to maximize acquisition of the most limiting resource (Bloom, Chapin, & Mooney, 1985;Chapin, Schulze, & Mooney, 1990). This means that both background ecosystem nutrient stoichiometry (e.g., N:P ratio) and the threefold increase in anthropogenic nitrogen deposition (Peñuelas, Sardans, Rivas-ubach, & Janssens, 2012;Reay, Dentener, Smith, Grace, & Feely, 2008) may alter seasonality of vegetation activity and phenology due to variations in nutrient availability (Piao et al., 2019). In particular, as anthropogenic N deposition (Reay et al., 2008) is far in excess of P deposition and mineral weathering (Mahowald et al., 2008;Peñuelas et al., 2012), the observed general shift toward a higher N:P ratio is expected to exert impacts on macro-scale ecosystem functioning and properties (Fay et al., 2015;Janssens et al., 2010). Previous studies show that changing N or P availability can change the timing of leaf unfolding and extension (Yang, Zavišić, Pena, & Polle, 2016) as well as the growing season length (Wang & Tang, 2019). Plants under different soil nutrients stoichiometry differ in uptake and transport of nutrients (Yang et al., 2016), which would further affect canopy structural development (Jing et al., 2017) and photosynthetic rate (Yang et al., 2016). As a result, changes in N:P ratio can alter the periodicity and magnitude of seasonal changes in both structural and functional properties of ecosystems (i.e., canopy structure and photosynthetic capacity). Understanding the links between nutrient availability and seasonal patterns of photosynthetic capacity (e.g., phenology of GPP max ) is thus crucial to understand temporal dry-down period, and seasonal amplitude, were extracted from time series of GPP max and used to represent the seasonality of vegetation activity. The seasonal amplitude of GPP max was higher for NT and NPT than CT, which was attributed to changes in structure and physiology induced by fertilization. PTDs were mainly driven by rainfall and exhibited no significant differences among treatments during the green-up period.
Yet, both fertilized sites senesced earlier during the dry-down period (17-19 days), which was more pronounced in the NT due to larger evapotranspiration and water usage. Fertilization also resulted in a faster increase in GPP max during the green-up period and a sharper decline in GPP max during the dry-down period, with less prominent decline response in NPT. Overall, we demonstrated seasonality of vegetation activity was altered after fertilization and the importance of nutrient-water interaction in such water-limited ecosystems. With the projected warming-drying trend, the positive effects of N fertilization induced by N deposition on GPP max may be counteracted by an earlier and faster dry-down in particular in areas where the N:P ratio increases, with potential impact on the carbon cycle of water-limited ecosystems.

K E Y W O R D S
GPP capacity, nitrogen deposition, nutrient imbalance, phenology, random forest, semi-arid, structure and physiology, tree-grass ecosystem, water availability patterns in biosphere-atmosphere carbon exchange Piao et al., 2019).
Few previous studies have considered nutrient effects on phenology and seasonality of vegetation activity, mostly focusing on leaf-level or individual plants at plot scale (Burner, Brauer, Snider, Harrington, & Moore, 2014;Peñuelas et al., 2013;Serrano, Filella, & Peñuelas, 2000). Those studies are valuable to explore the response of plants to different nutrient availability. However, due to limited sampling frequency, they normally cannot be used to quantify changes in canopy development and vegetation activity at high temporal resolution (e.g., daily). More importantly, they are limited in considering the integrated ecosystem response. For instance, they are limited to quantify the effects of interactions between the ecosystem compartments, for example, the interspecific competition among plants and feedbacks between plants, soil, and meteorology.
Few studies on nutrient and water interaction have been conducted in the semi-arid tree-grass ecosystem (Nomura & Kikuzawa, 2003;Xia & Wan, 2013;Yin, Zheng, Cao, Song, & Yu, 2016). One particularly understudied system is located in Mediterranean climate zones, where typical tree-grass systems have complex structure and changing limiting resources throughout the year, with water limitations in the summer dry season, nutrient limitations in the rainy seasons (Moreno, 2008;Morris, Nair, Moreno, Schrumpf, & Migliavacca, 2019;Nair et al., 2019;. These complex interactions imply that determining the interlinkage between nutrients (i.e., N and P) and water availability in shaping vegetation phenology and its activity in such ecosystems is a major challenge (Piao et al., 2019). Besides, water variability in such seasonally dry systems is important (Forkel et al., 2015;Jung et al., 2017;Zeng, Mariotti, & Wetzel, 2005) as it is the main environmental driver of carbon fluxes, which dominate the interannual variability of global carbon fluxes (Ahlström et al., 2015;Poulter et al., 2014).
In addition, the variability of vegetation activity can be attributed to changes in vegetation structure (e.g., LAI or biomass) and physiological factors (e.g., light use efficiency [LUE];Hu et al., 2018;Richardson, Hollinger, Aber, Ollinger, & Braswell, 2007;Wu et al., 2012). Decoupling the contribution of structure and physiology to vegetation activity at high resoultion (i.e., daily) is helpful to elucidate the effect of nutrient availability on the seasonal dynamics of these ecosystems.
In this study, we used data from a large-scale (24 ha) nutrient (N and P) manipulation experiment (MANIP: https://www. bgc-jena.mpg.de/bgi/index.php/Resea rch/Manip) conducted in a Mediterranean tree-grass ecosystem characterized by a high seasonal dynamic. The experiment was designed to provide an excess of N (mimicking the effects of 10 years of atmospheric N deposition), resulting in imbalanced N:P ratios in the footprint of one eddy covariance (EC) tower (NT) and relieve the imbalance by adding N and P together at the footprint of another tower (NPT). In addition, one tower site was kept as control (CT) without any nutrient addition (El-Madany et al., 2018). A P only treatment was not included because a trial fertilization on the herbaceous layer and the N:P ratio of vegetation show that the ecosystem is rather N-limited than P-limited and it does not respond to P fertilization alone Nair et al., 2019). We jointly used eddy flux observation, satellite data, and digital repeat photography Luo et al., 2018) to define metrics of the seasonality, namely phenological transition dates (PTDs), the start of critical phenological events such as the start of green-up or drydown; amplitudes of the seasonal cycle and slopes in the green-up and dry-down periods as rates of change (Luo et al., 2018). With these high-resolution data and extracted metrics, combining with other measured environmental factors, we aim to answer the following questions ( Figure 1): 1. How do nutrient availability and stoichiometry influence the seasonality of canopy development and photosynthetic capacity (GPP max )?
2. How do water availability and its interaction with nutrients drive the observed differences in GPP max ?
3. Are the observed changes in GPP max after fertilization due to changes in vegetation structure or due to physiological changes?
A schematic representation of the three questions addressed in this study is reported in Figure 1. Considering the characteristics of the Mediterranean climate (rainy autumn-winter, and warm and dry summer) and shape of the time series of canopy structure and photosynthetic capacity (two peaks in late autumn and spring, respectively), we adopted the concept of "hydrological years" as described in Luo et al. (2018)-a year defined from 1 September to 31 August. In this study, we mainly focus on the green-up at the autumn-winter period (green-up), the peak of GPP max (spring peak), and dry-down in the late spring-summer period (dry-down) as those periods are most important and sensitive to detect the ecosystem response to environmental changes (Luo et al., 2018).

| Description of experimental sites and instruments set-up
The three experimental sites are located in Majadas de Tiétar, Cáceres, Spain (39°56′24.68″N, 5°46′28.70″W). It is a Mediterranean tree-grass ecosystem, composed of a dominant herbaceous layer and low density (~20%; El-Madany et al., 2018) evergreen broadleaf oak trees (Quercus ilex L.). The annual mean air temperature of Majadas de Tiétar is 16.7°C. The annual precipitation approximates 650 mm and falls typically from November to May with a very dry summer (Luo et al., 2018;Perez-Priego et al., 2017). An initial dose of phosphorus fertilizer (50 kg P/ ha as triple superphosphate, Ca(H 2 PO 4 ) 2 ) was applied at NPT in successively added into NT (20 kg N/ha) and NPT (10 kg P/ha and 20 kg N/ha), respectively. The total doses of N were approximately 10 times higher than the current N deposition rate in this area . The P addition in the NP site was calibrated to compensate for the N imbalance and maintain the original N:P stoichiometry of the ecosystem's herbaceous layer .
The fertilization scheme was defined on the basis of a fully factorial small-scale fertilization experiment on the herbaceous layer Perez-Priego et al., 2015), which showed an N limitation at the experimental site and negligible response to P alone. This was also confirmed by the vegetation N:P stoichiometry at the site (both herbaceous layer and trees) that indicated an N limitation before fertilization (data not shown here). Moreover, it was logistically difficult to include a fourth treatment without overlap of the EC footprint climatology, given the large-scale fertilization applied (24 ha). For these reasons, we did not include in the experimental design a plot fertilized only with P.
An EC system was installed at 15 m height to measure the carbon, water, and energy fluxes at each site. Each EC system consisted of a three-dimensional sonic anemometer (R3-50; Gill LTD) and an infrared gas analyzer (LI-7200; LI-COR Bioscience) to measure mixing ratios of CO 2 and H 2 O. Additional vertical CO 2 and H 2 O concentration profiles were measured at seven levels (0.1, 0.5, 1.0, 2.0, 5.0, 9.0, and 15 m above ground with an LI-840; LI-COR Biosciences). Meteorological variables such as air temperature (T a ), shortwave incoming radiation (SW Rad ), incoming photosynthetically active radiation (PAR), vapor pressure deficit (VPD), and precipitation (Prec) were also measured at each site. Soil water F I G U R E 1 Illustration of research questions in this study. (1) How do nutrient availability and stoichiometry influence the daily maximum gross primary productivity (GPP max )? (2) How do the crucial environmental factors, specifically water availability, interact with the nutrient availability? (3) Are changes in structure and changes in physiology contributing to the variation in GPP max and at which periods are they contributing? In this schematic plot, the "+" stands for the increase in properties or positive contribution, while "−" represents the decrease in properties or negative contribution. "Yes" means supporting the hypotheses, while "No" represents that against the hypotheses. The brown circles are our original hypotheses content (SWC) was measured at 0.05, 0.1, and 0.2 m below ground were also selected to track the greenness changes in the experimental sites.
We summarized the main variables used in this study that were calculated from the above instruments and their roles in this study to better introduce the workflow in following method parts (Table 1).
Note ecosystem fluxes and properties before fertilization were compared between the three sites Nair et al., 2019). Specifically, pre-treatment measurements both of soil (i.e., soil texture, soil C, N, and P) and leaf N and P of the herbaceous layer indicated no statistically significant difference among treatments . Additionally, El-Madany et al. (2018) evaluated the differences in carbon, energy, and water fluxes among NT, NPT, and CT before applying fertilization and demonstrated that they were not significantly different.
2.2 | Data for the characterization of photosynthetic capacity and canopy structure 2.2.1 | EC data processing and flux partitioning to GPP EC data were collected at 20 Hz and then processed using EddyPro 6.2. The calculated CO 2 fluxes were then quality checked (Mauder & Foken, 2011;Rebmann et al., 2005) and corrected by adding storage fluxes (integrated CO 2 fluxes using seven levels of CO 2 profiles) to compute NEE. The friction velocity (u*) thresholds were detected for each year and tower and u* filtering was conducted following Papale et al. (2006).
The time series of NEE were gap filled using the marginal distribution sampling method (Reichstein et al., 2005). The portioning of NEE into GPP was conducted using nighttime-based methods (Reichstein et al., 2005). The calculation of fluxes was made according to El-Madany et al. (2018) and using the R package REddyProc v1.1.6 (Wutzler et al., 2018).

Research question Roles
Processed GPP max (daily maximum gross primary productivity) Note: The numbers 1, 2, 3 indicate that the variables were used for the analysis targeting for the scientific questions 1, 2, 3, respectively (shown in Section 1).

| Calculation of VIs from PhenoCam and Landsat 8
Digital numbers (DNs) of Red (R DN ), Green (G DN ), and Blue (B DN ) were extracted from each PhenoCam image and averaged over different regions of interest (ROIs, Figure S2). The green chromatic coordinate (GCC) was computed as in Equation (1; Richardson et al., 2007): Green chromatic coordinate was calculated individually from ROIs of trees and grasses and also at ecosystem scale considering grass and tree relative fractions within the tower footprint weighted by 0.8 and 0.2, respectively (Luo et al., 2018).
Atmospherically corrected surface reflectance products from Red (R) reflectance bands were used to calculate normalized difference vegetation index (NDVI; Equation 2) according to Tucker (1979).
Normalized difference vegetation index was further averaged over each selected area to calculate the 16-day NDVI time series in each EC site.

| Processing GCC, NDVI, and GPP max time series
From half-hourly GCC values, we applied a series of steps to calculate the time series of daily GCC (Luo et al., 2018), which allowed us to remove the outliers and retrieve the robust time series: 1. We discarded GCC values when PAR was below 600 µmol m −2 s −1 to filter out GCC measured during adverse meteorological conditions such as clouds or rain (Filippa et al., 2016;Migliavacca et al., 2011).
2. We extracted the 90th percentile of the GCC value from a 3-day moving window and assigned as the GCC value of the day at the center of the moving window (Sonnentag et al., 2012).
3. Daily GCC was further checked and the outliers were removed using a spline-based algorithm (Migliavacca et al., 2011;Richardson et al., 2018).
The daily maximum GPP (hereafter refers to GPP max ) was computed using the same procedure as the one used for GCC and using GPP data between 10 a.m. and 2 p.m. UTC.
To obtain daily NDVI data, the 16-day time series of NDVI were interpolated to daily using the spline method in the zoo R package (Zeileis & Grothendieck, 2005).
The resulting daily GCC, NDVI, and GPP max time series were used to investigate the fertilization effects on canopy structure (GCC, NDVI) and photosynthetic capacity (GPP max ). Direct in situ measurements of LAI for the herbaceous layer (Melendo-Vega et al., 2018) and indirect measurements for the trees (García et al., 2015) were used to calculate ecosystem-scale LAI and compared to Landsat NDVI and GCC ( Figure S3). The comparison between Landsat NDVI and observed LAI ( Figure S5) indicated a good agreement.

| Extraction of seasonality metrics
Three metrics of seasonal development of vegetation activity (hereafter referred as "seasonality metrics") were extracted to determine the changes from time series of photosynthetic capacity (GPP max ) and canopy structure (GCC and Landsat NDVI): PTDs, for example, start of the season (SOS) and end of season; slopes during both the green-up and dry-down periods (slopes); and the amplitudes of the time series in both autumn-winter peak and in spring (amplitudes). Since the time resolution of GPP max and GCC is comparable (daily), we extracted PTDs from these two time series. In contrast, although Landsat NDVI was extrapolated to daily resolution, the revisiting time of Landsat (16 days) makes it difficult to accurately describe the PTDs. Therefore, we decided to use the PTDs only from GCC and GPP. While Phenocams are robust for the determination of timing, they can be problematic for the comparison of the absolute values due to spectral calibration between individual cameras (but see Richardson, 2019). In addition, as cameras are fixed in one direction, their FOV would have some mismatches with EC footprint climatology. Therefore, we chose to use the interpolated daily Landsat NDVI to evaluate the absolute differences of amplitudes in canopy structure across sites. We also calculated slopes from NDVI time series in the green-up and dry-down phase and amplitudes of the seasonal cycle.
For the extraction of seasonality metrics time, series were first smoothed to better capture the seasonal dynamics of GCC/NDVI or GPP (Filippa et al., 2016;Migliavacca et al., 2011). PTDs, slopes, and amplitudes were extracted from the smoothed time series using the procedures described in Luo et al. (2018) and the Supporting Information (Note S1, Table S1, and Figure S4). Uncertainty of seasonality metrics was assessed by running the extraction steps repeatedly (100 times). The 100 times series were constructed by summing original data and random noise generated using the residuals between the spline-smoothed fitting and the observed data (Filippa et al., 2016).

| Water availability indexes
The time series of direct measured SWC well represent the variation of water availability in a small soil area. However, they have some  half-hourly WAI and CSWI at t (WAI t and CSWI t ) are calculated according to Equations (3) and (4), respectively: TA B L E 2 Differences of phenological transition dates (PTDs), slopes, and amplitudes for photosynthetic capacity (GPP max ) and canopy structure (GCC/NDVI) between different fertilized treatments during green-up and dry-down periods a,b

Green-up
Dry-down The differences and corresponding standard errors between treatments were shown in the table (mean ± SD). It is noted that it means minuend has faster-decreasing speed if the difference between the two treatments is negative. c All the amplitudes during green-up (Autumn-Winter) and dry-down (Spring-Summer) were used to conduct the Wilcoxon test between every two treatments. Hence, the results for the amplitudes refer to the whole season.

F I G U R E 4 (a)
Relationship between start of the season (SOS)50 and the date when daily precipitation first exceeds 5 mm (DatePreciover5) or cumulative precipitation reaches 50 mm (DateCumPreci50) in the green-up period. SOS 50 was extracted from daily maximum gross primary productivity (GPP max ) or green chromatic coordinates (GCC). For the definition of SOS 50 , readers can refer to Table S1. where R t and ET t are water recharge and evapotranspiration, whereas S t , P t , and S max are surface water storage, Prec, and maximum allowed storage (S max is set to 5 mm according to the setting of Nelson et al., 2018). Readers can refer to the details of calculating WAI and CSWI in Note S2.
We used the CSWI and WAI for different purposes. The CSWI was used to study the differences in water availability between different treatments as it was calculated using measured ET in the different treatments. In contrast, we used WAI rather than CSWI as one of the environmental drivers for the random forest (RF) analysis in Section 2.5 because it is based only on meteorological data and does not used measured fluxes (i.e., ET) as input driving variables as it would be with the CSWI.
After first calculating the half-hourly WAI and CSWI, then the data during the mid-day period (10:00-14:00) were selected to be consistent with the periods of GPP max and GCC/NDVI and averaged at the daily, weekly, and monthly scale to represent soil water availability for the three treatments.
To investigate how water availability influences the fertilization effects on seasonality metrics among different treatments, we calculated mean WAI during green-up and dry-down periods in each hydrological year. WAI at green-up was averaged from October to December, while WAI at dry-down was averaged during the period of May-July. Higher WAI stands for more water availability in the ecosystems.

| Disentangling the contributions from structure and physiology
The contribution of climatic and structural drivers to changes in GPP max could be inferred from the time series of meteorological data and LAI (here represented by NDVI). On the contrary, the quantification of the variability in physiological response is more complex. To disentangle the contribution of structure and especially physiology to the changes of GPP max after fertilization, we designed two different analyses: (a) we calculated the LUE, that is, the ratio between GPP max and APAR F I G U R E 5 Monthly mean of differences between different treatments for (a) evapotranspiration (ET), (b) conservative surface water index (CSWI) and cumulative differences for (c) ET or CSWI at the same period in each year (March-November). Difference between treatments of nitrogen fertilized and control (NT-CT) and the difference between treatments of nitrogen and phosphorus fertilized (NPT) and CT (NPT-CT) are in green, and blue, respectively. The black vertical dotted lines represent the dates when phosphorus and nitrogen fertilizers were applied. The ET data in July 2014 were removed because of a big gap in the measurements in CT due to a lightning strike we used a RF analysis to more quantitatively disentangle the contribution from structure and physiology to the changes of GPP max after fertilization. Besides, RF can take different confounding factors (e.g., differences between treatments in regard to soil water availability) into consideration, which could not be considered in the LUE analysis.

| Light use efficiency
Light use efficiency is a physiology-related property (Reichstein, Bahn, Mahecha, Kattge, & Baldocchi, 2014) and daily LUE was calculated using Equation 5 (Monteith,1972): where a fraction of absorbed photosynthetically active radiation (fAPAR) was computed by adopting a linear relationship between satellite NDVI and fAPAR in a semi-arid environment (Equation 6; Fensholt, Sandholt, & Rasmussen, 2004). Daily maximum PAR (PAR max ) is the PAR measured for the period retained as GPP max that is described in Section 2.2.3.
A significant linear regression and high correlation (r = .91) between measured green LAI and Landsat NDVI indicate that NDVI could be used to represent the green biomass absorbing radiation ( Figure S5).
Daily maximum gross primary productivity, Landsat NDVI, and LUE were also averaged at different periods: that is, green-up (October-December), spring peak (March-April), dry-down (May-July), and for the whole hydrological year (H-year).

| RF analysis
Three main steps were applied for RF analysis (Figure 2): 1. We trained RFs using daily climate (SW Rad , T a , VPD, and day length [DL]), soil WAI, and Landsat NDVI data to predict daily GPP max (i.e., ecosystem response). RFs were trained for both unfertilized (i.e., for CT) and fertilized conditions (for NT and NPT after fertilization). To evaluate the performance of RFs, we fAPAR = 1.51 × NDVI − 0.40, F I G U R E 6 (a-c) Mean and standard deviation (SD) of the difference in seasonality metrics (phonological transition dates [PTDs], slopes, and amplitudes) between treatments during the green-up (a-c) and dry-down (d-f) period. Red, green, and blue stand for the differences in metrics between Control (CT) and nitrogen fertilized (NT), nitrogen and phosphorus fertilized (NPT) and CT, and NPT and NT, respectively. The gray and orange dots represent the mean difference for the same hydrological year during green-up and dry-down periods, respectively. The vertical lines stand for SD of differences at the same hydrological year. WAI, water availability index (a) (c) conducted K-fold cross-validation (K = 10) for each treatment.
Specifically, the data were randomly split into 10 subsets. Then, one subset was used for validation and the remaining subsets for training the RF. This process was repeated until each of subset had served as the validation set. The root mean squared error, R-squared (R 2 ), and mean absolute error in cross-validation was computed for each treatment. For unfertilized condition, 1631 data points were used for cross-validation, while 1,396 data points were used for cross-validation for N and NP fertilized condition, respectively. The very good performances of RFs in cross-validation indicated their applicability (Table S2) 3. An analysis of the residuals of GPP max calculated between four combinations of the runs conducted at step 2 ( Figure 2) was used to factor out the role of changes in physiology and structure. For example, to evaluate the structural and physiological contribution to changes in the NT treatment, we used the predicted GPP max using (1) the RF trained at the CT and the NDVI of the CT tower; (2) the RF trained at the CT and the NDVI of the NT tower; and (3) the RF trained at the NT tower and the NDVI from the CT. If the residuals of GPP max between (2) and (1) were higher than residuals between (3) and (1), then we concluded that changes in structure (i.e., NDVI) were the most important driver of contributors to changes in GPP max . A more detailed description of methodology and procedures is found in Note S3 and Figure S6.

| Statistical analysis
We explored whether the fertilization effects were significant via conducting difference tests for the seasonality metrics derived from TA B L E 3 Mean of GPP max , Landsat NDVI and light use efficiency (LUE) for different treatments after fertilization at different periods a,b Green-up

Spring peak
Hydro Wilcoxon test was conducted to investigate whether significant differences of GPP max , NDVI, and LUE exist between treatments at different periods (for details, please refer to Section 2.6). The first sign in the brackets indicates whether there is a statistical difference between fertilized sites (NT or NPT) and unfertilized site (CT). The second one indicates whether there is a statistical difference between NT and NPT. p values are as follows: *p < .05, ns for p ≥ .05. The results at green-up and spring peak are highlighted with light green and dark green, respectively, for comparison. the time series of GPP max , NDVI, and GCC between different treatments. The time series of LUE in different treatments were also compared for the same reason. The difference tests applied are described below: 1. Seasonality metrics or LUE before fertilization were subtracted from the ones after fertilization. We used the differences (∆) calculated between pre-fertilization and post-fertilization (Equation 5) to conduct statistical tests to remove the differences before fertilization among treatments (even marginally).
In Equation (7), X represents metrics or LUE, while t 1 and t 0 stand for the post-fertilization and pre-fertilization, respectively.
1. The differences between every two treatments (Δ X T , i.e., NT − CT, NPT − CT, and NPT − NT) were then calculated and a Wilcoxon test was applied to investigate whether the differences are significant or not among different treatments.
Please refer to Note S4 for the details of difference tests for different metrics or LUE.
In Equation (8), Δ X T stands for the differences of ∆ calculated in step 1, between every two treatments. T1 and T2 stand for two different treatments of three, respectively.

| Fertilization effects on seasonality of photosynthetic capacity and canopy structure
Daily maximum gross primary productivity of the three treatments was similar before the fertilization (before November 2014), and the differences between fertilized (NT and NPT) and unfertilized sites (CT) started to be significant after N fertilization in March 2015 (Figure 3). The differences between treatments also showed large interannual variability. For instance, in the hydrological year 2016 (i.e., H2016: September 2016-August 2017, where Prec was 65 mm above the mean Prec; 12% increase compared to mean), there were large differences between treatments. In the hydrological year 2017 (H2017), when the rainfall was quite low in the autumn and winter seasons (229 mm less rainfall, i.e., 53% decrease) compared to the historical mean for the same period), the differences of GPP max between treatments became smaller.
Differences among treatments were also observed for metrics derived from GPP max (Table 2; Figure S8). For PTDs, there were no differences between treatments at green-up, while in both NT and NPT senescence started earlier (17-19 days) than at CT in the dry-down period. Both NT and NPT GPP max showed faster green-up (with slopes being 0.10 ± 0.03 (mean ± SD) and 0.14 ± 0.05 µmol m −2 s −1 day −1 higher than CT, respectively) and faster dry-down (0.12 ± 0.02 and 0.07 ± 0.05 µmol m −2 s −1 day −1 lower than CT, respectively). Within fertilized treatments, NPT had a significantly faster increase in GPP max in the green-up period, and a slower decrease in GPP max at dry-down compared to NT. Regarding the amplitudes, NT and NPT had significantly higher amplitudes compared to CT, while there was no significant difference between NT and NPT ( Table 2).
The canopy structure also differed among treatments ( Table 2). The PTDs extracted from GCC showed no significant differences between treatments both during green-up and dry-down periods. Fertilized sites (NT and NPT) generally had faster changes of NDVI compared to CT. However, there were no significant differences between NT and CT at green-up and between NPT and NT in the dry-down period. Similar to GPP max , amplitudes of NDVI in NT and NPT were significantly higher than in CT.

F I G U R E 8
Contribution from the structure (struct-) and physiological factors (physio-) to the changes in normalized daily maximum gross primary productivity (GPP max ) after fertilization from random forest analysis. The mean difference between nitrogen fertilized (NT) and Control (CT) and the mean difference between nitrogen and phosphorus fertilized (NPT) and CT are displayed. The error bars represent standard errors. Significant differences were determined by the Wilcoxon test. p-values are as follows: *.01 ≤ p < .05, ***p < .001. Note that struct-contribution is adjusted by subtracting the differences between NT − NPT and CT before fertilization F I G U R E 9 Contribution from changes in (a) structure (struct-) and (b) physiological factors (physio-) to the changes in normalized GPP max after fertilization at different periods from random forest analysis. Mean difference and standard error between fertilized sites (nitrogen fertilized: NT or nitrogen and phosphorus fertilized: NPT) and Control (CT) at periods of green-up (October-December), dry-down (May-July), spring peak (March-April), and the whole hydrological year (Hydro-year) were calculated and presented. The differences between NT − NPT and CT are all statistically significant different. The statistical differences between NT and NPT were also tested and the signs of statistical test are shown. p-Values are as follows: ns: p ≥ .05, ***p < .001. Note that struct-contribution is adjusted by subtracting the differences between NT − NPT and CT before fertilization

| Relationship between Prec and PTDs during green-up period
The relationship between Prec and the PTD SOS 50 (the timing when reaching 50% of the amplitude of GPP max , details in Table S1) extracted from GCC or GPP was investigated (Figure 4; Figure S9). SOS 50 was more strongly related to the date when cumulative Prec reached 50 mm (Date CumPreci50 ) than to the date of Prec onset, that is, when daily Prec first exceeds 5 mm (Date Preciover5 ; Figure 4). For instance, in H2016, even though one large rainfall event occurred quite early, the onset of the green-up occurred late and was more related to Date CumPerci50 (Figure 4b).

| Differences of ET and soil water availability between different treatments
The monthly mean of differences of ET and soil water availability between two fertilized sites (NT and NPT) and control, as well as their cumulative differences, were investigated to understand how fertilization affected water availability ( Figure 5). The results showed that NT had larger cumulative ET compared to CT, and the differences between them were largest during the spring (65.7 mm) and early summer (49.1 mm). The calculated monthly CSWI shows the NT treatment had lower water variability compared to CT and the largest difference (−0.71 mm/day) occurs in the summer period ( Figure 5b). The cumulative differences of ET and CSWI (Figure 5c) between NT and CT were higher after fertilization, which is consistent with the trend observed in SWC data ( Figure S12).
Compared to the differences between NT and CT (0.51 mm/day), the monthly mean differences of ET between NPT and CT (0.06 mm/ day) were much smaller ( Figure 5a). Likewise, from monthly mean CSWI ( Figure 5), the differences between NPT and CT were much smaller compared to the difference between NT and CT (e.g., −1.4 mm/month between NPT and CT, while −10.1 mm/month between NT and CT).
This pattern was also observed from the measured SWC; SWC in NPT was similar to that in CT ( Figure S12). The cumulative difference of ET showed NPT used a slightly larger amount of water compared to CT (difference of 18.9 mm/year) for ET, in contrast with a much larger difference between NT and CT (137.2 mm/year). The cumulative CSWI also indicated the difference in water availability between CT and NPT was much smaller compared to the difference between CT and NT ( Figure 5).

| Interannual differences of seasonality metrics between treatments
The differences in seasonality metrics between treatments for different water availability levels are displayed in Figure 6. With the exception of the differences of PTDs during the green-up period, the differences of metrics between treatments were larger with the increase in water availability (higher WAI). At green-up, with WAI increased from 24.3 (in H2017) to 55.3 (in H2015), the mean difference of slopes between treatments increased from −0.008 ± 0.019 to 0.037 ± 0.046 µmol m −2 s −1 day −1 . Likewise, the mean difference of amplitudes between treatments increased from −0.005 ± 0.030 to 0.07 ± 0.033 µmol m −2 s −1 . In the dry-down period, the WAI of different hydrological years changed within a small range (from 88.8 at H2016 to 93.9 at H2015). However, we can clearly observe that the absolute differences of seasonality metrics between treatments increased with the increase in WAI. The mean difference for slopes increased from −0.027 ± 0.033 (H2016) to −0.090 ± 0.042 µmol m −2 s −1 day −1 (H2015), while the mean absolute differences for amplitudes increased from 0.014 ± 0.042 (H2016) to 0.098 ± 0.050 µmol m −2 s −1 (H2015). Even with larger variation and more uncertainty, there was a similar trend for mean difference for PTDs between treatments in the dry-down period (ranges from −0.4 ± 17.3 to −5.6 ± 6.0 day).

| Comparison of GPP max , NDVI, and LUE between treatments at different periods
In general, the fertilized treatments (NT and NPT) had higher GPP max in all periods (Table 3). They were especially significantly higher during the green-up (the differences are 2.4 ± 0.14 [mean ± SE] µmol m −2 s −1 between NT and CT, and 1.9 ± 0.11 µmol m −2 s −1 between NPT and CT, respectively) and the spring peak (the differences are 2.9 ± 0.13 µmol m −2 s −1 between NT and CT, and 2.3 ± 0.15 µmol m −2 s −1 between NPT and CT, respectively). Within fertilized treatments, NT had a slightly higher GPP max but was not significantly different from NPT. Likewise, the NDVI in NT and NPT was significantly higher than NDVI in CT (difference of NDVI is 0.05 ± 0.001 between NT and CT, and 0.06 ± 0.001 between NPT and CT, respectively). The NDVI of NPT was slightly higher than in NT with the highest NDVI observed at the spring peak (0.72 ± 0.005).
After fertilization, LUE increased significantly both at NT and NPT compared to CT (Figure 7; Table 3), and this increase was most prominent in the green-up period (mean ± SE: 0.0020 ± 0.0001) when compared to NT and CT, and 0.0016 ± 0.0001 when compared NPT and CT, respectively; Figure 7). Moreover, apart from H2017 (with very dry autumn and winter), the LUE was highest in the green-up period even if a lower NDVI was observed at this time compared to spring peak (Figure 7).

| Results of RF analysis
The results of the RF analysis showed that most of the contribution to the changes of GPP max in the fertilized treatments can be attributed to changes in structure (Figure 8). In the NT, both structure and physiology contributed to changes in GPP max , while in NPT structure was the main driver.
Contributions to the variation of GPP max from physiology and structure varied during the year (Figure 9; Figure S14).
Consistent with LUE analysis, the RF analysis showed the most prominent contribution from physiology occurred at green-up, while structure both contributed largely at green-up and spring peak (Figure 9). In the dry-down period, contributions from the structure (positive contribution) and physiology (negative contribution) to the increase in GPP max were divergent, which resulted in the lowest total contribution to GPP max (Figure 9; Figure S14).

| D ISCUSS I ON
The overall fertilization effects on GPP max and the contribution from structure and physiology to changes of GPP max are summarized as the schematic plot in Figure 10.

| Water and nutrients are co-mediating the seasonality of vegetation activity
In this part, we discuss the factors driving the photosynthetic capacity as described in Figure 10, and the interactions between nutrients and water availability. We focus on the phenology in the green-up, dry-down periods, and the variation of GPP max amplitudes and the interannual variability. Specifically, the timing of green-up is solely controlled by Prec; the speed of green-up and amplitudes are modulated by the combination of water and nutrients; the earlier and faster senescence in fertilized sites are related to the water limitation and species changes; higher water availability amplifies the nutrients' effects on GPP max .

| Onset of green-up period: Prec is the controller of timing of green-up
Precipitation determines the timing of green-up and the onset of photosynthetic activity ( Figure 4). As the greenness of evergreen oak F I G U R E 1 0 General illustration of fertilization effects on photosynthetic capacity (daily maximum gross primary productivity [GPP max ]). Control (CT), nitrogen-fertilized (NT), and nitrogen-and phosphorus-fertilized treatments (NPT) are in red, dark green, and light blue, respectively. Three different categories of seasonality metrics (phenological transition dates [PTDs], slopes, and amplitudes) were used to evaluate the fertilization effects. The symbol "*" indicates statistically significant differences between treatments. Note that the symbol "<" for comparison of slopes during the dry-down indicated the faster speed of senescence (more negative value). The periods of green-up, spring peak, and dry-down were indicated by light green, dark green, and light brown, respectively. Physiological factors (Physio) mainly contribute to the increase in GPP max in the green-up period, while structure (Struct) contributes to the increase in GPP max both at green-up and spring peak. EOS, start of the season; SOS, end of the season trees remains relatively constant during the green-up period, the timing of green-up at the ecosystem scale is therefore mainly determined by the herbaceous layer (Luo et al., 2018). Some studies conducted on semi-arid regions in North America, Africa, and Europe (Gong, Fanselow, Dittert, Taube, & Lin, 2015;Jolly & Running, 2004;Kurc & Benton, 2010) showed that the onset of green-up is initialized by the Prec, which is in agreement with our results ( Figure S11, green-up is initialized by certain amount of rainfall), while other studies suggest a widespread pre-rainfall vegetation green-up in Africa (Adole, Dash, & Atkinson, 2018). After warm and dry summer, soil moisture especially that in the top soil layers, drops to levels unsuitable for the growth of herbaceous plants (Jolly & Running, 2004;Xu, Medvigy, Powers, Becknell, & Guan, 2016). With the onset of rainfall in autumn, soil moisture increases and grasses start to invigorate once water transport through the roots could sustain the turgor for biological activities in leaves (Condit et al., 2000;Hsiao, 1973). However, we show that the onset of the vegetation occurs not immediately after the first major Prec event, but once cumulative Prec reaches a certain amount, namely c.a. 50 mm (Figure 4). This is possible because both the frequency and intensity of Prec can influence the growth of vegetation (Rishmawi, Prince, & Xue, 2016;Svoray & Karnieli, 2011) through changing soil water availability. One single Prec event could not always provide a water pulse sufficient to trigger the start of green-up or ecosystem productivity (Reynolds, Kemp, Ogle, & Fernández, 2004). By comparing our 4-year record of daily rainfall and SWC in the green-up period ( Figure S11), we found that vegetation becomes active when the SWC in the top 20 cm is around or exceeds 20%.

| Increased speed of green-up and amplitudes due to larger resources usage or higher resource use efficiency
N and NP fertilization result in a faster increase in GPP max during the green-up period and higher amplitudes of GPP max in fertilized sites (Table 2; Figure S9). With increasing Prec and soil water availability in early autumn ( Figure S11), more organic and inorganic nutrients are accessible to plants (Agehara & Warncke, 2005;Zhang & Wienhold, 2002). Previous leaf biophysical and biochemical measurements indicate that leaves rapidly expand and pigments linearly increase in the green-up period (Croft, Chen, Froelich, Chen, & Staebler, 2015;Yang, Tang, & Mustard, 2014). With higher availability of nutrients and observed higher root biomass and root length density in fertilized sites , nutrients like N can be absorbed and used for leaves. In the study site, increasing N in leaf pools and increases in specific leaf area were observed . These increases cause a rise in the maximum carboxylation rate (V cmax ) and chlorophyll .
As a result, the photosynthetic capacity would increase (Fleischer et al., 2013;Tatarko, 2017), which likely contributes to the faster increase in GPP max and biomass (as demonstrated by NDVI) in the fertilized treatments in the green-up period (Table 2). At the same time, with higher total amount of nitrogen, fertilized sites were observed having higher LAI, which contributes to the notably higher amplitude of GPP max and NDVI in NT and NPT (Table 2). This is consistent with the current understanding of nutrient addition that results in increased biomass (Tilman, 1987) and productivity LeBauer & Treseder, 2008).
During the growing season, soil water availability supports the growth of plants (Grossiord et al., 2017;Meza, Montes, Bravo-Martínez, Serrano-Ortiz, & Kowalski, 2018). Compared to CT, more water usage in NT results in larger LE fluxes (i.e., higher ET; Figure 5). Accompanied with high nutrient availability, higher GPP capacity occurs in NT, which conforms the phenomenon observed in previous studies that higher availability of water and nitrogen induce the higher ecosystem productivity (Harpole, Potts, & Suding, 2007;Niu et al., 2009;Wang et al., 2012). In contrast, NPT uses a similar amount of water than CT (Figure S12; Figure 5) while sustaining similar GPP max and even slightly higher greenness than NT (Table 2). This is might due to enhanced water use efficiency observed at the NPT site, which is supported by observation on stable carbon isotope signature (δ 13 C) measured in a fully factorial N and P manipulation of the herbaceous layer of the site (Martini et al., 2019). Therefore, the faster green-up and increased amplitudes in the fertilized sites could either result from larger resource utilization (NT) or from improved resource use efficiency (NPT).

| Earlier and faster senescence in fertilized sites related to larger water usage and species changes
During the dry-down period (May-July), SWC drops dramatically on the topsoil layers due to increasing ET with the rise of air temperature and scarce rainfall (Battista et al., 2018;Luo et al., 2018;Montaldo, Albertson, & Mancini, 2008). Under these circumstances, faster usage of water in the NT accelerates the decrease in SWC on the topsoil ( Figure S11) and furtherly hampers the photosynthesis once the SWC drops below certain amounts .
Moreover, fertilization can induce changes in species diversity (Isbell et al., 2013;Soons et al., 2017) and N fertilization is likely to select early-senescing species (Wang & Tang, 2019). Several species surveys were conducted within the PhenoCam FOV between the peak of the season and the dry-down period in 2017-2018. The results show, with approaching the end of dry-down, the abundance of forbs was lower in the fertilized sites compared to unfertilized sites in 2018 (this phenomenon was also observed in the first year of a factorial nutrient fertilization experiment conducted in the same area; Martini et al., 2019;Migliavacca et al., 2017;Perez-Priego et al., 2015) and this was most significant in the NPT treatment ( Figure S13). As the forbs are generally the latest senesced during the dry-down (G. Moreno, personal communication), the decrease in abundance of forbs in fertilized sites in the dry-down period might lead to earlier senescence and yellowing in fertilized sites (Table 2

| Interannual variability of fertilization effect is modulated by water availability
The effect of the different treatments on GPP max seasonal course and seasonality metrics (PTDs, slopes, and amplitudes) is enhanced in years with higher water availability (Figures 3 and 6).
Apart from the amount of N added, the response of biomass and productivity of the ecosystem to exogenous N can also depend on the intensity of rainfall Lee, Manning, Rist, Power, & Marsh, 2010). With increasing amounts of Prec in rainfall events, N-induced stimulation of GPP max can be the main contributor to the increase in total GPP . This is because the low availability of water might prevent added N from becoming biologically available (Agehara & Warncke, 2005;Lee et al., 2010). Hence, with increased soil water availability, plants can utilize more nutrients to synthesize biomass (Lee et al., 2010), and enlarged differences between fertilized and unfertilized sites are observed ( Figure 5). Effects of amounts of Prec on ecosystem function should be further explored to identify the relationship between Prec and ecosystem response and its interaction with nutrient availability more clearly (Harper, Blair, Fay, Knapp, & Carlisle, 2005;Huxman et al., 2004;Rishmawi et al., 2016).

| Temporal variability
During the green-up, the changes of GPP max values are mainly attributed to the changes in the physiology and canopy structure. In contrast, the changes of GPP max in the spring peak are mainly attributed to the changes in structure (Figure 9). The analysis of LUE also confirms the results from the RF modeling. The largest increases in LUE (physiologically related index) in fertilized sites took place during the green-up period (Figure 7). Lower LUE but higher GPP max at the peak of the season (Figure 7; Table 3) indicates the larger contribution from the structure (i.e., NDVI, as a proxy of LAI and biomass) at this period ( Figure S5). In the green-up period, with more nutrients available at fertilized sites, LAI and chlorophyll increase simultaneously, hence both structure and physiology contribute to the increase in GPP max . At spring peak, leaves are matured and the total concentration of chlorophyll stays constant or even declines gradually (Gitelson, Viña, Ciganda, Rundquist, & Arkebauer, 2005;Yang et al., 2014). In contrast, LAI generally stabilizes along the whole peak period (Croft et al., 2015). Therefore, compared to the green-up period, the contribution from physiology to the increase in GPP max is decreasing and the main contribution can be ascribed to LAI. With absolute higher LAI after fertilization (Albaugh, Allen, Dougherty, Kress, & King, 1998;Gong et al., 2015), GPP max increased most at fertilized sites in spring peak compared to other periods (Table 2).

| Variability between treatments
Increased GPP max is mainly driven by structural changes both in NT and NPT while physiology contributes to the additional increase in GPP max , in particular in NT (Figure 7). The increase in N availability in fertilized sites might promote more synthesis of chlorophyll (Fleischer et al., 2013;Schlemmer, Francis, Shanahan, & Schepers, 2005) that stimulates the photosynthetic activities and increases GPP max (Figure 3) in NT and NPT. However, the different availability of phosphorus in these two treatments might exert different ecosystem responses. The nitrogen-only addition in NT results in higher N:P ratio in plants , which might change the ratio of RNA to proteins in plants (Loladze & Elser, 2011;Peñuelas et al., 2013), potentially producing more photosynthesisrelated proteins (Matzek & Vitousek, 2009), which then contributes to more increase in GPP max compared to NPTcom (especially at peak of season when they have the largest difference). A recent review related to molecular mechanism of N-P interplay in plants (Hu & Chu, 2019) points that high availability of N and P like in NPT can moderately stimulate the growth of plants, and the N-and P-related gene expression. However, only with high availability of N while simultaneously lacking P, plants tend to strongly activate the nitrate responsive and phosphate starvation-induced genes, which intensifies the utilization of nitrate and limiting phosphate in plants (Hu & Chu, 2019). This may explain the higher contribution of physiology to NT in the green-up period. On the other hand, in the dry-down period, earlier and faster senescence in the fertilized sites ( Figure 10) consequently results in fewer physiological contributions to the increase in GPP max compared to CT (Figure 9).

| Implication of this study for global change ecology
Here we show that the interaction of availability of nutrients and water influences canopy structure and photosynthetic capacity.
Water availability not only has an important influence on the timing of green-up but also modulates the effects of nutrients on the magnitudes of photosynthetic capacity response. At the same time, N addition also tends to deplete water faster in the spring peak and dry-down periods, which advances and accelerates the dry-down.
This response is alleviated with N:P ratio recovering after adding P into ecosystem (Table 2). These interactions between nutrients and water availability further draw more attention to the compensation effects of nitrogen addition on carbon budgets. Namely, the positive feedbacks of ecosystem productivity to nutrients addition at green-up and peak of the season may be offset by the decrease in carbon uptake during the dry-down period in water-limited ecosystems. With projected continuous global warming (Diffenbaugh, Swain, & Touma, 2015;Intergovernmental Panel on Climate Change, 2013) and the increase in regional drought (Mazdiyasni & AghaKouchak, 2015;Zhou, Zhang, Park Williams, & Gentine, 2019), compensation effects are expected to be intensified when growing season length is shortened with decreasing water availability (Ma, Huete, Moran, Ponce-Campos, & Eamus, 2015). We therefore advocate for more attention on changes of seasonality of canopy development and ecosystem carbon fluxes, how changes in plants water usage interacts with nutrients availability and N:P stoichiometry in seasonally dry ecosystems. Moreno acknowledges financial support from the grant agreement IB16185 of the Regional Government of Extremadura.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.