Nitrogen and Phosphorus Effect on Sun-Induced Fluorescence and Gross Primary Productivity in Mediterranean Grassland

: Sun-Induced ﬂuorescence at 760 nm (F 760 ) is increasingly being used to predict gross primary production (GPP) through light use e ﬃ ciency (LUE) modeling, even though the mechanistic processes that link the two are not well understood. We analyzed the e ﬀ ect of nitrogen (N) and phosphorous (P) availability on the processes that link GPP and F 760 in a Mediterranean grassland manipulated with nutrient addition. To do so, we used a combination of process-based modeling with Soil-Canopy Observation of Photosynthesis and Energy (SCOPE), and statistical analyses such as path modeling. With this study, we uncover the mechanisms that link the fertilization-driven changes in canopy nitrogen concentration (N%) to the observed changes in F 760 and GPP. N addition changed plant community structure and increased canopy chlorophyll content, which jointly led to changes in photosynthetic active radiation (APAR), ultimately a ﬀ ecting both GPP and F 760 . Changes in the abundance of graminoids, (%graminoids) driven by N addition led to changes in structural properties of the canopy such as leaf angle distribution, that ultimately inﬂuenced observed F 760 by controlling the escape probability of F 760 (Fesc). In particular, we found a change in GPP–F 760 relationship between the ﬁrst and the second year of the experiment that was largely driven by the e ﬀ ect of plant type composition on Fesc, whose best predictor is %graminoids. The P addition led to a statistically signiﬁcant increase on light use e ﬃ ciency of ﬂuorescence emission (LUE f ), in particular in plots also with N addition, consistent with leaf level studies. The N addition induced changes in the biophysical properties of the canopy that led to a trade-o ﬀ between surface temperature (Ts), which decreased, and F 760 at leaf scale (F 760leaf,fw ), which increased. We found that Ts is an important predictor of the light use e ﬃ ciency of photosynthesis, indicating the importance of Ts in LUE modeling approaches to predict GPP. air pressure (p, hPa), short wave downwelling radiation (Rin, m − 2 long wave downwelling radiation (Rli, W m − 2 ), air temperature (Ta, ◦ C), wind speed (u, m s − 1 ), soil moisture content (SMC, %), leaf area index (LAI m m − canopy height (h, m), a b cm


Introduction
An accurate estimation of gross primary production (GPP) by terrestrial ecosystems is crucial to understanding the variability of the global carbon (C) cycle [1]. One of the most common ways to estimate GPP relies on the use of light use efficiency (LUE) models (Equation (1)). In the LUE framework [2], estimates of GPP are based on three variables: (i) the fraction of photosynthetically active radiation (fAPAR) absorbed by the vegetation; (ii) the actual light use efficiency of photosynthesis (LUE p ), i.e., the conversion efficiency of absorbed radiation to fixed carbon; and (iii) incident photosynthetically active radiation (PAR).
The development and retrieval methods in passive sensing of sun-induced chlorophyll fluorescence (SIF), i.e., the radiation emitted by plants upon sun exposure, opens new possibilities to estimate GPP using remotely sensed data [3][4][5]. In the last decade, several studies have shown that sun-induced fluorescence at 760 nm retrieved from top-of-canopy (TOC) measurements (F 760 ) can track changes in APAR and LUE p , and therefore can be directly linked to GPP from leaves [6], ecosystem, [7][8][9][10] to regional and global scale [3,[11][12][13].
Although the mechanistic link between GPP and F 760 is not completely understood, recent advances in the field have contributed to explain the process under various conditions [14,15]. The reason F 760 and GPP correlate is that both processes start with the absorption of light by a chlorophyll molecule. Once the photon is captured by the antenna and reaches the reaction center of the photosystem II, the chlorophyll molecule can return to the ground state through photochemical quenching (PQ), through the non-photochemical quenching of the excited state (NPQ), as the photon is dissipated non-radiatively as heat [16], or it can be re-emitted as a photon of fluorescence [17]. Fluorescence emission cannot be physiologically regulated, and its quantum yield depends on the efficiency of PQ and NPQ [17]. The mechanisms regulating the partitioning of absorbed photosynthetically active radiation (APAR) into the different pathways is therefore fundamental to grasping the GPP-F 760 connection [18,19]. F 760 is usually described with a similar approach to the Monteith's LUE framework, as shown in Equation (2): where F 760 is equal to the product of fAPAR, PAR, the light use efficiency of fluorescence emission at 760 nm (LUE f ), and the escape probability of chlorophyll fluorescence at 760 nm (Fesc) [20]. Equations (1) and (2) can be combined into Equation (3), which shows that the only variables that control the relationship between GPP and F 760 are LUE p , LUE f and Fesc: Multiple factors can influence the different terms in Equation (3), and eventually GPP-F 760 relationship [5,8]. Among these, the ones that require more attention because they are not fully understood are: (i) leaf nutrient content, in particular nitrogen (N) and phosphorous (P); and (ii) canopy structural parameters such as leaf area index (LAI) and leaf angle distribution (LAD), which in grasslands are often related to the community structure of the canopy [8,21]. Quantifying the effect of nutrients and canopy structure on the partitioning of absorbed radiation and on LUE p , LUE f , and Fesc is the first step to shed light on GPP and F 760 changes under different nutrient availability.
Canopy N concentration (hereafter N%, N mass per gram of leaves of the whole canopy) is often related to the nutritional condition where the plant grows. Nitrogen is a fundamental constituent of leaves that is typically associated with higher LAI, and positively correlated with the amount of chlorophyll a and b (Cab) [22]. Higher LAI and Cab increase APAR, but at the same time should reduce Fesc due to higher absorption and scattering of emitted fluorescence [14]. Nitrogen is also positively related to the amount of ribulose-l,5-bisphosphate carboxylase and oxygenase (Rubisco) protein content [23,24], and thus the maximum carboxylation rates (Vcmax), which is a key determinant of the maximum photosynthetic rates, and therefore GPP [25]. Therefore, nitrogen can influence the partitioning of APAR into PQ, NPQ, and fluorescence emission [15], but different studies, mainly at leaf level, showed contrasting results [14,26]. Moreover, there is a lack of studies that investigate at canopy scale how LUE p , LUE f , and Fesc are modulated under varying nitrogen availability [14]. Canopy phosphorous concentration (hereafter P%) is another critical element for photosynthesis, being involved in the synthesis of Adenosine triphosphate (ATP) [27]. Leaf-level studies with active fluorescence measurements showed that P% deficient plants have lower chlorophyll fluorescence emission efficiency [28]. However, we are not aware of canopy level studies showing the effect of P% on F 760 and LUE f .
Canopy structural variables, such as LAI and LAD, influence the radiative transfer of incoming radiance and emitted SIF within the canopy [19]. LAD can vary on daily and seasonal bases and is strongly influenced by species composition and plant functional forms [29]. LAI and LAD can have a major influence on the sun/shaded leaf ratio through the canopy. This ratio has the potential to directly influence the level of NPQ in the canopy [30] (higher in sunlit, lower in shaded leaves) and therefore could indirectly influence the LUE f . Canopy structure, through absorption and scattering of the fluorescence emitted by the leaves, has a significant influence on observed F 760 , determining Fesc, the probability of F 760 to escape the canopy [31]. Absorption by chlorophyll is higher in the red region, whereas multiple scattering in the far-red region increases the probability of absorption by soil and woody elements. It has been shown recently with modeling studies that TOC observed F 760 (canopy scale) is only a fraction of the F 760 emitted at leaf scale (F 760leaf ) [32]. The decoupling between F 760leaf and F 760 , mainly mediated by Fesc, can have implications for the GPP-F 760 relationships. Recently, new methods to estimate Fesc are being developed, potentially allowing to downscaling the F 760 signal at the leaf level [31,33]. Finally, other variables such as soil moisture or surface temperature (Ts) also have the potential to impact the GPP-F 760 relationship. Heat and water stress have been proven to affect photorespiration, but not the PQ in Mediterranean species [34], thus decoupling photochemistry from F 760 [18]. Ts, in particular, contains information on both the activation of NPQ mechanisms and other processes related to stomatal closure and sensible heat losses [35]. Therefore, surface temperature might also help to better characterize the seasonal variations of LUE p and therefore to better predict GPP, in particular under stress conditions [35,36]. Figure 1 illustrates a theoretical framework that sums up current knowledge and our hypothesis regarding the interlinks between GPP and F 760 and their relationship with canopy structural parameters and leaf traits of vegetation. In Figure 1, solid colored lines represent the energy partitioning at both leaf and canopy level and dotted lines represent the hypothesized relationships.
All factors illustrated in Figure 1 play a role in determining GPP, F 760 , and their relationship. However, the strength of these influences, and whether leaf nutrient content and canopy structure influence the GPP-F 760 relationship directly (through LUE p , LUE f and Fesc) or occur indirectly (mediated by APAR or by a third variable), is not clear. In this study, we aimed to fill the gap in understanding on how nutrients and canopy structure control LUE p , LUE f and Fesc, and we investigated the mechanisms that drive GPP and F 760 in a nutrient manipulation experiment. We asked the following questions: How do the treatments (N, NP, and P) influence LUE p , LUE f , and Fesc?
What are the drivers of the light use efficiency equations terms (LUE p , LUE f , Fesc) that relate GPP and F 760 ?
What are the direct and indirect effects of nutrients (in particular N%) and canopy structure on GPP and F 760 ?
To answer these questions, we used GPP, F 760 , and additional data on vegetation properties from a nutrient manipulation experiment in Mediterranean grassland with addition of N, P and N and P together (NP). The aim of the fertilization was to induce a changed in both plant nutrient content and structural traits (through changes in LAD mediated by plant community and LAI) within the ecosystem. Figure 1. Energy partitioning at the leaf and canopy level representing the processes involved in the photosynthetic light use efficiency model (GPP = APAR x LUE p ) and fluorescence light use efficiency model (F 760 = APAR *LUE f * Fesc) are represented with solid arrows. Dotted arrows represent the hypothesized relationship between leaf traits, canopy structure and the various processes related to the allocation of energy and transfer of SIF within the canopy. Photosynthetic active radiation (PAR); absorbed (by vegetation) photosynthetic active radiation (APAR); PAR absorbed by chlorophyll a and b molecules (APAR green ), represented as the green bar in the equations on both sides of the figure; gross primary production (GPP); sun-induced fluorescence emitted by all leaves at 760 nm (F 760leaf ); sun-induced fluorescence at 760 nm observed at top of canopy (F 760 ); nitrogen concentration on a mass basis (N%); chlorophyll a and b on a mass basis (Cab); leaf mass per area (LMA); maximum carboxylation rate (Vcmax); leaf area index (LAI); leaf angle distribution (LAD).

Experimental Site
The study was conducted in a Mediterranean savannah located in Spain (39 • 56 24.68 N, 5 • 45 50.27 W; Majadas de Tietar, Caceres) characterized by a continental Mediterranean climate, with temperate winters and warm dry summers: mean annual temperature of 16.7 • C and annual precipitation of~650 mm distributed mainly between September and May [37].
The herbaceous layer is dominated by annual C3 species of the three main functional plant forms, namely grasses, forbs and legumes, which are green and active from October to end of May [38].
The site is managed as a typical wood pasture (Iberian Dehesa) with low intensity grazing by cows (~0.3 cows ha −1 ) [37].

Nutrient Manipulation Experiment, Gross Primary Production and Ancillary Data
A nutrient manipulation experiment focused on the herbaceous layer was established in early spring 2014 and 2015. The set-up consisted of four 20 m × 20 m width randomized blocks. Within each block we established four plots (9 m × 9 m) with 2 m of buffer between treatments ( Figure S1). We established four treatments (for details, see [37]): control (C) with no fertilization, N addition with one application of 100 kg N ha −1 as potassium nitrate (KNO 3 ) and ammonium nitrate (NH 4 NO 3 ), P addition with 50 kg P ha −1 as monopotassium phosphate (KH 2 PO 4 ), and nitrogen-phosphorous (NP) addition, with 100 kg N ha −1 and 50 kg P ha −1 as NH 4 NO 3 and KH 2 PO 4 , respectively.
Carbon Dioxide (CO 2 ) fluxes between the herbaceous layer and the atmosphere were measured in 32 collars of 60 cm × 60 cm for each field campaign around noon local solar time (Table 1). At each collar, GPP (µmol CO 2 m −2 s −1 ) was estimated as the difference between net ecosystem CO 2 exchange (NEE) measured with transparent chambers and ecosystem respiration (Reco) measured with opaque chambers. Measures CO 2 and water vapor mole fractions (W) were collected at 1 Hz by means of an infrared gas analyzer (IRGA LI-840, Lincoln, NE, USA) connected to the chambers. The flux calculations and corrections were conducted using the self-developed R package "RespChamberProc" (https://github.com/bgctw/RespChamberProc). Air temperature (Ta, • C) was measured with a thermistor probe (Campbell Scientific, Logan, UT, USA). Soil moisture content (%) at 5 cm depth was determined with an impedance soil moisture probe (Theta Probe ML2x, Delta-T Devices, Cambridge, UK). Vapor pressure deficit (VPD, hPa) was computed using relative humidity and Ta. Incident PAR (µmol m −2 s −1 ) was measured with a quantum sensor (Li-190, Li-Cor, Lincoln, NE, USA) mounted outside of the chamber. Surface temperature (Ts, • C) was measured with infrared thermometer installed in the chambers (Tc, IRTS-P, Apogee, UT, USA). The meteorological conditions for each field campaign are reported in Table 1. Destructive sampling of the vegetation in four parcels (0.25 m × 0.25 m each) within each plot was conducted to estimate LAI and green to dry biomass ratio [37]. The abundance of each functional group such as fraction of graminoids (%graminoids), forbs (%forbs), and legumes (%legumes) was determined. The Shannon biodiversity index (H) among plant functional types was determined as in [39]. N% and P% in plant tissues were determined as described in [37]. Carbon isotopic signature (δ 13 C) for the vegetation was determined from dried samples using a DeltaPlus isotope ratio mass spectrometer (Thermo Fisher, Bremen, Germany) coupled via a ConFlowIII open-split to an elemental analyzer (Carlo Erba 1100 CE analyzer; Thermo Fisher Scientific, Rodano, Italy). δ 13 C was calculated using the measured ratio between 13 C and 12 C in the sample and in a calibrated in-house-standard (Acetanilide: −30.06 ± 0.05% ) as in [40,41] (Equation (4) and Figure S2): where 13R sample and 13R standard are 13 C/ 12 C ratio of the sample and of the standard, respectively.

Transpiration Estimates
Two independent estimates of transpiration (expresses as latent heat fluxes, LE) were obtained: one from upscaling the δ 13 C measurements (LE ISO ) and the other from the runs of SCOPE optimized at the experimental site [42] to obtain the LE of canopy component (LE canopy,inv ).
LE ISO was calculated from δ 13 C, GPP and VPD according to Equation (5) [43], and then the units were converted from mmol H20 m −2 s-1 to W m −2 .
where VPD mean is the mean daytime VPD computed over the period between the beginning of the growing season and the plant sampling dates for the isotope measurements, and intrinsic water use efficiency (WUE i ) was calculated as following: where Ca is the CO 2 mole fraction in ambient air, b' is the mean fractionation during carboxylation and internal transfer (−27% ), a is the fractionation during diffusion through stomata (4.4% ) and ∆ lin is the community weighted mean of δ 13 C. Figure S3a,b displays LE ISO and LE canopy,inv , respectively, and Figure S3c shows the scatterplot of the two estimates. The two independent estimates have a good relationship, with Pearson correlation coefficient (r) of 0.701 and slope of 0.809. In Figure S3a there are no significant differences among treatments for each campaign in 2014 or 2015 in LE ISO . According to the ANOVA test, the LE canopy,inv shows significant differences in Campaign 2 in 2014 (F 3 , 11 = 11.4, p = 0.01) and the Tukey HSD post hoc-test identifies the P treatment as significantly different from the C treatment (p = 0.012). In addition, in 2015, in Campaign 7, there is a significant difference (F 3 , 10 = 5.47, p = 0.017) and the Tukey post-hoc identifies a significant decrease for N and P treatments in comparison with the control (p = 0.016, p = 0.042, respectively).

Field Spectroscopy, Retrieval of Sun-Induced Fluorescence and Biophysical Properties
TOC spectral radiances were collected under clear-sky conditions immediately before flux measurements at each collar [8,37]. The sampling strategy was designed to minimize the differences in solar zenith angle (SZA) between measurements, confirmed by the ANOVA test, which reports non-significant differences in SZA between treatments in each campaign (p = 0.43, p = 0.41, p = 0.33, p = 0.65, p = 0.99, p = 0.99, and p = 0.57 for Campaigns 1-7, respectively). The ranges of SZA for the spectral measurements are reported in Table 1. Two portable spectrometers (HR4000, OceanOptics, USA) were used to estimate chlorophyll fluorescence at the O 2 A band (i.e., F 760 ,) and reflectance in the spectral range 400-1000 nm. The measurements protocol was the following: We first measured the incident solar irradiance by nadir observations of a leveled calibrated standard reflectance panel (Spectralon, LabSphere, USA). We then acquired five measurements of TOC spectral radiances from nadir at 110 cm above the targeted area using bare fiber optics of 25 • of field of view (about 43 cm diameter at the ground, Figure S4). F 760 was estimated by exploiting the spectral fitting method [6]. The spectral interval used for F 760 was set to 759.00-767.76 nm. Albedo 400-900 was calculated from TOC spectral radiances as shown in Equation (7), assuming a Lambertian behavior of the reflected radiance.
where L r is the reflected radiance and E 400-900 is the Irradiance. fAPAR was estimated in three different ways: (i) fAPAR SCOPE was simulated by the process based SCOPE model [44]. (ii) fAPAR RENDVI was based on the established relationship between measured fAPAR and the red edge NDVI (RENDVI) found in maize, soybean and grasslands [45] (Equation (8)).
where RENDVI is calculated as shown in Equation (9): where R NIR and R RE are reflectance factors in spectral bands 770-800 nm and 700-710 nm, respectively. (iii) APAR Li&Moreau1996 was based on subtracting the integral (between 400 and 700 nm) of the incoming PAR (PAR inc ) from the integral (between 400 and 700 nm) of the reflected PAR (PAR refl ) measured by the spectrometers [7,46] and then multiplying by the proportion of canopy absorption (RAPAR) [47] (Equation (10)).
where RAPAR is calculated as: The fAPAR formulations are quite consistent with each other ( Figure S5), and therefore hereafter we use fAPAR RENDVI .

SCOPE Model Simulations
Forward and inverse simulations with the SCOPE model were conducted to assess the robustness of fAPAR, Fesc, and LE ISO derived from field observations. The forward runs model was parameterized using the structural and functional traits derived from the field sampling as well as meteorological and chamber data. Vapor pressure deficit (VPD, hPa), air pressure (p, hPa), short wave downwelling radiation (Rin, W m −2 ), long wave downwelling radiation (Rli, W m −2 ), air temperature (Ta, • C), wind speed (u, m s −1 ), soil moisture content (SMC, %), leaf area index (LAI m 2 m −2 ), canopy height (h, m), chlorophyll a and b content (Cab, µg cm −2 ), dry matter content (Cdm, g cm −2 ), maximum carboxylation rate (Vcmax, µmol m −2 s −1 ) and the parameters to characterize the leaf angle distribution (LAD), respectively, LIDFa and LIDFb, were used to parameterize the model run. SCOPE meteorological drivers were measured along with chamber measurements for the majority; in the case they were not available with the chambers, such as Rin, Rli, p, VPD, wind speed, atmospheric CO 2 concentration (Ca, ppm), and atmospheric O 2 concentration (Oa, ppm), they were derived by linearly interpolating two consecutive measurements around the chambers measurement time collected at the nearby eddy covariance flux tower at 10 min of temporal resolution. Canopy height was estimated in the field with a meter stick in five positions within the measurement collar. Additional parameters such as leaf equivalent water thickness, leaf width, Ball-Berry stomatal conductance parameter and dark respiration rate at 25 • C as fraction of Vcmax were obtained from the literature for C3 grasses [8]. The SZA at the time of the collection of the spectral measurements was used as model input. Soil reflectance spectra were collected in a dedicated field campaign in April 2015 and used for all the runs. Leaf angle distribution was parameterized in SCOPE as in [8] by assuming grasses to be erectophile, forbs spherical and legumes planophiles.
The accuracy of F 760 and GPP simulated with SCOPE (F 760FW and GPP FW , respectively) was evaluated by root mean-squared error (RMSE), slope, intercept, and the determination coefficient (R 2 ) of the linear regression between observed and modeled data ( Figure S6).
Inverse runs of SCOPE against reflectance, F 760 , GPP and thermal radiance, as described in [42], were carried out to obtain LE canopy,inv and Fesc (Fesc inv ).

Calculation of the Light Use Efficiency of Photosynthesis (LUE p ), Light use Efficiency of Fluorescence Emission (LUE f ) and Escape Probability of F 760 (Fesc)
For each plot and campaign, LUE p , LUE f and Fesc were computed. LUE p was calculated as in Equation (12): where GPP is the one measured with the chambers and APAR was calculated as in Equation (13): LUE f was computed as in Equation (14): where F 760 is the TOC fluorescence retrieved and Fesc fw is the escape probability calculated from forward runs of SCOPE and APAR radiance (mW m −2 nm −1 sr −1 ) is calculated from APAR (µmol m −2 s −1 ) as shown in Equation (15).
where 4.6 represent the conversion factor from µmol m −2 s −1 to W m-2 for radiation from 400 to 700 nm [48] and wl is the wavelength interval (300 nm), and π is used to transform irradiance to radiance.
We computed Fesc and F 760leaf in three alternative ways to evaluate their consistency: (i) Combination of forward runs of SCOPE and measured F 760 (Fesc fw ) as shown in Equation (16): where F 760leaf,FW and F 760leaf,fw are fluorescence emitted by all leaves at 760 nm as calculated by the forward SCOPE run (hemispherical and directional, respectively).
(ii) An empirical estimate of Fesc (Fesc emp ) computed according to Equation (17) [33]: NIR V was calculated as in Equation (18), where NIR T is the reflectance at 858 nm.
(iii) An estimation of Fesc using data from a SCOPE inversion (Fesc inv ) (Equation (20)). Fesc inv was obtained from inversion of SCOPE against reflectance, F 760 , GPP and thermal radiance, as described in [42] and was calculated as in Equation (20).
where F 760INV and F 760leaf,INV are the top-of canopy sun-induced fluorescence at 760 nm and sun-induced fluorescence emitted by all leaves at 760 nm as obtained from SCOPE inversion.
The three alternative Fesc and F leaf computed (F 760leaf,fw , F 760leaf,emp , and F 760leaf,inv ) were compared against each other ( Figure S7). The analysis presented below were conducted with all the different estimates of Fesc to evaluate the effect on the results presented. Hereafter, we report only the results obtained with Fesc fw and F 760leaf,fw .

Statistical Analysis
Our statistical analysis consisted of three parts. First, to answer Research Question (i), group differences among treatments were analyzed with Analysis of Variance (ANOVA) and differences among groups were tested with Tukey Honest Significant differences (HSD) post-hoc test. In the case of violation of the assumption of homoscedasticity of residuals, the ANOVA with the Welch's correction [49] and post-hoc analysis with Games-Howell test [50] were used. In addition, an analysis of Covariance (ANCOVA) was used to test if the relationship between GPP and F 760 (canopy scale) and F 760leaf,fw (leaf level) is changing with the treatment and in time.
Second, we addressed Research Question (ii) with the relative importance analysis with "lmg" (Lindeman, Merenda and Gold), a popular approach for quantifying the individual contributions of multiple regressors, assuming linear relationships, as implemented in the R package "relaimpo" [51]. Standard errors were computed by means of bootstrapping (n = 1000 realizations). Independent variables (i.e., predictors) used in the relative importance analysis are N%, %graminoids, %legumes, Ts, LAI, Shannon Biodiversity Index (H) and soil moisture. Additional relative importance analyses were carried out with the surface-air temperature (Ts -Ta) in place of Ts ( Figure S8), as Ts −Ta could be a good proxy for water stress [52].
Third, to answer Research Question (iii), a path analysis was used. The path analysis assumes linearity among variables and the effects are considered additive and not multiplicative. The structural model is based on expected relationships hypothesized and its model structure is shown in Figure S9. The user specifies the model structure, and the method outputs estimates of the path coefficients. The analysis was conducted with the R package "lavaan" [53]. The individual links among variables were evaluated by means of the p-value and standardized coefficient (β). It should be noted that in the analysis we used Ts in place of the reflectance based indexes because: (i) Ts contains information on NPQ [54]; and (ii) Ts is independent from the measurements used to estimate F 760 .
Chi-squared (χ 2 ), comparative fit index (CFI), standardized root mean square of residual (SRMR) and Root Mean Square Error of Approximation (RMSEA) were computed to evaluate the overall accuracy of the models. The standard error of β and of the model fit indices were obtained from bootstrapping the dataset (n =100 realizations). Additionally, to assess the stability of the individual paths across treatments and the robustness of the original model, we made intervention on the dataset by removing from the dataset one treatment and evaluating the impact on the individual β coefficients (Figures S10-S13).

Description of Fertilization Effects on Fluxes, Optical Data, and Vegetation Characteristics
The effect of the fertilization treatment on GPP, LUE p , F 760 , LUE f and Fesc fw is shown in Figure 2. All these variables show a wide variation in time (campaign) and with treatment. GPP is higher in the N and NP treatments in 2014 and more substantially in 2015 during Campaigns 5 (F 3 , 18 = 15.6, p < 0.01) and 6 (F 3 , 26 = 13.1, p < 0.01). LUE p in the N treatment is significantly different from the C treatment only during Campaign 6 (F 3,26 = 2.7, p < 0.05).
F 760 shows a significant increase during Campaign 2 for the NP treatment (F 3,11 = 5.9, p < 0.05) and during Campaigns 5 (for N and NP) (F 3,18 = 13.2, p < 0.01) and 6 (for N,NP, and P) ( . Data are divided among campaigns. Bar graphs represent means and error bars represent 1 standard error. Group differences in (a-h) were analyzed with ANOVA test and individual differences among groups were evaluated with Tukey HSD post hoc test. Group differences in (i,l) were analyzed with ANOVA with the Welch correction and individual differences among groups were evaluated with the Games-Howell post hoc test. "*" refers to a significant difference from the control treatment with p value < 0.05 and "**" refers to a significant difference from the control treatment with p value < 0.01. refers to a significant difference from the control treatment with p value < 0.05 and "**" refers to a significant difference from the control treatment with p value < 0.01.  . Data are divided among campaigns. Bar graphs represent means and error bars represent 1 standard error. Group differences in (e-h) were analyzed with ANOVA test and individual differences among groups were evaluated with Tukey HSD post hoc test. Group differences in (a,b,i,l) were analyzed with ANOVA with the Welch correction and individual differences among groups were evaluated with the Games-Howell post hoc test. "*" refers to a significant difference from the control treatment with p value < 0.05 and "**" refers to a significant difference from the control treatment with p value < 0.01.
All treatments show a significant increase in Albedo 400-900 during Campaigns 5 (F 3,18 = 29.3, p < 0.01) and 6 (F 3,26 = 13.6, p < 0.01) in 2015, but no significant treatment-induced changes in Albedo 400-900 are observed in 2014. Ts shows significant differences in Campaign 6 for the N, NP and P treatments (F 3,26 = 13.5, p < 0.01). LE ISO follows the phenological cycle with lower values in 2015 ( Figure S3a). There are differences in LE ISO among treatments (such as the increase during Campaign 2 of 2014 for N and NP), but these appeared not significant according to the ANOVA. LE ISO estimates are consistent also with independent simulations with SCOPE ( Figure S3c).

Temporal Variability of GPP-F 760 and GPP-F 760leaf.fw Relationship among Treatments
The results of the ANCOVA show that, in 2014, the intercept of the C treatment is significantly different from the other treatments for both F 760 (as shown in previous studies [8,37] and F 760leaf,fw (p < 0.05 and p < 0.05, respectively) ( Figure 4 and Table S1).  In 2015 the intercept is different for the C Treatment for both F760 and F760leaf,fw (p < 0.01 for both) and for the NP treatment with p<0.05 for both F760 and F760leaf,fw. In 2015 for the N treatment there is no significant interaction between F760 and Treatment (Figure 4c), but there is a significant interaction between F760leaf,fw and the N treatment (p < 0.05) (Figure 4d), with significant differences of the relationship GPP-F760leaf,fw. There is no significant effect of the year on the relationship between GPP-F760. For each treatment: p = 0.706, p = 0.323, p = 0.927 and p = 0.992 N, P and NP and C, respectively. Instead when substituting F760 with F760leaf,fw the effect of the year is not significant in the In 2015, the intercept is different for the C treatment for both F 760 and F 760leaf,fw (p < 0.01 for both) and for the NP treatment with p < 0.05 for both F 760 and F 760leaf,fw . In 2015, for the N treatment, there is no significant interaction between F 760 and treatment (Figure 4c), but there is a significant interaction between F 760leaf,fw and the N treatment (p < 0.05) (Figure 4d), with significant differences of the GPP-F 760leaf,fw relationship. There is no significant effect of the year on the GPP-F 760 relationship. For each treatment, p = 0.706, p = 0.323, p = 0.927 and p = 0.992 for N, P and NP and C, respectively. Instead, when substituting F 760 with F 760leaf,fw , the effect of the year is not significant in the treatments C and P (p = 0.659 and p=0.742), but is significant for the NP treatment with p < 0.05, and barely not significant for the N treatment with p = 0.057.
The results of the relative importance analysis for GPP, F 760 , and F 760leaf.fw show the importance of LAI that controls the seasonality of canopy structure and APAR ( Figure S15).
When substituting as predictor Ts with Ts -Ta, we found slightly better results than Ts alone when predicting GPP, F 760 , and F 760leaf,fw ( Figure S8). However, including Ts -Ta does not improve the overall prediction, as the contribution to R 2 of LAI decreases, but the total R 2 remains similar. When predicting LUE p , LUE f , and Fesc, Ts -Ta is a worse predictor of LUE p than Ts (R 2 = 0.28 ± 0.05).
The main predictor of Fesc is clearly %graminoids (R 2 = 0.52 ± 0.03), followed by soil moisture (R 2 = 0.03 ± 0.04) and Canopy N% (R 2 = 0.02 ± 0.02), the latter contributing only marginally. Results of the relative importance analysis for GPP, F760, and F760leaf.fw show the importance of LAI that controls the seasonality of canopy structure and APAR ( Figure S15).
When substituting as predictor Ts with Ts-Ta we find slightly better results than Ts alone when predicting GPP, F760, and F760leaf,fw ( Figure S8). However, including Ts-Ta does not improve the overall prediction, as the contribution to R 2 of LAI decreases, but the total R 2 remains similar. When predicting LUEp, LUEf, and Fesc, Ts-Ta is a worse predictor of LUEp than Ts (R 2 = 0.28 ± 0.05).  Figure 6 shows the output of the path analysis. The results of the final models are displayed as graphs. The overall model fit is evaluated: χ 2 = 129 ± 23, CFI = 0.901 ± 0.03, SRMR = 0.07 ± 0.02 and RMSEA= 0.19 ± 0.02. CFI and SRMR show excellent fit according to Hu & Bentler [55]. In contrast, the RMSEA is higher than expected. RMSEA is part of the parsimony-adjusted fit indexes, which reward low model complexity. Our goal is however to represent a holistic model that includes all the relevant processes and we do not use the path analysis a posteriori as a mean of model selection. Additionally, according to [56], "RMSEA over-rejects true models for 'small' n (n < 250)", which might be the cause of our RMSEA value, as our sample size is 133. Additionally, according to Figure 6. Path analysis displays the role of canopy nitrogen content (N%) and relative graminoids abundance (%graminoids) on the energy partitioning at the leaf and canopy level. Photosynthetic active radiation (PAR); Absorbed by vegetation photosynthetic active radiation (APAR), Fluorescence emission by all leaves at 760 nm calculated by forward runs of SCOPE (F760leaf,fw); gross primary production (GPP), Surface temperature (Ts) and observed fluorescence at 760 nm (F760). The strength of the relationship among variables is expressed by the standardized coefficient (β) of the path analysis. Each standardized coefficient has a standard error obtained from bootstrapping (n=100 times). The width of the arrows is proportional to their standardized coefficient (β). Colored lines (both solid or dotted) represent direct relationships between variables, whereas gray double-headed arrows represent the covariance among variables. Solid and dotted lines indicate significant (p < 0.05) and non-significant relationships, respectively. The width of the arrows is proportional to their standardized coefficient (β). The different colors are introduced to increase Figure 6. Path analysis displays the role of canopy nitrogen content (Canopy N) and relative graminoids abundance (%graminoids) on the energy partitioning at the leaf and canopy level. Photosynthetic active radiation (PAR); absorbed by vegetation photosynthetic active radiation (APAR); fluorescence emission by all leaves at 760 nm calculated by forward runs of SCOPE (F 760leaf,fw ); gross primary production (GPP); surface temperature (Ts); and observed fluorescence at 760 nm (F 760 ). The strength of the relationship among variables is expressed by the standardized coefficient (β) of the path analysis. Each standardized coefficient has a standard error obtained from bootstrapping (n = 100 times). The width of the arrows is proportional to their standardized coefficient (β). Colored lines (both solid or dotted) represent direct relationships between variables, whereas gray double-headed arrows represent the covariance among variables. Solid and dotted lines indicate significant (p < 0.05) and non-significant relationships, respectively. The width of the arrows is proportional to their standardized coefficient (β). The different colors are introduced to increase readability of the standardized path coefficients. The fit by the overall model is measured by means of Chi-squared (χ2), comparative fit index (CFI) and standardized root mean square of residual (SRMR). Figure 6 shows the clear effect of the %graminoids on F 760 . The N and NP treatments significantly affect N% with β of 0.44 ± 0.07 and 0.47 ± 0.08, respectively. N and NP treatments also affect significantly %graminoids with β of −0.27 ± 0.1 and −0.21 ± 0.09, respectively. N% has a significant relationship with four variables: APAR, Ts, GPP, and F 760leaf,fw with β of 0.37 ± 0.05, −0.37 ± 0.06, 0.12 ± 0.03 and 0.10 ± 0.04, respectively. %graminoids significantly affects APAR and F 760 with β of 0.27 ± 0.09 and −0.17 ± 0.02, respectively. The path between %graminoids and Ts is however not significant. APAR significantly influences GPP, F 760leaf,fw and Ts with β of 0.87 ± 0.02, 0.77 ± 0.03 and −0.25 ± 0.06. Finally, F 760leaf,fw and Ts have a significant covariance with β of −0.17 ± 0.04. F 760leaf,fw and GPP have a significant covariance with β of 0.07 ± 0.02 and so do GPP and Ts with β of −0.18 ± 0.03.

Mechanisms behind the Treatment Effect on GPP and F 760 at Leaf and Canopy Scale
Alternative models using different estimates of F 760leaf were tested and we found that the same paths are selected as significant, and the magnitude of the β coefficients are almost unchanged ( Figure S16). This suggests that the path analysis model is not strongly dependent by the estimation type of the fluorescence emission. The results of the intervention removing treatments show that the vast majority of the paths remain constant and significant. The only difference can be seen when removing the NP treatment ( Figure S11), where the links between canopy N and GPP and canopy N and F 760leaf,fw become non-significant.

Discussion
In the following section, we first discuss the treatment effects (N, NP, and P) on the LUE equation terms, second the predictors of LUE p , LUE f and Fesc fw , and third how the nutrient fertilization affects GPP and F 760 through changes in N%, plant community and canopy structure.

Treatment Effect on LUE p , LUE f , Fesc fw
The relative stability among treatments of LUE p , which is significantly different for the N treatment only in Campaign 6 and shows an increase of NP in Campaign 5 in 2015, suggests that our Mediterranean grasslands is quite constrained in its photosynthetic efficiency, and that any nutrient induced changes in GPP ( Figure 2) are mostly modulated by changes in structural parameters such as fAPAR.
The increase in LUE f in the NP treatment compared to N alone suggests a co-limitation of nitrogen and phosphorus on fluorescence efficiency. The role of P on the functional modulation of fluorescence efficiency at canopy scale has not yet been shown in the literature. However, a series of studies at leaf level showed a positive relationship between photochemical quenching and P in leaves as well as an effect of P on active fluorescence measurements [57]; these support the differences in LUE f observed in our study. Our study suggests that P, and in particular the co-limitation N and P, might have an important role on determining F 760 but is not conclusive on the mechanism, and more research is needed to understand the mechanism and also to support the current efforts to include P in terrestrial biosphere and photosynthesis models [27]. The fact that the magnitude of increase of Fesc fw is very similar in N and NP treatments support the idea that N addition is the main factor regulating canopy structure (Fesc fw and APAR). Other works show that N addition strongly impacts canopy structural parameters such as LAI and plant height in a short-grass prairie [58], although there are no studies focused on the effect of N and NP on Fesc.
Overall, the ecosystem responded in the first year to the fertilization, mainly in a functional way (higher LUE f ), whereas, in the second year of fertilization, we observed structurally mediated increase in GPP and F 760 (through higher APAR and Fesc fw ) (Figures 2l and 3d). The structurally mediated changes in 2015, driven by a decrease in abundance of erectophiles plants as the graminoids in the N containing treatments, caused a change in slope in the GPP-F 760 relationship for the N and NP treatment (Figure 4c), which is almost significantly different from the C for F 760 , but significantly different from the C for F 760leaf,fw in the NP treatment (Figure 4d).
The N treatment has proven to affect plant functioning and canopy structure (APAR and Fesc fw ), while P has only a marginal role on the LUE f . For this reason, in the next paragraphs, more attention is paid to the role of N%, together with meteorology and canopy structure, as driver of in LUE p , LUE f and Fesc fw , as well as GPP and F 760 .

Predictors of the Terms of the Light Use Efficiency Equation
Understanding the causes of variability of the parameters of LUE equations (LUE p , LUE f , and Fesc fw is fundamental to exploit remote sensing information such as F 760 for modeling spatio-temporal patterns of GPP [20]. We show that Ts is the main predictor of LUE p , and together with %graminoids is one of the two main predictors of LUE f . Ts is a good indicator of water stress and strongly related to phenology and green fraction of vegetation [59,60], which ultimately relates to temporal variability of LUE p . However, the fact that variables normalized by APAR such as LUE p and LUE f are driven by Ts indicates that it is not only a seasonal effect but rather physiological. In fact, Ts contains also information related to the activation of the xanthophyll cycle responsible for NPQ processes ( Figure S17) that ultimately is related to LUE p and LUE f [18]. Finally, many variables that have the potential to influence LUE p , such as photorespiration and chlororespiration, are influenced by leaf temperature [61], potentially explaining why Ts is being selected. Our results reinforce the idea that Ts should be used as additional input of LUE models aimed at the prediction of GPP [62].
The %graminoids is by far the best predictor of Fesc fw , independently by the method used for the calculation of Fesc. Graminoids are mainly erectophiles [29], because of this particular LAD, most of the fluorescence is emitted laterally and therefore scattered by the vegetation [8]. In this work, we tested different formulations of Fesc fw with consistent results, in particular between the model-based (Fesc fw ) and the data-driven (Fesc emp ) estimates. The fact that %graminoids is a good proxy for the effect of structure on F 760 and Fesc also opens interesting perspective to use F 760 as well as Fesc to assess taxonomic diversity, when diversity is somehow represented by changes in canopy architecture [63].
N% is an additional predictor selected for LUE f and LUE p , although the additional explained variance seems marginal ( Figure 5). N% is positively related to Vcmax [24,64], and under light saturated conditions a higher Vcmax leads to an increase of LUE p and, to less extent, to increase of LUE f [65]. As hypothesized, from this analysis, it appears that the effect of N% on F 760 and LUE equation terms is not direct and, in Section 4.3, we discuss the relationships between N%, canopy structure, and the observed variables.

Mechanisms behind the Treatment Effect on GPP and F 760 at Leaf and Canopy Scale
The effect of canopy structure on F 760 manifests itself mainly through variation in APAR and Fesc fw (Figures 6 and 2i,l, respectively). With the path analysis, we conclude that %graminoids positively influences APAR that leads to an increase of F 760leaf,fw indirectly. Moreover, %graminoids negatively influence Fesc fw . The changes of canopy structure mediated by changes in plant community at plot level are the most important factors controlling the pathway between F 760leaf,fw and F 760 , and ultimately GPP and F 760 .
By analyzing the relationships between different components measured in the manipulative experiment presented here, we were able to disentangle the pathways, mostly unknown [14,20], through which N% influence the different components of the LUE equations. Our results show that the largest effect of N% on fluorescence emission is not direct, but rather mediated by APAR and Ts (Figure 6), which in turn affect F 760leaf,fw .
There are two indirect ways in which N% affects F 760leaf,fw : (i) Higher N% in the green fraction of the vegetation is associated to an increase of photosynthetic pigments and in particular Cab in leaves [64] and in the canopy [22], which ultimately has a positive effect on APAR [15,66]. Increase in APAR causes higher fluorescence emission at leaf and canopy level ( Figure 6) [67]. There are contrasting results in the literature regarding the effect of N% on fluorescence and all the studies conducted at the leaf level [14,15,26]. Our study at canopy level supports the findings in [15] that at varying levels of N available APAR modulates F 760leaf,fw and F 760 , and its relationship with GPP. (ii) N% influences positively F 760leaf,fw through Ts. N% has a negative effect on Ts and F 760leaf,fw exhibits a negative relationship with Ts. The first hypothesized mechanism is related with the observed increased in Albedo 400-900 (Figure 3e,f) associated with the higher N%. The effect of N% on albedo, despite being quite debated in the literature [68,69], has been demonstrated both at canopy scale [70,71] and at leaf level [72] and has to do with the increase in near infra-red (NIR) reflectance that is larger than the decrease of the reflectance in the visible region due to higher Cab and light absorption. Therefore, the increase of Albedo 400-900 with increasing N% results in less available energy in the canopy, which eventually leads to a decrease of Ts if other conditions such as soil moisture and VPD are similar [69,72]. The second has to do with the modulation of transpiration due to the fertilization (Figure 3g,h), which cools down the canopy, as the leaf surfaces lose heat when water evaporates through the stomata. Our estimate of LE ISO show an increase in N and NP treatments during the peak of the growing season, but it is not significant (Figure S3a,b) and lower than the changes in in Albedo 400-900 for N, NP and P, in particular in 2015 (Figure 3c,d). Given the strong response of GPP in the N and NP treatments in 2015 (Figure 2b), the mild change in LE ISO (Figure S3a,b) suggests an increase of water use efficiency, which is backed by δ 13 C measurements, which show a significant increase in the N and NP treatment of Campaign 6 ( Figure S2) (where less negative values correspond to higher WUE [73]). Therefore, we can conclude that, although transpiration might be involved in the regulation of Ts at the peak of the season, biophysical variables such as Albedo 400-900 are much more affected by N% and contribute to reduce Ts.
Given that a large amount of N is invested in Rubisco protein [23], N can impact directly the carboxylation rates. The direct link between carboxylation rates and F 760leaf is not yet clear [74]. However, we found a direct, though weak, relationship between N% and F 760leaf,fw ( Figure 6) that is likely mediated by the ceiling effect mechanism described in the literature in an elevated CO 2 manipulation experiment [19,65], but not yet observed in nutrient manipulation experiments.

Conclusions
This study analyzed and explained the underlying mechanism responsible for the changes in gross primary productivity (GPP) and sun-induced fluorescence at 760 nm (F 760 ), and their relationship, due to a nutrient fertilization with nitrogen (N), phosphorous (P), and the combination of the two nutrients (NP). The nitrogen additions (N and NP) had an effect mainly through changes in absorbed photosynthetically active radiation (APAR) and escape probability of fluorescence (Fesc fw ). Changes in APAR are directly related to changes in GPP and F 760 and are due to the combination of changes in canopy chlorophyll content and in species composition that modifies the canopy structure. Changes in Fesc fw are mainly due to the changes in the abundance of erectophile vegetation with N addition. In the treatment with the addition of N, forbs (non-erectophile) increased while graminoids (erectophile) decreased, which ultimately led to changes in leaf angle distribution and modified the F 760 observed in particular in 2015. This has an effect on GPP-F 760 relationship both across treatments and from year to year. Phosphorous addition had a significant effect on the light use efficiency of fluorescence, in particular when combined with high nitrogen availability. This result points toward the need of better understanding the thus far neglected role of phosphorous on modulating sun-induced fluorescence.
With a path analysis, we also revealed that N% not only affects F 760 indirectly through APAR and Fesc fw , but also is tightly related with surface temperature (Ts). The negative relationship between N% and Ts is biophysically mediated by higher albedo observed after the fertilization, and only marginally physiological mediated by increase in transpiration. We also found a trade-off between F 760 and Ts (likely mediated by the non-photochemical quenching mechanisms), indicating the importance of measuring simultaneously these two quantities. We finally found that Ts is also the main predictor of the light use efficiency of photosynthesis, which is a fundamental parameter to improve the predictability of GPP. In conclusion, our results show that both nutrient availability and their indirect effect on biodiversity are fundamental drivers of sun-induced fluorescence, and its relationship with gross primary productivity. Our results also reveal the interlink among fluorescence, surface temperature and GPP, and support the importance of tandem missions such as the FLuorescence EXplorer (FLEX) and Sentinel-3, providing concomitant estimates of sun-induced fluorescence, vegetation related spectral indices, and land surface temperature.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/21/2562/s1, Figure S1: Aerial photograph of the experimental site (SMANIE). Figure S2: Group differences among treatment of carbon isotopic signature (δ 13 C). Figure S3: The two transpiration estimates and the albedo 400-900 across treatments. Figure S4: Schematic of the radiometric and chamber footprint. Figure S5: Scatterplot of the two APAR estimations. Figure S6: Scatterplot of modeled vs. observed GPP and F 760 . Figure S6: Relationship between F 760leaf from forward runs of SCOPE, inverse runs and empirical estimates. Figure S7: Scatterplot of GPP-F 760 at leaf and canopy scale across treatments. Figure S8: Relative importance analysis of GPP, F 760 , F 760leaf,fw , F 760leaf,inv , LUE p , LUE f , and Fesc fw with Ts −Ta instead of Ts. Figure S9: Set of equations that represent the model structure for the path analysis. Figure S10: Path analysis without the nitrogen treatment. Figure S11: Path analysis without the nitrogen and phosphorus treatment. Figure S12: Path analysis without the phosphorus treatment. Figure S13: Path analysis without the control treatment. Figure S14. Bar graph representing differences among treatments of %graminoids, %Forbs and %Legumes. Table S1: Evaluation of the relationship between GPP and F 760 and between GPP and F 760leaf,fw among different treatments. Figure S15: Relative importance analysis of GPP, F 760 , F 760leaf,fw , F 760leaf,inv , LUE p , LUE f , and Fesc fw . Figure S16: Path analysis with fluorescence emission at 760 nm calculated from SCOPE inversion. Figure S17: Scatterplot of Ts and PRI . Table S1: Evaluation of the relationship between Gross Primary Production (GPP) and Fluorescence at 760 nm (F 760 ) and between GPP and Fluorescence at emission level at 760 nm (F 760leaf,fw ) among different treatments.