1. Introduction
Silver birch is one of the pioneer tree species found in the boreal forests. Within the global forest area, the boreal forest comprises 30%. The economic and ecological importance of silver birch paved way to research in their genotypic variation. The recent advances in genome sequencing of silver birch contribute to a deep understanding about the evolution of this tree species and thereby characterizing the genomic adaptation [
1]. The wide range of intraspecific variation or genetic diversity is an important component that enables tree species to adapt to changing environmental and biotic conditions. Leaves are primary organs for the photosynthesis, transpiration, and acclimation of plants to variable and changing environmental conditions through various physiological, chemical, and structural modifications [
2,
3,
4]. Silver birch genotypes have been shown to vary widely based on their leaf morphological and physiological traits [
5,
6]. Leaf spectral reflectance is an important leaf trait that varies based on the leaf surface and internal properties and their interaction with the environment.
The leaf spectral reflectance is affected by changes in the leaf structural, physiological, and biochemical traits. The reflectance in the visible range is mainly determined by the absorption properties of leaf pigments such as chlorophyll, carotenoids, and anthocyanins [
7,
8]. Reflectance on the leaf surface is also influenced by the leaf morphological features such as the trichomes and the density of epicuticular waxes [
9] that may allow less light to penetrate through the leaves in the visible and near-infrared range [
9,
10]. The leaf internal structures such as the cell shape and aerial interspaces play a role in the spectral range of 700–1300 nm, where the absorption by water is relatively weak [
9,
11]). Reflectance in the shortwave infrared region is strongly influenced by the water content, with the main absorbance bands at 1450 nm, 1940 nm, and 2500 nm [
11]. In addition to water, the shortwave infrared reflectance is influenced by many biochemical compounds such as cellulose, lignin, protein, and nitrogen. These chemical constituents correlate with the structural features such as specific leaf area (SLA).
Optical remote sensing technology provides a possibility for a non-invasive and fast method of capturing a large amount of information with high spatial resolution compared with the traditional methods. Species differentiation by hyperspectral remote sensing, analyzing the reflectance spectra has shown great promise in diverse ecosystems [
12]. Variation in the leaf spectral reflectance can be used to investigate biodiversity across different plant ecosystems [
3]. Analysis of reflectance spectrum has been utilized in mapping the leaf biochemistry [
13,
14,
15] and the native and non-native forest tree species in rainforests [
13]. In Finnish forests, leaf reflectance properties in the visible and near-infrared region as well as baseline variability of photochemical reflectance index (PRI) have been shown to differ greatly among three tree species: silver birch (
Betula pendula Roth.), Scots pine (
Pinus sylvestris L.) and Norway spruce (
Picea abies (L.) Karst) [
15,
16]. Hyperspectral imaging has been utilized in detecting the intraspecific variance of boreal tree species for seasonal variation [
17], within-canopy variation [
16,
17,
18], and the leaf sides, i.e., the adaxial and abaxial side of the leaves [
15,
17]. However, the applicability of hyperspectral imaging in studying genotypic variation of boreal tree species, especially in silver birch, has hardly been assessed. A study by Madritch et al. (2014) reported intraspecific genetic diversity in trembling aspen (
Populus tremuloides) forests at canopy level reflectance measured with satellite-based remote sensing. The aspen genotypes can be distinguished from each other based on their spectral features, and this variation is further reflected on soil traits and below-ground processes. However, to our knowledge, other tree species have not been so far studied in respect to the genotypic differences in their reflectance spectra.
This study focuses on silver birch that has a wide range of intraspecific variation and thus has a large natural range covering the entirety of Eurasia [
19,
20]. Intraspecific variation in silver birch tree morphology and physiology, herbivory, and leaf litter decomposition is well-documented [
5,
21,
22,
23]. This variation in leaf secondary metabolites is clear in silver birch and affects the herbivory and abiotic interactions [
24,
25,
26,
27]. Leaf surface secondary metabolites differed clearly among the genotypes and the chemical profile differed among the provenances of silver birch in our previous study [
28]. Differences among the genotypes within a population in chlorophyll content, leaf area (LA), and SLA were reported in Possen et al. (2014), whereas differences among the provenances were unclear in chlorophyll content, LA, and SLA [
6]. This work provides a detailed study of the variance in reflectance properties among provenances of origin (populations) and the genotypes (individuals) for a deciduous tree species, silver birch (
Betula pendula Roth), in a common garden for two consecutive years (2015 and 2016). In this study, the leaf spectral reflectance of silver birch in the wavelength range of visible and near-infrared (VNIR, 400–1000 nm) and shortwave infrared (SWIR, 1000–2500 nm) were measured. The objective of the present study is (1) to examine the genetic variation in the leaf spectral reflectance among and within three Finnish silver birch provenances and (2) to identify the key spectral wavelengths that influence the discrimination of the leaves based on the provenance and genotype. It is also hypothesized that the leaf traits measured in this study could be associated with the reflectance features.
2. Materials and Methods
2.1. Sampling
This study was conducted in a botanical garden in Joensuu, Finland (62°35′ N, 29°46′ E) (see Heimonen et al. (2015) for common garden experimental details) on 10 June 2015 (161 doy) and 21 of July 2016 (203 doy) and the sampling was done between 12:00 and 14:00 h Finnish time. Micropropagated plantlets of 26 genotypes from six provenances were planted in five blocks in botanical garden on July 2010. Silver birch trees originating from three provenances of different latitudes: 60° N (southern Finland, Loppi), 62° N (central Finland, Vehmersalmi), and 66° N (northern Finland, Rovaniemi) were selected for this study. Four genotypes from each latitude were sampled: L1, L6, L14, L15 (60° N, Loppi), V1, V4, V5, V14 (62° N, Vehmersalmi), and R3, R8, R11, and R15 (66° N, Rovaniemi). Details of the common garden experimental design and setup are described in Heimonen et al. (2015). Samples were collected from all of the five experimental blocks, one leaf from each genotype in each block (n = 9–11 for genotype, n = 39–41 for provenance) for 2015 and from three experimental blocks, three leaves from each genotype in each block (n = 3 × 3 for genotype, n = 12 × 3 for provenance) for 2016.
2.2. Data Acquisition and Ecophysiological Measurements
Mature leaves from the southern side of the tree from the upper third of the canopy were chosen for the measurements. Before the sample collection, the chlorophyll content indices were measured in the field with a chlorophyll meter, Dualex 4 Scientific (Dx4) (FORCE-A, Paris, France). The chlorophyll content was measured as the proportion of transmittance at 850 to 710 nm [
29]. Then, the measured leaves were detached from the tree, collected in separate plastic bags, kept on ice in a portable cooler and moved to the laboratory. The fresh weight (FW, g) was measured before the reflectance measurements. The reflectance measurements were performed within three hours of the leaf collection in the laboratory. The leaves were oven dried overnight (16 h) for 40 °C after the reflectance measurements and the dry weight (DW, g) was measured the next day. Subsequently, the water content (WC [%]) was calculated as follows, WC = [(FW − DW)/FW] × 100). The leaf area was calculated from the images (0.4 mm resolution) acquired during reflectance measurements in the visible near-infrared range with Evince v.2.7 hyperspectral imaging software (Prediktera, Umeå, Sweden). Specific leaf area (SLA) was calculated by dividing leaf area (cm
2) by leaf dry weight.
2.3. Reflectance Spectral Measurement Set-Up
The reflectance spectra of leaves were measured using a pushbroom hyperspectral imaging system by Specim (Spectral Imaging Ltd., Oulu, Finland). The system includes a visible and near-infrared (VNIR) camera combined with an imaging spectrograph (ImSpector V10E) ranging from 400–1000 nm and a short wavelength infrared (SWIR) camera with a spectrograph (ImSpector N25E) that covers the spectral range of 1000–2500 nm. The illumination to the samples was provided uniformly by eight tungsten halogen lamps of 35 W in an orientation of 45 degrees with the leaf surface. Two cameras encased side by side (back and front) within the frame with the light panels connected to the left and right were attached to a conveyor belt, moving the cameras and light panels. The leaves were placed at 48.5 cm and 47 cm from the objectives of the VNIR and SWIR on a fixed table under the camera. The image focus was adjusted for both cameras. The exposure time was set at 8.1 ms and 1.2 ms for the acquisition of the leaf images by the VNIR and SWIR cameras, respectively. The VNIR camera has a spatial resolution of 1032 pixels and 240 spectral channels (FWHM = 3.5 nm) and the resulting image had a spatial resolution of 0.4 mm per pixel. The SWIR camera has a spatial resolution of 320 pixels and 256 spectral channels (FWHM = 12 nm). The resulting image had a spatial resolution of 1.2 mm per pixel.
Leaves were arranged in a four-times-six grid (24 leaves) on a plywood plate (40 × 60 cm) with the adaxial side facing the camera. The plywood was painted homogeneously with Maston 100 series Matt black to reduce reflectance from the wood (as in Deepak et al., 2019). The leaves were arranged in random order in five batches for 2015. For 2016, the leaves were imaged in a random order in five batches with a similar arrangement for all batches, except that the last batch included only 12 leaves for 2016. The leaves were covered by a high transmission glass plate (Pilkington Optiwhite) with high transparency from 400 nm to 2500 nm (as in Deepak et al., 2019) to press the leaves flat against the surface of the plywood. Each batch in 2016 was imaged three times in a sequence and three times after switching the places of the leaves in random order to eliminate possible position-related effects. The average of all six scans was used for the reflectance of each leaf.
2.4. Spectral Correction and Preprocessing
The SpectralonTM (Labsphere, North Sutton, NH, USA) standard white reference image (close to 99% of reflectance from 400 nm to 2500 nm) was captured with the same imaging configuration as for the leaves to compensate the non-uniformity of the illumination. The dark noise image was acquired by closing the mechanical shutter of the cameras to remove non-uniformity distribution due to the system. The instrument measurement values of the images were converted to the reflectance value with Evince v.2.7 hyperspectral imaging software (Prediktera, Umeå, Sweden). Spikes caused by sensor faults were replaced by the median of the intensity value of neighboring pixels.
The reflectance of leaves samples was calculated from the raw images as follows:
where
is the spectral reflectance image,
is the raw spectral image,
is the dark image,
is the white-reference-acquired image, and
is the reflectance spectra of the white-reference sample in the range from 400 nm to 2500 nm. The leaves were segmented from the background with thresholding. The mean reflectance was extracted for each leaf by averaging the pixels values of each wavelength, excluding the leaf margins.
2.5. Statistical Analysis
SIMCA-P+ 14.1 (Umetrics, Umeå, Sweden), IBM SPSS Statistics version 21 (IBM Corp, Armonk, NY, USA) and R software version 3.5.0 (R Core Team, 2018, Vienna, Austria) were used for the statistical analysis. Spectral ranges from 386 to 480 nm, 979 to 1009 nm, and 2313 to 2538 nm were removed before preprocessing due to excess noise. Before multivariate analysis, first derivative transformation with five points smoothening (Savitzky–Golay filter, quadratic polynomial order) and mean centering were applied to the spectra. Principal component analysis (PCA), an unsupervised method, was used for visualization of possible genotype- and provenance-related patterns and for the dimensionality reduction of the reflectance dataset for both years separately. Since the hyperspectral dataset contains highly correlated neighboring bands, PCA can be effectively used [
30]. In order to differentiate between the groups within the reflectance dataset, the discriminant analysis of principal components (DAPC) were employed [
31] as implemented in the adegenet package [
32] of the R program (R Development Core Team, 2014). Model validation was done using cross-validation by retaining the optimal number of PCs. The number of PCs associated with the lowest RMSE values was chosen as the optimal number of PC’s, and thus four PCs were retained. The contribution plot estimated the contribution of each wavelength to the DAPC model. Apart from DAPC, partial least squares discriminant analysis (PLS-DA) was also used to differentiate the groups, and analysis of variance testing of cross-validated predictive residuals of the model (CV-ANOVA) was used for cross-validation of the model. Variable importance in projection (VIP) estimated the importance of each variable in the PLS-DA model. Variation in the leaf traits and ten individual wavelengths due to provenance and the genotype within each provenance were analyzed for both years separately with linear mixed models (IBM SPSS Statistics version 21, R Core Team 2014) with the provenance as fixed factor and genotype within the provenance as a random factor. In addition, marginal R2 (R2m, variance explained only by fixed factor) and conditional R2 (R2c, variance explained by both fixed and random factors) [
33] for the leaf traits were calculated using the r.squaredGLMM function in the MuMIn package [
34]. Bar charts for the leaf traits and spectral indices were prepared using the R package ggpubr [
35]. Relationships among the chlorophyll and their reflectance indices were studied by means of Pearson correlation analysis.
Two normalized difference vegetation indices (NDVIs), NDVI
670,780: (R
780 − R
670)/(R
780 + R
670) [
36] and NDVI
705,750: (R
750 − R
705)/(R
750 + R
705) [
37] were used in this study. Two NDVIs were used to compare a red edge-based index to chlorophyll absorbance peak-based index. Our previous study on silver birch canopy layer variation [
18] found variation in the performance of these two indices. Apart from this, chlorophyll reflectance index CRI
550,672,708: R
672/(R
550 × R
708) [
38] was also included in this study.
3. Results
The average leaf reflectance spectra separated the central provenance in the green hump and showed clear difference among the three provenances in the NIR plateau, continuing up to 1800 nm in the SWIR range (
Figure 1a). In the green hump, leaves from the central provenance had higher reflectance than leaves from the southern or northern provenance, whereas in the NIR plateau, the leaf reflectance among the provenances followed a latitudinal order, increasing from south to north. In the first derivative reflectance spectra, differences among the provenances were found at the green hump, red edge, and at the two minima in the SWIR range (
Figure 1b). The central provenance differed from the others at the red edge, with the red edge inflection point (REIP) at 702 nm for 2015 and 712 nm for 2016, whereas for the southern and the northern provenances it was at 704 nm for 2015 and 714 nm for 2016 (
Figure 1b). At two minima of the first derivative spectra (1389 nm and 1874 nm), the provenances followed a latitudinal order with the highest values for the southern provenance.
The magnitude of variation in the reflectance of 2015 and 2016 is shown in
Figure 1. The green hump and the NIR plateau had high reflectance in 2015 compared with 2016, whereas the red edge showed lower reflection in 2015 (
Figure 1a,b). In the first derivative spectra, the difference among the provenances was clearer in 2016 than in 2015 (
Figure 1c,d).
Furthermore, the mean leaf reflectance spectra of the green hump, red edge, and NIR around 750 nm for individual genotypes for both years are shown in
Supplementary Figure S3. The differences in the first derivative reflectance spectra, among the genotypes at the green hump, red edge, and at the two minima in the SWIR range are shown in
Supplementary Figure S4. Provenance showed significant difference in all of the wavelengths analyzed with the linear mixed model except at 1874 nm in 2016. However, genotype was significant only at 704 nm in 2015 and 712 nm in 2016. Both years showed significant difference among the genotypes at 750 nm.
In the PCA analysis for 2015, the provenances formed overlapping clusters on the first and the second PCs, explaining 59% and 16% of the variance, respectively (
Figure 2a). For 2016, there was a variance of 61% and 14% for the first and the second axis (
Figure 2b). The data was further analyzed with PLS-DA (
Figure 2c,d), explaining 59% and 12% of variance for 2015 and 61% and 13% variance for 2016. The PLS-DA model was cross-validated with CV-ANOVA, which showed a significant
p value (
p < 0.0001) for the provenance in both years.
The most important wavelengths that influenced the discrimination between the provenances were identified as 517 nm, 519 nm, 694 nm, 696 nm, 722 nm, 725 nm, 727 nm, 1383 nm, 1389 nm, 1862 nm, 1868 nm, and 1887 nm (
Figure 3) with PLS-DA. Wavebands influential for the separation of the southern and northern provenances were on the first discriminant axis at 504–527 nm, 683–750 nm, 1364–1408 nm, and 1862–1874 nm in 2015 (
Figure 3a) and at 507–579 nm, 686–753 nm, 1370–1389 nm, and 1855–1868 nm in 2016 (
Figure 3b). The wavebands that influenced the separation of the central provenance were on the second axis at 507–532 nm, 683–748 nm, 1370–1414 nm, and 1862–1906 nm in 2015 (
Figure 3c) and 509–574 nm, 686–750 nm, 1351–1421 nm, and 1849–1893 nm in 2016 (
Figure 3d). The most influential waveband in both the first axis and the second axis was 694 nm in 2015 and 696 nm in 2016 nm.
DAPC and PLS-DA were used to explore genotypic differences in reflectance spectra. After cross-validation, four PCs were retained and the DAPC space formed by the first two axes showed clustering of the genotypes within their respective provenances for both years (
Figure 4). The analysis with PLS-DA (
Supplementary Figure S1) followed by cross-validation with CV-ANOVA (
p > 0.05) indicated a non-significant model, even though the patterning indicated some grouping. The PLS-DA was further applied for the discrimination of genotype within each provenance separately (
Supplementary Figure S2). The genotypes from Loppi showed a clear separation of L1 and L15 from L14 and L16 with CV-ANOVA (
p < 0.01) and the genotypes from Vehmersalmi showed a clear separation of V1 from the other genotypes with
p < 0.01. The Rovaniemi genotypes did not show significance (CV-ANOVA,
p = 0.088 for 2015 and
p = 0.771 for 2016) for the PLS-DA model.
The discrimination of genotypes in the DAPC model was most influenced by wavelengths at 696 nm, 720 nm, 1389 nm, 1887 nm, and 1893 nm (
Figure 5). Wavelengths at 694–727 nm, 1395–1408 nm, and 1868–1906 nm were influential in the separation of the southern and northern genotypes on the first discriminant axis in 2015 data, whereas in 2016, it was 694–725 nm, 1363–1414 nm, and 1862–1893 nm (
Figure 5a,b). The wavelengths at 519–532 nm, 686–709 nm, and 1874 nm-1906 influenced the separation of the central genotypes on the second axis in 2015 data, and in 2016, it was 689–707 nm, 857–860 nm, 876–886 nm, and 1874–1906 nm (
Figure 5c,d).
Provenance affected the chlorophyll content, NDVI 705, 750, and CRI (
Figure 6 and
Figure 7,
Table 1) in both 2015 and 2016. The contents of all the three leaf traits were the highest in the northern provenance and lowest in the central provenance in both years. Variance explained by provenance was relatively low for 2015 in comparison with 2016. Variance explained by provenance was more prominent compared with the total variance explained by both genotype and provenance (
Table 1). Provenance explained most of the variance for the model with the R
2c (variance explained by both provenance and genotype within the provenance) explaining only 0 to 19% of variance.
The CI correlated with the NDVI705,750 (r = 0.84, p < 0.001) and CRI550,672,708 (r = 0.79, p < 0.001), whereas NDVI670,780 (r = 0.19, p = 0.28) showed no correlation in 2015. The CI correlated with the NDVI705,750 (r = 0.450, p < 0.001) and CRI550,672,708 (r = 0.368, p < 0.001) also in 2016.
To explore the association of the measured ecophysiological traits and reflectance indices to the DAPC result, the genotype averages of the traits on the DAPC space for 2015 and 2016 were plotted. The only traits to show a clear pattern related to the discrimination of the genotypes were related to the chlorophyll content of the leaves (CI, NDVI
705,750, CRI
550,672,708) (
Figure 8, shown for CI). The differences among the provenances and genotypes in leaf spectral reflectance thus followed the differences in chlorophyll content in both study years.
4. Discussion
In this study, differences in the spectral reflectance among the silver birch provenances and genotypes in two study years, 2015 and 2016, were demonstrated. Among provenances, the average reflectance spectra and the first derivative spectra revealed differences around 550 nm, at the red edge, in the NIR plateau, and at two first derivative spectrum minima in the SWIR. The leaf reflectance showed a shift in the red edge inflection point among the provenances. The shift in REIP had more prominent yearly variation than provenance-related variation. Year-to-year variation in the REIP shift is probably due to seasonal variation, as analysis was performed on 10 June for 2015 and 21 July for 2016.
The provenances were not well separated from each other in the principal component analysis, describing the overall variation present in the reflectance spectra. Further analysis with PLS-DA obtained a significant separation of provenances for the both years, 2015 and 2016. Analysis with PLS-DA showed a similar pattern for the separation genotypes for both years, but the model was not significant. However, further analysis with PLS-DA separately for provenances showed significant separation for some of the genotypes indicating within the provenance genotypic variation for both years. A clear separation of L1 and L15 from L6 and L14 in the Loppi provenance were obtained for both years. In the Vehmersalmi provenance, V1 genotype formed a strong group separate from other the genotypes during both years. In a discriminant analysis of principal components performed for separation of genotypes, the provenances were distinguished from each other more clearly, following a similar pattern in both years. Individual genotypes were mostly overlapping with each other, even though some groupings were consistent between the two study years. The provenance-related variation in the leaf reflectance appeared thus much stronger than genotypic variation. This runs contrary to our previous study demonstrating a higher intraspecific (or genotypic) variation of leaf surface secondary metabolites as compared to provenance-related variation in the same genetic materials [
28]. In silver birch, clear genotypic variation has been found in secondary chemistry and herbivore resistance [
24,
26,
39,
40,
41], and a high degree of phenotypic plasticity in the growth among genotypes was also found [
41].
Among the ecophysiological traits and spectral reflectance indices, the only traits showing significant provenance variation were related to the chlorophyll content of the leaves, dualex chlorophyll index (CI), chlorophyll reflectance index (CRI), and the red edge-based NDVI. This is in line with the PLS-DA VIP result that showed the strong influence of wavelengths at the red edge. By contrast, Tenkanen et al. (2019) found clear differences among Finnish silver birch provenances in height growth increment, stomatal conductance, and water use efficiency (WUE), and less strongly for photosynthetic rate, but no differences in chlorophyll content index. These datas were gathered from the same common garden site; however, Tenkanen et al. (2019) studied younger trees (in 2013). The higher chlorophyll content (CI), as well as higher values for indices CRI and red edge-based NDVI in 2016 than in 2015, were probably due to the seasonal effect, since the 2016 sampling was later in the season. It is known that chlorophyll content of silver birch leaves increases in early June (Tenkanen et al., 2019).
The PLS-DA VIP plot indicated that also wavelength bands at 1389 nm, 1395 nm, 1868 nm, and 1887 nm in the SWIR were influential in distinguishing the provenances for 2015 and 1383 nm, 1389 nm, 1862 nm, and 1868 nm for 2016. The wavelengths do not correspond exactly with any of the absorption bands reported for plant chemical constituents in the literature [
42,
43], but two wavebands in the same region (1377–1395 nm and 1887 nm) were found to be influential in separating silver birch canopy layers in our previous study [
18]. These wavebands were quite close to the main water absorbance bands in the SWIR, and since the water content affects reflectance widely in the SWIR region, this might have influenced the results. However, there was no significant difference in the water content of the leaves among the provenances. There was no yearly variation in the influential wavelengths for the separation of the provenances.
In accordance with the observed trait differences among the provenances, plotting the genotypic averages of the measured traits and spectral indices in the DAPC result space exhibited clear patterns of increasing magnitude only for traits associated with chlorophyll content. Remarkably, two common variants of the NDVI, based either on 670 and 780 nm or 705 and 750 nm, showed very different results. The NDVI
670,780 did not correlate significantly with chlorophyll index estimated with Dualex; it also did not differ among provenances and failed to show any tendency for a pattern in the DAPC space. By contrast, the red edge-based NDVI
705,750 not only differed significantly among the provenances, but also showed a strong correlation with the Dualex chlorophyll estimate which exhibited a genotype and provenance-related pattern in the DAPC space of the leaf spectral reflectance. The performance of spectral indices used for chlorophyll estimation, or strongly associated with it, have been found to vary widely among different leaf types deciduous trees [
44]. Chlorophyll reflectance indices related to red edge position have been successfully used to estimate chlorophyll content in various studies [
45,
46]. In our study on silver birch canopy layer reflectance [
18], the red edge-based indices also performed better than indices based on chlorophyll peak absorbance at about 680 nm. Correlations of various reflectance indices with chlorophyll content or other plant traits have been widely studied. The chlorophyll content is strongly influenced by the plant species, environmental conditions, and possible experimental treatments, and the results of this study suggest that the differences in the chlorophyll content lead to differences in the leaf reflectance spectrum that are consistent in the data of two consequent years. The canopy variation was clearly reported differently for variants of NDVI indices for same species (silver birch) in Deepak et al. (2019). To our knowledge, this is the first study to associate foliar chlorophyll content and the chlorophyll-related spectral indices for the discrimination of provenances and genotypes of a naturally occurring tree species.