Seasonal Patterns in Spectral Irradiance and Leaf UV-A Absorbance Under Forest Canopies

Plants commonly respond to UV radiation through the accumulation of flavonoids and related phenolic compounds which potentially ameliorate UV-damage to crucial internal structures. However, the seasonal dynamics of leaf flavonoids corresponding to epidermal UV absorbance is highly variable in nature, and it remains uncertain how environmental factors combine to govern flavonoid accumulation and degradation. We studied leaf UV-A absorbance of species composing the understorey plant community throughout two growing seasons under five adjacent tree canopies in southern Finland. We compared the relationship between leaf flavonol index (Iflav—repeatedly measured with an optical leaf clip Dualex) and measured spectral irradiance, understorey and canopy phenology, air temperature and snowpack variables, whole leaf flavonoid extracts, and leaf age. Strong seasonal patterns and stand-related differences were apparent in Iflav of both understorey plant communities and individual species, including divergent trends in Iflav during spring and autumn. Comparing the heterogeneity of the understorey light environment and its spectral composition in looking for potential drivers of seasonal changes in Iflav, we found that unweighted UV-A irradiance, or the effective UV dose calculated according to the biological spectral weighting function (BSWF) for plant growth (PG action spectrum), in understorey shade had a strong relationship with Iflav. Furthermore, understorey species seemed to adjust Iflav to low background diffuse irradiance rather than infrequent high direct-beam irradiance in sunflecks during summer, since leaves produced during or after canopy closure had low Iflav. In conclusion, we found the level of epidermal flavonoids in forest understorey species to be plastic, adjusting to climatic conditions, and differing according to species' leaf retention strategy and new leaf production, all of which contribute to the seasonal trends in leaf flavonoids found within forest stands.


Repeated optical leaf trait measurements from understorey plants
Modelling the change in time with flavonol index, Iflav was complex as repeated measurements of the same leaves on a plant would be expected to be autocorrelated, and similar local light conditions at the measurement points could produce high spatial correlations i.e. within stands/measurement points/species. Consequently, to the sake of simplicity to visually compare trends, we decided use loess-based fits (obtained with R function "loess") to data from different stands with 95% confidence intervals (CI). It was difficult to closely fit the trends in Iflav using a simple function. We tried several approaches, for instance fitting 2nd and 3 rd order polynomial functions to trends in Iflav with DOY from different stands. Although all the terms were significant at the 5 % level, R²values were low ranging from 0.31 to 0.54 between different stands thus indicating a poor fit of the model (A2 Figure S1). Furthermore, visual inspection of the residuals plotted against fitted values showed problems with the model fit of these polynomials. Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.72015 0.01119 64.348 < 2e-16 *** poly(DOY, 3)1 -6.29844 0.29736 -21.181 < 2e-16 *** poly(DOY, 3)2 2.77382 0.29736 9.328 < 2e-16 *** poly(DOY, 3)3 2.06010 0.29736 6.928 9.68e-12 *** ---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.  Figure S1. A 3 rd order polynomial function fitted to Iflav data from different stands in 2015.
To investigate whether temperature was an important driver of changes in adaxial Iflav during the consecutive years 2015-2016, we chose minimum daily air temperature out of the colinear weather-related variables as explanatory variable in the model. This decision was based on high correlation with Iflav and expected response of Iflav to low temperature spells. Furthermore, we examined the relationship between minimum air temperature and abaxial Iflav separately for 2016.
We found additive models to give a better fit than linear regression models for both adaxial and abaxial Iflav (A2 Figure S2). The model selection was based on a combination of visual inspection (normality, outliers, residuals plotted against fitted values & residuals plotted against explanatory variables), comparison AIC/BIC and t -and F-statistics when suitable. Out of different models explaining changes in adaxial Iflav, we found the best model to be GAMM with correlations structure (R function corCompSymm) and variance structure (R function VarIdent) allowing for different residual spread in different stands (A2 Figure  S3). However, for the abaxial Iflav models we did not find a satisfactory fit by simply adding a variance structure (R function VarComb structure to allow different spread over values of explanatory variables, as well as in different stands which did not converge) and/or correlation structure (R function corCompSymm or CorAr1) to GAMM. Instead we tested complicated GAMMs including global and/or group-level smoothers with differing smoothness with or without shared trend and/or interaction term. However, the resulting model had still low R 2 -values, and low edf values, although model selection process and visual inspection indicated better fit for smoothing models (A2 Figure S4). Figure S2. Standardised residuals of linear regression model for adaxial (above) and abaxial (below) Iflav with minimum daily air temperature as explanatory variable, plotted against fitted values showing problems i.e. pattern (above) and differences in spread of residuals (below) related to the model.   Figure S4. Actual (filled circles) and predicted (empty circles) values from the final model for changes in abaxial Iflav with minimum daily air temperature as explanatory variable.

