**Impacts of Climate Change on Tree Physiology and Responses of Forest Ecosystems**

Editor

**Mariangela Fotelli**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Mariangela Fotelli Hellenic Agricultural Organization Dimitra Forest Research Institute Vassilika, Thessaloniki Greece

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Forests* (ISSN 1999-4907) (available at: https://www.mdpi.com/journal/forests/special issues/Tree Physiology).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-2750-5 (Hbk) ISBN 978-3-0365-2751-2 (PDF)**

Cover image courtesy of Mariangela Fotelli

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**

## **About the Editor**

**Mariangela Fotelli** graduated with distinction from the Aristotle University of Thessaloniki, Greece, and received her PhD with honors in 2002 from the Albert Ludwig University of Freiburg, Germany, working on nitrogen and water balance in European beech natural regeneration. In the following years, she participated in several international research projects and carried out post-doctoral studies on the carbon metabolism of N-fixing legumes at the Agricultural University of Athens and on the nitrogen ecophysiology, water and carbon status of forest species under climatic variability at the Democritus University of Thrace and the Forest Research Institute of Thessaloniki, Greece. She has been a researcher and head of the lab of Forest Ecophysiology at the Forest Research Institute of the Hellenic Agricultural Organization Dimitra since 2017.

### *Editorial* **Impacts of Climate Change on Tree Physiology and Responses of Forest Ecosystems**

**Mariangela N. Fotelli**

Hellenic Agricultural Organization Dimitra, Forest Research Institute, Vassilika, 57006 Thessaloniki, Greece<EMAIL_ADDRESS>
#### **1. Tree Adaptive Responses to Climatic Variability**

In a changing climate, forest trees have to deal with a range of altered environmental conditions. Within this Special Issue, researchers aimed at revealing the responses of different forest species to various globally changing abiotic parameters, such as increasing atmospheric CO2, elevated temperatures, and limited water availability.

In an experiment with open-top chambers, Kim et al. [1] studied the physiological and transcriptome responses of Poplar hybrids to elevated CO2, being one of the main components of global change. They identified differentially expressed genes and alterations in chlorophylls, starch, and soluble sugars in response to increasing CO2. The impact of elevated CO2 was also assessed by Wang et al. [2], who reported photosynthetic declines in *Phoebe bournei* Hemsl., an endemic tree threatened by habitat loss in China. However, they concluded that this negative impact of increased CO2 may be alleviated by the application of ammonium or nitrate fertilization. The effect of drier growth conditions on the tree radial profile of sap flow was assessed in a coastal Mediterranean Aleppo pine ecosystem [3]. A steeper decline in sap flux density with increasing sapwood depth was observed in drier periods. The limiting role of water availability in such xerothermic environments was further supported by the fact that sap flow responded positively to a favorable microclimate only when water availability exceeded a certain threshold. The effect of rising temperature was studied by Sambayo-Maldonado et al. [4] who assessed the thermal time and cardinal temperatures of *Cedrela odorata* in order to estimate changes in the germination of red cedar under different climate change scenarios. It was found that as the temperature increases in line with the predictions for 2050 and 2100, the germination of the soil seed bank will accelerate, which may have implications for the future distribution of this important forest species in tropical environments.

Common gardens and provenance trials are an experimental approach offering longterm information regarding the potential of different forest species and origins to adapt to future climate. In order to assess the suitability of different Mediterranean oak species for potential assisted migration schemes in central Europe, Bantis et al. [5], established a four common garden experiment and evaluated the early field performance of oak seedlings of different species and provenances. Based on survival rates, morphological, and physiological features, they concluded that *Quercus pubescens* may be a suitable candidate forest tree for poor and dry soils under the expected summer heat and drought conditions in Central Europe. In the same context, a provenance trials approach was used to evaluate nine European beech (*Fagus sylvatica* L.) provenances under contrasting climates [6]. An acclimation potential of provenances transferred to warmer and drier conditions was indicated by the development of smaller leaves with decreasing stomatal density and length of guard cells. However, the higher phenotypic plasticity in stomatal traits was also associated with increasing mortality under more arid conditions.

**Citation:** Fotelli, M.N. Impacts of Climate Change on Tree Physiology and Responses of Forest Ecosystems. *Forests* **2021**, *12*, 1728. https:// doi.org/10.3390/f12121728

Received: 5 December 2021 Accepted: 6 December 2021 Published: 8 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **2. Responses to Changing Climate at the Ecosystem Level**

Climate change is also expected to induce alternations in the distribution range of forest species, species' replacement and losses of forest land, with important economic and ecological consequences.

In this context, Zindros et al. [7] monitored the changes in tree line over a 40-year period in the core of Olympus Mt. National Park in Greece and reported a warming-related upward shift in the tree line. Apart from spatial shifts, temporal shifts in forest ecosystems ecophysiological and phenological responses are also related to climate change. Marquardt et al. [8], compared the sensitivity and temporal stability of climate–growth relationships in *Pinus arizonica*, and *P. ponderosa* var. *brachyptera* consisting of two needle-types, in Arizona, USA, in order to characterize their growth responses to changing climate. They concluded that the species responded to summer precipitation pre-1950 but are now responding to spring precipitation post-1950, representing a shift in allocation of tree resources from maximizing biomass in the summer to reproduction spring. In the same lines, Ogaya and Peñuelas [9] found that a small but continuous decrease in soil water availability for 21 years in a Mediterranean holm oak (*Quercus ilex* L.) forest resulted in the gradual replacement of holm oak by another evergreen shrub (*Phyllirea latifolia* L.), which is more tolerant to hot and dry conditions.

#### **3. Conclusions**

The above-mentioned and other forest trees' and ecosystem responses to rapidly changing climate are reviewed by Kijowska-Oberc et al. [10]. In their review, they brought together literature reports on a range of adaptive responses to both forest abiotic threats (such as drought and heat) and biotic threats (such as pathogens and insects) and they discuss the role of epigenetics, phenology, and phenotypic plasticity for tree adaptations to climate change.

The studies included in this Special Issue illustrate the clear effects of climate change on forest trees' function, evidenced by different responses on genes expression, biochemical, and physiological levels. However, it is also indicated that different forest species from a broad range of climate zones and ecosystem types are able to develop adaptive mechanisms to cope with the current climatic variability. On the ecosystem level this may lead to distribution shifts and changes in vegetation composition. Further research is needed to identify if measures, such as assisted migration or others, may facilitate the transition to climate change resilient forest ecosystems.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

### *Article* **Physiological and Transcriptome Responses to Elevated CO2 Concentration in** *Populus*

**Tae-Lim Kim 1, Hoyong Chung 2, Karpagam Veerappan 2, Wi Young Lee 1, Danbe Park 1 and Hyemin Lim 1,\***

**Abstract:** Global climate change is heavily affected by an increase in CO2. As one of several efforts to cope with this, research on poplar, a representative, fast growing, and model organism in plants, is actively underway. The effects of elevated atmospheric CO2 on the metabolism, growth, and transcriptome of poplar were investigated to predict productivity in an environment where CO2 concentrations are increasing. Poplar trees were grown at ambient (400 ppm) or elevated CO2 concentrations (1.4× ambient, 560 ppm, and 1.8× ambient, 720 ppm) for 16 weeks in open-top chambers (OTCs). We analyzed the differences in the transcriptomes of *Populus alba* × *Populus glandulosa* clone "Clivus" and *Populus euramericana* clone "I-476" using high-throughput sequencing techniques and elucidated the functions of the differentially expressed genes (DEGs) using various functional annotation methods. About 272,355 contigs and 207,063 unigenes were obtained from transcriptome assembly with the Trinity assembly package. Common DEGs were identified which were consistently regulated in both the elevated CO2 concentrations. In Clivus 29, common DEGs were found, and most of these correspond to cell wall proteins, especially hydroxyproline-rich glycoproteins (HRGP), or related to fatty acid metabolism. Concomitantly, in I-476, 25 were identified, and they were related to heat shock protein (HSP) chaperone family, photosynthesis, nitrogen metabolism, and carbon metabolism. In addition, carbohydrate contents, including starch and total soluble sugar, were significantly increased in response to elevated CO2. These data should be useful for future gene discovery, molecular studies, and tree improvement strategies for the upcoming increased-CO2 environments.

**Keywords:** carbon dioxide; climate change; open-top chamber; RNA sequencing; gene expression analysis; *Populus*

#### **1. Introduction**

The CO2 concentration in the atmosphere is rapidly increasing due to industrialization and deforestation, causing global warming and abnormal weather conditions. The reduction in the CO2 concentration in the atmosphere through photosynthesis highlights the importance of forests (IPCC, 2013). According to the IPCC Fourth Assessment Report (IPCC 2007), the atmospheric CO2 concentration will increase from 389 ppm in 2005 to 550–700 ppm by 2050 and to 650–1200 ppm by 2100, which will cause global warming of 1.1–6.4 ◦C by the end of this century [1]. The concentration of CO2 in the atmosphere directly determines the ratio of plant photosynthesis, and indirectly affects plant productivity and fitness, and thus could serve as a form of selective pressure, though there is little evidence to support this contention. Carbon dioxide is a major component in photosynthesis that converts solar energy into energy stored in carbohydrates, thereby controlling plant yield. An elevated CO2 concentration raises the photosynthetic rate and plant's net production or biomass [2], but a rising temperature increases the ratio of photorespiratory loss of carbon to photosynthetic gain, thereby having an opposite effect [3]. The increase in

Veerappan, K.; Lee, W.Y.; Park, D.; Lim, H. Physiological and Transcriptome Responses to Elevated CO2 Concentration in *Populus*. *Forests* **2021**, *12*, 980. https://doi.org/ 10.3390/f12080980

**Citation:** Kim, T.-L.; Chung, H.;

Academic Editors: Mariangela Fotelli and Cate Macinnis-Ng

Received: 30 March 2021 Accepted: 21 July 2021 Published: 23 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the concentration of atmospheric CO2 is likely to have a significant effect on the photosynthesis, metabolism, and development of plants [4,5]. In plants, increased primary carbon sources accelerate their metabolism and growth, especially when other growth sources are abundant [6,7].

*Populus* is a significant model for studying the effects of abiotic stresses on trees. In addition, *Populus* is now recognized as a model tree [8], enabling genomic resources to be deployed to answer questions of ecological and evolutionary significance on plant response and adaptation to climate change. The *Populus* species is an excellent model for examining drought stress responses, which affect not only survival but also biomass accumulation [9,10]. Additionally, *Populus* is an economically important wood tree, and in recent years, there has been increasing interest in studying its genotype, transcriptome and drought response mechanisms [11,12]. *Populus* is the internationally accepted model system for physiological and molecular studies in woody plants, in part due to the availability of the complete genome sequence of *Populus trichocarpa* [13]. The availability of the poplar genome plays an important role in understanding the molecular processes of growth, metabolism, and stress responses to environmental changes. The *Populus alba* × *Populus glandulosa* "Clivus" clone is a hybrid species made by NIFoS (National Institute of Forest Science) in Korea. It is a good resource for genetic engineering as it is sterile [14]. The *Populus euramericana* "I-476" clone is an interspecific hybrid produced from the cross of *Populus nigra* and *Populus deltoides*. Many *P. euramericana* clones have been commercialized, used in forestry production, and to promote ecosystem stability [15].

Next-generation sequencing technologies have been used to elucidate the molecular bases of poplar physiological and developmental mechanisms [16,17], and the responses of the *Populus* transcriptome to both biotic and abiotic stresses [18–20]. In particular, RNA-Seq analysis have been an important breakthrough for sensitive, quantitative, annotationindependent, and high-throughput analyses [21].

We aimed to study physiological and transcriptional changes in poplar trees grown under ambient (400 ppm) and elevated (560 and 720 ppm) atmospheric CO2, and the physiological parameters are correlated with changes in transcript abundance.

#### **2. Materials and Methods**

#### *2.1. Plant Materials and Growth Conditions*

Poplar seedlings *(Populus alba* × *Populus glandulosa* hybrid "Clivus" clone and *Populus euramericana* "I-476" clone) were grown in pots containing appropriate soil moisture in the open-top chambers (OTCs) system (Figure 1). These *Populus* clones were used as solitary maternal plants. Rooted cuttings from these maternal *Populus* clones were cultivated in pots in a greenhouse for 4 weeks for acclimation and then transferred to OTCs. Ten plants were grown and 8 plants were used in the experiments. The experiment was conducted for 16 weeks in the National Institute of Forest Science in Suwon, Korea (37◦15 04 N, 136◦57 59 E), under natural environmental conditions [22]. Three treatment levels of CO2 concentration were applied to the OTCs: ambient (×1.0, ~400 ppm), ×1.4 (~560 ppm), and ×1.8 (~720 ppm). Although the air temperature inside was 1.2–2.0 ◦C higher than that of the outside, the temperature differences among the OTCs were less than 0.2 ◦C.

#### *2.2. Measurement of the Chlorophyll Content*

The chlorophyll content was determined according to the method of Sibley et al. (1996) [23]. From all treatments, 0.1 g fresh samples were taken in triplicate, homogenized thoroughly with dimethyl formamide (DMF) and centrifuged at 14,000× *g* for 10 min at 4 ◦C. The supernatant was used as the chlorophyll source. The chlorophyll levels were determined by reading the supernatant absorbances at 647 nm and 664 nm with a Biospectrometer (Eppendorf, Hamburg, Germany). The chlorophyll contents and their means were calculated as follows for each plant and treatment:

• Carotenoids = (1000A470 − 1.82Chl a − 85.02Chl b)/198. (A, absorbance; pigment concentration in mg/g fresh weight (FW)).

**Figure 1.** A schematic representation of the experimental design and the phenotypes of poplar clones. (**a**) A photograph of the whole open-top chambers (OTCs) facility. (**b**) The conceptual model. (**c**,**d**) The effects of elevated CO2 on the shoot growth of the poplar Clivus and I-476 clones. The values are the means ± SDs (*n* = 10) Different lowercase letters indicate significant differences (Tukey's HSD, *p* < 0.05 for (**d**)).

#### *2.3. Extraction and Measurement of Starch and Soluble Sugar*

To extract metabolites, approximately 0.1 g poplar leaves from the plants grown for 16 weeks were ground in liquid nitrogen to a fine powder, and then the pulverized tissues were extracted twice with 80% (*v/v*) ethanol at 80 ◦C. To analyze starch content, the resulting sediments from aqueous ethanol extractions were re-suspended in distilled H2O and enzymatically digested to glucose according to the method described by Walters et al. (2004) [24]. The sugar concentrations were determined enzymatically by a method described by Stitt et al. (1989) [25] using a Biospectrometer (Eppendorf, Hamburg, Germany).

Total soluble sugars were extracted from leaf tissues by 80% ethanol by a modified method of Irigoyen et al. (1992) [26] as follows. After fresh weight determination, the leaves were homogenized by grinding in liquid nitrogen with mortar and pestle. Then, 2 mL of 80% (*v/v*) ethanol was added and the sample was vortexed for 1 h. After centrifugation at 14,000× *g* for 10 min, the supernatant was collected. The supernatants were added with chloroform and completely mixed. After centrifugation at 14,000× *g* for 10 min, 50 μL of supernatant was reacted with 4.95 mL of freshly prepared anthrone reagent (500 mg anthrone + 50 mL 72% H2SO4) at 100 ◦C for 15 min. After cooling on ice, the total soluble sugar content was determined at 620 nm by a Biospectrometer (Eppendorf, Hamburg, Germany) using glucose as the standard.

#### *2.4. Measurement of Malondialdehyde (MDA), Proline and H2O2*

From each treatment, approximately 0.1 g of poplar leaves were harvested in triplicate and extracted with a buffer consisting of 20% TCA (*w/v*) and 0.5% thiobarbituric acid (TBA) (*w/v*), followed by warming at 95 ◦C for 30 min. The reaction was terminated by placing the mixture on ice for 30 min and then centrifuging it at 14,000× *g* for 10 min. The absorbance of the supernatant was read at 532 nm using a Biospectrometer (Eppendorf, Hamburg, Germany). The MDA content was derived according to the method of Heath and Packer (1968) [27].

Proline was extracted from a sample of 0.5 g fresh leaf samples in 3% (*w/v*) aqueous sulfosalicylic acid and estimated using the ninhydrin reagent according to the method of Bates et al. (1973) [28]. The absorbance of the fraction with toluene separated from the liquid phase was read at a wavelength of 520 nm. Proline concentration was determined using a calibration curve and expressed as mmol proline g−1 FW.

Hydrogen peroxide was measured spectrophotometrically after reaction with KI. The reaction mixture consisted of 0.5 mL 0.1% trichloroacetic acid (TCA), leaf extract supernatant, 0.5 mL of 100 mM K-phosphate buffer and 2 mL reagent (1M KI *w/v* in fresh distilled H2O). The blank probe consisted of 0.1% TCA, K-phosphate buffer, and KI reagent in the absence of leaf extract. The reaction was developed for 1 h in darkness and absorbance was measured at 390 nm. The amount of hydrogen peroxide was calculated using a standard curve prepared with known concentrations of H2O2 [29].

#### *2.5. RNA Isolation, Library Preparation and qRT-PCR Analysis*

Samples obtained from the leaves of eight plants in each CO2 treatment were frozen in liquid nitrogen, grinded, and stored at −80 ◦C. RNAs from 8 plants for each treatment were extracted, and equal amounts of RNA were pooled to obtain single RNA samples; RNA pools were used for RNA sequencing. cDNAs for qRT-PCR were synthesized from all 48 samples. Total RNA was isolated using a RNeasy plant mini kit (Qiagen, Hilden, Germany). Approximately 2 μg of RNA from each tissue was used to construct cDNA libraries for sequencing according to the Illumina TruSeq RNA Sample Preparation protocol. In short, the workflow included isolation of poly-adenylated RNA molecules using poly-T oligo-attached magnetic beads, enzymatic RNA fragmentation, cDNA synthesis, ligation of bar-coded adapters, and PCR amplification. The libraries were sequenced using an Illumina HiSeq 4000 platform with 101 paired-end sequences at the Macrogen (Korea).

For real-time quantitative RT-PCR analysis, first-strand cDNA was synthesized from 1 μg of DNase-treated total RNA using RNA to cDNA EcoDryTM Premix (TaKaRa, Shiga, Japan). All reactions were performed using the IQtm SYBR Green Supermix (BIO-RAD, Hercules, CA, USA) and carried out in a CFX96 Touch Real-Time PCR Detection System (BIO-RAD, Hercules, CA, USA) according to the manufacturers' instructions. The genespecific primers used for the quantitative RT-PCR were found in Supplementary Table S1. The reaction cycle was: 1 cycle of 95 ◦C for 30 s, followed by 40 cycles of 95 ◦C for 5 s, and 60 ◦C for 34 s. Relative quantification was performed to calculate expression levels of target genes in different treatments using the 2−ΔΔCt methods [30]. In the reaction plate, each sample was measured in triplicate. The expression level of *ACTIN1* was used for the normalization of quantitative real-time PCR results [31].

#### *2.6. RNA Sequencing Analysis*

The sequenced raw reads from *Populus* young leaf tissues were processed with the following RNA sequencing pipeline. The quality of the RNA-Seq reads from all of the six tissues was checked using FastQC. The Trimmomatic software (v0.0.14) was used to process raw reads to remove adaptor sequences and low quality reads [32]. The obtained high quality reads were mapped into *Populus alba* reference genome (GCA\_005239225.1\_ ASM523922v1\_genomic.fna) from NCBI (https://www.ncbi.nlm.nih.gov/genome/13203, accessed on 15 January 2020) using Hisat2 [33]. Further, SAMtools was implemented to convert SAM files from BAM files [34]. Then, we used Feature Counts [35] to estimate the

uniquely mapped gene counts by using *Populus alba* genome annotation (GCA\_0052392 25.1\_ASM523922v1\_ genomic.gff).

#### *2.7. Differentially Expressed Gene (DEG) Analysis*

The generated gene count files from each *Populus* tissue (Clivus\_ OTC1 (C\_01), Clivus\_OTC2 (C\_02), Clivus\_OTC3 (C\_03), I-476\_OTC1 (I\_01), I-476\_OTC2 (I\_02), and I-476\_OTC3 (I\_03)) were used for differential expression genes analysis using DEseq2 version 1.20 package in R analysis environment [36]. Further, Fisher's exact test and likelihood ratio test methods were implemented to perform differential expression analysis of digital gene expression data following a binomial distribution. This experiment lacks biological replicates and log2 fold change ≥2 was used as a threshold to identify the DEGs. Further, only DEGs with a gene symbol were considered for downstream analysis. Common DEGs in both clones from all comparison (400 vs. 560, 400 vs. 720, and 560 vs. 720) were extracted using Venny 2.1.0 for further understanding.

#### *2.8. Annotation of DEGs*

The Omicsbox (Blast2GO) program (Biobam, Purchased version, Valencia, Spain) was used to perform functional annotation analysis. In detail, we performed a BLASTx programbased homology search for the *Populus alba* gene sequences against the Arabidopsis protein database (NCBI Arabidopsis protein sequences https://ftp.ncbi.nlm.nih.gov/genomes/ all/GCF/000/001/735/GCF\_000001735.4\_TAIR10.1/GCF\_000001735.4\_TAIR10.1\_protein. faa.gz, accessed on 3 February 2020) using a cutoff E-value < 10−5 and a maximum number of allowed hits were fixed at 10 per query. The alignment results with the smallest E-value were considered to select best top hits. Further, GO (Gene Ontology) terms were identified for all genes which are associated from BLASTx (NCBI, https: //blast.ncbi.nlm.nih.gov/, accessed on 11 February 2021) search results with Omicsbox (Biobam, Spain, https://www.biobam.com/omicsbox/, accessed on 11 February 2021) and WEGO (Web Gene Ontology Annotation Plot) (http://wego.genomics.cn, accessed on 4 October 2020) used for plotting [37]. The corresponding Arabidopsis gene symbol, along with DEG specific fold values, were used to perform gene set enrichment analysis by WebGestalt (http://www.webgestalt.org/, accessed on 2 May 2020) [38]. Cluster analysis was performed using ClustVis (https://biit.cs.ut.ee/clustvis/, accessed on 11 February 2021) [39].

#### *2.9. Statistical Analysis*

Analyses were carried out using a one-way ANOVA with multiple comparisons using Dunnett's T3 or Tukey's HSD. *p*-values < 0.05 were considered significant. Values are presented as means with SD.

#### **3. Results**

#### *3.1. Growth and Physiological Changes in Response to Elevated CO2*

We observed the phenotypes and measured the physiological changes of the two poplars' leaves at the seedling stage. The height of the Clivus grown with an elevated CO2 level of 720 ppm was significantly reduced as a response to elevated CO2 (Figure 1). However, except this instance, there were no significant physiological differences observed between ambient or elevated CO2 poplar plants. The leaves' chlorophyll content was correlated with photosynthetic activity, and the chlorophyll levels were affected by elevated CO2. The chlorophyll A content of I-476 poplar was significantly reduced at elevated CO2 levels (560 and 720 ppm). Furthermore, in I-476, reduced chlorophyll: carotenoid content was also observed at elevated CO2 (720 ppm) compared to the controls (Figure 2). In addition, the total chlorophyll content of I-476 poplar was also significantly reduced at elevated CO2 levels. In both Clivus and I-476, there was a change in carotenoids at an elevated CO2 level of 560 ppm. Except for the above-mentioned changes, we found no

significant changes in the levels of photosynthetic pigments in poplar grown in any CO2 conditions.

**Figure 2.** Effects of CO2 treatment on photosynthetic pigments in poplar clones. (**a**) Chlorophyll A. (**b**) Chlorophyll B. (**c**) Carotenoids. (**d**) Total chlorophyll/carotenoids. (**e**) Total chlorophyll. (**f**) Chlorophyll A/B. Values are means ± SDs of three independent measurements. Different lowercase letters indicate significant differences (Tukey's HSD, *p <* 0.05 for (**a**,**c**,**e**); and Welch's ANOVA with Dunnett's T3, *p <* 0.05 for (**d**)).

To investigate the effects of elevated CO2 on the carbon metabolism and stressresponsive elements, the MDA, H2O2, proline, and soluble sugar contents were examined. The MDA content was measured to determine the degree of tissue damage caused by stress [40]. The change in H2O2 content is a good indicator of the state of the ability to remove free radicals in oxidation stress [40]. Elevated CO2 increased the concentrations of MDA, proline, total soluble sugar, and starch in leaves of the Clivus, but did not alter the concentration of H2O2. In contrast, elevated CO2 significantly increased the concentration of H2O2, proline, total soluble sugar, and starch in leaves of the I-476. Only the level of MDA at 520 ppm decreased in I-476. In particular, elevated CO2 tended to increase the concentration of starch in both the Clivus and I-476 (Figure 3).

**Figure 3.** Effects of CO2 treatment on the metabolites and carbohydrate contents. (**a**) MDA. (**b**) H2O2. (**c**) Proline. (**d**) Total soluble sugar. (**e**) Starch. Different lowercase letters indicate significant differences (Tukey's HSD, *p* < 0.05 for I-476 in (**a**–**e**); and the same for Clivus, except Welch's ANOVA with Dunnett's T3, *p* < 0.05 for Clivus in (**a**)).

#### *3.2. RNA Sequencing Analysis*

RNA-Seq analysis was performed in order to study the differential gene expression patterns exerted in response to elevated CO2 treatments in two poplar species (Clivus and I-476). After the adapter trimming and quality check, the clean reads were mapped to the reference genome *Populus alba*. The mapping percentage varied significantly between the species; Clivus (*P. alba* × *P. glandulosa* hybrid) showed higher percentage matching than the I-476 (*P. euramericana*), as expected. The detailed report on raw read number, data in GB, clean read number and clean read rate for each comparison are detailed in Table 1. The quality score, Q30 clean base rate, was approximately 94% in all comparisons (Table 1).

**Table 1.** A summary of the RNA sequencing results.

#### *3.3. Differentially Expressed Genes (DEGs)*

The DEG analysis was used as an exploratory data to identify the genes responsive to elevated CO2 treatment (Supplementary Figures S1 and S2). Elevated CO2 treatment displayed a change in expressed genes in both the clones (Figure 4). Clivus clone showed a greater number of DEGs than I-476 in all comparisons (400 vs. 560, 400 vs. 720, and 560 vs. 720). The higher percentage of mapping observed in Clivus clone to *Populus alba* reference genome might be the reason for the distinct number of DEGs observed between clones. DEGs that were consistently expressed genes in both elevated CO2 treatment could be candidate responsive genes for elevated CO2 treatment in poplar clones. Therefore, we extracted common DEGs from all comparisons. A total of 29 common genes were observed in Clivus, of which two were upregulated and 24 downregulated (Figure 4). Most of these DEGs belong to the cell wall protein family, especially hydroxyproline-rich glycoproteins (HRGP) [41,42], and the fatty acid metabolism pathway (Figure 5a). HRGP are proline-rich cell wall proteins that have a wide range of functions in signal transduction cascades, such as plant development and stress tolerance. As observed in Figure 5a, the last three genes showed increased or decreased expression in comparison 400 vs. 560, while increasing the CO2 concentration (400 vs. 720) did not show concomitant increase or decrease. Similarly, in I-476, there were seven genes upregulated and eight genes downregulated in both the elevated CO2 treatment compared to ambient CO2 concentration, and the last 10 genes showed varied expression pattern when CO2 increased. Common DEGs obtained in the I-476 clone belonged to the HSP chaperone family and metabolic pathways, especially photosynthesis, nitrogen metabolism, and carbon metabolism (Figure 5b).

**Figure 4.** Comparison of differentially expressed genes. (**a**) Clivus (*Populus alba* × *P. glandulosa*, a poplar clone hybrid); and (**b**) I-476 (*Populus euramericana*). Upregulated-log2 fold change (≥2); Downregulated-log2 fold change (≤−2). Each species was treated with three different concentrations of CO2. C-400, Clivus CO2 400 ppm; C-560, Clivus CO2 560 ppm; C-720, Clivus CO2 720 ppm; I-400, I-476 CO2 400 ppm; I-560, I-476 CO2 560 ppm; I-720, I-476 CO2 720 ppm.

#### *3.4. Functional Classification of Two Populus Species*

The GO study displayed a similar pattern of enrichment in both the poplar species. Under the biological process categories, the greater number of DEGs falls under response to stress, responsive to stimulus and external stimulus sub-classes. The CO2 treatment has triggered the responsive stimulus in the leaves to a greater extent, which was well observed in GO enrichment analysis (Figure 6). Additionally, HGRP and other membrane-associated

DEGs found in elevated CO2 treated leaves correlated with the plant response in protecting its cell wall.

The KEGG enrichment analysis on DEGs resulted in unanimity pathways in all comparisons, irrespective of species differentiation. Carbon metabolism, starch sucrose metabolism, glutathione metabolism, and fatty acid metabolism are the most common enriched pathways in top 10 lists (Figure 7).

**Figure 5.** Heat map showing the gene expression pattern of common DEGs in (**a**) Clivus and (**b**) I-476. Corresponding log2 fold change is displayed inside the box. Each species was treated with three different concentrations of CO2. C-400, Clivus CO2 400 ppm; C-560, Clivus CO2 560 ppm; C-720, Clivus CO2 720 ppm; I-400, I-476 CO2 400 ppm; I-560, I-476 CO2 560 ppm; I-720, I-476 CO2 720 ppm.

**Figure 6.** *Cont.*

**Figure 6.** Gene Ontology (GO) enrichment analysis of DEGs. (**a**,**b**) Clivus; (**c**,**d**) I-476. Grey and red represent the upregulated and downregulated DEGs, respectively.

**Figure 7.** KEGG pathway enrichment analysis. (**a**,**b**) Clivus; (**c**,**d**), I-476.

#### *3.5. Validation of RNA-Seq Results by qRT-PCR*

We performed a qRT-PCR for the top 10 most significant DEGs, including D5086\_0000 265220, D5086\_0000080910, D5086\_0000065160, D5086\_0000242970, D5086\_0000065170, D5086\_0000289410, D5086\_0000064330, D5086\_0000042670, D5086\_0000250770 and D5086\_<PHONE_NUMBER> in Clivus (Figure 8a–e,k–o); and D5086\_0000176630, D5086\_0000043470, D5086\_<PHONE_NUMBER>, D5086\_0000215870, D5086\_0000038460, D5086\_0000212360, D5086\_0000328260, D5086\_0000281500, D5086\_0000305900, and D5086\_0000277080 in I-476 (Figure 8f–j,p–t). To verify the differences in gene expression among OTC1, OTC2, and OTC3, 20 genes were selected for qRT-PCR. In qRT-PCR assay, most of the genes showed similar expression patterns as were observed in RNA-Seq data. We found that common DEGs, D5086\_0000080910, D5086\_0000065160, D5086\_0000242970, and D5086\_0000065170, were downregulated, while D5086\_0000265220 was upregulated in Clivus species in both OTC2 and OTC3 compared with OTC1. The common DEGs of I-476, except the upregulated D5086\_0000038460, D5086\_0000176630, D5086\_0000043470, D5086\_0000080180, and D5086\_0000215870, were all downregulated in both OTC2 and OTC3 compared with OTC1. We analyzed the sucrose related genes of in both Clivus and I-476; compared with OTC1, both OTC2 and OTC3 showed upregulation of DEGs D5086\_0000064330, D5086\_0000042670, D5086\_0000250770, and D5086\_0000252880, whereas D5086\_0000289410 was downregulated in Clivus. In the case of I-476, D5086\_0000212360, D5086\_0000328260, D5086\_0000281500, D5086\_0000305900, and D5086\_0000277080 were downregulated in both OTC2 and OTC3 compared with OTC1. Among the above results, the decreased expression levels of D5086\_0000042670 and D5086\_000305900 in OTC2 and OTC3 were opposite to each other when compared with RNA-Seq data.

**Figure 8.** Validation of differential expressions of 20 genes by quantitative real-time PCR (qRT-PCR). qRT-PCR analysis data for the top 10 most significantly differentially expressed genes (DEGs) in Clivus and I-476, respectively. qRT-PCR data were

analyzed using the 2−ΔΔCt method with Actin1 gene as an internal control. Three biological replicates were performed for each sample. Error bars show mean standard error (*n* = 3). (**a**–**e**) Common Clivus genes; (**f**–**j**) common I-476 genes; (**k**–**o**) sucrose genes of Clivus; and (**p**–**t**) sucrose genes of I-476. Different lowercase letters indicate significant differences. (**a**, DREB2A-interacting protein 2/E3 ubiquitin protein ligase DRIP2-like; **b**, ABI-1-like 1/probable protein ABIL5; **c**, fatty acid desaturase 2/delta(12)-acyl-lipid-desaturase-like; **d**, myb domain protein 9/transcription factor MYB53-like; **e**, fatty acid desaturase 2/delta(12)-acyl-lipid-desaturase-like; **f**, HSP20-like chaperones superfamily protein/23.6 kDa heat shock protein, mitochondrial-like; **g**, galactinol synthase 1/galactinol synthase 2-like; **h**, Ribulose bisphosphate carboxylase (small chain) family/ribulose bisphosphate carboxylase small chain clone 512-like protein; **i**, HSP20-like chaperones superfamily protein/22.0 kDa class IV heat shock protein-like; **j**, glucuronoxylan 4-O-methyltransferase-like protein (DUF579); **k**, Glycogen/starch synthases, ADP-glucose type; **l**, glycosyl hydrolase 9C3; **m**, glycosyl hydrolase 9B13; **n**, glycosyl hydrolase 9B18; **o**, Glycosyl hydrolases family 31 protein; **p**, glycosyl hydrolase 9C3; **q**, Glycosyl hydrolase superfamily protein; **r**, plant glycogenin-like starch initiation protein 1; **s**, Glycosyl hydrolase superfamily protein; and **t**, plant glycogenin-like starch initiation protein 3) (Tukey's HSD, *p* < 0.05 for (**e**,**f**,**h**,**i**,**k**,**l**,**r,t**); and Welch's ANOVA with Dunnett's T, *p* < 0.05 for (**a**–**d**,**g**,**m**–**q**,**s**).

#### **4. Discussion**

As a rapid increase in CO2 is expected within decades, studies on the reaction of plants to high concentrations of CO2 have been steadily progressing. As the initial growing environmental conditions of the seedlings are very important for development and survival, a response study to high concentrations of CO2 was conducted. The physiological response of poplar to elevated CO2 varies greatly depending on the species, growth period, additional nutrients, and temperature [7,43–46]. We observed changes in Clivus and I-476 poplar seedlings in OTCs without any additional factor intervention. However, OTCs have an intrinsic constraint in that the temperature (by 1.2–2 ◦C) is higher than outside due to their limited airflow. Growth deterioration due to spatial constraints of OTCs did not affect growth at the seedling level. Under these conditions, only Clivus clone at an elevated CO2 level of 720 ppm showed a decrease in growth (Figure 1). In the previous study, impaired plant growth was observed after long-term exposure to a high CO2 concentration and insufficient nitrogen by altering the primary metabolism [47].

Previous studies have shown that the chlorophyll content may decrease or increase in response to changes in the external environment [48]. In our result, it was confirmed that levels of chlorophyll A and carotenoids showed significant changes (Figure 2). It has been reported that the content of chlorophyll is highly correlated with the concentration of CO2 in the atmosphere and the photosynthetic ability [49,50]. An elevated CO2 concentration initially increases the photosynthetic rates of plants, but over time, the rate of photosynthesis is reduced due to the feedback of photosynthetic products and reduction in Rubisco content [51]. Referring to the above, our result showed that treatment with a high concentration of CO2 for 16 weeks passed the initial stage and reached the mid-term or long-term reaction. In particular, a decrease in chlorophyll A, which is involved in carbon fixation, is suggested by the feedback of accumulated soluble sugar and starch [50,52]. In a recent report, it was shown that in *Populus* grown at high CO2 concentrations, the content of carotenoids slightly decreased or remained unchanged [53]. This was in accordance with the reduction in the level of carotenoids in our experiment (Figure 2). The nitrogen content of the leaves can decrease and photosynthetic products do not move to the sink organs and accumulate as starch, reducing the rate of photosynthesis [51]. This photosynthetic acclimation is also consistent with the large increase in starch content at elevated CO2 concentrations in our results (Figure 3). In addition, we observed a decrease in the expression of glycosyl hydrolase genes, which is related to the degradation of starch, and is also consistent with the result of increasing starch content [51,54]. It is well known that an increase in the soluble sugar content in cells inhibits photosynthetic genes [45,55]. In both Clivus and I-476 seedlings we found an increase in the soluble sugar at elevated CO2 concentrations, and this was associated with a decrease in the level of chlorophyll A in both poplars (Figures 2 and 3). The presence of many DEGs in sucrose and starch

metabolic pathways and the related physiological changes observed have confirmed the influential changes in carbon partitioning affecting metabolic changes in poplar under high CO2 concentration.

Plants grown at elevated CO2 concentrations have been reported to be more tolerant to drought stress, as stomata can be kept small due to high CO2 pressure [43]. We confirmed that the overall level of proline increased at elevated CO2, which is thought to be helpful in improving drought tolerance by affecting the intracellular osmotic pressure regulation of proline (Figure 3). In this study, Clivus showed greater changes in carbon partitioning and stress indicators than I-476. These observations shown in Clivus were operating more pronouncedly at mildly (520 ppm) elevated CO2 concentrations. As reported previously, H2O2 is produced primarily by plant cells during photosynthesis and photorespiration, and affected by environmental stresses [31]. In our results, on the contrary, in I-476, H2O2 was increased even though photosynthetic pigment was decreased at an elevated CO2 concentration. The alterations in the stress-related indicators at elevated CO2 concentrations might be involved in overlapping stress and elevated CO2 regulatory networks [56].

Here, we have studied the genetic response to elevated CO2 on two poplar clones. The transcriptome profiling of poplar clones in elevated CO2 was utilized as an exploratory data to find the gene response pattern. Due to the lack of statistical testing, transcriptome data from this study has been interpreted with caution. We extracted common DEGs in both clones that are responsive for the elevated CO2 treatments. The common DEGs from Clivus pooled mainly in HGRP and fatty acid metabolism, whereas for I-476 they were mainly HSPs and from diverse metabolic pathways. This correlates with the experimental data, where in Clivus there is increased MDA production, reflecting the transcriptome analysis showing more DEGs related to fatty acid metabolism. Lesser DEGs related to MDA production were seen in I-476, which explains the lower MDA levels in this clone (Figure 4). The HGRP-related DEGs correlates with the free sugar and proline increase at elevated CO2 in Clivus [57]. Small HSP maintain ROS homeostasis and cooperate with other stressrelated genes in response to diverse biotic and abiotic stresses to protect plant cells and tissues [58]. In I-476, DEG analysis showed an evident HSP chaperone-based defense mechanism involving scavenging of H2O2. By treating changes in CO2 concentration, the stress indicators such as MDA, proline, H2O2, and expression of HSP-related genes were changed. In Clivus and I-476, elevated CO2 concentrations significantly decreased the expression of sucrose synthase (SUS) genes. This reduction in expression level is related to the fact that the regulation of photosynthetic amount was due to photosynthetic acclimation and influenced the synthesis of sucrose, the first product of photosynthesis. The plant glycogenin-like starch initiation protein (PGSIP)-like gene's expression was also decreased in both Clivus and I-476. This seems to suggest that the decrease in starch synthesis genes is related to the reduction in photosynthetic pigments and feedback regulation of accumulated starch in leaves. Although, the transcriptome analysis throws light on increased CO2 concentration induced gene expression changes in these poplar clones, the absence of biological replicates in transcriptome analysis warrants further studies for increased specificity.

#### **5. Conclusions**

To summarize, our results present a transcriptome dataset and physiological analysis, which contributes to the understanding of gene expression profiling in response to elevated CO2. By comparing the transcriptome data of three different concentrations of CO2, we found that many DEGs were identified as unique to Clivus and I-476. In addition, this study presented significant physiological alterations with respect to chlorophyll A and carotenoids. Additionally, we reported alterations in starch and soluble sugar. Several synthesis and degradation-related sucrose and starch genes showed altered expression dynamics in our transcriptome and qRT-PCR analysis. The transcriptome data and physiological response of poplars at elevated CO2 presented here could serve as a resource for further research on poplar trees in different CO2 environments.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/f12080980/s1, Figure S1: MA plot of DEGs, Figure S2: Volcano plot of DEGs, Table S1: Primer of qRT-PCR in Clivus and I-476.

**Author Contributions:** Conceptualization, H.L. and T.-L.K.; methodology, H.L.; software, H.C.; validation, T.-L.K. and D.P.; formal analysis, K.V.; investigation, T.-L.K.; resources, W.Y.L.; data curation, H.C.; writing—Original draft preparation, T.-L.K.; writing—Review and editing, H.L.; visualization, T.-L.K.; and supervision, H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Data is contained within the article or Supplementary Materials.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Field Performances of Mediterranean Oaks in Replicate Common Gardens for Future Reforestation under Climate Change in Central and Southern Europe: First Results from a Four-Year Study**

**Filippos Bantis 1,\*, Julia Graap 1, Elena Früchtenicht 1, Filippo Bussotti 2, Kalliopi Radoglou 3 and Wolfgang Brüggemann 1,4**

**Abstract:** Climate change imposes severe stress on European forests, with forest degradation already visible in several parts of Europe. Thus adaptation of forestry applications in Mediterranean areas and central Europe is necessary. Proactive forestry management may include the planting of Mediterranean oak species in oak-bearing Central European regions. Five replicate common gardens of Greek and Italian provenances of *Quercus ilex*, *Q. pubescens* and *Q. frainetto* seedlings (210 each per plantation) were established in Central Italy, NE Greece (two) and Southern Germany (two, including *Q. robur*) to assess their performance under different climate conditions. Climate and soil data of the plantation sites are given and seedling establishment was monitored for survival and morphological parameters. After 3 years (2019) survival rates were satisfactory in the German and Italian sites, whereas the Greek sites exerted extremely harsh conditions for the seedlings, including extreme frost and drought events. In Germany, seedlings suffered extreme heat and drought periods in 2018 and 2019 but responded well. Provenances were ranked for each country for their performance after plantation. In Greece and Italy, *Q. pubescens* was the best performing species. In Germany, *Q. pubescens* and *Q. robur* performed best. We suggest that Greek or Italian provenances of *Q. pubescens* may be effectively used for future forestation purposes in Central Europe. For the establishment of *Quercus* plantations in Northern Greece, irrigation appears to be a crucial factor in seedling establishment.

**Keywords:** *Quercus*; morphology evaluation; survival rate; extreme frost; heat and drought

#### **1. Introduction**

Climate change scenarios bring forward particularly challenging effects on European forests leading to the necessary development of new strategies for the maintenance of functional forests. Extreme temperatures and severe drought events are predicted to increase in frequency, like it has previously been recorded in 2018 and 2019 across Europe [1,2]. Environmental changes can be so abrupt that natural replacement of forest tree species with better adapted species may not occur at sufficient speed to prevent severe forest degradation. Southern regions may lose large areas of their forest cover altogether, except for their higher mountain ranges. The total net financial loss of the value of European forest land from these shifts is estimated to be up to hundreds of billions Euro in the period to 2100 [3], while ecological implications are also substantial.

**Citation:** Bantis, F.; Graap, J.; Früchtenicht, E.; Bussotti, F.; Radoglou, K.; Brüggemann, W. Field Performances of Mediterranean Oaks in Replicate Common Gardens for Future Reforestation under Climate Change in Central and Southern Europe: First Results from a Four-Year Study. *Forests* **2021**, *12*, 678. https://doi.org/10.3390/f12060678

Academic Editor: Carlos Gonzalez-Benecke

Received: 2 March 2021 Accepted: 24 May 2021 Published: 26 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Proactive forestry management may include the planting of Mediterranean oak species in oak (*Quercus robur* L., *Q. petraea* Matt.) bearing Central European regions [4,5]. This approach aims at the establishment of (at least) seed donor and (at best) timber plantations of future climate-tolerant species within existing oak stands to ensure the continuous presence of closely related oak species, thus maintaining oak-dominated ecosystem functions including oak-associated biodiversity. Such a concept is currently realized e.g., in the South Hesse Oak Project (SHOP [5]). However, since Central Europe and low/medium elevation forests in the Mediterranean will face–in principle–the same problems in the forthcoming century, this approach may be extended to forest management strategies in mountain regions in Mediterranean countries as well.

One way of testing the response of species or populations under certain environmental conditions is by establishing common gardens in several locations of ecological interest [6]. Common gardens are considered essential in determining the potential of selected species and/or provenances [7] to be implemented for 'assisted migration' due to climate change [8–10].

For reforestation purposes, different strategies can be applied. One option is to start the process by active seed dispersal, usually leading to healthy plants with welldeveloped root systems, especially in tap rooting species like most oaks [11]. However, buried seeds are often prone to herbivory in the field (e.g., Pausas et al. [12]). Therefore, foresters, especially in Germany, often rely on pre-growth of seedlings in tree nurseries and transplanting established containered seedlings into the field sites. This strategy minimizes the need for seed material, but implies additional working steps and costs. To further minimize potential losses, various methods of fast and easy selection of vigorous seedlings have been suggested (for review see e.g., Mattson [13]). For example, in the case of the Mediterranean oak species *Quercus ilex* and *Q. coccifera*, Tsakaldimi et al. [14] tested the use of morphological parameters like height, leaf area or stem diameter as easily measurable parameters for the prediction of seedling survival in the first year after outplanting. Since Mattsson [13] also suggested the use of chlorophyll content and fluorescence as potential, easily measurable predicitive parameters, we attempted to verify their usefulness for the prediction of first-year survival in the field under different climate conditions.

The objectives of the present study were to evaluate the survival and growth responses of seedlings from three oak species established on experimental plantations exposed to different environmental conditions including climate, water availability and soil properties. Specifically, *Quercus ilex* L., *Quercus frainetto* Ten., and *Quercus pubescens* Willd. from Greek and Italian provenances were grown in replicate common gardens in Southern Germany, Central Italy, and Northern Greece. The German plantations also include a local oak species, *Quercus robur* L. for comparative purposes.

The specific questions addressed were:

#### **2. Materials and Methods**

#### *2.1. Plant Material*

In autumn 2015, seeds of three common Mediterranean oak species, *Q. ilex*, *Q. frainetto* and *Q. pubescens* were collected in Central Italy and Northeastern Greece (Supplementary Material Table S1). Seeds from all species were shipped by air mail to commercial tree nurseries in Germany, Italy and Greece, and sown in plastic trays (Quickpot 15T/16 or 24/18). In each country, the trays were filled with different substrate according to the

regional practices in order to achieve optimum growth of oaks. Specifically, in Germany trays were filled with a mixture of humus, volcanic lapilli and different organic fibres (FE Typ Bio Blumenerde # 03010, HAWITA Co., Vechta, Germany), in Italy trays were filled with a 40/40/20 mixture of peat, loam and volcanic lapilli along with 25 g·L−1 pelleted humified manure and2gL−1 Osmocote Exact fertilizer, and in Greece trays were filled with a 70/30 mixture of peat and black peat along with 12 kg m−3 clay. In Italy and Greece, the emerged seedlings overwintered outdoors, while in Germany they overwintered in a nonheated, but frost-free glasshouse. These procedures are common practice for the first winter with pot-seeded material in the different countries (in Germany, other frost-protection strategies are also applied). In Germany, 1-year-old *Q. robur* seedlings (provenance 81707 "Oberrheingraben", supplied by Darmstädter Forstbaumschulen, Darmstadt, Germany) were added to the ensemble in spring 2016. During 2016, all seedlings were grown outdoors with regular watering, while in Germany seedlings were grown under nettings in autumn. In autumn 2016, emerged seedlings were measured and scored with respect to growth and physiological parameters. In spring 2017, 420 (Germany and Greece) and 210 (Italy) seedlings per species and provenance with the highest score (cf. paragraph "Plant fitness parameter scoring and ranking" below) were transferred to the plantation sites. In each site, 210 seedlings per species and provenance were established, except for Italian *Q. frainetto* in the Greek plantations, where only 170 seedlings were available due to lower seedling emergence rates. In Greece, seedlings were planted on 24 February–4 March 2017, while in Italy and Germany seedlings were planted between 15 April and 1 May 2017.

#### *2.2. Common Garden Experiments*

Two plantations were established in Germany, two in Greece, and one in Italy. In Germany, one plantation site is located in a managed oak and pine forest, on fluvial sand, with about 2 m deep groundwater table, in Schwanheim, Frankfurt ("SWA", 50◦04 12.6 N, 08◦33 42.2 E, 114 m a.s.l.). SWA is a typical oak site where oak forests exist for at least 500 years. The other plantation is in Riedberg, Frankfurt ("RIE", 50◦10 12.8 N, 08◦37 53.2 E, 130 m a.s.l.) on loamy soil, where planted seedlings were irrigated when necessary, as defined by monitoring weather conditions and potential drought symptoms on the trees. RIE is a cleared agricultural field, natural vegetation here would be (mixed) beech forest. In Greece, one plantation was established in Olympiada, Chalkidiki ("OLY", 40◦36 33.6 N, 23◦45 05.0 E, 48 m a.s.l.) in natural *Q. ilex* and *Q. pubescens* forest stands with loamy soil, where water was solely provided by precipitation. The particular location of the OLY plantation is a former grazing site. The other plantation is in Stratoniki, Chalkidiki ("STR", 40◦31 05.5 N, 23◦47 19.2 E, 248 m a.s.l.) in a mixed forest in the *Q. pubescens* zone, with loamy soil, where water was also solely provided by precipitation. STR is a cleared forest site. In Italy, the plantation is established at Sant' Anatolia di Narco, Umbria, in the Central Apennine ("SAN", 42◦44 18.4 N, 12◦50 17.6 E, 290 m a.s.l.), surrounded by mixed forest including *Q. ilex* and *Q. pubescens*, on calcareous soil. SAN is a former tree nursery site. During 2017, the seedlings here were also irrigated when necessary. Moreover, weed growth was intensive due to precipitations, thus the sites were cleared 2 to 3 times per year.

Plantation schemes in the common gardens are shown in supplementary Figure S1. In each plantation, 10 plots per species and provenance were established consisting of 21 seedlings, in a 3/5/5/5/3 arrangement called "grouped" scheme (Figure S1D in supplementary material) [15–17]. In the case of Italian *Q. frainetto* in the Greek plantation, only 17 seedlings per plot were established due to the abovementioned limitation of good seedlings. This particular scheme is well suited for enhancing tree quality in oaks, as well as more economic compared to row plantation [15]. To avoid clines in the plantations, the selected trees from the evaluation/scoring procedure were distributed randomly across the plantation in a RCBD design. Except for SAN, the original individual tree IDs from the evaluation procedure in the nursery were recorded for each tree during the planting procedure, thus allowing comparison of evaluation scoring data from 2016 with data measured on trees in later years.

#### *2.3. Climate Data*

At RIE and SWA, climate data were recorded on-site with iMetos sm SMT280 weather stations (Pessl Instruments, Weiz, Austria). For SAN, data in Val Nerina, Casteldilago (coordinates 42◦34 58.8 N, 12◦46 30 E at 270 m a.s.l.) were recorded with a Davis Vantage Pro2 weather station (Davis Wetterstationen, Neumarkt, Germany). For OLY and STR, data were recorded at the closest by weather stations of Hellas Gold SA, at 40◦36 09.45 N, 23◦45 01.39 E, 28 m.a.s.l. and at 40◦31 03.42 N, 23◦47 47.96" E, 208 m a.s.l., respectively. Due to technical problems with recording of Tmax at STR, the Tmax values from the closest weather station at Arnaia (ARN: 40◦29 30.88 N, 23◦36 08.49 E, 565 m a.s.l.) were used in the STR/ARN diagram in Figure 1. From DOY 146 on in 2017, additional data were collected in the OLY and STR sites with (unshaded) HOBO onset 8K data loggers placed at ca. 1.5 m above ground to address the real heat stress which the seedlings experienced. Precipitation data reported for OLY and STR were all taken from the Hellas Gold weather stations.

**Figure 1.** Climate diagrams for SWA, SAN and OLY. Data were averaged for 1981–2010 (SWA: Frankfurt Airport, climate station 1420 of DWD (2018)), 2012–2018 (SAN) and 2008–2017 (OLY), respectively. Black lines represent the average temperature, while blue lines represent the average precipitation.

#### *2.4. Soil Analysis*

5–10 representative samples of 200 mL soil each, of the top 20 cm after removal of the top organic layer were randomly collected from various parts of the respective site. Soil pH was measured in 10 mM CaCl2 extracts according to DIN 19 865, part 1 (1977). Granulometric distribution was measured according to DIN 19 683 part 1 and 2 (1973). After rattle sieving with 2 mm sieves, fine soil was subjected to H2O2 treatment (for destruction of organic matter) and dispersed in 0.4 N Na4P2O7. After wet filtration for particle sizes of coarse (2000–630 μm), medium (630–200 μm) and fine sand (200–63 μm), smaller fractions were analyzed by sedimentation analysis according to Köhn [18]. Inorganic and organic (humus) carbon contents were analyzed according to DIN 19 684, part 2 (1977) with gas analyzers LECO RC-412 and LECO EC-12 (for acidic samples) from LECO Instruments (Mönchengladbach, Germany), respectively. Organic matter was calculated from carbon content of the organic fraction by multiplication with a factor of 1.72 [19]. Initial water contents and field capacity were determined gravimetrically in fresh and water saturated samples, respectively, before and after drying for three days at 105 ◦C [20].

#### *2.5. Morphological Attributes*

In summers of 2016, 2017, 2018 and 2019, seedling heights ("*h*" in the formulas given below) were measured with rulers from soil surface to the top of the longest stem (branch), i.e., in general to the top bud. Root collar diameter was measured with an electronic Vernier caliper at 5 mm above soil level. In 2016, also the numbers of emerged leaves per seedling ("*leaf no*") were counted for fitness determination.

#### *2.6. Relative Chlorophyll Content and Chlorophyll Fluorescence*

In general, relative chlorophyll contents were measured in summers of 2016, 2017, 2018 and 2019 with a SPAD 502 Plus chlorophyll meter (Konica Minolta Sensing Co., Munich, Germany) on vital, fully developed leaves emerged in the respective year (important for *Q. ilex*) and given as dimensionless numbers ("*SPAD*" in the formulas below). For comparison with data obtained in Greece in 2016 with an OPTI-Sciences (Hudson, NH, USA) CCM-200 chlorophyll meter, a comparative measuring row of 152 leaves from all oak species covering the observed range of values was performed and used for determination of a calibration coefficient (to obtain "calculated SPAD values" from the CCM-200 values).

Performance index (*PIabs* in the formulas below) is a parameter summarizing the effects of light trapping, quantum efficiency of reduction of QA and efficiency of electron transport from QA- to the intersystem carriers of the electron chain [21,22]. *PIabs* was determined in mid- to late summer, i.e., after full leaf development with a Pocket PEA chlorophyll fluorometer (Hansatech, King's Lynn, UK) either at night or in the early morning after dawn, on leaves pre-darkened for a minimum of 20 min with leaf clips, and data were analyzed with PEA Plus <IP_ADDRESS> (Hansatech) software. A detailed analysis of the results of these analyses is provided in a separate paper [2].

#### *2.7. Plant Fitness Parameter Scoring and Ranking*

To compare plant fitness, the morphological and physiological parameters of the emerged seedlings in summer 2016 were used to perform a ranking of the individual seedlings for selection of the best-suited individuals for planting. For this purpose, an arbitrary fitness indicator (*FI*) was calculated for each individual seedling *i* according to (1):

$$FI = 0.25 \ast \frac{h\_i}{h\_{\max}} + 0.25 \ast \frac{leaf \,no\_i}{leaf \,no\_{\max}} + 0.25 \ast \frac{SPAD\_i}{SPAD\_{\max}} + 0.25 \ast \frac{PIabs\_i}{PIabs\_{\max}}\tag{1}$$

with "i" indicating each individual seedling and "max" indicating the highest value observed for this parameter in the respective population (species/provenance). For the plantations, the 420 (Italy: 210) best performing seedlings were then selected. In Greece, for *Q. frainetto* from Italy only 340 seedlings were available due to low germination rates.

Each individual tree was labeled and could be traced over time in the common garden plantations. To assess the potential predictive value of the 2016 *FI* values for seedling survival and development, we grouped the individual seedlings into "survivors" and "nonsurvivors" (in 2017) and compared their respective 2016 *FI* values by t-test. Furthermore, within the surviving seedlings, the rank of each individual obtained in 2016 was compared to its rank in 2019 by Spearman's correlation analysis. The 2019 ranking values of the individual were obtained omitting leaf number and SPAD, thus weighing height and *PIabs* with a factor of 0.5 each. To compare population (species/provenance) performance at the different sites, we took survival rates into account. We calculated, for each site (or, in the case of Germany, for both sites combined, see Table 1), mean values of individual parameters ("x") for each population k (i.e., species/provenance), ("xk"), then identified the maximum mean values of all population ("xmax") and calculated an arbitrary Relative Population Fitness (*RPF*) score (for each site) according to (2):

$$RPF = 0.40 \ast \frac{survival\_k}{survival\_{max}} + 0.30 \ast \frac{h\_k}{h\_{max}} + 0.30 \ast \frac{PIabs\_k}{PIabs\_{max}} \tag{2}$$

**Table 1.** Relative fitness of populations (*RPF*, cf. Equation (2)) in the different countries and common garden sites. Data were obtained in summer 2019.

In addition, the root collar diameter in 2016 was compared (by *t*-test) between the abovementioned "survivors" and "non-survivors" (in 2017) in order to check the parameter's potential as an indicator of second-year survival as reported by Tsakaldimi et al. [14]. Moreover, the potential of seedling height used alone as a predictive parameter for the German and the Greek sites was tested.

#### *2.8. Statistical Analysis*

Statistical analysis was conducted using PRISM 2.0 software (GraphPad Software, San Diego, CA, USA). Specifically, Kruskal-Wallis test with Dunn's post-hoc test at a significance level of *p* = 0.05 was conducted to identify potential significant differences between populations (i.e., site/species/provenance combinations).

#### **3. Results**

#### *3.1. Soil Parameters*

The soil parameters are given in Table 2. The German sites are very different in soil properties, with more calcareous and loamy conditions at RIE and an acidic, sandy soil in the recently partially cleared site SWA with its very high organic matter in the topsoil. The Italian site is typically calcareous and characterized by good water holding capacity. Together with its situation in the river valley, it provides very good conditions for seedling growth. Both Greek sites show slightly acidic, loamy soil with good water holding capacity,

but since they both are far away from recent superficial water sources, they both rely on precipitation for seedling water supply.

**Table 2.** Soil parameters in the plantation sites. Data are means ± SD of n = 5–6 determinations (except for SWA: means ± SE, n = 2).

#### *3.2. Climate Data and Extreme Climate Events*

While the young seedlings were under controlled watering conditions during their first year of growth, the trees in SWA, OLY and STR were dependent on precipitation after planting. In RIE and SAN, watering supply was installed and used, if the weather conditions indicated drought situations. Monthly temperature and precipitation patterns at the common gardens (2017, 2018 and 2019) are shown in Table S2, a few of which are noteworthy. In January 2017, the potted seedlings were exposed to extreme frost (subzero temperature down to −10 ◦C for several days) before establishing the plantations in Greece. In July and August 2018 and 2019, extremely high temperatures in combination with prolonged drought were encountered by the seedlings in Germany. Especially in SWA temperatures were higher than RIE and even exceeded 40 ◦C for a few days.

#### *3.3. Seedling Survival and Growth*

Initially, 68% of seeds germinated on average and gave rise to seedlings. The survival rates in the third year after planting (2019) are given in Table 3. Seedlings grown in Germany and Italy exhibited relatively high survival rates, i.e., 59–99%. In Greece, environmental conditions were particularly challenging for the seedlings, especially for *Q. ilex* with few individuals surviving. In general, seedlings in STR had greater survival rates compared to OLY plantation. In all five plantations, *Q. pubescens* revealed the highest survival rates among the three oak species. Moreover, the drier (sandy soil) German plantation site (SWA), and the warmer and drier (lowland) Greek plantation site (OLY) exhibited higher seedling mortality compared to RIE and STR, respectively.

**Table 3.** Survival rates (%) of seedlings after three years in the five plantation sites (August 2019).

Figure 2 depicts the initial size and the height growth of the seedlings during the first three years in the common gardens. As a general pattern, seedlings in Germany were tallest except for the Italian *Q. ilex*, which showed the highest initial size in the Italian nursery (Gubbio, Italy). Statistical analyses after three years in the plantation sites revealed that seedling height (Table S3) was greatest for *Q. ilex* and Italian *Q. pubescens* in RIE and SAN, while *Q. frainetto* seedlings were taller mainly in the German plantations (RIE and SWA) and secondarily in SAN. In the Greek plantations, seedling height of all species was significantly lower showing practically no height growth.

**Figure 2.** Development of seedling height in the different plantations sites. Greek provenances: blue, Italian provenances: green, *Q. robur*: red. Box plots with 10–90% whiskers.

In general, root collar diameter followed the same trend with height. Height × root collar diameter correlation analysis in 2019 revealed moderately positive (Greek *Q. frainetto* and both *Q. pubescens*) to strong positive (Italian *Q. frainetto* and both *Q. ilex*) relationships between the two parameters (Figure 3).

**Figure 3.** Pearson's correlation analyses between seedling height and root collar diameter in 2019 for every species × provenances tested. The "r" value represents Pearson's correlation coefficient. In all species × provenances, *p*-value was <0.0001.

#### *3.4. Relative Chlorophyll Content*

Table 4 depicts the SPAD values of each species established in the common gardens in the nursery and during three years after planting. *Q. ilex* provenances in general showed the lowest SPAD values in all years in STR and high values in SAN and the German sites. SPAD values of both *Q. frainetto* provenances were generally greater in SWA compared to the rest of the plantation sites. Moreover, both Greek and Italian *Q. pubescens* revealed significantly greater SPAD values mainly in SWA and secondarily in SAN. In general, SPAD values in SWA and RIE showed a tendency to increase after establishment in the field, whereas seedlings at SAN revealed the opposite response.

**Table 4.** Relative chlorophyll content (SPAD) of each species × provenance from the nursery (2016) until after three years in the different plantations sites (2017–2019). Mean values (± SE) of the different plantation sites for each species/provenance within a line for each year followed by different letters are significantly different

#### *3.5. Quality Scoring of Seedlings as Predictors for Survival and Growth and Population Fitness*

To analyze, whether root collar diameter measured in the first year to predict secondyear survival of oak seedlings, the length data from 2016 of the perished (2017) and the live seedlings in 2017 were arranged into separate groups for all seedlings planted in Germany (favourable conditions) and in Greece (unfavourable conditions), respectively, and yielded statistically significant differences were tested, as suggested by Tsakaldimi et al. [14]. No significant differences were observed between survived (2017) and perished seedlings in SWA and STR. However, survived *Q. frainetto* from both provenances and *Q. pubescens* IT in RIE, and *Q. frainetto* IT in OLY developed significantly thicker stems (in 2016) compared to the respective perished seedlings (Supplemental Material Table S4). To further elucidate potential differences between the survivors and the perished trees at the sites, where stem diameter alone revealed no indicative value, a more complex parameter was also developed, taking into account both morphological and physiological fitness (i.e., *FI*, a function of tree length, leaf development, chlorophyll content and efficiency of the photosynthetic apparatus (Equation (1)). However, in no combination of provenance and country of plantation the 2016 *FI* values of the subpopulations (live in 2017 versus dead in 2017) yielded significant differences, meaning that the 2016 *FI* values could not be used as predictors for seedling survival (Supplemental Material Table S5).

In addition, the relationship between the *FI* values of individual trees in 2016 and 2019, respectively, was studied, for the German plantations to assess the predictive value of seedling measurements for future growth (Supplemental Material Table S6). In the German plantations, where 72–99% of the seedlings survived into 2019 in all provenances, the Spearman correlation coefficients between the *RPF* in 2016 and 2019 were weak (i.e., <0.21) and calculated *RPF* obtained in 2016 were considered of little predictive value for the future performance of individual seedlings. If data for seedling height in 2016 were compared to those of 2019, in all cases weak correlations were found both in the German and in the Greek plantations (Pearson's r < 0.3).

In Table 1, all populations per site have been evaluated for their relative vitality, using survival rates, morphometric and physiological data obtained in 2019. Data for STR had to be omitted due to road blocking by storm-felled trees in 2019 which prevented accessibility at night for *PIabs* measurements. However, the site was accessible at daytime and data for survival and height at STR are given in Table 3 and Figure 2, respectively. The best performing population in Germany was the Italian *Q. pubescens* followed by *Q. robur* and the Greek *Q. pubescens*, all in RIE, while *Q. frainetto* showed the lowest *RPF*. In Germany, RIE populations showed higher vitality scores compared to SWA. In Italy, *Q. pubescens* showed the best overall performance with *Q. ilex* following. In Greece, *Q. pubescens* showed relatively higher vitality scores compared to the other oak species.

#### **4. Discussion and Conclusions**

The present article provides an early evaluation of three oak species' establishment in four common garden studies, including nursery stage as well as field performance during the first three years after field plantation. The seedling production process followed the protocols of local nurseries. It was decided to not impose a standardized procedure on the nurseries, since each of them had its own year-long experience with the local oak provenances under the local climatic conditions. The same argument holds for the different plantation procedures. It was not the aim of the study to experiment with the local standardized growth techniques, since this would have introduced additional technical problems in procedure (e.g., shipping of plantation soil from one country to another, missing availability of frost-free greenhouses in Gubbio and Olympiada) and site management. The authors were aware that the different treatments at different sites would prevent a more comprehensive comparison of seedling performance with solely the different macroclimates (and soil conditions) as variables. However, the cooperation with the local nurseries and their local experiences with the growth of oaks reflect the situation in real-life transfer of genotypes from one country to the other.

The potential to include a German population of *Q. pubescens* was also examined. In German nurseries, no *Q. pubescens* seeds or seedlings are available from German populations, since these are mostly very small, often include hybrids with the local species, and, to our knowledge, are thus usually not included in certified seedling material lists. It was decided against self-collecting such seeds due to the genetic uncertainties as well as due to the fact that it was not possible to predict if there would be sufficient seeds available in 2016. Since the aim of the study was not translocating German material to the Mediterranean countries, the German control material was confined to the local species at the plantation sites (i.e., *Q. robur*).

#### *4.1. Seedling Survival and Establishment*

Seedling emergence was good for all collected seed populations, giving rise to sufficient plant material for an initial selection step, except for *Q. frainetto* of Italian origin at the Greek nursery. Transportation of the seeds in perforated plastic bags by airmail may have led to partial anoxia in the center of the bags, thus reducing germination rates, but in general, the outcome of seedlings was satisfactory. After their respective plantation in winter 2016/17, survival rates were assessed in summer 2019, and except for the Greek plantations, seedlings had well established themselves in the Common Gardens (Table 3). Initially it was intended to irrigate the seedlings in one of the Greek sites but this ultimately turned out to be impossible due to technical problems, i.e., the large distance in the finally available plantation site from available water. At RIE and SAN, facilities for watering of the seedlings had been established beforehand, so during severe soil drought in summer, these seedlings could be watered. Here, as well as in SWA, where no watering occurred, survival rates of all seedlings were satisfactory. In Greece, seedlings encountered a period of extreme frost in early January 2017 before establishing the plantations, and the evergreen *Q. ilex* was severely affected. Low seedling vitality has also been described for other broadleaves after freezing [23] and in summer-dry, calcareous Mediterranean habitats (e.g., Pausas et al. [12]: only ca. 15% survival of *Q. ilex* ssp. *ballota* after 3 years). During planting at the beginning of March, the seedlings still looked alive with most of the leaves still green, but turned out to be lethally damaged when the growth period started. For *Q. ilex* from its northernmost population in Italy, North of Lake Garda, it is known that frost tolerance (50% survival) of 1-year-old seedlings is around −10 (leaves) to −15 ◦C (sprout cambium) [24]. Such temperatures have, to our knowledge, never been recorded in the area of Olympiada, where the seedling material stems from, nor in Latium, Italy, close to the sea, where the Italian seeds have been collected. Since temperatures at soil level, where the seedlings overwintered outdoors, may also have been well below the −9.7 ◦C recorded at the Olympiada weather station (i.e., at 2 m height), we therefore conclude, that stems, roots and leaves of most *Q. ilex* seedlings were irreversibly damaged already before the planting. In contrast to the seedlings, the predominant mature *Q. ilex/Q. pubescens* forest in the surroundings of the plantation at OLY was not visibly affected by the frost event (unpublished field observations in May 2017), in accordance with the findings at Lake Garda, where adult trees revealed frost tolerance of −25 ◦C (stem cambium, 50% survival, [24]). The deciduous species, *Q. frainetto* and *Q. pubescens*, which occur naturally at higher elevations in Chalkidiki, northern Greece, obviously developed better frost tolerance of stem and cambium (and even roots). Although *Q. frainetto* seedling survival rates in the field are better when planted earlier and allowed to develop a better root system (i.e., Dec/Jan, [25]), in the current experiment we opted for a later planting date because of the severe frost in January 2017. However, in March/early April 2017 precipitation was low, especially at OLY (Table S2), and many seedlings died of drought due to weakly developed root systems before sufficient rainfalls occurred in mid-April. *Q. pubescens* showed better survival rates compared to *Q. frainetto* under identical conditions (Table 1) confirming similar findings in Northern Greece [26]. *Q. frainetto* has been found to be the least drought-tolerant oak among four Mediterranean species including *Q. pubescens* and *Q. ilex* [27], and this may explain its low survival rate at OLY and STR after the spring drought

immediately after plantation. At SAN and RIE, in contrast, seedlings were watered when necessary and therefore mortality was much lower. In SWA, although seedlings relied on precipitation, mortality was also low. Here, strong rainfalls occurred in mid-April 2017 before the flushing of the trees, which is usually observed one month later in Germany than in Greece.

Concerning the predictive value of morphological and physiological properties (chl content, *PIabs*) of 1-year old seedlings for their future survival and development in the field, our results did not meet our expectations. Bayala et al. [28], on five tropical tree species, and Tsakaldimi et al. [14] on five Mediterranean tree species, including *Q. ilex* and *Q. coccifera*, pointed to stem diameter as a good indicator for second-year survival after outplanting. We examined this proposed indicator and the results were variable among the common gardens. In RIE, where seedlings showed the greatest survival, two-year old seedlings of *Q. frainetto* (both provenances) and *Q. pubescens* IT were successfully grouped to survived or perished depending on their root collar diameter, thus confirming the results of Tsakaldimi et al. [14] (Table S4). In general, root collar diameter of *Q. frainetto* showed a greater potential to predict second-year survival but differences were not always significant. While the morphological features height and root collar diameter were well correlated in older seedlings, confirming data of Pausas et al. [12] on *Q. ilex* ssp. *ballota* (Figure 2), the *FI* values in the subgroups of established (survivors into 2017) and perished seedlings (dead in 2017) did not differ significantly (Table S5). Secondly, the Spearman r values between the *FI* values in 2016 (using seedling height, number of leaves, chl content and *PIabs*) and 2019 (using seedling height and *PIabs*) were considered too low to be of practical value (cf. Supplementary Material Table S6: r was either not significant or even negative in all but one cases). Quite similarly, if measurements are confined to height the predictive value is also weak, with Pearson coefficients for the respective 2016 and 2019 values being negative or below 0.17 in all cases, thus height alone is neither considered suitable for preselection for plantation purposes.

In conclusion, of the parameters tested, only in good growing conditions (i.e., RIE), root collar diameter of 2nd year seedlings had a moderate predictive value for seedling survival in three out of seven species/provenances, while other non-destructive parameters (included in the *FI* values) appear to be of no practical value.

#### *4.2. Cumulative Growth Conditions at the Common Garden Sites as Reflected in Seedling Performance*

Height and root collar growth during the first three years after planting (measured in summers 2017, 2018 and 2019) can be used as first predictors for the suitability of the respective oak provenances for establishment at the different sites. Generally, both parameters were relatively high at SAN for most species/provenances over the three years, reflecting the good nutrient conditions (cations) in the basic soil and the sufficient water supply. In the German sites, *Q. pubescens* and *Q. robur* revealed no height growth during the first year, pointing to the costs of acclimation (i.e., root development) after planting. However, in the following years the Italian *Q. pubescens* in RIE developed strongly and reached comparable height with the SAN homologue most probably due to a well-developed root system. In the first year at the German sites, *Q. frainetto* and *Q. ilex* seedlings at SWA in general showed a stronger height growth than at RIE, possibly because of the lower light availability at SWA (thinned pine umbrella plantation with 50 trees/ha [5]) as compared to the open setup at RIE. However, after three years both *Q. ilex* provenances as well as *Q. robur* were significantly promoted by the relatively favorable conditions in RIE (i.e., irrigation) compared to SWA. After three years at the Greek sites, especially *Q. pubescens* revealed better height growth at STR than at OLY, again possibly due to partial shading (cf. Supplementary Material Table S3).

The quite consistent higher chl contents of *Q. frainetto*, *Q. pubescens* and *Q. robur* at SWA vs. RIE can be explained by partial shading by pine trees at SWA, whereas seedlings at RIE grew in full sunlight. At SAN, seedlings were also partially shaded by weeds during the period of measurements which caused early senescence of *Q. frainetto* and subsequently low chl content values. Leaf senescence depends on the carbon fixation efficiency which is modulated by plant photoreceptors such as phytochrome A [29]. STR and OLY seedlings generally showed similar (low) chl content over the three-year period as compared to the other sites. At STR, the high irradiation, at least compared to SWA and RIE, may have caused the lower chl contents. Additionally, OLY seedlings were exposed to full sunlight, combined with high temperature and a long drought period. This combination accelerated photoinhibition [2]. Frequent photoinhibition, in turn, may cause photodamage to pigments. Similarly, Cotrozzi et al. [30] reported leaf injuries (i.e., yellowing) induced by drought stress in seedlings of *Q. ilex*, *Q. pubescens* and *Q. cerris*. In the German plantations, lower SPAD values were observed in 2017, 2018 and 2019 (outdoor grown seedlings) vs. data from 2016 (greenhouse-grown seedlings).

The advantage which the seedlings in Germany had by initial greenhouse growth over the Italian seedlings at the time of planting (compare seedling heights in 2016 in Figure 2) was (mostly) lost in the course of outdoor growth in 2019. This is partly due to prolonged extreme heat and drought periods in 2018 and 2019 (Table S2) which showed similar intensity with the heat and drought in Europe in 2003 [31]. The most severe growth conditions existed at the Greek sites, nearly extinguishing the *Q. ilex* populations by the very unusual extreme frost in January 2017 and affecting growth in the deciduous species e.g., by spring drought.

#### *4.3. Preliminary Conclusions for (Re)forestation Strategies*

In summary, we can conclude that:

Presently, differences between Italian and Greek populations are small, and it will be necessary to monitor seedling development for more years to decide, whether the Italian populations show better performance in the long run. At OLY and STR, growth conditions are extremely hard for young oaks. To regenerate *Q. ilex* forest (especially at the lower, dryer site OLY), several attempts in successive years may be necessary to avoid frost damage, if extreme situations like in January 2017 are to occur in future years again. It is further recommended to irrigate new plantations of rooted seedlings in this region during the first one or two years [33]. In Germany, *Q. pubescens* and *Q. robur* showed the highest relative vitality scores (Table 1) in both sites, despite the fact that the two species showed minimal height growth during the first year after planting. *Q. frainetto* showed inferior quality compared to the other oaks. Consequently, *Q. pubescens* may be considered a potential future forest tree on poor and dry soils for the expected future summer heat and drought conditions in Central Europe, as were observed in Germany in 2018 and 2019. Whether Italian or Greek provenances are better suited for planting of Mediterranean species in southwest Germany should be addressed by future research in the common gardens.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/f12060678/s1, TableS1: Oak species, collection sites and nursery sites. Table S2: Average minimum and maximum monthly temperatures and total monthly precipitation recorded in the common gardens in the course of 2017, 2018 and 2019. Data in RIE, SWA were recorded on site, data in SAN at the climate station in Casteldilago, and data for OLY and STR close to the sites at the climate stations of Hellas Gold. Table S3: Statistical analysis of height (in cm) in the common gardens, measured in summer 2019. Values followed by different letters within a row are significantly different at *P* < 5%, as calculated by Kruskal-Wallis one-way-ANOVA tests with Dunn´s post-hoc test. Means of *Q. robur* were analyzed by t-test at *P* < 5%. Table S4: Statistical analysis of root collar diameter (mm) in the common gardens, measured in summer 2016, to examine the parameters potential as indicator of second-year survival after outplanting. Values followed by different letters within a row are significantly different at *P* < 5%, as calculated by t-tests. Table S5: Mean values of fitness indicators (*FI*) of the plants in 2016, separated into those who established live in 2017 and those who perished in the first year after plantation (dead in 2017). Table S6: Spearman correlation analyses between Relative Population Fitness (*RPF*) ranking positions in 2019 vs. 2016 at the German sites. *RPF* were calculated as described in Materials and methods. Pearson´s r value was also calculated for the height data in 2019 vs. 2016 at the German and Greek sites. Significance values for the correlations are given at *P* < 5% (\*), 1% (\*\*), 0,1% (\*\*\*). Figure S1: Plantation schemes at RIE, SWA, and SAN (1A), OLY (1B), and STR (1C). Within each plot, individual trees were planted according to 1D; in the case of *Q. frainetto* from Italy at OLY and STR, trees number 02, 09, 13 and 20 were omitted.

**Author Contributions:** Conceptualization, W.B.; methodology, W.B., K.R. and F.B. (Filippo Bussotti); investigation, F.B. (Filippos Bantis), J.G. and E.F.; data curation, W.B., E.F. and F.B. (Filippos Bantis); writing—original draft preparation, F.B. (Filippos Bantis) and W.B.; writing—review and editing, F.B. (Filippos Bantis), J.G., E.F., K.R., F.B. (Filippo Bussotti) and W.B.; supervision, project administration and funding acquisition, W.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financed by grant No. 01DS15014 by the German Federal Ministry of Education and Research (BMBF) to W.B.

**Data Availability Statement:** The data presented in this study are available within the article and its supplementary materials.

**Acknowledgments:** Special thanks are due to the City of Frankfurt/Main (Germany), to Greek Nurseries SA (Olympiada, Greece), Hellas Gold SA (Athens, Greece) and to Agenzia Forestale Regionale (Perugia, Umbria, Italy) for providing the respective plantation sites and help with establishment and maintenance of the plantations as well as for providing climate data (OLY and STR). Thanks are due to Martina Pollastrini (Univ. Firenze), Vera Holland, Sonja Ströll, Nathalie Reininger and Lisa Schäfer (Univ. Frankfurt) for help with the data acquisition, and also to Mariangela Fotelli for commenting on the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Climate Change Effects in a Mediterranean Forest Following 21 Consecutive Years of Experimental Drought**

**Romà Ogaya 1,2,\* and Josep Peñuelas 1,2**

**Abstract:** *Research Highlights:* A small, long-term decrease in the water availability in a Mediterranean holm oak forest elicited strong effects on tree stem growth, mortality, and species composition, which led to changes in the ecosystem function and service provision. *Background and Objectives:* Many forest ecosystems are increasingly challenged by stress conditions under climate change. These new environmental constraints may drive changes in species distribution and ecosystem function. *Materials and Methods:* An evergreen Mediterranean holm oak (*Quercus ilex* L.) forest was subjected to 21 consecutive years of experimental drought (performing 30% of rainfall exclusion resulted in a 15% decrease in soil moisture). The effects of the annual climatic conditions and the experimental drought on a tree and shrub basal area increment were studied, with a focus on the two most dominant species (*Q. ilex* and the tall shrub *Phillyrea latifolia* L.). *Results:* Stem growth decreased and tree mortality increased under the experimental drought conditions and in hot and dry years. These effects differed between the two dominant species: the basal area of *Q. ilex* (the current, supradominant species) was dependent on water availability and climatic conditions, whereas *P. latifolia* was more tolerant to drought and experienced increased growth rates in plots where *Q. ilex* decay rates were high. *Conclusions:* Our findings reveal that small changes in water availability drive changes in species growth, composition, and distribution, as demonstrated by the continuous and ongoing replacement of the current supradominant *Q. ilex* by the subdominant *P. latifolia*, which is better adapted to tolerate hot and dry environments. The consequences of these ecological transformations for ecosystem function and service provision to human society are discussed.

**Keywords:** basal area increment; climate change; forest dieback; Mediterranean forest; stem growth; water availability

#### **1. Introduction**

Climate change is increasing the levels of environmental stress conditions for many forest ecosystems around the world due to higher rates of evapotranspiration and increases in the intensity and duration of extreme climate events [1,2]. Some water-limited forest ecosystems, such as Mediterranean forests, may be particularly vulnerable to the expected decrease in water availability and may experience associated decreases in stem growth [3,4], greater levels of canopy defoliation [5], and increases in tree mortality [6,7]. While rising concentrations of atmospheric carbon dioxide (CO2) are associated with higher growth rates and productivity in forests [8], it is possible that decreases in stem growth may occur in water-limited ecosystems, if the effects of water scarcity exceed those of fertilization from the increasing levels of CO2.

*Quercus ilex* L. is an evergreen tree species and it is the most abundant and representative of the Mediterranean forest [9]. *Phillyrea latifolia* L. is a tall shrub evergreen species, and it is also very common in the Mediterranean forest. Both species are well adapted to resist a drought period (usually the summer season) due to their stomatal closure around midday to avoid excessive water loss by transpiration [9]. Despite these two species having

**Citation:** Ogaya, R.; Peñuelas, J. Climate Change Effects in a Mediterranean Forest Following 21 Consecutive Years of Experimental Drought. *Forests* **2021**, *12*, 306. https://doi.org/10.3390/f12030306

Academic Editor: Mariangela Fotelli

Received: 14 January 2021 Accepted: 1 March 2021 Published: 6 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

experienced similar distributions and adaptations for coping with drought conditions, several studies have shown that *Q. ilex* is more sensitive to dry and hot conditions than other co-occurring tall shrub species, such as *P. latifolia* [3,4,10,11]. For example: within a shared forest, *Q. ilex* was found to suffer from high rates of stem mortality during periods of extreme drought, which had little impact on *P. latifolia* [10,12]; while wood carbon allocation decreased in *Q. ilex* during dry years and remained stable in *P. latifolia* [11]; and fruit production correlated with annual rainfall in *Q. ilex*, but not in *P. latifolia* [11].

Mediterranean vegetation is characterized by the physiological adaptations that allow its tolerance to, and recovery from, episodes of dry environmental conditions; however, if extreme drought events exceed the tipping points of climate tolerance, some sensitive species may fail to survive, whereby species diversity, composition, and associated ecosystem function become permanently altered, leading to perturbations in the ecosystem service provision [13]. *Q. ilex* is the current dominant species of the Mediterranean forest, and is widely distributed across the whole Mediterranean Basin in South Europe and North Africa, and from Portugal to Syria. Under hotter and drier environmental conditions, a decrease in the abundance of *Q. ilex* in Mediterranean forests may be expected, along with a reduction in forest capacity to mitigate the effects of climate change through the sequestration of atmospheric CO2; there may also be a progressive replacement of *Q. ilex* by other tall shrub species, such as *P. latifolia*, that transforms the structure and ecosystem function of the forest.

In water-limited ecosystems, the principal effects of hydrological processes on vegetation are mediated through soil water content [14]; thus, experimental methods that affect soil moisture effectively simulate the likely impacts of decreases in water availability on plants. Furthermore, long-term climate manipulation experiments in natural ecosystems allow the prediction of the effects of environmental change on ecosystem function and development [15,16] and may highlight the dampening effect on the size of treatments that has been widely reported in global-change experiments [17].

Here, our main objective was to determine the effects of climate change-mediated drought on the development of the Mediterranean forest, paying special attention to the development of the different species and their possible divergent evolution. This study was conducted as a long-term climate manipulation experiment that simulates a 15% decrease in soil moisture, as projected by general circulation models for the Mediterranean region [18].

#### **2. Materials and Methods**

#### *2.1. Study Site*

The study was conducted in a Mediterranean holm oak forest located on a south facing slope (25%) of the Prades Mountains, Catalonia, in north-eastern Spain (41◦20 38" N, 1◦2 0" E, 950 m a.s.l.), in which the vegetation comprises dense multi-stem, coppiced woody species (16.616 stems ha−1), created for charcoal production; no coppicing activity has taken place in the forest during the preceding 70 years. The dominant species are *Q. ilex* and *P. latifolia* (both species represent more than 80% of the total forest basal area). There are also other evergreen shrub species such as *Arbutus unedo* L., *Erica arborea* L., and *Juniperus oxycedrus* L.; deciduous tree species include the occasional *Sorbus torminalis* L. and *Acer monspessulanum* L. The soils comprise a Dystric Cambisol, and the range in soil depth (35 to 100 cm) results in contrasting levels of water and nutrient availability, reflected by two forest-type structures [19]: a tall canopy forest (8–10 m high) dominated by *Q. ilex*, and a low canopy forest (4–6 m high), with a greater abundance of *P. latifolia* than in tall canopy (Table A1).

#### *2.2. Experimental Design*

Eight 15 × 10 m plots were delimited at the same altitude along the slope (four in the tall canopy stand and another four in the low canopy stand) in 1999. Four of these plots (two each in the tall and low canopy stands) were randomly selected as drought plots, and received an experimental partial rainfall exclusion; the four remaining plots were left untreated as a control. In the drought plots, six plastic strips, comprising four rigid PVC tubes connected with transparent PVC, were installed under the canopy along the longitudinal slope, at about 1 m above the soil surface and extending to 2 m outside the plot margin to reduce edge-effects; the strips covered 30% of the plot surface area for the interception of rain throughfall and its removal downslope. Soil moisture was measured monthly using a time domain reflectometer (Delta-T Devices, Cambridge, UK) placed in the upper 15 cm of the soil profile of four randomly selected loci per plot, along with the leaf water potential (measured seasonally, in two leaves from two different trees per species, collected in the upper layer of the canopy and fully exposed to solar radiation) at the same loci using a Scholander-type chamber pressure (PMS Instrument Company, Albany, OR, USA). An automatic meteorological station installed at the study site monitored temperature and rainfall during the 1999–2019 study period.

#### *2.3. Basal Area Determination*

In 1999, all of the living tree and shrub stems in the plots were marked at 50 cm above ground level using permanent paint, and their circumference (diameter > 2 cm) was measured at the painted mark annually in winter. The basal area of each stem was determined from the stem circumference measurements, assuming each stem section to be a circle; the annual basal area increment was calculated from the difference in the basal area measured in two consecutive years. For each plot, species, and year, we determined the basal area increment (BAI), basal area lost by stem mortality (BAL), and total basal area increase (TBAI) which was determined as the basal area increment minus basal area lost by stem mortality.

#### *2.4. Data Analysis*

Seasonal differences in soil moisture and leaf water potential were tested using repeated measures analysis of variance (ANOVA). January, February and March were combined in "winter", April, May and June were combined in "spring", July, August and September were combined in "summer", and October, November and December were combined in "autumn". The relationships between the BAI, BAL, and TBAI per species and canopy stand (dependent variables) and the annual mean temperature, rainfall, and SPEI climatic indices [3,4,20] (independent variables) were tested using linear and polynomial regression analyses. We used general linear modeling to test the determinant factors of the annual BAI, BAL, and TBAI of *Q. ilex* and *P. latifolia.* The species, treatment and canopy stand were considered independent factors, and the annual mean temperature and rainfall were covariates. Finally, we conducted the same models without the species factor, and including all tree and shrub species present in the plots named "overall species". All analyses were performed using Statistica 10.0 (StatSoft Inc., 2010, Tulsa, OK, USA).

#### **3. Results**

#### *3.1. Soil Moisture and Leaf Water Potential*

The meteorological conditions of the study site were typically Mediterranean: summers were hot and dry, winters were relatively cold, and rainfall events tended to be concentrated in spring and autumn (Figure 1). Accordingly, both soil moisture and leaf water potential levels showed a clear seasonality; their values decreased in summer, recovered in autumn, and were greatest in spring and winter when evapotranspiration rates were lower. Soil moisture levels were 15% lower in the drought plots than in the control plots (*p* < 0.05) (Figure 2). The decrease in leaf water potential values in summer was stronger in *P. latifolia* than in *Q. ilex* (*p* < 0.05), whereas in autumn, *Q. ilex* maintained a higher water potential than *P. latifolia*; in winter and spring, the levels were similar in both species (Figure 2). The drought treatment did not exert any significant additional effect on the leaf water potential levels.

**Figure 1.** Average meteorological conditions at the study site during the experimental period (1999–2019). Red dots indicate mean monthly air temperature and grey columns indicate mean monthly rainfall. Vertical bars in the grey columns and red dots indicate the standard error of the mean (*n* = 21 years).

**Figure 2.** Average seasonal leaf water potential (upper panel) and soil moisture (lower panel) during the study period (1999–2019). Vertical bars indicate the standard error of the mean (*n* = 21 years).

#### *3.2. Climate Effects on Stem Growth*

BAI was positively related to the annual rainfall in the two species separately (*Q. ilex*: *p* < 0.05; *P. latifolia*: *p* < 0.01) and for all of the tree and shrub species (*p* <0.01) (Figure 3); there was a trend towards a positive relationship between the BAL and the mean annual temperature in *Q. ilex* (*p* = 0.06) and all of the tree and shrub species (*p* = 0.06), but not in *P. latifolia* (Figure 4). As a consequence of lower stem growth in dry years and greater stem mortality in hot years, the TBAI was positively related to the annual rainfall in *Q. ilex* (*p* = 0.05), *P. latifolia* (*p* < 0.01), and all of the tree and shrub species (*p* < 0.01), but these relationships were not linear; they tended to saturate at high precipitation values. On the other hand, there was a negative relationship between the TBAI and the mean annual temperature in *Q. ilex* (*p* = 0.05) and all of the tree and shrub species (*p* < 0.05); in *P. latifolia*, there was a trend towards a negative relationship between the TBAI and the mean annual temperature (*p* = 0.06). Also, these relationships were not linear; they tended to increase at very low temperature values, but strongly decrease at high temperatures. Thes effects of climate on stem growth differed between the two canopy stands, where the mean annual temperature tended to be more determinant for the TBAI in the tall canopy stand for *Q. ilex* (*p* < 0.05) and all of the tree and shrub species together (*p* < 0.05), whereas the annual rainfall was more determinant for the TBAI in the low canopy stand for *Q. ilex* (*p* = 0.05), *P. latifolia* (*p* < 0.01) and all of the tree and shrub species together (*p* < 0.05) (Figure 5). However, we did not obtain significative results between the BAI measurements and SPEI climatic indices.

**Figure 3.** Linear regression analysis of annual basal area increment (BAI) and annual rainfall in *Q. ilex*, *P. latifolia*, and all tree and shrub species. *n* = 21 years.

**Figure 4.** Linear regression analysis of annual basal area lost by stem mortality (BAL) and mean annual temperature in *Q. ilex*, *P. latifolia*, and all tree and shrub species. *n* = 21 years.

**Figure 5.** Polynomial regression analysis of total annual basal area increase (TBAI) and mean annual temperature in the tall canopy stand and total annual rainfall in the low canopy stand, for *Q. ilex*, *P. latifolia*, and all of the tree and shrub species. *n* = 21 years.

#### *3.3. Drought Effects on Stem Growth*

Between 1999 and 2019, regardless of treatment, the BAI in *Q. ilex* was greater in the tall canopy stand than in the low canopy stand, while in *P. latifolia* it was lower in the tall canopy stand than in the low canopy stand (species × canopy interaction: *F* = 114.1, *p* < 0.01), reflecting the recorded basal area of the species in 1999 and 2019 (Table A1); there was no difference in the BAI between the tall and low canopy stands when all of the tree and shrub species were considered together (Figure 6). While there were no effects of the drought on the BAI in *Q. ilex*, *P. latifolia*, or all of the tree and shrub species, the BAL was greater in the drought plots in *Q. ilex* and *P. latifolia* (*F* = 8.4, *p* < 0.01), and all of the tree and shrub species (*F* = 8.1, *p* < 0.01), with no overall BAL canopy type differences in *Q. ilex* or all of the tree and shrub species. In P. latifolia, the BAL was greater in the low canopy stand (Figure 6); however, the BAL was always higher in *Q. ilex* than in *P. latifolia* (*F* = 33.9, *p* < 0.01). Across the study period, the TBAI reflected clear differences between the two species: *Q. ilex* experienced a higher TBAI in the tall canopy stand, whereas *P. latifolia* experienced higher TBAI in low canopy stand (species x canopy interaction: *F* = 29.1, *p* < 0.01); there was also a clear decrease in the TBAI levels induced by drought (*F* = 5.7, *p* = 0.02), especially in *Q. ilex* (species x treatment interaction: *F* = 5.1, *p* = 0.02). We found that the maximum TBAI in *Q. ilex* was recorded in non-restricted rainfall (control) plots located in the tall canopy stand, where there was a slight positive association between the TBAI and drought treatment in the tall canopy stand and a slight negative association between the TBAI and non-restricted rainfall in the low canopy stand; the lowest TBAI was recorded in drought plots of the low canopy stand (Figure 7). When all of the tree and

shrub species were considered together, there was a decrease in the TBAI levels in drought plots (*F* = 12.9, *p* < 0.01) and there was no difference in the canopy stand type (Figure 7).

**Figure 6.** Drought treatment and canopy stand type effects on basal area increment (BAI) and basal area lost by stem mortality (BAL) in *Q. ilex*, *P. latifolia*, and all of the tree and shrub species during the study period (1999–2019). Vertical bars indicate the standard error of the mean (*n* = 21 years). Different letters indicate drought treatment differences between and within the canopy stand types (*p* < 0.05).

**Figure 7.** Drought treatment and canopy stand type effects on total basal area increment (TBAI) in *Q. ilex*, *P. latifolia*, and all of the tree and shrub species during the study period (1999–2019). Vertical bars indicate the standard error of the mean (*n* = 21 years). Different letters indicate drought treatment differences between and within the canopy stand types (*p* < 0.05).

Overall, the TBAI of the most representative species of the study forest followed three main response patterns to changes in soil moisture: species with a greater TBAI under non-drought conditions, such as *Q. ilex*, A. unedo, and J. oxycedrus; species with a similar TBAI, regardless of drought conditions, such as *P. latifolia*, E. arborea, and Rhamnus alaternus; and species with decreases in the TBAI, regardless of drought conditions, albeit with larger decreases under drought conditions, such as S. torminalis (Figure 8).

**Figure 8.** Total basal area increment (TBAI) of the most representative forest tree and shrub species under full (control) and partial (drought) rainfall conditions during the study period (1999–2019). Vertical and horizontal bars indicate the standard error of the mean (*n* = 21 years). Species with a higher TBAI in control plots than in drought plots are depicted with blue circles, species with a similar TBAI in both control and drought plots are depicted with green circles, and the species with a negative TBAI in both the control and drought plots are depicted with a red circle. Q.i. means *Quercus ilex*, A.u. means *Arbutus unedo*, J.o. means *Juniperus oxycedrus*, P.l. means *Phillyrea latifolia*, E.a. means *Erica arborea*, R.a. means *Rhamnus alaternus*, and S.t. means *Sorbus torminalis*.

#### **4. Discussion**

As observed previously at this experimental study site [3,4,12,20], we found that the tree and shrub stem growth was constrained by water availability and tree mortality increased in dry and hot years. These effects were particularly marked in the dominant species, *Q. ilex*, whereas in other species, such as *P. latifolia*, that are more tolerant of lower levels of water availability, there were no effects of experimental drought and there was positive stem growth, especially where the *Q. ilex* abundance was declining in the droughted low canopy stand. Many studies report an increase in basal area increment with tree size [21], and a different climate sensitivity depending on tree size [22], but the *Q. ilex* stem size tended to be greater than that of *P. latifolia*. Despite this, the BAI was lower in *Q. ilex* than in *P. latifolia*, and in the low canopy stand where the tree size of both of the species was similar, climate sensitivity was higher in *Q. ilex* than in *P. latifolia*. These two species share a strategy to avoid excessive water loss due to transpiration; they close their stomata when they are subjected to a high vapor pressure deficit and low water availability. However, previous findings from this long-term experiment revealed only a small decrease in net photosynthetic CO2 uptake induced by the drought treatment in the two species; this small reduction in the photosynthetic rate is insufficient to explain the considerable decrease in the BAI observed in *Q. ilex* from the drought plots [23]. Changes in the leaf water potential and other physiological parameters may explain the lower BAI in droughted *Q. ilex* in this study; for example, lower production of new leaves in *Q. ilex* canopies has been recorded under drought conditions than under control conditions, with no such treatment differences in the *P. latifolia* leaf demography [24].

We found that the experimental drought increased stem mortality in all of the tree and shrub species, but this effect was stronger in *Q. ilex* than in *P. latifolia*. Increases in both the air temperature and dry conditions are considered to be the most important drivers of tree mortality due to hydraulic failure and carbon starvation [6,7,25]. In our study, mortality increased markedly with the mean annual temperature in *Q. ilex*, whereas *P. latifolia* was more tolerant of higher temperatures, as indicated by the similar mortality rates across the study years [12]. Despite their similar physiological strategies to increase tolerance of dry conditions, detailed xylem architecture observations revealed wider xylem vessels in *Q. ilex* than in *P. latifolia* [26] that allow the transport of greater amounts of sap under well-watered conditions, but are more susceptible to embolism under drought conditions.

Shifts in species distributions induced by climate change are expected to occur globally in several biomes [27,28]. Across the duration of this long-term experiment, the total basal area of *Q. ilex* decreased exponentially with an increase in air temperature and decrease in rainfall, while a previous study showed the contrasting presence of abundant *P. latifolia* seedlings and young individuals and scarce *Q. ilex* seedlings in this forest [29], indicating that the replacement of moribund *Q. ilex* may only be possible from resprouting shoots. These results indicate that the warming air temperatures and decreasing water availability that are predicted for the Mediterranean Basin [1] will lead to a receding *Q. ilex* distribution and progressive substitution by shrub species better adapted to drier conditions, particularly under extreme arid conditions, such as those in the low canopy stand. By the first time, we detected in our experimental site a clear substitution of *Q. ilex* by other shrub species, as it was clearly observed in the low canopy stand.

Across the study period, the total basal area of the drought plots did not increase, and this reduction in plant productivity was accompanied by other important changes in the ecosystem function of the forest. For example: leaf material from the drought treatment contained greater concentrations of nitrogen (N)−15 than that from the control, due to greater ecosystem N losses under the experimental drought conditions [30]; bacteria diversity has been shown to be particularly sensitive to water availability and fungi diversity was more tolerant to seasonal and experimental changes in soil moisture [31]; and drought conditions led to decreases in soil phosphatase activity [32] and altered trophic interactions, such as increases in leaf sugar and phenol concentrations to resist drier conditions, that boosted levels of insect herbivory [19].

In order to increase ecosystem contributions to climate change mitigation, considerable efforts have been made to preserve species composition and the associated ecosystem function and service provision [33]. For Mediterranean forests, such as that studied here, with high levels of stem density and competition and an associated low BAI [34], the reduction of stem density by selective thinning may represent an appropriate management option to improve the BAI [12,35].

#### **5. Conclusions**

The results from this long-term study demonstrated a progressive replacement of the dominant species, *Q. ilex*, by another shrub species, *P. latifolia*, that was better adapted to the hotter and drier climatic conditions that are predicted to impact this Mediterranean region over the coming decades. The slight decrease in water availability, as simulated by the drought treatment, that reduced soil moisture by 15% exacerbated this transition in relative species abundance and contributed to a decrease in stem growth and an increase in tree mortality in *Q. ilex*. The experimental drought also elicited additional, ecosystem level impacts, including losses in N, changes in the metabolic activity of soil organisms, and the capacity to resist attack by plant pathogens. Importantly, some ecosystem services of this forest type, such as the amelioration of climate change impacts as a carbon sink, are expected to decline further with the ongoing and continued changes in the climate.

**Author Contributions:** Conceptualization, experimental design, statistical analyses, and discussion of the results: R.O. and J.P.; writing: R.O.; and, supervision: J.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was financially supported by the Spanish government project PID2019- 110521GB-I00 ELEMENTAL SHIFT, the European Research Council Synergy grant ERC-2013-SyG-2013-610028 IMBALANCE-P, and the Catalan government grant SGR-2017-1005.

**Acknowledgments:** This study was undertaken with the approval of Ester Trullols, director of the PNIN de Poblet Natural Park.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** Mean (±SEM) basal area (m2 ha−1) of *Q. ilex*, *P. latifolia*, and all tree and shrub species in the plot types at the start and the end of the study period (1999–2919) (*n* = 2 plots).

#### **References**

### *Article* **Sap Flow in Aleppo Pine in Greece in Relation to Sapwood Radial Gradient, Temporal and Climatic Variability**

**Evangelia Korakaki 1 and Mariangela N. Fotelli 2,\***

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +30-2310-461172 (ext. 238)

**Abstract:** Research Highlights: The radial gradient of sap flux density (Js) and the effects of climatic factors on sap flow of Aleppo pine were assessed at different time scales in an eastern Mediterranean ecosystem to improve our understanding of the species water balance. Background and Objectives: Aleppo pine's sap flow radial profile and responses to environmental parameters in the eastern Mediterranean were, to our best knowledge, originating to date from more arid planted forests. Information from natural forests in this region was lacking. Our objectives were to (a) determine the species' radial variability in Js on a diurnal and seasonal basis and under different climatic conditions, (b) scale up to tree sap flow taking into account the radial profile of Js and (c) determine the responses of Aleppo pine's sap flow over the year to climatic variability. Materials and Methods: Js was monitored in Aleppo pine in a natural forest in northern Greece with Granier's method using sensors at three sapwood depths (21, 51, and 81 mm) during two periods differing in climatic conditions, particularly in soil water availability. Results: Js was the highest at 21 mm sapwood depth, and it declined with increasing depth. A steeper gradient of Js in deep sapwood was observed under drier conditions. The same patterns of radial variability in Js were maintained throughout the year, but the contribution of inner sapwood to sap flow was the highest in autumn when the lower seasonal Js was recorded in both study periods. Not taking into account the radial gradient of Js in the studied Aleppo pine would result in a c. 20.2–27.7 % overestimation of total sap flow on a sapwood basis (Qs), irrespective of climatic conditions. On a diurnal and seasonal basis, VPD was the strongest determinant of sap flux density, while at a larger temporal scale, the effect of soil water content was evident. At SWC > 20% sap flow responded positively to increasing solar radiation and VPD, indicating the decisive role of water availability in the studied region. Moreover, in drier days with VPD > 0.7 KPa, SWC controlled the variation of sap flow. Conclusions: There is a considerable radial variability in Js of the studied Aleppo pine and a considerable fluctuation of sap flow with environmental dynamics that should be taken into account when addressing the species water balance.

**Keywords:** sap flux; radial profile; sapwood depth; Aleppo pine; diurnal variation; seasonal variation; climate

#### **1. Introduction**

Sap flux density and sap flow measurements are widely used in ecophysiological studies to investigate aspects of plant–water interactions such as canopy conductance [1–5], whole-plant hydraulic conductance [3,6], and leaf stomatal conductance [3,6–9]. To estimate these variables from sap flow monitoring, sap flux density (Js) measurements need to be upscaled to the tree level. In several studies [10–17], sap flow is estimated from sap flux density measurements with sensors installed in the outer part of the xylem, assuming uniform sap flux density across the sapwood. Berdanier et al. [18] reported that out of 122 reviewed articles (published between 2013 and 2015) that scaled

**Citation:** Korakaki, E.; Fotelli, M.N. Sap Flow in Aleppo Pine in Greece in Relation to Sapwood Radial Gradient, Temporal and Climatic Variability. *Forests* **2021**, *12*, 2. https://dx.doi.org/10.3390/ f12010002

Received: 23 November 2020 Accepted: 18 December 2020 Published: 22 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

up sap flux observations to whole-tree or stand levels, 58% assumed homogenous flow throughout the xylem.

However, neglecting radial variability may lead to systematic errors in the total sap flow [19–25]. Cermak and Nadezhdina [ ˇ 20] showed that upscaling sap flux density, measured in the outer half of the sapwood of different species, led to overestimation of total tree sap flow (by 10–140%), whereas when measured in deep inner sapwood layers, total tree sap flow was underestimated (by 40–80%). Similarly, Delzon et al. [22] showed that ignoring the variability in radial profile of sap flow in maritime pine trees led to an overestimation of transpiration and stand water use that increased with increasing stem diameter from 5% for the young stands to 47% for the oldest stands. To overcome these uncertainties, the radial variability in sap flux density in different trees has been determined [19,26–30] or correction factors accounting for the spatial variation of sap flow with increasing sapwood radius have been developed [23,31]. Therefore, taking into consideration radial variability is highly important when upscaling sap flux density measurements to accurately estimate tree and stand level water balance, particularly for tracheid-bearing conifers that tend to have deep functional sapwood [18,32].

Variation in the sap flux radial profile in conifer stems, as well as variation in sap flow on the tree level, may be associated with time-dependent changes in environmental conditions. A number of studies have reported that changes in radial patterns were related to environmental variables, such as vapor pressure deficit and soil water availability. Phillips et al. [33] indicated that nighttime recharge of stem water storage increased with increasing soil moisture depletion. Cermak and Nadezhdina [ ˇ 20] showed that, under long-term drought, deeper layers of the sapwood contributed more to water transport. Similarly, Ford et al. [27] found that the radial profile of sap flux density was related to daily changes in evaporative demand, and that during periods of high transpiration, the inner xylem contributed more to total tree sap flow. Hernandez-Santana et al. [34] observed that Js increased at all explored xylem depths with increasing VPD from 8:00 to 11:00 GMT, but with higher VPD (till 15:00 GMT), Js continued increasing only in the deepest xylem, whereas the outer one decreased. In the study of Sanchez-Costa et al. [31], *P. halepensis* showed higher sensitivity of Js to VPD, than the rest diffuse-porous species examined and a higher soil moisture threshold for Js reductions, consistent with the species strict stomatal control under drought. Finally, Cohen et al. [35] supported that the depth of conducting sapwood in *P. halepensis* trees was larger in trees growing under a semi-arid climate than under a sub-humid climate. Thus, assessment of the climatic controls in different time scales on sap flow patterns is important in our effort to understand and predict Aleppo pine's responses to climate variability.

Aleppo pine is widely distributed in coastal and low elevation areas of the Mediterranean basin [36,37], being predominately threatened by the prevailing climate change [38], particularly at its eastern limits [39]. Thus, assessing this species' water balance in relation to environmental factors is crucial. Still, to our knowledge, limited information is available about the radial patterns of sap flux and its seasonal and diurnal fluctuation, in relation to climate, in *P. halepensis* at eastern Mediterranean ecosystems. Mostly, if not solely, the relevant literature refers to Aleppo pine plantations in Israel, such as the works of Schiller and Cohen [40] who first reported on the radial profile of sap flux in this species and Cohen et al. [35] who studied the radial sap flux gradients of Aleppo pines under different climatic conditions.