Foliage Biophysical Trait Prediction from Laboratory Spectra in Norway Spruce Is More Affected by Needle Age Than by Site Soil Conditions

: Scaling leaf-level optical signals to the canopy level is essential for airborne and satellite-based forest monitoring. In evergreen trees, biophysical and optical traits may change as foliage ages. This study aims to evaluate the effect of age in Norway spruce needle on biophysical trait-prediction based on laboratory leaf-level spectra. Mature Norway spruce trees were sampled at forest stands in ten headwater catchments with different soil properties. Foliage biophysical traits (pigments, phenolics, lignin, cellulose, leaf mass per area, water, and nitrogen content) were assessed for three needle-age classes. Complementary samples for needle reﬂectance and transmittance were measured using an integrating sphere. Partial least square regression (PLSR) models were constructed for predicting needle biophysical traits from reﬂectance—separating needle age classes and assessing all age classes together. The ten study sites differed in soil properties rather than in needle biophysical traits. Optical properties consistently varied among age classes; however, variation related to the soil conditions was less pronounced. The predictive power of PLSR models was needle-age dependent for all studied traits. The following traits were predicted with moderate accuracy: needle pigments, phenolics, leaf mass per area and water content. PLSR models always performed better if all needle age classes were included (rather than individual age classes separately). This also applied to needle-age independent traits (water and lignin). Thus, we recommend including not only current but also older needle traits as a ground truth for evergreen conifers with long needle lifespan. L.H., F.O.; L.H., writing—review L.H., J.A.;