Comparison of optically measured Iflav and extracted flavonoids
In this section we tested the consistency of trends derived from the optical leaf-clip method i.e. Iflav compared with whole leaf flavonoid extracts. We used the mean absorbance of leaf extracts from UV region (UV-B, UV-A and combination of UV-B + UV-A), since we found a weak relationship between absorbance of leaf extracts at 375 nm and Iflav. All these regions were co-linear (A2 Figure S5) and fitting them all in the same model was not possible, hence we made separate models for each explanatory variable. The relationship in some measurements was non-linear and fitting a species-specific linear model showed clear patterns in visual inspection of residuals plotted against fitted values, indicating a violation of homogeneity. We fitted GLS with different variance structures (allowing for heterogeneity of residuals among different measurement times, but assuming homogeneity within one), which indicated a need for further extension of temporal correlation structure in some models. If species-specific GLS did not pass model selection process, we tried more complex models, namely GAMM. Nevertheless, some low edf values indicating a linear relationship were found in GAMM, but in most cases the relationship was not linear for all samples within a species (also see results section for possible explanation).
A2 Figure S5. Collinearity between explanatory variables: mean absorbance of leaf extracts calculated over different wavelength regions (UV-B, UV-A, combination of UV-B & UV-A) and absorbance at 375 nm. Final variable is the explanatory variable i.e. optically measured Iflav.
The best species-specific models chosen in the model selection process: Aegopodium podagraria: The linear regression model had clear problems when residuals were plotted against fitted values (A2 Figure S6). Although seemingly the extended GLS gave a good fit in this species, the predicted values were all negative. Since we kept iterating between finding a good model in visual inspection and did not find a solution within these models, we tried more complex GAMM. The selected model allowed smoothers to change over different measurement times (A2 Figure S7). Throughout the model selection process, models with mean absorbance over UV-B region as explanatory variable gave the best fit.

Anemone nemorosa:
Both the linear regression model and the extended GLS showed clear problems in visual inspection of the residuals plotted against fitted values (A2 Figure S8). Tests indicated a better fit with smoothing models over linear models. Throughout the model selection process, the best GAMMs had absorbance over UV-A region or some over UV-B and UV-A region combined as explanatory variable (A2 Figure S9). Compared to the other species, a variance structure allowing for different spread of residuals over the values of explanatory variable and in different stands was better than a variance structure allowing for different spread over different measurement times.

Convallaria majalis:
The linear regression model showed clear problems in visual inspection (A2 Figure S10). Furthermore, all linear regressions had particularly low R 2 -values in this species, indicating a poor fit of the model. GLS with extensions taking into account heterogeneity, did not improve model performance, thus we moved on to using GAMM. Throughout the model selection process, the models with UV-B region or with UV-B and UV-A region combined as explanatory variables were found best suited for C. majalis data (A2 Figure S11). The selected model in this species was not as good as for other species, which might have been related to lower Iflav values in general.    Figure S11. Actual (filled circles) and predicted (empty circles) values from the final model of C. majalis for changes in adaxial Iflav with mean absorbance over UV-A region as explanatory variable.

Hepatica nobilis:
The linear regression model showed problems in the model with visual inspection (A2 Figure  S12). However, the best model was found with GLS including a variance structure allowing for different residual spread over UV-A region and for different measurement times. Throughout, best models for H. nobilis were found with UV-A region as explanatory variable (A2 Figure S13). Figure S12.  Figure S13. Actual (filled circles) and predicted (empty circles) values from the final model of H. nobilis for changes in adaxial Iflav with mean absorbance over UV-A region as explanatory variable.

Oxalis acetosella:
The linear regression model showed patterns in visual inspection (A2 Figure S14) and GLS extended with different variance structures did not fully fix problems related to different spreads over different UV regions. Smoothing models indicated a better fit and thus we used GAMM (although 3 measurement times had edf values of 1). For some reason linear regression models and GLS gave a better fit with models using absorbance over UV-B and UV-A region as explanatory variable, but smoothing models using absorbance over UV-B as explanatory variable were considered best (A2 Figure S15). Visual inspection of the final model was adequate.

Relating understorey Iflav and understorey spectral irradiance
In this section we compared the obtained Iflav values with corresponding spectral irradiance measurements. However, a direct comparison of these data was not possible, since our observational unit was different for Iflav measurements (plant individual) compared to understorey irradiance measurements (measurement point). Instead we used the mean Iflav calculated for each measurement point to match our observational units. We considered this approach to be adequate for our objectives to test whether any spectral component aligned well with the coarser changes observed with Iflav. However, this approach unfortunately excluded for instance the possibility to investigate whether different species had more sensitivity for any spectral components. Since these data were compared with the mean Iflav from the same measurement point, all different understorey positions i.e. shade, sunfleck or position where radiation was transmitted through canopy i.e. leaf position, were compared with the same mean Iflav. Including all understorey positions into the same model was not functional, since excluding any of these positions if found non-significant would not be possible. Hence, we made separate models for spectral irradiance from each understorey position to investigate which of these understorey positions (or rather the spectral irradiance) seemed to best explain changes in Iflav.
The linear regression models made for different spectral components as explanatory variables (A2 Figure S16), either with or without variable stand in the model showed unwanted patterns in visual inspection of the residuals plotted against fitted values. Furthermore, some models also showed patterns when residuals were plotted against explanatory variables. Figure S16. Relationship between mean Iflav and unweighted UV-B irradiance measured in different understorey positions with linear regression trendline (blue) and with loess fit (red). The grey band is 95% confidence interval. Adding stands as variable to the models improved the model fit.

A2
Understorey shade: The selected model for spectral irradiance measured in understorey shade had effective UV dose calculated according to biological spectral weighting function (BSWF) for plant growth (PG, Flint & Caldwell 2003) as explanatory variable. However, the model with UV-A irradiance as explanatory variable was almost similar based on AIC/BIC values. Furthermore, throughout the model selection process best models had PG or UV-A as explanatory variable (A2 Figure S17). However, the difference between these selected models and models using other spectral regions (e.g. PAR) as explanatory variable was not indisputable (R 2 -values: 0.86 vs 0.80). The selected model with spectral irradiance from understorey shade as explanatory variable was better in terms of AIC/BIC, R 2 -value and visual inspection than the models with spectral irradiance from other understorey positions as explanatory variable.    Understorey sunflecks: The selected model for spectral irradiance measured in understorey sunflecks had UV-A irradiation as explanatory variable (A2 Figure S18). Throughout the model selection process, the best models had UV-A irradiation or effective UV dose calculated according to BSWF for plant growth (PG, Flint & Caldwell, 2003) as explanatory variable. However, as seen from A2 Figure S18, the model predicts opposite trend for the Picea abies stand compared to all the other stands, whereby highest UV-A irradiance results in lowest Iflav values which makes the model biologically less reliable and reduces its explanatory value.   Figure S18. Actual (filled circles) and predicted (empty circles) values from the final model for changes in mean adaxial Iflav with unweighted UV-A as explanatory variable. Spectral irradiance is measured from understorey sunflecks.
Understorey leaf position: The selected model for spectral irradiance measured in understorey leaf position had effective UV dose calculated according to BSWF for flavonoid accumulation (FLAV, Ibdah et al., 2002) as explanatory variable (A2 Figure S19). However, the final model seemed to explain changes in mean Iflav poorly compared to spectral irradiance from other understorey positions.