Introduction
The physiological status of Norway spruce forests in Central Europe is currently affected by several abiotic and biotic stress factors. Some of them are connected with ongoing climate change such as increasing drought [1], especially during the summer, which, in combination with high temperatures, negatively affects spruce growth, particularly at low-elevation stands [2][3][4]. Natural abiotic disturbances including drought, temperature extremes, and windstorms can be followed by biotic ones: bark beetle attacks [5][6][7] or fungal infections [8] being the most common and currently most devastating agents in Norway spruce productive forests in Central Europe (e.g., summarized in [9]). Central European spruce forests also suffer from remarkable soil acidification, a heritage from past heavy air pollution loads resulting in large-scale Norway spruce forest decline [10][11][12].
In the 1990s, a network of small Norway spruce dominated stream catchments known as GEOMON (http://www.geology.cz/geomon/english) was established in the Czech Republic to focus on biogeochemical monitoring of input-output element budgets [13,14]. The concentrations of the main 1980s pollutant, sulphur dioxide, have decreased considerably since the 1990s [15]; the atmospheric acid deposition was reduced and led to the recovery of stream water chemistry at small forested catchments [16]. Although there are apparent signs of ongoing soil recovery (e.g., increasing pH), other parameters (low base cation to aluminum ratios) still indicate unfavorable conditions in the spruce rooting zone [17,18]. The GEOMON network provides a rich variety of data sources related to ecosystem functioning, focusing mainly on soil-stream water chemistry [16,19,20], but recently also supporting the development of remote sensing solutions for forest recovery and biomass mapping [21,22]. This study further expands the breadth of GEOMON datasets by analyzing variability in Norway spruce optical properties and traits related to soil properties of different headwater catchments.
Nowadays, spatio-temporal changes in forest cover and forest health can be monitored via remote sensing optical methods [10,11,[23][24][25][26][27]. Scaling of leaf-level properties to the canopy level is an essential component for the development of airborne and satellite-based tools for vegetation monitoring. Leaf optical properties reflect biochemical composition and structure [28][29][30] and thus provide information about the physiological status of an individual tree or forest stand [25]. Furthermore, canopy optical signals have recently been used for the estimation of soil properties [31,32]. Some attempts have been made using Eucalyptus grandis canopy spectra for soil carbon estimation [33] or spectral chlorophyll indices of Picea rubens and Abies balsamea to predict the dissolved organic C and N of forest soils [34]; though, to our knowledge, such attempts usually had a case-study character.
Coniferous species of the Pinaceae family, such as Norway spruce (Picea abies (L.) H. Karst.) and Scots Pine (Pinus sylvestris L.), are very common and widespread in European forests. Therefore, when using leaf optical properties for forest health assessments, the challenges brought on by conifers' leaf longevity, needle structure, and shoot architecture for vegetation spectroscopy and remote sensing applications should be taken into account [35]. To contribute to these topical issues, in the present study we focused on the effect of needle age on its biophysical and optical properties in Norway spruce. It is well known that conifer needles change their biochemical composition with age, usually accumulating photosynthetic pigments [36,37], secondary metabolites (e.g., phenolic compounds [38,39]) and dry mass in general [40]. On one hand, some mobile elements (e.g., nitrogen, phosphorus, magnesium) can be retrieved from current biochemical compounds and redistributed from older to younger needles in case of soil nutrient scarcity. On the other hand, less mobile nutrients which play preferentially structural roles (e.g., calcium) accumulate with needle age and thus cannot be efficiently re-distributed [41][42][43]. The structural and biochemical changes during needle lifespan are, to some extent, reflected in changes in needle optical properties [35]. In addition to needle age, other factors cause within-crown variability in needle biophysical and optical properties; the vertical light gradient being the most influential [36,44,45] and the cardinal direction branch orientation almost negligible [46,47]. However, in the present study, we focused primarily on the variability in needle traits imposed by needle age in the crown part most related to canopy remote sensing and thus analyzed only the upper sunlit branches in a productive part of a crown.
There is not a consensus in the remote sensing community on how to deal with needles developed in different seasons (needle age classes, NACs) while taking biophysical and leaf level optical measurements as ground truth. Some remote sensing and spectroscopic studies use only one needle age class, usually the youngest one, i.e., the needles developed in the current vegetation season (NAC1, e.g., [48,49]); however, in other cases, more needle age classes are analyzed, e.g., [10,50,51]. On one hand, current year needles (NAC1) and one-year old needles (NAC2) are the most photosynthetically active in the canopy of several conifers [40,[52][53][54] and, therefore, they should well reflect trees' physiological status. On the other hand, the remote optical signal acquired at the canopy level is derived from more needle age classes than only the current and previous-year needles (NAC1 and NAC2). Moreover, with increasing spatial resolution of UAV-borne, airborne, and even satellite acquired data to a sub-canopy scale, knowledge of age dependency of the biophysical and optical needle traits is becoming more important.
The present study aims to provide deeper insight into needle-age related variability of foliage biophysical and optical traits in conifers with long-living needles. The main goal was to evaluate the effect of Norway spruce needle age on biophysical traits predictions based on laboratory leaf-level spectra and to make recommendations on how to deal with needle age classes while taking ground truth. We hypothesize that age-dependence of needle traits influence the power of prediction models. In line with this, we first examined how the measured needle biophysical traits and optical properties change with needle aging, and second, how they vary between sites (characterized by soil properties). Finally, we evaluated the power of partial least square regression (PLSR) models of needle biophysical traits, including different needle age classes separately and together. The secondary goal was to evaluate if separation of Norway spruce-forested catchments is driven by needle biophysical traits or by site conditions (i.e., soil properties).

Description of Study Sites
The forested headwater catchments (median size of 62 ha) dominated by Norway spruce (Picea abies L. Karst) were selected for the study (catchment characteristics are summarized in Table 1, location, and sampling design in Figure 1). All sites are part of the GEOMON network (established in 1994) that primarily focuses on biogeochemical monitoring of various ecosystem elements, atmospheric deposition, and output by surface runoff [13,14,16]. Most of the catchments are covered by managed Norway spruce monoculture over 40 years old; the only exceptions are UDL and UHL sites, which were completely deforested due to heavy acid depositions loads in the 1980s; therefore, the spruce stands there are younger (21-40 years

Soil Sampling and Analyses
Soil sampling was conducted at GEOMON catchments in 2015. At each catchment, 5-15 sampling pits were excavated, forest floor and mineral soil were sampled and processed as described in [55]. Only the results of pooled samples of the fermentation and humus horizons are used in the present study, reported as FH horizon soil properties. Exchangeable cations in soils were analyzed in 0.1 M BaCl 2 extracts by flame atomic absorption spectrophotometry (FAAS, AAnalyst Perkin Elmer 100, Perkin Elmer, Waltham, MA, USA). Furthermore, available Ca, Mg, K (FAAS) and P (spectrophotometry) were determined in Mehlich extracts. Total content of base cations (FAAS) and P (spectrophotometry) were analyzed after acid digestion of ash after sample combustion (350 • C). Exchangeable acidity was determined by titration of the 0.1 M BaCl 2 extracts. Cation exchange capacity (CEC) was calculated as the sum of exchangeable base cations and exchangeable acidity. Base saturation (BS) was determined as the fraction of CEC associated with base cations. The following parameters of FH horizon were used for data analysis: exchangeable acidity-Al+H (meq kg −1 ); base saturation-BS [%]; pH in water and in KCl; exchangeable base cations (BaCl 2 ) of Ca ex , K ex , Mg ex (mg kg −1 ); available Ca avail , K avail , Mg avail and P avail (Mehlich) (mg kg −1 ); and total Ca tot , K tot , Mg tot and P tot (total digest) (mg kg −1 ) and C:N ratio (g g −1 ). Soil organic C and N were measured simultaneously using a Thermo Scientific Flash 2000. For descriptive statistics of selected FH horizon properties see Supplementary Table S1.

Needle Sampling and Analyses
An extensive field campaign at ten GEOMON sites was conducted during 2-23 August 2017. The needle sampling was conducted during sunny days with minimum clouds and was synchronized with acquisition of airborne hyperspectral data (not used in the present study). Actual weather conditions (daily mean temperature, daily mean relative humidity and daily precipitation) are shown in Supplementary Figure S2. In total, 123 Norway spruce trees were sampled, 9-15 trees per site, always with 3 trees adjacent to a soil sampling pit ( Figure 1b). None of the sampled trees exhibited any visual damage at the time of sampling and held at least 6-8 needle age classes, which corresponds to a healthy tree. One branch from the sun exposed production crown part was collected from each tree by a tree climber and transported to a laboratory. Samples for three needle age classes were taken for laboratory spectroscopy measurements and assessment of needle biophysical traits ( Figure 1c): current needle age class (NAC1) represents needles developed in 2017 growing season, thus, less than one year old needles at the time of sampling; previous year needle age class (NAC2) represents needles developed in 2016 growing season-almost two years old needles at the time of sampling and a mixed sample of 4-6 years old needles (NAC4) developed in 2014 and earlier. On one hand, NAC1 and NAC2 represented young photosynthetically active foliage [40,[52][53][54]. On the other hand, NAC4 represented rather older, less active needles. Three-year-old needles (NAC3) were not sampled to keep the number of samples reasonable for processing time. Needle sampling and spectroscopic measurements were always conducted not later than 48 h after the branch cutting.

Needle Spectra Acquisition and Preprocessing
Needle directional-hemispherical reflectance factor (DHRF) and directionalhemispherical transmittance factor (DHTF) between 350 and 2500 nm were measured using an ASD FieldSpec 4 Wide-Res portable spectroradiometer (Malvern Panalytical, Malvern, UK) attached to an ASD RTS-3ZC integrating sphere (Malvern Panalytical, Malvern, UK). The measurements of needle optical properties were done according to protocols used in previous studies [56,57]. Immediately before the spectra measurement, needles were detached from shoots and arranged to specifically designed carriers, which make it possible to place needles close to each other with small gaps to cover the field of view of the integrating sphere [56,57]. The geometry between the sample, light source and spectroradiometer is, therefore, given for reflectance and transmittance measurements by the configuration of sample ports in the integrating sphere. Distance between a needle sample and a light source varies for the reflectance and transmittance configuration, which causes different portions of a needle sample to be illuminated. The illumination zone for the reflectance and transmittance of the sample was captured by taking a digital photograph of a blank paper placed in the illuminated sample port. The proportion of gaps (gap fraction, GF) between needles within the field of view of the integrating sphere was computed from scans of needles in the carriers by applying a mask of the actual light beam position in the sample port. Equations (1) and (2) show the computation of final DHRF and DHTF: R sample and T sample correspond to reflectance and transmittance, respectively, of needles measured in digital number (DN) values, R stray is stray light in DN values, R is is reflectance of the inner sphere wall, WR is a white reference reflectance measured in DN values and GF R and GF T is computed gap fraction for the reflectance and transmittance sample, respectively. All spectra were smoothed out using a filter based on weighted linear least squares and a 2nd degree polynomial model and a filter window width of 41 nm (smooth "rloess" function in Matlab R2018b). Spectra below 450 and above 2000 nm were too noisy and, therefore, were not used for principal component analysis (PCA) and PLSR modeling. Unfortunately, transmittance measurements of needle leaf objects are often affected by gaps between needles placed in the sample port, which produces negative transmittance values in the visible wavelengths (e.g., [58][59][60]). This was also the case of our samples and those samples were corrected by applying an offset computed from the minimum transmittance value in the red region. Therefore, needle transmittance spectra were used only to illustrate differences among NAC1, NAC2 and NAC4 and were not used for trait modeling.

Assessment of Needle Biophysical Traits
Photosynthetic pigments (chlorophyll a and b-Cab; total carotenoids-Car) were extracted in dark conditions, at 4 • C for seven days with N,N-dimethylformamide [61] and their content was determined spectrophotometrically. Pigment concentrations were calculated according to the equations reported by [62] and related to the unit of dry matter. Another set of needle samples was used for fresh and dry weight assessment and computation of water content as the equivalent of mass reduction after needle drying. After homogenizing dry needles, total soluble phenolic compounds were extracted in 80% methanol and determined spectrophotometrically using a Folin-Ciocalteau phenol reagent and gallic acid as a standard [63]. Lignin content was assessed by thioglycolate solubilization [64] and then determined spectrophotometrically at a wavelength of 280 nm using hydrolytic lignin (Aldrich Chemical Company, USA; [8672-93-3]) as a standard. Cellulose content was analyzed after Loader et al. [65] from pellets after lignin extraction, which was repeatedly washed with 17% NaOH, dried and related to the original dry weight. Oven-dried needles were homogenized and used for determination of C and N contents and their mass ratio. Total C and N were determined simultaneously using a Thermo Scientific Flash 2000 analyzer (Thermo Fisher Scientific, Waltham, MA, USA).
Needles used for the optical properties measurements were removed from the carrier, scanned for projection area and their hemi-surface area was calculated according to [36]. Needles were afterward dried at 80 • C to a constant weight and the leaf mass per area (LMA, g cm −2 ) was determined. LMA was used for conversion of other needle traits (Cab, Car, water soluble phenolic compounds, lignin, and cellulose) contents to area-based units.

Statistical Analyses and PLSR Modeling
Variation in leaf traits among the catchments and due to needle age classes were tested by two-way ANOVA with needle age and catchment as fixed factors. Differences were considered significant at α = 0.05. ANOVA was performed by NCSS 9 software (NCSS, LCC Kaysville, UT, USA).
Principal component analysis (PCA) was applied on needle traits and catchment soil properties. Three separate models were created. Each model included needle traits from one needle age class (NAC1, NAC2 and NAC4) and FH horizon soil properties. Needle traits were averaged for three sampled trees adjacent to the individual soil pit. In total, 41 samples for each variable entered the model.
Visualization of leaf optical properties and the differences among NAC1, NAC2 and NAC4 were performed in Matlab R2018b. Differences among needle age classes in spectral range 450-2000 nm were tested by linear mixed-effect model, with nested tree effect. Linear mixed-effect model was performed in R (4.0.3), using the package "nlme" [66].
For the purpose of separating sites and needle age classes from needle spectra, PCA was applied also on needle reflectance and transmittance spectra to reduce dimensionality and reveal if there is any structure in needle spectra related to needle age or catchment. Principal Component Analysis was based on the NIPALS algorithm [67], the model was cross-validated with random 20 segments. PCA was conducted for needle reflectance and transmittance separately, including spectra from all catchments and needle age classes.
Prediction models for selected needle traits were established, using PLSR analysis with the same NIPALS algorithm as the PCA. For each model, the following input variables were used: the measured needle spectra (DHRF between 400 and 2000 nm) represented independent variables X and the needle traits were then dependent variables Y. Using the measured reflectance spectra, the individual PLSR models were computed for numerous needle traits.
The PLSR modeling included the standard calibration and the validation stages [68,69]. To assess the model performance, a random-cross-validation was conducted. The models were evaluated based on the standard measures [70,71], such as: the lowest number of factors (a measure of model robustness), the lowest root mean square error of prediction (RMSEP), and the coefficient of determination (R 2 ). Both PCA and PLSR modeling were performed in Unscrambler 11 (CAMO Analytics, Oslo, Norway).

Needle Biophysical Traits-Age Effects
Needle pigments (Cab and Car) including their ratio Car:Cab, phenolic compounds ( Figure 2) and dry mass (LMA- Figure 3a) on an area basis showed a very strong effect of needle age, accumulating the compounds in older needles. This pattern was significant and consistent at all GEOMON sites. Needle nitrogen content was age-dependent too ( Figure  3e,f). The N content decreased, usually in NAC 4 needles, probably due to dilution by accumulating structural compounds or due to translocation to younger tissues. The only exception in this pattern was the PLB site showing the accumulation of N in older needles instead. The C:N ratio increased usually towards older needles and confirmed the nitrogen dilution effect by accumulating dry mass. There was a group of needle traits independent of needle age: water, cellulose, and lignin contents (Figure 3a-d).   Table 1). (a) Chlorophyll a+b content; (b) total carotenoids content; (c) the ratio of total carotenoids to chlorophylls (Car:Cab); (d) soluble phenolics content in needles at ten GEOMON sites (ANE, CER, LIZ, LYS, NAZ, PLB, POM, SAL, UDL and UHL). NAC-needle age class: 1-current year needles; 2-previous year needles; 4-four year and older needles. Two-way ANOVA, effect of needle age, site, and their interaction (NAC x Site): p-values; *-significant at α = 0.05.

Variation among the GEOMON Catchments in Needle Traits and Soil Properties
Principal component analysis (PCA) showed that the GEOMON catchments differ mainly by their soil properties; the needle traits showed less influence in the model ( Figure 4). Score plots of PCA models based on needle traits for three different needle age classes showed a very similar separation of GEOMON catchments. Roughly, we could distinguish three groups of catchments based mainly on their exchangeable acidity in FH horizon, needle nitrogen content, base saturation, and calcium and magnesium concentration in FH horizon (Figure 4a,c,e). The first group with high exchangeable acidity, high N needle content and low Ca and Mg included three catchments: UHL, CER, LYS. Four catchments-LIZ, UDL, POM, SAL-localized scores position closer to the coordinates' center, exhibiting rather moderate values of exchangeable acidity, N needle content and Ca and Mg in FH horizon. Finally, catchments ANE, NAZ and PLB formed the third group, showing higher C:N ratio, higher Ca and Mg content, higher base saturation and low exchangeable acidity in the FH horizon.    Table 1). For each model, the needle traits from one needle age class were included: (a,b) NAC1; (c,d) NAC2; (e,f) NAC4. Needle traits: Carotenoids-total carotenoids content; Chlorophyll-total chlorophyll content; Lignin-lignin content; Phenolics-content of total soluble phenolics; Water-water content; LMA-leaf mass per area; Car:Cab-total carotenoids to total chlorophyll ratio; N% needle-needle nitrogen content related to dry mass. Soil properties: AlH_ex-total exchangeable acidity; pH_KCl, pH_H2O-pH assessed in KCl and water, respectively; BS-base saturation; Ca_ex, K_ex, Mg_ex-exchangeable base cations (BaCl2); Ca, K, Mg, P available Ca, K, Mg, P (Mehlich); Ca_tot, K_tot, Mg_tot and P_tot-total Ca, K, Mg, P (total digest); C/N-C:N ratio in FH horizon.  Table 1). For each model, the needle traits from one needle age class were included: (a,b) NAC1; (c,d) NAC2; (e,f) NAC4. Needle traits: Carotenoids-total carotenoids content; Chlorophyll-total chlorophyll content; Lignin-lignin content; Phenolics-content of total soluble phenolics; Water-water content; LMA-leaf mass per area; Car:Cab-total carotenoids to total chlorophyll ratio; N% needle-needle nitrogen content related to dry mass. Soil properties: AlH_ex-total exchangeable acidity; pH_KCl, pH_H2O-pH assessed in KCl and water, respectively; BS-base saturation; Ca_ex, K_ex, Mg_ex-exchangeable base cations (BaCl2); Ca, K, Mg, P available Ca, K, Mg, P (Mehlich); Ca_tot, K_tot, Mg_tot and P_tot-total Ca, K, Mg, P (total digest); C/N-C:N ratio in FH horizon. PC1, driven particularly by exchangeable acidity, total phosphorus content, pH, base saturation, and content of base cations (their exchangeable and total pools) accounted for 28% (±1%) of total variation for NAC1, NAC2 and NAC4 (Figure 4). The PC2 accounted for 19, 18 and 16% of variation for NAC1, NAC2 and NAC4, respectively, driven mainly by the needle traits (Figure 4b,d,f). The most influential needle traits in PC2 were LMA, lignin, phenolics, water and cellulose, consistently for all models regardless of needle age. Position of chlorophyll, total carotenoids and their ratio in loading plots varied among models for individual needle age classes: In the model based on needle traits of NAC1, both pigments contributed to both PC1 and PC2 (Figure 4b). In the model based on needle traits of NAC2, carotenoid and chlorophyll content contributed preferentially to PC1 in contrast to their ratio contributing only to PC2 (Figure 4d). Finally, the model based on needle traits of NAC4 exhibited carotenoid and chlorophyll contribution mainly to PC2 with other needle traits (Figure 4f).
The first two components explained 47% of total variation if modeled for NAC1 and NAC2; 43% if modeled for NAC4. In total, 7 components were included in the PCA models and together they explained 86, 86 and 85% of total data variation if modeled from needle traits for NAC1, NAC2 and NAC4, respectively. The percentage of variation explained by PCs 3-7 was very similar for models of all needle age classes.

Needle Optical Properties
Needle optical properties-reflectance and transmittance (DHRF and DHTF)-for individual needle age classes averaged for all GEOMON sites are shown in Figure 5. Both reflectance and transmittance significantly differed among needle age classes except for some spectral intervals ( Figure 5). On one hand, reflectance did not show significant differences among NACs across two broader spectral intervals: 1138-1354 nm and 1679-1737 nm. On the other hand, three narrow spectral bands were equal in needle transmittance (450-487 nm, 669-683 nm, and 1888-2000 nm). Both reflectance and transmittance decreased with needle age in the green part of the visible (VIS) spectrum and similarly in the near infrared (NIR) spectral interval (800-1000 nm). For transmittance, the age-dependent differences in NIR were more apparent and pronounced in the shortwave infrared (SWIR) wavelengths. Needle-age effect on reflectance and transmittance were consistent across all studied GEOMON catchments (see Supplementary Figure S1 for details).

Needle Optical Properties
Needle optical properties-reflectance and transmittance (DHRF and DHTF)-for individual needle age classes averaged for all GEOMON sites are shown in Figure 5. Both reflectance and transmittance significantly differed among needle age classes except for some spectral intervals ( Figure 5). On one hand, reflectance did not show significant differences among NACs across two broader spectral intervals: 1138-1354 nm and 1679-1737 nm. On the other hand, three narrow spectral bands were equal in needle transmittance (450-487 nm, 669-683 nm, and 1888-2000 nm). Both reflectance and transmittance decreased with needle age in the green part of the visible (VIS) spectrum and similarly in the near infrared (NIR) spectral interval (800-1000 nm). For transmittance, the age-dependent differences in NIR were more apparent and pronounced in the shortwave infrared (SWIR) wavelengths. Needle-age effect on reflectance and transmittance were consistent across all studied GEOMON catchments (see Supplementary Figure S1 for details).

Needle Optical Properties-Age and Site Effects
PCA of both reflectance and transmittance spectra showed a separation trend among needle age classes (Figure 6a,c), highlighting the difference between NAC1 and NAC4. In the case of DHRF, the separation of needle spectra due to needle age was weaker than for DHTF. The stronger separation trend among needle age in DHTF confirmed the results of the ANOVA, showing only minor spectral intervals without significant differences among needle age classes (Figure 5b). The first two components PC1 and PC2 accounted for 94 and 95% of overall variation in DHRF and DHTF, respectively. In contrast to needle age, the effect of the catchment on needle spectra separation was almost negligible (Figure  6b,d). shaded area-standard deviation; n = 123 per needle age class. Differences among needle age classes were tested by linear mixed-effect model. p-values (right y-axis) indicate significant (α = 0.05; dark green) and not significant (light green) differences between needle age classes for each wavelength.

Needle Optical Properties-Age and Site Effects
PCA of both reflectance and transmittance spectra showed a separation trend among needle age classes (Figure 6a,c), highlighting the difference between NAC1 and NAC4. In the case of DHRF, the separation of needle spectra due to needle age was weaker than for DHTF. The stronger separation trend among needle age in DHTF confirmed the results of the ANOVA, showing only minor spectral intervals without significant differences among needle age classes (Figure 5b). The first two components PC1 and PC2 accounted for 94 and 95% of overall variation in DHRF and DHTF, respectively. In contrast to needle age, the effect of the catchment on needle spectra separation was almost negligible (Figure 6b,d).

Needle Traits PLSR Modeling-NACs Effect
Results of cross-validation of PLSR models for selected needle traits are presented in Table 2. Models always performed better if all needle age classes were included: they showed higher coefficient of determination (R 2 ), lower root mean square error (RMSE) in absolute values of given trait and lower RMSE expressed as a percentage of the training dataset range (RMSE%), and higher residual prediction deviation (RPD) than models for individual needle age classes. The model for needle specific mass (LMA) showed the highest predictive power (R 2 = 0.79, RMSE% 9% and RPD 2.18). The second-best predicted needle traits were soluble phenolics content (R 2 = 0.73, RMSE% 11% and RPD 1.93), total carotenoid content (R 2 = 0.72, RMSE% 11% and RPD 1.91) and chlorophyll content (R 2 = 0.71, RMSE% 11% and RPD 1.86). The model for water content showed acceptable results (R 2 = 0.60, RMSE% 11% and RPD 1.58), but the models for lignin, chlorophyll to carotenoids ratio (Car:Cab), cellulose, nitrogen, and C:N ratio were rather unsuccessful (R 2 = 0.38; 0.31; 0.15; 0.23 and 0.23, respectively).

Needle Traits PLSR Modeling-NACs Effect
Results of cross-validation of PLSR models for selected needle traits are presented in Table 2. Models always performed better if all needle age classes were included: they showed higher coefficient of determination (R 2 ), lower root mean square error (RMSE) in absolute values of given trait and lower RMSE expressed as a percentage of the training dataset range (RMSE%), and higher residual prediction deviation (RPD) than models for individual needle age classes. The model for needle specific mass (LMA) showed the highest predictive power (R 2 = 0.79, RMSE% 9% and RPD 2.18). The second-best predicted needle traits were soluble phenolics content (R 2 = 0.73, RMSE% 11% and RPD 1.93), total carotenoid content (R 2 = 0.72, RMSE% 11% and RPD 1.91) and chlorophyll content (R 2 = 0.71, RMSE% 11% and RPD 1.86). The model for water content showed acceptable results (R 2 = 0.60, RMSE% 11% and RPD 1.58), but the models for lignin, chlorophyll to carotenoids ratio (Car:Cab), cellulose, nitrogen, and C:N ratio were rather unsuccessful (R 2 = 0.38; 0.31; 0.15; 0.23 and 0.23, respectively). Figure 7 shows validation of the best performing PLSR models for individual needle age classes (Figure 7 a,c,e,g) and all NACs pooled together (Figure 7 b,d,f,h). LMA was predicted very well using NAC2 and NAC1 reflectance with model parameters close to the model including all NACs. In the case of phenolic compound models, younger needles (NAC1 and NAC2) had better performance in comparison to NAC4. The model for chlorophyll content showed the best parameters (R 2 , RMSE, RMSE% and RPD) for NAC1 in comparison to older needles, NAC2 and NAC4. For water content, the NAC1 model showed the best prediction power and, in some parameters (R 2 and RPD), even outperformed the model based on all NACs together. Carotenoids showed the best performing model for NAC4 (regarding R 2 , RMSE% and RPD) in comparison to NAC1 and NAC2. Table 2. Parameters of partial least square regression (PLSR) models' cross-validation of selected needle traits from needle spectra (450-2000 nm). Parameters for separated models for NAC1, NAC2 and NAC4 are presented, then parameters for the model of all NACs pooled are shown. R 2 -coefficient of determination; RPD-the ratio of performance to deviation; RMSE-root mean square error of validation dataset (in absolute values of given trait); RMSE%-RMSE expressed as a percentage of the training dataset range. The best performing model according to the highest R 2 is highlighted in bold. Several needle traits such as chlorophyll or water content have more specific absorption features than others (e.g., LMA). We used the PLS regression coefficients in the spectral region 450-2000 nm to better explain model performance by highlighting the wavelengths contributing most to the model (Figure 8). A nearly identical pattern in the PLS regression coefficients was found for chlorophyll and carotenoid content (Figure 8a). Significant model contributions were found in the green region (negative coefficients with a peak around 535-545 nm), the red region (positive coefficients between 620 and 700 nm with a peak around 690 nm), the red-edge region (negative coefficients with a peak around 725 nm) and two more regions in SWIR (positive coefficients with a peak around 1485 nm and negative coefficients with a peak around 1713 nm). Significant model contribution of the red region is attributed to chlorophyll absorption maxima [30]. In accordance with other studies [30,72], the other sharp negative peak in regression coefficients for chlorophyll model was found in the red edge region (725 nm). Other significant wavelengths were detected by the PLSR coefficients besides the SWIR peak at 1485 nm. These were located at 945, 1080, and 1930 nm and are related to the water content absorptions [29], showing also the intercorrelation between these two biophysical traits. Surprisingly, the absorption placed at very 1713 nm was found to be very important for modeling the chlorophyll content. This absorption correlates with leaf dry matter and has been used for estimations of lignin content [29,73,74]. For modeling needle water content (Figure 8b), the two most pronounced negative peaks in PLS regression coefficients were identified in red edge (734 nm) and SWIR (1380 nm) wavelengths, which are the water diagnostic absorption features [29]. Additionally, another important spectral region with positive coefficients was placed between 1040 and 1140 nm. Furthermore, contributions from spectral regions around 1280, 1480 and 1900 nm were found and the weaker correlation between water content and the absorption placed at 1713 nm, related to lignin, was also identified. LMA (Figure 8b) shows the highest correlation with the red region with a peak placed at 626 nm, and the red edge region with the maxima placed at 742 nm. Furthermore, we could observe significant contributions of wavelengths around 1100, 1300 and 1490 nm, typically associated with water absorption. The LMA also showed a significant correlation with the 1713 nm lignin absorption feature.

Needle Biophysical Traits-Age Effects
First, we discuss the effect of needle aging on their biophysical traits. As expected, most of the studied traits showed a significant age-dependence: LMA, photosynthetic pigments and phenolic compound contents gradually increased from current to older needles (Figures 2 and 3a). Similarly, the ratio of total carotenoids to chlorophylls increased in older needles (NAC4), corresponding with importance of accumulation of protective carotenoids [75]. However, the differences between NAC1 and NAC2 were not as profound and, for many sites, the Car:Cab ratio in NAC2 was lower than for NAC1. The accumulation of photosynthetic pigments with needle age appears to be a general phenomenon in conifers, as it was shown for Norway spruce (e.g., [76]) and other species, e.g., Pinus pinaster [54]. Accumulation of photosynthetic pigments in older needles may be important for maintenance of high photochemical efficiency in darker microenvironments, as shown in Abies alba [40]. Nitrogen needle content decreased from current to older needles ( Figure  3e) in most studied catchments. This is probably due to N recycling from older to young foliage [77][78][79] together with the dilution effect of accumulating dry matter (increase in LMA, Figure 3a) and also the C:N ratio increase in the majority of the studied catchments For phenolic and lignin content (Figure 8c), we found the same patterns of regression coefficients previously described for the biophysical traits above. When we compare the wavelength position and how pronounced the peaks are, we can see that lignin has the stronger correlation with its diagnostic absorption at 1713 nm and correlated also with the red edge peak placed at 730 nm, while phenolics show stronger correlation with the low reflectance-absorption placed at 540 nm and the red edge peak shifted to longer wavelengths (750 nm). The correlation with the absorption at 1713 nm is weaker when compared to lignin.

Needle Biophysical Traits-Age Effects
First, we discuss the effect of needle aging on their biophysical traits. As expected, most of the studied traits showed a significant age-dependence: LMA, photosynthetic pigments and phenolic compound contents gradually increased from current to older needles (Figures 2 and 3a). Similarly, the ratio of total carotenoids to chlorophylls increased in older needles (NAC4), corresponding with importance of accumulation of protective carotenoids [75]. However, the differences between NAC1 and NAC2 were not as profound and, for many sites, the Car:Cab ratio in NAC2 was lower than for NAC1. The accumulation of photosynthetic pigments with needle age appears to be a general phenomenon in conifers, as it was shown for Norway spruce (e.g., [76]) and other species, e.g., Pinus pinaster [54]. Accumulation of photosynthetic pigments in older needles may be important for maintenance of high photochemical efficiency in darker microenvironments, as shown in Abies alba [40]. Nitrogen needle content decreased from current to older needles ( Figure  3e) in most studied catchments. This is probably due to N recycling from older to young foliage [77][78][79] together with the dilution effect of accumulating dry matter (increase in LMA, Figure 3a) and also the C:N ratio increase in the majority of the studied catchments (Figure 3f). Two catchments where the N and C:N ratio did not follow the same pattern as observed at the other catchments were PLB and UHL. PLB is the only site with Mgrich serpentinite bedrock and UHL was one of two completely deforested catchments in the 1990s due to heavy acid deposition loads. Both PLB and UHL belong to the most N saturated catchments [16], which was manifested by high dissolved and gaseous N losses. Therefore, at PLB and UHL catchments, tree growth is not N-limited and the nitrogen translocation from older to young needles is not necessary.
Higher contents of soluble phenolic compounds in older Norway spruce needles in comparison to the current ones were described before [38,47,76,80]. The enhanced accumulation of phenolics, either soluble phenolics or lignin in needles, is considered a non-specific stress indicator [81]. This phenomenon is probably connected with defense mechanisms against fungal and herbivore attacks [38].
Structural compounds (cellulose and lignin) and water content were independent of needle age in the present study, although some age-dependent trends were described in the literature even for these traits in various conifer species. Lignin, cellulose, and C content decreased with needle age in Picea abies × obovata from northern taiga [82]. The similar cellulose content decrease from current to older needles was observed in Norway spruce growing on a base cation-poor crystalline bedrock, with LYS from the current study being one of the forest stands in focus [83]. Probably due to great variance in cellulose needle content (Figure 3b), there was no general trend related to needle age, which is consistent also with [83]. Needle water content trends varied among the studied sites and the age effect was not statistically significant (Figure 3d). Some catchments, however, showed a trend of lower water content in older needles (LIZ, LYS, NAZ, PLB), which is in accordance with the decrease in water content for older needles observed for Abies alba [40] and Pinus contorta [84].

Needle Traits and Soil Properties-Site Effects
The ten GEOMON catchments significantly differed from their soil properties (Supplementary Table S1), particularly by exchangeable acidity, and base cations concentrations, and less from their needle biophysical traits. Despite the sulphur deposition-the main acidifying factor, which generally decreased at the GEOMON catchments between 1994 and 2014 by ca 70-80% [16]-, the differences among the catchments given by the parent bedrock or forest management became more apparent. On one hand, high base cations content, higher base saturation, and low exchangeable acidity in NAZ and PLB catchments ( Figure 4), were observed in FH horizon, which corresponded to Mg-rich serpentinite bedrock (PLB, [85] and Ca-Mg rich amphibolite bedrock (NAZ [20]). On the other hand, UHL and LYS exhibited high exchangeable acidity, and low base cations due to Ca-poor granite bedrock [86]. The UDL catchment was seriously affected by air-pollution, which caused forest decline in the past, and forest management interventions, particularly the application of dolomitic limestone between 1980 and 2007 [87], were most pronounced there. As a result, UDL catchment, naturally underlined by a Ca poor bedrock (orthogneiss) exhibited relatively low exchangeable acidity and high base saturation in the forest floor. This was further accompanied by high N soil concentration and low soil C:N ratio.
Some needle traits, however, also contributed to among-site variability in PC1. These were chlorophyll and carotenoid contents mainly in young needles (NAC2 and NAC1). Surprisingly, photosynthetic pigment content correlated negatively with Mg and Ca soil contents and positively with total soil potassium, phosphorus, and exchangeable acidity. These relationships almost disappeared for older needles, where P and K resorption to younger tissues may be the cause [41]. Regardless of the needle age, LMA, and nonpigment compounds (lignin, phenolics, water and cellulose), showed less dependence on soil properties, only being influential in PC2. In the PCA models for current-and previous-year needles, soluble phenolic needle content negatively correlated with N needle content, which is consistent with the finding that N fertilization leads to a decrease in various carbon-based secondary metabolites in current Norway spruce needles, such as flavonoids and tannins [38]. There was also an apparent negative correlation between needle N content and organic soil C:N ratio, which diminished in older needles. A similar relationship was observed by Aitkenhead-Peterson et al. [83] for C:N in Norway spruce forest floors. This further highlights the employment of foliar N as a good predictor of site trophic status.
To conclude, the larger part of variation among the ten studied sites can be described by intrinsic soil properties. Needle photosynthetic pigments (chlorophyll and carotenoids) and N content are correlated with soil properties, particularly in NAC1 and NAC2 needles. Needle structural traits (LMA, lignin, cellulose) and water are rather independent of soil properties.

Variation in Needle Optical Properties-Age and Site Effects
Needle optical properties varied according to needle age consistently in all ten studied sites ( Figure 5 and Supplementary Figure S1). Reflectance in the VIS region, particularly in green wavelengths, decreased with needle age, which corresponds well with chlorophyll accumulation. In NIR between 800 and 1000 nm, reflectance decreased with needle age, but in the shortwave infrared wavelengths there was almost no difference among age classes ( Figure 4). Transmittance in the VIS region behaved the same as reflectance, showing a decrease with needle age. In NIR and SWIR, transmittance exhibited significant differences among needle age classes with the highest values for the NAC1. The more pronounced agedependence in transmittance in the NIR and SWIR ranges (in comparison to reflectance) resulted in considerably better needle age class spectra separation by PCA (Figure 6a,c). The age effect on needle optical properties was similar, as shown by Hovi et al. [88] for nine boreal coniferous species. Although Hovi et al. [88] focused only on NAC1 and NAC2, we showed that the optical properties changed similarly with further needle aging. The decreasing transmittance in NIR and SWIR wavelengths could be attributed to cell expansion and accumulation of structural compounds expressed as increasing LMA (Figure 3a), both of which result in the reduction in intercellular spaces in the mesophyll, as shown by Rock et al. [89], in Picea rubens and Tsuga canadensis. Our results, similar to those of [88], clearly show that the age-related differences in needle optical properties are wavelength dependent (Figure 4), which is important to be aware of if dealing with simple ratio or normalized difference vegetation indices. A study by Dengel [90] detected stress by various vegetation indices along an altitudinal gradient in subarctic Pinus sylvestris, and some of the routinely used indices (NDVI, chlNDVI and REP) were equal for the five youngest needle age classes. Needle age-effect on the optical properties at the leaf level has been covered by several studies (e.g., [88]) to date. However, much less attention has been paid to upscaling this phenomenon to the canopy level. A recent study by Wu et al. [91] conducted on Chinese fir showed that needle-age dependent seasonal changes in leaf optical properties resulted in canopy NIR reflectance trends across the season. Thus, the needle age effect should be considered with seasonal effects for better interpretation of satellite remote sensing data.
We discussed considerable variation in needle optical properties given by the needle age, which was consistent across all studied catchments (Supplementary Figure S1). Regarding differences in optical properties among individual GEOMON sites, PCA did not reveal any trend to catchment separation and sample scores were homogeneously scattered (Figure 6b,d). Lately, there have been efforts to use plant spectra as surrogates for important belowground processes [31,32] and canopy spectra appeared to correlate well with soil nitrogen and activity of selected microbial enzymes [32]. Although soil conditions at GEOMON catchments represent substantial gradient in soil properties, we were not able to find any site-related pattern in needle reflectance and transmittance. Recently, significant improvement from acidification was observed in the GEOMON headwater catchments network [16], which could result in the absence of differences in optical properties of studied needle age classes. In extreme cases, like post-mining sites with substrates rich in S and trace metal (As, Cu, Zn), the Pinus sylvestris needle spectral and biophysical properties can reflect the soil properties and separate the sites with different soil conditions [92]. However, in the present study, we did not focus directly on relationships or modeling of soil traits from needle spectra, which will be a subsequent step in exploiting a rich dataset from the GEOMON network.
In general, the measurement of optical properties of narrow, irregular leaves of conifers is challenging as the narrow needles do not fully cover the sample port in integrating spheres and the most common practice is to fix needles in a special carrier with some gaps between them [57,60,[93][94][95]. In particular, needle transmittance measurements are influenced by the gap effect (additional signal transmitted through gaps, signal multiple scattering at the edge of needles and light interaction with a carrier), which often results in erroneous negative transmittance values mainly in the VIS range (as shown, e.g., by Olascoaga et al. [60] for Picea abies and Pinus sylvestris, and as was also the case in this study). Therefore, results based on spectral transmittance should be interpreted with care and transmittance was not used for trait modeling. Nevertheless, the PCA analysis pointed to larger spectral variation within a single tree due to needle age effect than between stands growing in different environmental conditions.

Needle Traits PLSR Modeling-NACs Effect
We selected PLSR as the target regression method for modeling Norway spruce needle biophysical traits for following reasons: as a nonparametric approach, PLSR takes advantage of model training with the full spectrum of information, is designed to deal with spectral data multicollinearity, and includes a dimensionality reduction step [67,96]. Among the multitude of approaches in biophysical parameter retrieval, PLSR has recently been successfully applied in plenty of vegetation studies both on leaf level spectral data, e.g., [92,[97][98][99], canopy-level imaging spectroscopy, e.g., [72], and other studies mentioned in the review by the authors of [96].
The best performing PLSR models based on mixed NACs (Table 2, Figure 7) exhibited slightly lower R 2 than other recent studies: LMA R 2 = 0.79 LMA (present study) compared to 0.90 in [97]; phenolics R 2 = 0.73 (present study) compared to 0.83 in [98]; chlorophyll R 2 = 0.71 (present study) compared to 0.89 and 0.94 in [92,100], respectively; and water R 2 = 0.60 (present study) compared to 0.89 in [97]. In some cases, a lower model performance in the present study probably resulted from the needle-like leaf form of Norway spruce, which may cause difficulties in measuring optical properties at the needle level, as discussed above. The cited studies for model performance comparison [97,98] were performed on broad leaved species-annual crops and poplar cuttings, and Asian broadleaved forest tree species, respectively.
Our results clearly show that needle age has a significant effect on PLSR models' performance (Table 2). For almost all studied traits, the models based on mixed NACs performed the best regarding R 2 and RPD. This result is to be expected for needle traits, which are age-dependent themselves, such as chlorophyll, carotenoid and phenolic content and LMA, due to the extension of the training dataset range. However, the mixed NACs model performed the best also for lignin, a trait without significant age dependence. This finding supports using several needle age classes for functional traits' retrieval from the hyperspectral data.
If we focus on models for individual NACs, the results were not consistent-for chlorophyll a+b and water, the best performing models resulted for NAC1. Considering chlorophyll and water as the most usually monitored traits from the optical signal, it may not bring a serious bias if only current year needles are sampled for model calibration, as in case of study by [101] and many others. On the contrary, LMA, soluble phenolics and lignin showed the best performing models for NAC2 and carotenoids and Car:Cab ratio even for NAC4. The hyperspectral data contain information about a broad range of leaf traits connected to plant physiological performance [102], carbon-nitrogen status and source-sink balance [97]. Therefore, the best way to cover as many traits as possible is to use mixed age classes data. In some chemometric studies dealing with evergreen conifers, the biochemical traits were modeled by PLSR from spectra measured on dried and homogenized needles and the models were always trained on mixed age classes [99,103]. Similarly, Schlerf et al. [104] included NAC1 and NAC3 for chlorophyll and nitrogen retrieval in Norway spruce, both at the leaf and canopy level. However, neither of the studies evaluated the model performance for the NACs separately. One of the scarce studies focused on needle functional traits relation to its optical properties for individual age classes was conducted by Kováč et al. [51]. The green reflectance continuum removal-based optical index was correlated to the xanthophyll de-epoxidation state of Norway spruce foliage. The relationships for various NACs were compared and the linear regression between optical index and the xanthophyll de-epoxidation state showed higher R 2 for NAC2 (0.63) than for NAC1 and NAC2 mixed (0.61). However, needles from different vertical canopy positions were examined, which could influence such dynamic leaf traits, such as the xanthophyll de-epoxidation state. The recent study by Wu et al. [91] strongly supports our idea of always using more needle age classes for taking ground truth for evergreen conifers with a long lifespan foliage. The authors went even further and showed the differences in optical properties seasonal dynamics for new needles (NAC1) and mature needles (NAC2 and older), which had a serious impact on NIR canopy reflectance of evergreen coniferous forests [91].
To sum up, we have analyzed the variability in needle traits induced by needle age only at the top level of the tree crown. However, we are convinced that if we also included needles across the vertical crown profile, we would probably find even more pronounced differences due to needle acclimation to shade conditions [36,44,45]. Understanding of foliage traits' heterogeneity within a tree crown can help to better scale and interpret remotely sensed spectral information by upscaling (for instance, via 3D radiative transfer models).

Conclusions
This study provided deeper insight into the variability of needle biophysical traits and optical properties of the most common productive tree species in Central Europe-Norway spruce. Our database of needle traits was compiled from 123 spruce trees growing in different environmental conditions and was further complemented by information about sites' soil properties. We confirmed the age-dependence of most studied needle biophysical traits: pigments (chlorophyll a+b and carotenoid contents), leaf mass per area, water and phenolics content. We found out that variation in needle optical traits was mainly driven by the presence of different needle age classes within a single tree crown rather than different catchments environmental conditions. Needle reflectance in the VIS and NIR regions and transmittance in the VIS, NIR and SWIR regions decreased with increasing needle age. This also influenced the accuracy of PLSR prediction models of needle traits, which varied among the needle age classes. Traits that could be modeled with moderate accuracy were needle pigments (chlorophyll a+b and carotenoid contents), leaf mass per area, water and phenolics content. Traits with low prediction accuracies were lignin, cellulose, and nitrogen content. PLSR models always performed better if all needle age classes were included in comparison to single age-class based models, which applied also to needle-age independent traits (water and lignin). Variability among ten catchments could be largely attributed to different soil properties (mainly to exchangeable acidity, total phosphorus content, pH, base saturation, and content of base cations that formed the PC1 axis, explaining about 27% of the total variability) rather than to differences in needle biophysical traits (which formed the PC2 axis, explaining about 18% of the total variability). Due to the significant age-dependence of many needle traits, and also the effect of needle-age on PLSR models' performance, we recommend including not only the current but also older needle traits as ground truth for Norway spruce, and probably for other conifers bearing several needle age classes.

Data Availability Statement:
The data presented in this study are available on reasonable request from the corresponding author.