*Article* **Integration of the Physiology, Transcriptome and Proteome Reveals the Molecular Mechanism of Drought Tolerance in** *Cupressus gigantea*

**Pei Lei 1,2,†, Zhi Liu 1,3,†, Jianxin Li 1, Guangze Jin 1, Liping Xu 1, Ximei Ji 1, Xiyang Zhao 2, Lei Tao <sup>1</sup> and Fanjuan Meng 1,2,\***


**Abstract:** Drought stress can dramatically impair woody plant growth and restrict the geographical distribution of many tree species. To better understand the dynamics between the response and mechanism of *Cupressus gigantea* to drought and post-drought recovery, a comparative analysis was performed, relying on physiological measurements, RNA sequencing (RNA-Seq) and two-dimensional gel electrophoresis (2-DE) proteins. In this study, the analyses revealed that photosynthesis was seriously inhibited, while osmolyte contents, antioxidant enzyme activity and non-enzymatic antioxidant contents were all increased under drought stress in seedlings. Re-watering led to a recovery in most of the parameters analyzed, mainly the photosynthetic parameters and osmolyte contents. Transcriptomic and proteomic profiling suggested that most of the differentially expressed genes (DEGs) and differentially expressed proteins (DEPs) were specifically altered, and a few were consistently altered. Drought induced a common reduction in the level of DEGs and DEPs associated with photosynthesis. Notably, DEGs and DEPs involved in reactive oxygen species (ROS) scavenging, such as ascorbate oxidase and superoxide dismutase (SOD), showed an inverse pattern under desiccation. This study may improve our understanding of the underlying molecular mechanisms of drought resistance in *C. gigantea* and paves the way for more detailed molecular analysis of the candidate genes.

**Keywords:** *Cupressus gigantea*; oxidative stress; antioxidant enzymes; gene expression; drought responses

#### **1. Introduction**

Environmental stresses are known to adversely affect plants' growth and distribution [1]. Among these many adverse factors, drought is a serious detrimental environmental factor constraining seed germination, plant growth and the economic value of crops [2]. However, being sessile organisms, plants are unable to escape environmental stresses and thus have evolved sophisticated strategies to acclimate to dehydration, including morphological, anatomical, physiological and molecular adaptive strategies [3]. Accordingly, the initial response of plants is closely related to the reduction in water evaporation, and consequent reductions in photosynthesis, transpiration, stomatal closure and the accumulation of osmolytes [4]. Of these, the decline in photosynthetic process under drought stress is mainly attributed to stomatal closure, reductions in CO2 fixation and disturbances in photosynthetic electron transport (PET) [5]. PET is composed of multi-subunit protein complexes, such as Photosystems I and II (PSI, PSII), the cytochrome b6/f complex (Cytb6/f) and ATPase synthase, which absorb light energy and transduce solar energy into chemical

**Citation:** Lei, P.; Liu, Z.; Li, J.; Jin, G.; Xu, L.; Ji, X.; Zhao, X.; Tao, L.; Meng, F. Integration of the Physiology, Transcriptome and Proteome Reveals the Molecular Mechanism of Drought Tolerance in *Cupressus gigantea*. *Forests* **2022**, *13*, 401. https:// doi.org/10.3390/f13030401

Academic Editor: Romà Ogaya

Received: 28 January 2022 Accepted: 23 February 2022 Published: 1 March 2022

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

**Copyright:** © 2022 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/).

energy [6]. Among these, PSII functions as a water-plastoquinone oxidoreductase embedded within the thylakoid membrane, which contains a series of peripheral light-harvesting complexes (LHC), and is vital to the initiation of photosynthesis and electron transport [7]. The Cyt b6/f complex, a rate-limiting step in photosynthesis, is involved the linear electron transport from PSII to PSI, giving rise to the production of ATP and NADPH [8]. Accumulation of some osmolytes such as proline and soluble sugar can play a prominent roles in preventing membrane disintegration and enzyme inactivation under water deficit [9].

When a reduction in water evaporation is insufficient to mitigate the stress stimulus, plants mainly respond to dehydration by activating a defense system of enzymatic and non-enzymatic antioxidants to cope with the overaccumulation of reactive oxygen species (ROS) [10]. Superoxide dismutase (SOD) is the first line of plant ROS defense, catalyzing the conversion of oxygen ions (O2 ·−) to oxygen (O2) and hydrogen peroxide (H2O2), which is then eliminated by the coordinated action of peroxidase (POD), catalase (CAT), ascorbate peroxidase (APX) and the AsA-GSH cycle [11]. Additionally, some enzymes, including APX, glutathione reductase (GR), dehydroascorbate reductase (DHAR) and monodehydroascorbate reductase (MDHAR), maintain the balance of AsA-GSH [12]. In addition, numerous drought-related transcripts and proteins are induced, which involve signaling transduction, activation/regulation of transcription, antioxidant capacity and ROS scavengers [13]. Thus, joint transcriptomic and proteomic profiling is considered to be effective and acceptable for disentangling the sophisticated processes of plants' responses to water deficit at the molecular levels [14,15].

More recently, to understand such complex regulatory mechanisms, integrated multiomics approaches (e.g., transcriptomics, proteomics and metabonomics) have been applied because multi-omics techniques offer powerful ways to reveal the relationships between genotype and phenotype. For example, comparative proteomic and transcriptomic analyses have provided insight into the formation mechanisms of seed size in castor bean [16]. Similarly, comparative transcriptomic and proteomic analyses were performed to determine the effect of pigment content on drought resistance in a wheat mutant [17]. In maize, an integrated transcriptomic, proteomic and metabolomic analysis found that its UV-B stress response involves signal transduction and signal molecules [18]. Dry-farm plants are promising candidates for studies on drought-related genes, proteins and metabolites [19]. While candidate genes and proteins have been discovered for wheat, rice, and soybean under dehydration conditions, no such reports are yet available for *C. gigantea* [20–22].

Being a rare tree species that is remarkably drought-tolerant, *Cupressus gigantea* W.C. Chen et L.K. Fu has high ecological and medicinal values due to its extensive use in afforestation, traditional Tibetan medicine and the construction industry [23]. It was added to the Red List of Threatened Species of the International Union for Conservation of Nature (IUCN), given the slow natural regeneration of this species, the anthropogenic disturbances it faces and its geographic isolation [24]. To date, studies on *C. gigantea* have mainly focused on its genetic diversity and population structure [25], total protein extraction [26], photosynthetic capacity [27], determination of the complete chloroplast genome [28] and comprehensive transcriptome characterization [29]. Nevertheless, to our best knowledge, no study has comprehensively investigated the molecular mechanisms enabling *C. gigantea* to tolerate drought. Overall, the mining of genes and proteins related to drought is an indispensable step towards deciphering the adaptive mechanisms of *C. gigantea*.

In this study, a cross-disciplinary approach combining classical physiological measurements, RNA sequencing (RNA-Seq) and two-dimensional gel electrophoresis (2-DE) was carried out to gain insight into the responses of *C. gigantea* to drought stress. This study identified differentially expressed genes (DEGs) and differentially expressed proteins (DEPs) linked to the tree's drought stress responses, which are of great and timely importance for understanding the mechanisms by which *C. gigantea* ameliorates the effects of water deficit at the physiological and molecular levels.

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

#### *2.1. Plant Material and Drought Treatments*

In this study, we selected *C. gigantea* from the Tibet Agriculture and Animal Husbandry College (Linzhi, Tibet, China). Experiments were conducted at the Northeast Forestry University of (Harbin, Heilongjiang, China). Three-year-old seedlings were sown in plastic pots (25 cm in diameter and 25 cm high) containing a 1:1 mix of perlite and soil. The soil was humus soil, with a pH in H2O of 4.5–5.5, an organic matter content of 12.5 g/kg, and available nitrogen (N), phosphorous (P) and potassium (K) contents of 100, 40 and 90 mg/kg, respectively. All seedlings were grown in a greenhouse at a temperature of 17 ◦C during the night and 22 ◦C during the daytime with a relative humidity of 60% and a 12 h photoperiod. Three-year-old seedlings were exposed to the drought treatment (DTs), in which irrigation of the plants was suspended for 21 consecutive days (hereafter referred to as DT0, DT7, DT14 and DT21). At the end of the stress treatment, the seedlings were re-watered over a period of 7 days (hereafter referred to as RW). Well-watered control plant seedlings (CKs) irrigated every 7 days served as the control group (hereafter referred to as CK0, CK7, CK14, CK21 and CK28). Leaves were collected from DT and CK plants. Three biological replicates were used for the total RNA and protein extractions, with another 5 biological replicates measured for their photosynthesis and physiology parameters.

#### *2.2. Physiological Measurements*

To evaluate the physiological changes in *C. gigantea* under drought stress, the net photosynthetic rate (*P*n), stomatal conductance (*G*s), intercellular CO2 concentration (*C*i) and transpiration rate (*T*r) of the leaves were measured using a Li-6400 portable photosynthesis system (LI-COR Inc, Lincoln, NE, USA) at 9:00–11:00 a.m. The CO2 concentration, stimulation light intensity, and gas flow rate were set as 450 <sup>μ</sup>mol·mol<sup>−</sup>1, 1000 <sup>μ</sup>mol·m−2·s−<sup>1</sup> and <sup>500</sup> <sup>μ</sup>mol·s<sup>−</sup>1, respectively. Five plants were chosen from the control and drought-stressed groups, and each plant was measured 5 times.

#### *2.3. Biochemical Analysis*

The leaves were infiltrated with 0.1 mg/mL 3, 3 -diaminobenzidin (DAB) (Sigma, USA) in 50 mM Tris-acetate buffer (pH 5.0) for H2O2 staining or with 0.1 mg/mL nitro blue tetrazolium (NBT) (Sigma, USA) in 25 mM K-HEPES (pH 7.6) for superoxide staining, according to the method described by Chen [30].

Proline from the leaves was extracted by the method described by Zhang [31]. The leaves (0.3 g) were extracted in 4 mL of 1% sulfosalicylic acid. After centrifugation at 6000× *g* for 5 min, 1 mL of the supernatant was added to 1.5 mL of glacial acetic acid and 2 mL of a ninhydrin solution. The mixture was incubated at 100 ◦C for 30 min. After cooling, the contents were separated with 2 mL toluene, and the optical density was measured at 520 nm.

The soluble sugars of the leaves were analyzed by the anthrone method as described by Zhao [32]. Tissue oven-dried at 65 ◦C for 24 h (100 mg) was added to 10 mL of 80% ethanol and incubated in a thermostat water bath at 90 ◦C for 10 min, then the supernatant was collected. The pellet was extracted again as described above, and the supernatant was obtained and combined with the previous aliquot. After adding 3.5 mL of anthrone reagent to 0.1 mL of the supernatant, the mixture was heated in a thermostat water bath at 90 ◦C for 10 min. After the mixture had cooled to room temperature, the absorbance was measured at 620 nm.

The SOD (EC 1.15.1.1) activity of the leaves was detected with NBT [33]. To measure enzyme activity, leaves (0.3 g) were ground to a fine powder in liquid nitrogen and dissolved in 4 mL of potassium phosphate buffer (PBS) (50 mM, pH = 7.8). The assay mixture contained 50 mM PBS (pH = 7.8), 195 mM methionine, 0.3 mM ethylene diamine tetraacetic acid, 1.125 mM NBT, 70 μL of the extracting solution and 60 μM riboflavin. Enzyme activity was detected at 560 nm by a spectrophotometer.

POD (EC 1.11.1.7) and CAT (EC 1.11.1.6) activity was measured in leaves according to the method of Kosar et al. [34]. The activity of POD was assayed in a mix containing 50 μL of the extracting solution, 50 mM PBS, 14 μL guaiacol and 19 μL H2O2 (30%, *v/v*). Enzyme activity was measured at 470 nm. The activity of CAT was assayed in a mix containing 0.2 mL of the extracting solution, 50 mM PBS (50 mM, pH = 7.8), 500 μL water and 10 mM H2O2, and the decrease in absorbance at 240 nm was monitored for 3 min.

The contents of ascorbate and glutathione were estimated in leaves according to the protocol of Sun et al. [35]. The activity of ASA was assayed in a mix containing 0.1 mL of the extracting solution, 0.7 mL ddH2O, 10% TCA, 44% H3PO4, 4% 2,2 -dipyridyl and 3% FeCl3, then the mixture was heated in a thermostat water bath at 37 ◦C for 60 min. After the mixture had cooled to room temperature, the absorbance was measured at 525 nm. The activity of ASA + DHA was assayed in a mix containing 0.1 mL of the extracting solution, 0.5 mL ddH2O, 10 mM DTT, 10% TCA, 5% N-ethylmaleimide, 44% H3PO4, 4% 2,2 -dipyridyl and 3% FeCl3, then the mixture was heated in a thermostat water bath at 37 ◦C for 60 min. After the mixture had cooled to room temperature, the absorbance was measured at 525 nm. The activity of GSH was assayed in a mix containing 0.2 mL of the extracting solution, 150 mM (pH = 7.4) NaH2PO4 and 0.6 mM DTNB, then the mixture was heated in a thermostat water bath at 30 ◦C for 5 min. After the mixture had cooled to room temperature, the absorbance was measured at 412 nm. The activity of GSH + GSSG was assayed in a mix containing 0.2 mL of the extracting solution, 50 mM (pH = 7.4) PBS, 0.6 mM DTNB, 2 mM NADPH and 2 U GR, then the mixture was heated in a thermostat water bath at 25 ◦C for 10 min. After the mixture had cooled to room temperature, the absorbance was measured at 412 nm.

The activity of APX (EC 1.11.1.11), DHAR, GR and MDHAR was determined in leaves as described previously [36]. To measure the enzyme activity, leaves (0.3 g) were ground to a fine powder in liquid nitrogen and dissolved in 4 mL of a potassium phosphate buffer (PBS) (50 mM, pH = 7.8). To measure APX, the extraction solution was added to a mix containing 50 mM PBS (50 mM, pH = 7.8), 5 mM ASA and 20 mM H2O2, and the decrease in absorbance at 290 nm was monitored for 3 min. To measure DHAR, the extraction solution was added to a mix containing 50 mM PBS (50 mM, pH = 7.8), 0.5 mM DHA and 5 mM GSH, and the decrease in absorbance at 265 nm was monitored for 3 min. To measure GR, the extraction solution was added to a mix containing 50 mM PBS (50 mM, pH = 7.8), 2 mM NaDPH and 10 mM GSSG, and the decrease in absorbance at 340 nm was monitored for 3 min. To measure MDHAR, the extraction solution was added to a mix containing 50 mM PBS (pH = 7.8), 0.1 mM ASA and 0.55 U AAO, and the decrease in absorbance at 340 nm was monitored for 3 min.

#### *2.4. Protein Preparation and Quantitative Proteome Analysis*

In each sample, protein was extracted from the leaves using the TCA/acetone method [37]. Plant tissue (2 g) frozen by liquid nitrogen was ground into powder, then suspended in 4 mL of a cold extraction buffer (0.5 M Tris-HCl (SRL, India) pH = 7.5, 0.7 M sucrose (Himedia, India)) for protein solubilization. Samples were incubated overnight at 4 ◦C, after which an equal volume of phenol saturated with Tris-HCl (pH = 7.5) was added to each sample. The homogenate was centrifuged (5000× *g* at 4 ◦C for 30 min) and the clear white pellet was washed three times with ice-cold acetone (Himedia, Mumbai, India) at −20 ◦C, followed by centrifugation at 5000× *g* for 5 min. The pellet was then freeze-dried and stored at −80 ◦C. The protein powder was dissolved in a rehydration buffer (8 M urea (Sigma, St. Louis, MO, USA), 2% (*w*/*v*) CHAPS (Sigma, St. Louis, MO, USA), 50 mM DTT (Sigma, St. Louis, MO, USA), and 0.2% (*v*/*v*) Biolyte (Bio-Rad, Hercules, CA, USA)) for 1 h at 37 ◦C. After centrifugation, the protein concentration of the supernatant was determined by the Bradford assay (Bio-Rad, Hercules, CA, USA).

The 2-DE was carried out according to the methodology of Wang et al. [38]. First, the immobilized pH gradient (IPG) strips [pH 4–7, 13 cm, (GE Healthcare, Chicago, IL, USA)] were separated. Next, 1.5 mg of the protein sample was applied to each IPG strip and this was covered with 800 μL of mineral oil. The voltages used were as follows: 30 V for 13 h, 100 V for 1 h, 500 V for 1 h, 1000 V for 1 h and 3000 V for 1 h, and then 8000 V for 7 h. After focusing, the strips were equilibrated in Equilibration Buffer I (0.1 g DTT and 10 mL of a SDS balancing buffer) and in Equilibration Buffer II (0.15 g iodoacetamide (IAM) (Sigma, St. Louis, MO, USA) and 10 mL of a SDS balancing buffer). After 15 min, the strips were washed with a SDS balancing buffer. The equilibrated strips were analyzed by 12.5% sodium dodecyl sulfate–polyacrylamide gel electrophoresis. The electrophoresis was run on Electrophoreses Power Supply EPS 601 (GE Healthcare, Chicago, IL, USA) at a constant current of 2 W for 30 min and then modulated to 8 W. The gels were visualized by staining them overnight with colloidal Coomassie brilliant blue R-250. Stained gels were scanned on an Image Scanner II (GE Healthcare, Chicago, IL, USA), and the captured images were analyzed in Melanie 7.0 software (GeneBio, Geneva, Switzerland). Three independent replicated gels were analyzed and characterized. Protein spot detection was based on a fold-change of ≥2 or ≤0.5, for which a threshold of *p* ≤ 0.05 was used to distinguish the differentially expressed protein spots.

The protein spots were carefully excised from the gels. The gel spots were washed twice for 20 min with deionized water, then incubated and dehydrated with acetonitrile (ACN). Proteins were digested for 18 h at 37 ◦C in 10 μL of a trypsin solution (15 ng μL<sup>−</sup>1). Next, the supernatants were collected, and the gel spots were extracted twice with a 50 μL extraction buffer (50% ACN (Sigma, St. Louis, MO, USA) and 5% TFA) for 1 h at 37 ◦C. The extractions and the trypsin supernatant of the gel spots were combined and then vacuumdried. The peptides' mass spectra were detected via matrix-assisted laser desorption ionization time-of-flight/time-of-flight (MALDI-TOF/TOF) (ABI 4700, AB Systems, CA, USA) and the proteins were identified using the UniProt database (http://www.uniprot.org, accessed on 22 February 2020). The gene ontology of identified proteins was determined using the Blast2Go v2.3.6, with the Target P program used for their functional classification.

#### *2.5. RNA-Seq Library Construction and Transcriptome Analysis*

The RNA-Seq libraries' preparation and their sequencing were carried out in leaves by Novogene equipment (Beijing, China). Briefly, RNA degradation and contamination were monitored on 1% agarose gels, RNA purity was checked using a NanoPhotometer® spectrophotometer (Implen, Westlake Village, CA, USA) and RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA), with 1 μg RNA per sample used as the input material for RNA sample preparation. The sequencing libraries were generated using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA), following the manufacturer's recommendations, with unique index codes added to attribute the sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System, using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA), according to the manufacturer's instructions. After cluster generation, the libraries of the preparations were sequenced on an Illumina NovaSeq platform, from which 150-bp paired-end reads were generated. Both the reference genome and the gene model annotation files were downloaded from the genome website directly. Through use of Hisat2 v2.0.5 (Johns Hopkins University, Baltimore, MD, USA), the index of the reference genome was built, and paired-end clean reads were aligned to it. The clean data were deposited in NCBI Sequence Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra, accessed on 18 July 2021) under the accession number SRR15183944-SRR15183948.

Gene ontology (GO) enrichment analysis of the DEGs was implemented by the R package 'clusterProfiler', which corrected any gene length bias. GO terms with a corrected *p*-value of <0.05 were considered to be significantly enriched by the DEGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding the high-level functions and utilities of biological systems (http://www.genome.jp/kegg, accessed on 22 March 2020). The clusterProfiler package in R was used to statistically test for the enrichment of DEGs in the KEGG pathways.

#### *2.6. Quantitative Real-Time PCR (RT-PCR)*

Total RNA was extracted by using the OmniPlant RNA Kit (DNase I) (Cowin Biosciences, Beijin, China). A summary of the procedure is as follows: the samples were individually milled in a mortar with liquid nitrogen and then incubated with 500 μL of RLS Buffer. The samples were centrifuged at 12,000 rpm for 2 min at 4 ◦C, and a half volume of absolute ethanol was added to the supernatant after centrifuging, then 52 μL RNase-Free Water, 8 μL 10 × Reaction Buffer and 20 μL DNase I was then added to the spin column after incubating at room temperature for 15 min. The samples were centrifuged at 12,000 rpm for 2 min at 4 ◦C, then the precipitates were washed in RNA Wash Buffer II and dried. RNA was dissolved with RNase-Free Water, then 0.5 μg of total RNA was used to synthesize the first-strand cDNA by using the ReverTra Ace® qPCR RT Master Mix with the gDNA Remover (Toyobo, Osaka, Japan). Amplifications were performed on a Roche 480 PCR System (Roche, Rotkreuz, Switzerland) with a THUNDERBIRD® qPCR Mix (Toyobo, Osaka, Japan). The primer sequences used for the real-time qPCR are listed in Table S7.

#### *2.7. Statistical Analysis*

The data were subject to analysis of variance (ANOVA) (*p* < 0.05). Multiple comparisons of the treatment effects were analyzed using the least significant difference (LSD). All the tests were performed using SPSS 22.0 software (Windows, USA). All figures were developed using GraphPad Prism 8.2 (GraphPad Software, CA, USA) and Adobe Illustrator CC2018 (Adobe, CA, USA).

#### **3. Results**

#### *3.1. Drought Stress Induced Growth and Physiological Changes in C. gigantea*

To investigate the morphological and physiological response mechanisms of *C. gigantea* under drought stress, 3-year-old seedlings were planted in the same environment and then subjected to a drought treatment for 3 weeks, then re-watered for 1 week (Figure 1a). At the onset of stress, no significant differences were observed in plant growth. After 21 days, seedlings in the drought treatment had leaves and stems that were badly curled and wilted when compared with those of the control group. Interestingly, the droughtstressed seedlings were able to largely restore their viability during the re-watering period. Upon closer examination, the seedlings' *P*n, *G*s, and *T*r significantly decreased by 38.8%, 72.3% and 57.9%, respectively, whereas *C*i gradually increased by 20.9% after completing the drought treatment compared with the control (*p* < 0.05). Surprisingly, RW induced a rapid recovery in *C. gigantea*, reaching parameters similar to those recorded in the control group (Figure 1b–e), indicating that photosynthesis was inhibited by drought.

Physiological indexes such as histochemical staining, soluble sugar content, proline content and antioxidant enzymes are commonly used to evaluate the stress resistance capacity of plants in response to drought stress. Here, DAB and NBT histochemical staining analyses were conducted to measure the H2O2 and O2 − contents. Under normal conditions, there was no apparent difference between the treatment and control; however, during prolonged drought stress, the DAB and NBT staining intensities of the treatment plants' leaves were stronger than those of the control leaves (Figure 2a,b). Importantly, we found that the proline and soluble sugar contents were higher in the seedlings under the drought stress treatment than in the control (Figure 2c,d). Furthermore, sharp increases of 17.9%, 68.6% and 82.6%, respectively, were observed in SOD, POD and CAT activity (*p* < 0.05) after 21 days of the drought treatment (Figure 2e–g), indicating that *C. gigantea* promoted its ROS scavenging by modulating the activity of key antioxidant enzymes, such as SOD, POD and CAT. Similarly, non-enzymatic antioxidant contents, namely those of ascorbic acid (ASA), dehydroascorbic acid (DHA), reduced glutathione (GSH) and oxidized glutathione (GSSG), were starkly increased under the drought stress treatment relative to the control seedlings (*p* < 0.05, Figure 2h–k). The antioxidant enzyme activities of APX, DHAR, GR and MDHAR were further analyzed; this showed that they increased in content after

the drought treatment (*p* < 0.05), and their accumulation was highest at 21 days of the drought treatment (*p* < 0.05, Figure 2i–o). Subsequent to rehydration of the seedlings', the parameters of those in the treatment returned to similar levels to those of the control group. Overall, these results revealed that seedlings of *C. gigantea* enhanced their drought resistance by increasing their osmolyte content, activating the antioxidant system and maintaining ROS homeostasis.

**Figure 1.** The effect of drought stress on the morphological and photosynthesis traits of *C. gigantea*. (**a**) Photographs showing the phenotypes of the seedlings over time. (**b**–**e**) Line charts showing the physiological changes corresponding to *P*n, *G*s, *T*r, *C*i, respectively. Values are the means ± SD (*n* = 5). Columns headed by different letters denote means differing significantly from one another (*p* < 0.05).

**Figure 2.** The effect of drought stress on the physiological indexes of *C. gigantea*. (**a**,**b**) Photographs showing DAB and NBT histochemical staining of the seedlings over time. (**c**–**o**) Line charts showing the physiological changes corresponding to soluble sugar, proline, superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), ascorbic acid (ASA), dehydroascorbic acid (DHA), reduced glutathione (GSH), oxidized glutathione (GSSG), ascorbate peroxidase (APX), dehydroascorbate reductase (DHAR), glutathione reductase (GR) and monodehydroascorbate reductase (MDHAR), respectively. Values are the means ± SD (*n* = 3). Columns headed by different letters denote means differing significantly from one another (*p* < 0.05).

#### *3.2. Characterization and Functional Distribution of the C. gigantea Leaf Proteome*

To identify the drought-induced proteins of *C. gigantea*, we performed a 2-DE analysis coupled with MALDI-TOF/TOF at different time points (7, 14, 21, and 28 days) in seedlings under drought stress. Total protein extracts were separated in all samples, with a pI range of 4–7 and a molecular mass range of 14.4–116 kDa. All these protein spots were localized and detected on Coomassie brilliant blue-stained gels. Accordingly, 66 protein spots changed significantly in abundance between the drought-stressed and control seedlings of *C. gigantea* (Figure S1). Temporally, of these, 26 (7 days), 16 (14 days), 21 (21 days) and 3 (28 days) protein spots were found to be altered by drought stress.

In total, 66 DEPs were identified in the control vs. treatment comparison: 61 (i.e., 23 upregulated and 39 downregulated) in the 7DT vs. 7CK comparison, 62 (21 upregulated and 41 downregulated) in the 14DT vs. 14CK comparison, 66 (36 upregulated and 30 downregulated) in the 21DT vs. 21CK comparison and 64 (29 upregulated and 35 downregulated) in the 28RW vs. 28CK comparison (Figure 3a,b). In addition, there were four proteins

(Spots 6, 10, 11 and 17) whose expression was significantly upregulated by drought, but was upregulated/unchanged during the recovery period. Conversely, the expression of two proteins (Spots 36 and 39) was upregulated/unchanged by drought yet downregulated during the recovery phase of the seedlings (Figure S2).

**Figure 3.** Identification and functional classification analysis of the proteomics data from *C. gigantea* seedlings. (**a**) Bar graph showing the upregulated and downregulated differentially expressed proteins (DEPs) for each pairwise comparison. (**b**) Venn diagrams depicting the overlaps of upregulated and downregulated differentially expressed proteins (DEPs) across four comparisons. The total number of differentially expressed proteins (DEPs) is provided in parentheses. (**c**) Functional categorization of the differentially expressed proteins (DEPs) identified in *C. gigantea* in response to drought.

The experimental MW and pI values of the protein spots were then estimated and compared with the theoretical counterparts of the corresponding proteins. As Table S1 shows, most of the experimental values matched up well with the theoretical ones, indicating unambiguous identifications. According to their biological function, these identified proteins were classified into nine major groups: stress response (27.27%), photosynthesis (16.67%), energy (9.09%), metabolism (12.12%), photosynthesis electron transfer (4.55%), protein transport (3.03%), amino acid synthesis (10.61%), protein synthesis and processing (12.12%), and unknown (4.55%) under the drought stress treatment (Figure 3c).

#### *3.3. Characterization and Functional Distribution of C. gigantea's Leaf Transcriptome*

To further investigate the regulatory network of *C. gigantea*'s response to drought stress, cDNA libraries were built using leaves from seedlings under the two treatments: severe stress (21 days) and re-watering (28 days), and these were sequenced on an Illumina NovaSeq platform (Illumina Inc, San Diego, CA, USA). High-quality reads of the treatment and control samples were acquired after implementing data-cleaning procedures (Table S2). The length distribution of the assembled transcripts of *C. gigantea* was obtained (Figure S3a). To validate the accuracy and reproducibility of our RNA-Seq data, we selected 15 DEGs for a follow-up qRT-PCR; their gene expression levels as inferred by RNA-Seq were strongly correlated with those from the qRT-PCR (R2 = 0.92; Figure S3c).

Many DEGs were identified via two pairwise comparisons between the treatment and control (i.e., 21DT vs. 21CK, and 28RW vs. 28CK). Under severe stress, 215 genes were expressed differentially in 21DT vs. 21CK, including 36 upregulated genes and 84 downregulated genes, while 442 genes were expressed differentially under re-watering in 28R vs. 28CK, including 250 upregulated and 97 downregulated genes (*p* < 0.05, fold-change ≥ 2). Another 30 upregulated and 65 downregulated genes overlapped in the pairwise comparisons between 21DT and 21CK, and between 28RW and 28CK (Table S3, Figure S3b).

To illustrate the biological functions and pathways of the DEGs, GO annotation (at *p* < 0.05) was performed to analyze the drought-responsive genes in *C. gigantea* (Table S4), and 24 enriched GO terms were subsequently selected (Table S4, Figure 4). According to the GO classification graphic, "GO:0050896: response to stimulus", "GO:0016020: membrane" and "GO:0003824: catalytic activity" were highly enriched during both severe drought stress and re-watering of the seedlings, suggesting that the DEGs identified are likely involved in the modulation of ROS scavenging. All DEGs were further examined for KEGG pathway enrichment (Table S5, Figure 5). The top 30 enriched KEGG pathways of DEGs between the 21DT and 21CK, and between the 28RW and 28CK groups are shown in Figure 5a,b. The pathway of phototransduction, starch and sucrose metabolism, as well as plant hormone signal transduction, were also significantly enriched for DEGs between 21DT and 21CK, and between 28RW and 28CK. Accordingly, we performed a careful manual annotation of our list of DEGs to identify candidate genes that were directly involved in ROS scavenging and photosynthesis. Based on a combination of the GO and KEGG pathway results, we identified 12 such candidate DEGs (Table S6). According to our manual annotation, eight DEGs could be assigned generically to ROS scavenging, with another four DEGs assigned to photosynthesis; this suggested that these 12 candidate DEGs are probably key regulatory components of the drought stress response in *C. gigantea*.

**Figure 4.** Differentially expressed genes (DEGs) identified by Gene Ontology (GO) classification in *C. gigantea*.

**Figure 5.** Differentially expressed genes (DEGs) identified by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis in *C. gigantea*. (**a**) Top 30 KEGG enrichment results in 21DT vs. 21CK. (**b**) Top 30 KEGG enrichment results in 28RW vs. 28CK.

#### **4. Discussion**

Much accumulated empirical evidence has indicated that drought can alter cell membrane fluidity, disrupt protein complexes and inhibit photosynthesis. Therefore, a comprehensive study that combined physiological measurements, RNA-Seq and 2-DE analysis was conducted here to screen out some important candidate drought-induced changes in key physiological parameters, DEGs and DEPs. Our results revealed that the transcriptional regulation of proteins related to photosynthesis and the post-transcriptional regulation of genes involved in ROS scavenging might be the reasons for drought resistance in *C. gigantea.*

As it is an important secondary messenger, the rapid accumulation of ROS leads to rapid stomatal closure to limit photosynthesis under conditions of water deficit [39]. A recent study found that the *P*n, *G*s and *T*r of rubber trees decreased significantly under dehydration conditions [40]. In our study, similar results were detected in the droughtstressed seedlings (Figure 1b–e). Additionally, how a plant responds to stress not only entails changed photosynthesis but also the accumulation of osmolytes, such as soluble sugars and proline, which reflect the degree of injury [41]. We found that the content of soluble sugars and proline both enlarged as the drought stress continued (Figure 2c,d). Similar findings were reported recently for *Scutellaria baicalensis* Georgi [42].

In the present study, histochemical staining by both DAB and NBT together revealed a pronounced drought-induced accumulation of H2O2 and O2 ·− in *C. gigantea* leaves (Figure 2a,b). This result is consistent with earlier findings that drought promotes the generation of ROS [43]. Moreover, we demonstrated that the drought-stressed seedlings showed higher SOD, POD and CAT activity compared with the control group under drought stress (Figure 2e–g). This is consistent with another analysis of antioxidant enzymes in plants under drought conditions [44]. Similarly, drought stress significantly enhanced the contents of AsA, DHA, GSH and GSSG during drought, and these antioxidant molecules peaked under severe drought conditions (Figure 2h–o), which also is in line with prior research [45]. Compared with the control group, the characteristic withered and curved growth phenotype of the treatment seedlings suggested that *C. gigantea* likely relies on limited photosynthesis, the accumulation of osmolytes and ROS scavenging strategy to adapt to drought conditions.

Underlying many of these changes is the differential expression of genes in response to drought stress. Thus, both RNA-Seq and 2-DE were used to further explore the critical genes and pathways to obtain a better understanding of the molecular mechanisms underlying the response to drought. Here, we successfully quantified and compared 551 DEGs and 66 DEPs in the leaves of *C. gigantea* seedlings under drought stress. However, only a few genes were commonly regulated at the transcriptomic and proteomic levels, although this is consistent with a previous study [46], whose data indicated that post-transcriptional regulation plays a crucial role in the drought response of cassava plants. In the case of *C. gigantea*, the categories of overlapping DEGs and DEPs were most related to photosynthesis and ROS scavenging, suggesting their involvement in the drought response of this tree species.

Photosynthesis the most fundamental and intricate physiological process in plants, and it is highly susceptible to drought [47]. Consistent with previous studies of photosynthesisrelated genes' responses to drought, we found that genes involved in the Cyt b6/f complex and ATPase synthase were all downregulated after seedlings experienced the drought treatment [48]. Similarly, we found that many DEPs associated with PSII, LHC and ATP synthase were also downregulated, as inferred by the proteomic analysis (Figure S2). Similar responses in maize plants were reported recently [49]. Taken together, our results all implied that modulating PET is probably crucial to increasing the tolerance of *C. gigantea* to drought.

The ROS scavenging system, an important defense mechanism against drought, maintains the cellular homeostasis of ROS. Overproduction of ROS during drought not only affects the stress-related proteins but also modulates the transcription level of genes [50]. In

the current work, four of these enzymes, including GST, SOD, APX and POD, were identified. Among these genes, the *SOD* gene was downregulated in *C. gigantea* seedlings under drought stress. A significant reduction in the expression of the *APX* gene was obtained under drought stress conditions. Notably, our proteomics data showed that the abundance of APX, SOD [Cu-Zn], ASR and GST was also enhanced in seedlings of *C. gigantea* to cope with the imposed drought stress. These results are consistent with the observed activities of antioxidant enzymes. These proteins could also work together to maintain the dynamic balance of ROS levels and thereby avoid the oxidative damage incurred by plants. This type of coordination occurs in alkaligrass [51], rapeseed [52] and cassava [53], perhaps due to the fact that enzyme activity is regulated by post-translational modifications [54]. Finally, POD can also contribute to how plants respond to drought stress [5]. Nevertheless, in this study, drought stress led to downregulation of POD at both the protein and gene levels, but, interestingly, the activity of POD was significantly increased by drought stress. Evidently, to explain the changes in enzyme activity and the expression of genes and proteins of *C. gigantea* under drought stress more precisely, further in-depth research investigations and analyses are still needed.

#### **5. Conclusions**

The main aim of this study was to expand our knowledge about the mechanism underlying the impact of drought stress in seedlings of *C. gigantea*. An attempt was made to identify the critical genes and pathways by comparing the findings from transcriptomic and proteomic analytical methods. On the basis of results, here, we propose a schematic diagram for the regulation of drought resistance in *C. gigantea* (Figure 6). When plants are subjected to drought, ROS induces oxidative stress and activates the ROS scavenging system. Simultaneously, oxidative stress impairs photosynthesis and inhibits PET, ultimately inhibiting plant growth. Taken together, the omics data suggest that post-transcriptional regulation of PET and ROS scavenging plays a crucial role in the drought response of *C. gigantea* seedlings.

**Figure 6.** Schematic diagram of drought-induced changes in the physiological parameters, genes and proteins in *C. gigantea* during stress and recovery.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/f13030401/s1. Figure S1: 2-DE gel images of the control and treatment seedlings under drought extracted at different time points (7 d, 14 d, 21 d and recovery). The squares represent the drought-responsive proteins identified by MS. Figure S2: Expression of DEPs in the control and drought treatment at different time points (7 d, 14 d, 21 d and recovery). Figure S3: Characteristics of *C. gigantea* unigenes. (**a**) Length frequency distribution of all unigenes. (**b**) Venn diagram showing the number of DEGs. (**c**) Correlation analysis between the transcriptome results and the qRT-PCR results, and validation of the RNA-Seq results by qRT-PCR. Table S1: The corresponding induction factor (percent volume of the spot under stress conditions/percent volume of the spots under control conditions) of drought-responsive proteins of seedlings. Table S2: Summary of unigenes identified in *C. gigantea*. Table S3: Summary of DEGs identified in *C. gigantea*. FDR ≤ 0.001 log2 fold change ≤ −2 or ≥2. Table S4: List of GO enrichment terms of DEGs in *C. gigantea* during drought stress. Table S5: KEGG pathway enrichment analysis of DEGs in *C. gigantea* during drought stress. Table S6: DEGs and DEPs involved in ROS scavenging and photosynthesis in *C. gigantea*s during drought stress. Table S7: Primers of the genes used for qRT-PCR and the expression correlation between the RNA-seq and qRT-PCR analyses.

**Author Contributions:** P.L., Z.L., L.T. and F.M., conceived and designed the experiments; P.L., Z.L. and J.L. performed the experiments; P.L., Z.L. and L.X. analyzed or interpreted the data for the work; P.L., X.Z., L.T. and F.M. wrote the manuscript; P.L., X.J., L.X., G.J., X.Z. and F.M. ensure that questions related to the accuracy or integrity of any part of the work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project was supported by the Fundamental Research Funds for the Central Universities (2572021AW15).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Sequencing data have been deposited at the Short Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra, accessed on 18 July 2021) under accession numbers SAMN20236055 to SAMN20236058 (Bio-Project PRJNA746687). The proteome data can be found at: https://github.com/lppaper/datebase, accessed on 18 February 2022.

**Acknowledgments:** We thank the reviewers and editors for their work.

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

#### **References**


### *Article* **Water Retention Capacity of Leaf Litter According to Field Lysimetry**

**Taehyun Kim 1, Jungyoon Kim 1, Jeman Lee 1, Hyun Seok Kim 1,2,3, Juhan Park <sup>4</sup> and Sangjun Im 1,2,\***


**Abstract:** The water retention capacity of forest leaf litter was estimated through lysimeter measurements under field conditions. Six lysimeters were placed in *Pinus koraiensis* and *Quercus acutissima* forests and filled with the surrounding leaf litter to represent the effects of litter type on the water retention capacity. Two years of measurements for rainfall and litter weight have been conducted in all lysimeters at 30 min intervals. Field measurements showed that *P. koraiensis* litter stored more water during rainfall periods than did *Q. acutissima* litter. As a result, immediately after the cessation of rainfall, 1.82 mm and 3.00 mm of water were retained per unit mass of *Q. acutissima* and *P. koraiensis* litter, respectively. Following rainfall, after the gravitational flow had entirely drained, the remaining water adhered to the litter was estimated to be 1.66 ± 1.72 mm and 2.72 ± 2.82 mm per unit mass per rainfall event for *Q. acutissima* and *P. koraiensis* litter, respectively. During the study period, approximately 83.7% of incident rainfall drained into the uppermost soil layer below the *Q. acutissima* litter, whereas 84.5% of rainfall percolated through the *P. koraiensis* litter. The moisture depletion curves indicated that 50% of the water retained in the *Q. acutissima* and *P. koraiensis* litter was lost via evaporation within 27 h and 90 h after the cessation of rainfall, respectively. This study demonstrated the water retention storage of leaf litter and its contribution to the water balance over floor litter according to litter and rainfall characteristics. The results also proved that lysimetry is a reliable method to quantify the variation of litter moisture under natural conditions.

**Keywords:** litter lysimeter; water retention capacity; litter drainage; *Pinus koraiensis* litter; *Quercus acutissima* litter

#### **1. Introduction**

Forests occupy approximately 63% of the land area of the Republic of Korea. The national forest restoration project led to an increase in the volume of forests from 7.3 m3/ha in 1955 to 165.2 m3/ha in 2020 [1], thereby contributing to the development of a thick layer of leaf litter on the forest floors. Fallen leaves resting on the underlying mineral soil are the source of most plant nutrients and provide habitats for a great diversity of organisms. Hydrologically, leaf litter on the forest floor acts as a porous interface between mineral soil and ambient air and influences the processes of subsurface runoff, overland flow, and soil erosion [2,3]. Leaf litter also influences the physical properties of soils, thereby increasing the ability of water to infiltrate into the soil [2,4].

Rainfall partitioning occurring on the forest floor is a minor but fundamental component of the water balance in forest watersheds [4]. A small amount of rainwater is temporarily captured by the surface of leaf litter, and the retained water is lost via evaporation within a few hours or days after the rainfall. When the amount of intercepted water exceeds the threshold, known as the water retention capacity, water percolates into the

**Citation:** Kim, T.; Kim, J.; Lee, J.; Kim, H.S.; Park, J.; Im, S. Water Retention Capacity of Leaf Litter According to Field Lysimetry. *Forests* **2023**, *14*, 478. https://doi.org/10.3390/f14030478

Academic Editors: Yanhui Wang, Karl-Heinz Feger and Lulu Zhang

Received: 8 February 2023 Revised: 21 February 2023 Accepted: 24 February 2023 Published: 27 February 2023

**Copyright:** © 2023 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/).

soil along litter funnels under the influence of gravity [2,5]. These hydrological processes repeatedly occur during successive rainfall events.

Leaf litter intercepts rainfall much in the same way as canopy interception. The extent of rainfall retained in the litter layer is highly related to the litter's ability of water absorption. Where water retention storage is small, less water is available for evaporation, but more water reaches the upper soil layer for runoff. Water retention capacity also plays a major role in predicting the ignition and spread rate of forest fires [2–4].

The water retention capacity of litter differs among leaf types [6–8]. Owing to the large size and curved shape of the leaves, broadleaf litter has an advantage in storing water on the leaf surface during rainfall [9]. Conversely, the low porosity of the needle shape of coniferous litter contributes to high water retention within the litter as a consequence of the greater resistance of the packed litter to the percolation of water [10]. Considerably more water is required to completely saturate thicker needle litter [11] before it moves into the soil underneath.

Numerous attempts have been made to continuously measure the water retention storage of floor litter over the past couple of decades. According to Helvey and Patric [11], these methods may categorize into laboratory and field experiments. A detailed review of measuring techniques can be found in Gerrits and Savenije [12].

Several laboratory experiments have been conducted to analyze the water dynamics of leaf litter as it relates to retention and drainage capacity [9,10,13,14]. These experiments have provided a reliable estimation; however, accuracy is limited by several factors. First, the hydrologic response of leaf litter under simulated rainfall in a lab likely differs from that under real incidents of rainfall. Further, inconsistencies in data from natural rainfall events can be attributed to spatial and temporal variations, which would not occur in a lab setting. In addition, the relatively short duration and high intensity of artificial rainfall in laboratory studies cannot completely saturate the leaf accumulation, thereby leading to low measurements of water retention capacity [10], which may be inaccurate.

To overcome some of the drawbacks of laboratory experiments, field lysimeters have recently been employed in experimental forest hydrology. This device was designed to measure the weight of the litter within a certain time interval. Schaap and Bouten [15] developed a concrete lysimeter equipped with a load cell sensor. The present study is based on a field experiment conducted by Gerrits et al. [16], similar to that of Schaap and Bouten [15], to measure the weight of the litter layer. Lysimetric techniques are scalable and suitable for directly quantifying the water budget of a forest floor through changes in litter weight over the shortest intervals, but the challenges of using this tool include high manufacturing and operational costs [17].

Forests can retain a significant amount of rainfall via foliage, branches, stems, and leaf litter. Although the litter layer can store only a few millimeters of water, water absorption and depletion processes of litter can result in a considerable reduction of rainfall reaching to soil surface. However, the role of litter on rainfall interception and retention remains unclear and may even be disregarded in many hydrological studies because it is commonly considered a minor process and has practical difficulties in making accurate measurements [8,13]. Therefore, an accurate understanding of how rainfall is partitioned over the floor litter can provide a better perspective on rainfall interception and evaporation processes.

The objectives of this study were to measure the water stored in the litter layer under natural conditions and to examine the effects of litter type and rainfall amount on the water balance of the litter layer. In this study, six lysimeters were placed in deciduous and coniferous forests, and litter weight was measured during the experimental period. The reliable estimation of water retention capacity can contribute to understanding the hydrologic processes of the forest floor litter layer and can be utilized to predict fire ignition potential based on fuel moisture variation.

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

#### *2.1. Study Area*

This study was conducted between November 2015 and November 2017 in the Mt. Taehwa University Forest (TUF) of Seoul National University, located approximately 30 km southeast of the Seoul capital area (N37◦18 –20 , E127◦17 –20 ) (Figure 1). The TUF encompasses a 796-ha mountain forest at a latitude ranging from 150 m to 644 m above sea level. The climate of the TUF is characterized as temperate, with hot wet summers and cold dry winters. The mean monthly temperature ranges from −3.3 ◦C in January to 25.2 ◦C in August, according to the most recent 30-year data (1991–2020). The mean annual precipitation is 1308 mm, most of which falls in the summer season between June and September [18].

**Figure 1.** Locations of lysimeters in the Mt. Taehwa University Forest (TUF), Korea.

The TUF comprises 496 ha of natural deciduous forest and 300 ha of the coniferous plantation, with the remaining being composed of mixed forest. The dominant tree species in the deciduous forest include *Quercus acutissima*, *Quercus variabilis*, and *Quercus mongolica*, and the plantation forest comprises *Pinus koraiensis* and *Larix kaempferi*. The bedrock underneath the TUF is composed largely of weathered granite, and the soil texture is primarily sandy loam according to USDA classification. There has been no record of fire in the TUF over the last few decades. A detailed description of the TUF is provided by Im et al. [19].

#### *2.2. Data Collection*

Field lysimeters were used to measure the amount of water retained in the litter layer within each forest (Figure 2). The lysimeter comprises an aluminum container, a steel frame, a load cell, and a data storage device. The square litter container has a surface area of 1 m<sup>2</sup> with a permeable wire mesh bottom that allows excess water to drain easily into the underlying soil. Four load cells (BCL-5L, CAS, Korea) were embedded between the container and the steel frame to measure the container weight at 1 min intervals. The container weight, recorded in grams, was averaged over 30 measurements to calculate the variation in litter weight every 30 min.

Each lysimeter was placed in a gap area where other species or understory vegetation were completely absent and filled with the surrounding leaf litter. Three lysimeters, TBK1, TBK2, and TBK3, were located in *Q. acutissima* forest, while three, TCK1, TCK2, and TCK3, were placed in *P. koraiensis* plantation forest that was planted in 1964 (Figure 1). A description of each forest is presented in Table 1. The broadleaf litter lysimeters were set up in a 34-year-old *Q. acutissima* forest with an average tree diameter (n = 12) at breast height (DBH) of 24.1 ± 8.5 cm (mean ± std. dev) and an average tree height of 13.6 ± 3.3 m. The needle litter lysimeters were set up in a *P. koraiensis* plantation with an average DBH of 33.2 ± 4.5 cm and an average tree height of 18.4 ± 1.5 m.

**Figure 2.** Litter lysimeter (**a**) and schematic diagram of the lysimeter (**b**).

**Table 1.** Tree and leaf litter characteristics.


The physical characteristics of the leaf litter used in the lysimeter experiment are also presented in Table 1. The uneven distribution of leaf size and shape can influence the water retention capacity of the leaf litter [9,10]. Therefore, samples of litter in 10 locations near each lysimeter were collected in a zipper storage bag and immediately moved to the laboratory for analysis. Surface area parameters of the litter were extracted from scanned image data using the LeafArea R package [20]. Litter samples were put into a drying oven for 24 h to obtain the oven-dry weight. Similar to the characteristics of most *Q. acutissima*, individual leaves in the deciduous samples tended to be broad and round-tooth shaped, and they were, on average, 155.1 ± 22.0 mm long and 41.2 ± 7.5 mm wide at their widest point. The leaves in the *P. koraiensis* litter samples had a noticeably needle-shaped surface, with an average length of 84.9 ± 16.1 mm. The oven-dry weight of *Q. acutissima* litter was lower than that of *P. koreiensis* litter.

Prior to weighing the lysimeters, all lysimeter containers were filled with fallen leaves, which were obtained from random locations in the study area with as little disturbance to their accumulation as possible. Litter load in each lysimeter approximately corresponded to the litter accumulated on a 1 m × 1 m area adjacent to each lysimeter location. The undecomposed litters were carefully collected on the ground surface, avoiding leaves with obvious symptoms of pathogen or herbivore attack or with a decomposed entity. After litter collection, decomposed organic matter and twigs were manually removed prior to piling the litter in the lysimeter container.

The water retention capacity of leaf litter can be assessed per unit of dry weight. The oven-drying method is the simplest and most commonly adopted technique for determining the oven-dried weight of the litter; however, as this direct measurement requires destructive sampling, we decided not to utilize it in this study to maintain an intact litter structure as much as possible. Alternatively, the dry weight of litter samples was indirectly estimated by multiplying the air-dried litter weight by the ratio of the weights of air-dried and oven-dried litter, which were obtained from the 10 litter samples collected near each lysimeter. Table 2 shows the variation in litter weight within the forest, which may be attributed to the varying thickness and consolidation of litter accumulation.


**Table 2.** Characteristics of rainfall and leaf load at each lysimeter location.

<sup>1</sup> Weight of litter piled in the lysimeter container; <sup>2</sup> 1-h maximum intensity for a rainfall event.

Gross rainfall was measured at KoFlux towers that were located in an open field adjacent to each site (approximately 100 m from each site). The characteristics of the rainfall events between November 2015 and November 2017 are presented in Table 2. Some of the observations were missing values owing to equipment failures, power outages, or severe weather conditions. Therefore, 99 rainfall events were selected for TBK1, 70 for TBK2, and 85 for TBK3. In the coniferous forest, 91, 80, and 87 rainfall events were used for the analyses of TCK1, TCK2, and TCK3, respectively. Evaporation was derived from the eddy covariance flux measurements of the KoFlux towers. A more detailed explanation of the eddy covariance method can be found in Kang et al. [21].

#### *2.3. Estimation of Water Balance Parameters*

To understand the hydrologic behavior of leaf litter, the water retention capacity and free drainage were estimated from lysimeter measurements. The term retention capacity refers to the ability of the leaf litter surface to store water and is expressed by the maximum and minimum retention [9,10]. The maximum retention is the maximum volume of water held on the litter during a rainfall event, and this value is obtained just before the rainfall stops. The minimum retention is quantified as the amount of water stored in the leaf litter after drainage (or gravitational flow) completely ceases. This value represents the inherent ability of litter to retain water owing to surface tension and adhesive force [10,13].

The water retained in the litter layer was later depleted via evaporation after rainfall, which was characterized by a litter moisture depletion curve [10]. The depletion curve is the lower part of the falling limb of a litter moisture curve and provides the length of time required to achieve a certain weight of litter moisture, assuming no further rainfall. Litter moisture depletion is a drying process that can be expressed as a simple exponential form [10,22],

$$S\_{\mathbf{l}} = S\_{\text{mx}} e^{-\mathbf{a}t} = S\_{\text{mx}} k^{\mathbf{t}} \tag{1}$$

where *St* is the retention storage at time *t* (mm), *Smx* is the maximum retention storage after the cessation of rainfall (mm), *t* represents the elapsed time (h), and *k* is the depletion constant (=*e*−*a*).

The depletion constant *k* represents the time-dependent decline of litter moisture and determines the line of best fit through litter moisture measurements with time. The drying of litter did not occur in a uniform manner across the experiments. Thus, a representative *k* value for each litter type was obtained by averaging over all rainfall events.

As no drainage data were available in this study, drainage outflow from the litter layer was indirectly determined from the difference in litter weight over a duration of 30 min. During a rainfall event, evaporation typically occurs at a rate of 0.1 to 0.5 mm/h because the ambient air near the litter surface is closely saturated [23,24]. Therefore, evaporation loss can be disregarded in the calculation of free drainage during rainfall. During the rainfall period, the free drainage for a duration of 30 min can be calculated as follows:

$$D\_t = R\_t - \Delta S\_t,\text{ where }\Delta S\_t = S\_t - S\_{t-1} \tag{2}$$

where *Dt* and *Rt* are the amounts of water in the litter drainage and rainfall at time *t* (mm/30 min), respectively. Δ*St* implies the change in litter moisture for a duration of 30 min.

When the rain ceases, the evaporation process gradually recovers. Therefore, drainage can be estimated as follows:

$$D\_t = \Delta S\_t - E\_{t\prime} \text{ for } \Delta S\_t > E\_t \tag{3}$$

$$D\_t = 0, \qquad \text{for } \Delta S\_t \le E\_t \tag{4}$$

where *Et* is the amount of evaporative loss for a duration of 30 min (mm/30 min).

#### *2.4. Statistical Analysis*

Tests for normality must be checked before full statistical analysis can be conducted. Because the data didn't follow normal distribution due to high skewness and heterogeneous variance, differences among treatments were tested using the Kruskal–Wallis test. Bonferroni correction was also applied to minimize the alpha inflation when assessing the treatment effects. Two litter types were designed in this study to evaluate the effects of litter and rainfall on the water retention capacity with three lysimeters. Differences were considered statistically significant at *p* < 0.05. All statistical analyses were performed using the R statistical package, version 4.2.0 (22 April 2022) [25].

#### **3. Results**

#### *3.1. Water Retention Capacity of Leaf Litter*

The water retention capacity of the leaf litter was estimated from the litter lysimeter measurements and is presented in Table 3. For the entire period of measurement (Table 2), the maximum retention capacity (mean ± std. dev) of the *Q. acutissima* litter (1.69 ± 1.28 mm for TBK1, 1.46 ± 1.33 mm for TBK2, and 1.67 ± 1.57 mm for TBK3) was equivalent to 3.07 ± 2.23, 1.16 ± 1.06, and 1.23 ± 1.16 mm per unit litter mass (kg/m2), respectively. The maximum retention capacity in the *P. koraiensis* litter (1.75 ± 1.79 mm for TCK1, 1.90 ± 1.79 for TCK2, and 2.01 ± 1.99 mm for TCK3) was equivalent to 2.04 ± 2.08 mm, 4.13 ± 3.89 mm, and 3.46 ± 3.42 mm per unit litter mass (kg/m2), respectively, which was higher than that of the *Q. acutissima* leaf litter. This demonstrated that more rainfall was retained in needle litter than in broadleaf litter, regardless of rainfall amount, and this result was consistent with the finding of Li et al. [10]. Maximal water retention depends on the development of free drainage in the litter layer during rainfall events [9]. Needle litter forms an extremely dense accumulation that obstructs the dispersion and percolation of rainfall. However, in broadleaf litter such as that of *Q. acutissima*, biomat flow can be developed in the litter layer that allows water to move laterally and vertically [26].


**Table 3.** Water retention in leaf litter per rainfall event.

In contrast to the previous laboratory experiments, where gravitational water was readily drained within a short period (approximately 30 min) after rainfall cessation [9,10], this in situ study demonstrated that water continued to drain out of litter samples beyond 30 min after rainfall ended. Therefore, the lowest, nearly asymptotic water retention was considered in this study as the minimum retention [10]. The minimum retention varied, ranging from 1.35 ± 1.32 mm (TBK2) to 1.56 ± 1.47 mm (TBK3) in broadleaf litter and from 1.58 ± 1.75 mm (TCK1) to 1.97 ± 1.90 mm (TCK3) in needle litter. The average values of minimum retention were 1.66 mm and 2.72 mm per unit mass for *Q. acutissima* and *P. koraiensis* litters, respectively. No significant differences were observed in retention among the groups with different litter weights (*p* = 0.1912); however, a statistically significant difference was observed between litter types (*p* < 0.01).

According to previous studies [9,10], the water retention capacity of litter depends on litter mass, regardless of its thickness. The relationship between retention and litter mass is shown in Figure 3. The coefficients of determination for linear relationship were also presented in Figure 3, implying a measure of how well linear regression represents the measurements. As expected, both the maximum and minimum retention values increased linearly as litter mass increased. The maximum retention value increased considerably in the *P. koraiensis* litter, whereas it tended to change only slightly in the *Q. acutissima* litter. As shown in Figure 3, all lysimeter measurements presented a linear relationship between minimum retention and litter mass, with slope coefficients of 0.187 for *Q. acutissima* and 0.399 for *P. koraiensis* litter. This implies that a stronger adhesive force, which can resist the vertical movement of pore water, exists among the elements of the needle litter compared to that of the broadleaf litter. The difference in maximum water retention between litter types was small because it is controlled by rainfall intensity. Unlike the maximum retention values, the minimum retention values exhibited high variability between litter types. There were significant differences in both maximum and minimum retention storages between broadleaf and needle litters (*p* < 0.01).

**Figure 3.** Variation in water retention with litter amount.

As shown in Figure 4, the high variability in the minimum retention capacity within litter type was likely due to natural variability in rainfall intensity and duration. The minimum retention capacity as a percentage of rainfall could be as high as 80% during light rainfalls (<10 mm), and it became no more than 20% in cases where rainfall amount was >30 mm. The minimum retention capacity of broadleaf litter was higher than that of needle litter for light rainfall events (<20 mm) but could be low for heavy rainfalls. Broadleaf litter can easily capture water with large, curve-shaped surfaces when rainfall does not exceed the storage capacity [9,10]. Compared to *Q. acutissima* litter, *P. koraiensis* litter formed a relatively dense barrier layer owing to its smaller physical dimensions and lower porosity. Thus, the strong adhesion and surface tension of water molecules in the *P. koraiensis* layer caused a reduction in water movement through the litter layer for heavy rainfalls [10].

**Figure 4.** Variation in minimum water retention with rainfall and the corresponding envelope curve.

The percentage of rainfall retained by litter decreased as rainfall amount increased. The logarithmic relationship between retention capacity and rainfall for broadleaf (Equation (5)) and needle (Equation (6)) litter can be expressed as:

$$\text{Log}(S\_{\text{mm}}) = 3.461 - 0.043 \, R\_{\text{tot}} \, \text{(mm)} \left( R^2 = 0.600 \right) \tag{5}$$

and

$$\text{Log}(S\_{mn}) = 3.173 - 0.027 \, R\_{tot} \, \text{(mm)} \, \left( R^2 = 0.407 \right), \tag{6}$$

respectively, where log(*Smn*) is the natural log of the minimum water retention (%) and *Rtot* is the total amount of rainfall (mm). The derived relationships between water retention and rainfall amount were statistically significant (*p* < 0.01) and had residual standard errors of 1.159 and 0.604 for the broadleaf and needle litters, respectively.

For heavy rainfall with a long duration, the retention capacity reaches the potential value regardless of rainfall amount because the initial abstraction by litter is satisfied, and litter can be completely saturated. Figure 4 shows the upper boundary of the minimum retention value under natural conditions. The upper envelope curve (Figure 4) was fitted using the Aston curve [25]. The potential values of minimum retention per unit mass (kg/m2) can be expressed as follows,

$$S\_{mn,p} = 14.705 \left( 1 - e^{-0.081 R\_{\text{tot}}} \right) \text{ for needle filter} \tag{7}$$

$$S\_{mn,p} = 8.226 \left( 1 - e^{-0.173 R\_{tot}} \right) \text{ for broadband filter} \tag{8}$$

respectively, where *Smn*,*p* is the potential (maximum) values of minimum retention per unit mass of litter (mm/kg/m2).

The minimum retention value increased exponentially with light rainfall but approached asymptotic boundaries when rainfall exceeded 30 mm. This phenomenon was similar to the findings of Sato et al. [9], Li et al. [10], and others [27]. Light rainfall produced greater retention potential in broadleaf litter than in needle litter. The influence of rainfall

on retention capacity was not significant in heavy rainfall events, although the capacity slightly increased as rainfall increased. Figure 4 indicates that the upper envelope for minimum retention was 8.226 mm per unit mass for *Q. acutissima* litter and 14.705 mm for *P. koraiensis* litter.

#### *3.2. Litter Moisture Depletion Curve*

In this study, a moisture depletion curve was plotted for each litter type. As shown in Figure 5, a simple exponential relationship was observed in the water depletion curve of the leaf litter, implying that average *k* values for broadleaf and needle litters were 0.975 and 0.991, respectively.

**Figure 5.** Depletion curve of litter moisture after the cessation of rainfall.

The litter drainage led to the decline in litter moisture content immediately after rain stopped, and as time elapsed, water absorbed by the litter asymptotically approached the minimum retention value of the litter. Figure 5 shows that half of the water retained in the *Q. acutissima* litter was depleted within approximately 1.1 days (27 h) after the cessation of rainfall. Over 3.8 days (90 h), the remaining moisture was slowly extracted from the *Q. acutissima* litter layer, which then reached approximately 10% of its maximum retention value. Water depletion in *P. koraiensis* litter followed a similar pattern to that in the broadleaf litter, depleting 50% and 90% of its maximum retention value within 3.4 days (80 h) and 11.2 days (268 h), respectively.

#### *3.3. Water Balance Analysis for Rainfall Periods*

This study revealed an increase in litter weight due to water retention and a decrease due to evaporation loss. The water balance for rainfall periods was analyzed, and the total rainfall was partitioned into litter retention, evaporation loss, and free drainage, as presented in Table 4. Because there were no significant differences among lysimeters within litter type, the average values across three lysimeters within each litter type were calculated. Regardless of rainfall amount, approximately 83.7% of incoming rainfall drained from the *Q. acutissima* litter into the uppermost soil layer, whereas 84.5% of rainfall percolated through the *P. koraiensis* litter.


**Table 4.** Water balance analysis for rainfall periods.

<sup>1</sup> Rainfall events used for lysimeter measurements are included; <sup>2</sup> Parentheses indicate the portion of rainfall amount.

After litter drainage was complete, the remaining water retained in the litter layer was available for evaporation or retention. Evaporation losses during the depletion periods accounted for 5.2% of total rainfall in the broadleaf litter layer and 1.6% in the needle litter layer. Naturally, some retained water was not entirely lost via evaporation and contributed to the antecedent moisture of the litter for subsequent rainfall events. When lysimeter measurements were considered, the portion of total rainfall retained within the litter layer after the completion of litter drainage was 11.1% in the *Q. acutissima* litter and 14.0% in the *P. koraiensis* litter.

This was attributable to the physical properties of the litter, as reported in previous studies [10,28]. A horizontally oriented and densely compacted needle litter layer retains more water than curved and loosely compacted broadleaf litter layer. However, a relatively stronger adhesion and surface tension in needle litter can resist the evaporation of water from the litter surface and consequently enhance litter drainage, compared to broadleaf litter.

As shown in Figure 6, litter drainage varied with rainfall amount. For light rainfall of less than 10 mm, 25% and 32% of rainfall was intercepted or evaporated by the *Q. acutissima* and *P. koraiensis* litters, respectively, resulting in a decrease in the amount of rainfall reaching the mineral soil. The portion of rainfall that drained increased as the rainfall increased. These findings correspond with those of previous studies [10,23].

**Figure 6.** Variation in litter drainage with rainfall.

#### **4. Discussion**

The water balance of leaf litter was analyzed based on lysimeter measurements. In this study, total rainfall was partitioned into three components: retention, evaporation loss, and litter drainage. Water balance analysis revealed a proportionately greater amount of litter drainage due to lower water retention of leaf litter. The water retention values observed at the end of a rainfall period fell within the range of interception storage reported by Li et al. [10]. The amount of water retained is dependent on the initial moisture content of the litter prior to rainfall. Leaf litter that is very dry will retain a larger proportion of the rainfall arriving on the litter layer, and that which is close to the saturation point will retain very little additional moisture. The influence of antecedent moisture conditions on water dynamics in the litter layer was not quantitatively considered in this study.

Water retention per unit mass of leaf litter has been reported by previous studies based on laboratory experiments. Sato et al. [9] reported that the maximum water retention levels of *Lithocarpus edulis* (broadleaf litter) and *Cryptomeria japonica* (needle litter) were 1.56 mm and 1.59 mm, respectively, and the minimum retention of *L. edulis* and *C. japonica* were 1.53 mm and 0.81 mm, respectively. Li et al. [10] reported that the average maximum retention of broadleaf litter ranged from 3.96–6.56 mm for *Q. variabilis* and 5.26–6.25 mm for *Q. acutissima*; however, the values for needle litter ranged from 0.35–0.43 mm for *Abies holophylla* and 1.65–2.47 mm for *Pinus strobus*. Putuhena and Cordery [13] reported that the maximum and minimum water retention of pine litter was 1.25 mm and 0.97 mm, respectively. Compared to previous results measured in laboratories, the water retention values derived from the in situ lysimeter experiments in this study are within the same order of magnitude but exhibit wide variation due to the high variability of natural rainfall.

According to the lysimeter measurements in this study, water retention capacity was 10%–15% of the rainfall for storm periods. Other studies also published the water retention of gross rainfall or annual value for different litter types as follows: 8%–12% for mixed oak stands in India [29], 8%–16% for *Quercus petraea*, 12.1% for *Pinus patula*, 8.5% for *Eucalyptus grandis*, and 6.6% for *Acacia mearnsii* in South Africa [30]. Gerrits et al. [16] found litter interception to be as high as 22% in a beech forest and 18% in needle leaf litter in a cedar forest, while Helvey and Patric [11] found litter interception to be 15%–34% in a poplar stand in the USA. Comparison of the current study with these past studies is limited due to differences in climate conditions, tree species, litter mass, and methods of measurement.

Litter drainage was indirectly calculated as the difference between the amount of rainfall and the amount of water retained in the leaf litter. This concept is valid only if the evaporation loss is neglected for shorter durations. Even if this condition is not satisfied, it leads to a small discrepancy in almost all instances because the evaporation rate is too small to have a substantial impact on the drainage calculation for a very short period [24,31].

The gravimetric method is the most widely used technique for the moisture determination of litter [6,7], where litter samples are transported to the laboratory for experimental measurements. Although this is a very simple and accurate method, it disturbs the samples and requires 12–24 h to dry in an oven to weigh the litter sample under dry conditions. Therefore, the gravimetric method is not applicable to continuously measure the water partitioning in field conditions. A field lysimeter is the best solution to quantify litter interception and retention under natural conditions [12]. Because the experiment is executed in the original environment, litter samples are less disturbed compared to the gravimetric method.

Naturally, some biological and pathological changes could have occurred during the measurement period [32,33]. Decomposition can modify the physical characteristics of leaf litter, resulting in a change in its weight [34], which would influence the hydrologic processes in the litter by reducing its water retention capacity and resisting water penetration into the layer. In this study, litter decomposition was not considered because it was the intent of the experiment to obtain moisture measurements from relatively undisturbed litter samples. Litter-decomposing fungi are also known to have higher water resistance

(hydrophobicity) due to fungal mycelia [33], although they have rarely been observed in litter samples.

#### **5. Conclusions**

The water retention capacity of leaf litter was estimated from lysimeter measurements taken over a duration of 30 min. Regardless of rainfall characteristics, needle litter had a greater capacity to retain rainfall than broadleaf litter did. During the experiment, approximately 1.66 ± 1.72 mm per rainfall event was stored per unit mass (kg/m2) of broadleaf litter after rainfall, while 2.72 ± 2.82 mm of rainfall reaching needle litter was retained per unit mass (kg/m2). Approximately 83.7% of incident rainfall drained into the uppermost soil layer below the *Q. acutissima* litter, whereas 84.5% of rainfall percolated through the *P. koraiensis* litter. The depletion curves indicated that the water retained in the *Q. acutissima* litter was more easily lost via evaporation than the *P. koraiensis* litter.

The duration and intensity of rainfall are known to affect litter water retention capacity. Most previous studies have been conducted under limited conditions, where the short duration and/or high intensity of artificial rainfall might lead to an underestimation of the retention capacity of litter, likely due to the litter accumulation not being thoroughly wetted prior to the simulated rainfall. On the contrary, as this study was conducted under natural rainfall conditions, it resulted in a higher variation in retention capacity, likely due to the variation in rainfall. Practically, the upper envelope curves derived from lysimetry experiments can provide the upper boundary of the retention capacity, reflecting the effect of rainfall characteristics under natural conditions.

This lysimetry experiment has some limitations. First, evaporation from the litter layer has not been measured but can be derived in this study. It would require intensive field experiments where hydrological components such as rainfall, evaporation, percolation, and drainage are measured simultaneously. Second, litter decomposition as a result of fungal and soil faunal activity also occurs over time. Changes in the physical traits of leaf litter can affect hydrological function; however, this phenomenon was not considered in the current study. Therefore, long-term measurements are required to accurately understand the hydrological role of leaf litter, considering spatial and temporal variations in production, accumulation, and decomposition.

This study highlighted the litter retention and evaporation processes. However, relatively little is known on water dynamics occurring on the forest floor. Biomat flow likely follows preferential flow paths through the litter layer. Further research can be directed to the clear understanding of how biomat flow influences the residence time of water within the layer, time to litter saturation, and, in consequence, the water retention capacity of litter.

**Author Contributions:** Conceptualization, T.K. and S.I.; methodology, T.K. and S.I.; software, J.K. and J.L.; validation, T.K., J.K. and J.L.; investigation, J.K., J.P. and J.L.; resources, H.S.K.; data curation, T.K. and J.P.; writing—original draft preparation, T.K.; writing—review and editing, S.I. and H.S.K.; visualization, J.K. and J.L.; supervision, S.I.; project administration, S.I.; funding acquisition, S.I. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was carried out with the support of 'R&D Program for Forest Science Technology (Project No. 2021343C10-2323-CD01)' provided by Korea Forest Service (Korea Forestry Promotion Institute).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We acknowledge the Forest Ecophysiology Laboratory of Seoul National University, Korea, for providing the experimental data, and also thankful for the support of the Mt. Taehwa Seoul National University Forest.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Victoria Virano-Riquelme 1,\*, Karl-Heinz Feger <sup>1</sup> and Stefan Julich 1,2**


**Abstract:** The structure of forests in temperate climates has been changing to ensure the resilience of trees. This change affects the local water balance. Knowledge of soil hydraulic properties (SHP) is essential to assess the water cycle in ecosystems. There is little knowledge about the impact of tree species on SHP and the water balance. Based on a compilation of 539 related studies we aimed at identifying the effects of tree species and age on SHP in temperate climates. However, most studies concentrated on soil biogeochemical properties, whereas only 256 studies focused on SHP. The literature presents no standard methods for assessing SHP and there is no knowledge of their variations in forests. We present a systematic overview of the current state of knowledge on variations in SHP based on forest type in temperate climates. We identify the gaps and weaknesses in the literature and the difficulties of evaluating the reviewed studies. More studies following standardised methodologies are needed to create a robust database for each forest type and soil texture. It would improve the assessment of the forest water balance through calibrated plot/site-scale process models. Such a database does not yet exist, but it would greatly improve the management and development of future forest ecosystems.

**Keywords:** soil hydraulic properties; saturated hydraulic conductivity; bulk density; temperate forests

#### **1. Introduction**

In recent decades, forest ecosystems have shown a decline in productivity, exhibiting tree species with negative physiological processes undergoing high and abrupt mortality [1]. According to Allen et al. [2], the increasing tree death rate is caused by climate change with concomitant drought and heat as the main factors compromising the prosperity of the forests. The severity of droughts and the high stress in temperate forest ecosystems create conditions favouring pests and pathogens [3], creating an environment of high risk, especially in areas with low water availability [4].

In addition, future scenarios foresee increasing drought frequency, intensity, and duration as well as additional impacts such as wildfires and severe storms, leading to elevated tree mortality rates [1,2,4]. Long-term droughts, coupled with rapid tree decline, represent considerable challenges to the management and policy-making communities. In response to changing conditions, forests of the temperate zone are transforming, being managed to withstand these adverse constraints.

In the past, monocultures of needle-leaved trees (i.e., of Norway spruce—*Picea abies*) were widely used in Central Europe as a promising economic resource. However, in the context of climate change, spruce monocultural stands have weakened and are rapidly dying out. Therefore, current forest management considers the use of native broad-leaved tree species as mitigation to endure drought stress and preserve the vitality of trees. As reported by Pretzsch et al. [5], the process of transformation and adaptation from non-natural

**Citation:** Virano-Riquelme, V.; Feger, K.-H.; Julich, S. Variation in Hydraulic Properties of Forest Soils in Temperate Climate Zones. *Forests* **2022**, *13*, 1850. https://doi.org/ 10.3390/f13111850

Received: 1 September 2022 Accepted: 29 October 2022 Published: 4 November 2022

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

**Copyright:** © 2022 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/).

and widespread coniferous monocultures towards mixed forests including deciduous tree species will be of great ecological interest considering different rates and species.

This strategy plays a key role in overcoming the current climatic conditions and in the development of future forest ecosystems. The introduction of broad-leaved and hardwood species in needle-leaved stands should be performed carefully as it changes several soil properties (i.e., humus, litter deposition), likely affecting the soil hydraulic properties [6–8]. Therefore, the assessment of possible consequences on the water balance is essential for efficient forest productivity [9].

The evaluation of the water balance of forest stands in terms of different species can be achieved by using calibrated plot/site-scale process models such as LWF-Brook90 [10,11], which depend on soil hydraulic properties (SHP). SHP are essential to assess the water cycle in ecosystems. The evaluation of these parameters allows the description of several hydrological processes within the soil profile, such as water infiltration, surface runoff and water storage, the availability of water for root uptake, percolation, and the generation of surface runoff as a function of the soil management practices [12,13]. Furthermore, the soil hydraulic conductivity is particularly relevant since it allows the estimation of water movement towards the root zone and the recharge of underground aquifers, as well as the transport of minerals and nutrients into the deeper layers of the soil profile [14,15].

Some studies compare SHP between forested and non-forested sites, but only a few focus on the different management practices within forest ecosystems. In contrast to other land uses, studies reveal that forest soils usually show higher hydraulic conductivity [14,16,17], improved water retention [17], decreased runoff and increased soil moisture [17,18], with often lower density [19]. Moreover, the infiltration capacity and soil water retention in forested soils control the formation of surface runoff acting as a natural flood regulator [16,20,21]. However, most studies addressing the variation in soil properties as a result of forest management focus on aboveground characteristics including biomass [22–24], vegetation traits [25–27], and soil biochemical properties [28–32].

Articles assessing SHP in temperate forests are extremely scarce and present a great variety in findings. Differing values for soil bulk density, hydraulic conductivity, and the accumulation of soil organic carbon in relation to forest type and tree age have been observed [8,14,21,33–35]. In addition, the relief of the sites has a noticeable implication for the assessment of the SHP. On steep slopes, needle-leaved forests have revealed enhanced hydraulic conductivity in contrast to smoother surfaces [36]. However, under similar conditions, broad-leaved forests create conditions less prone to surface runoff generation compared to needle-leaved forests [34]. SHP may vary according to a vast combination of soil texture, tree species, root distribution, and the inherent complexity of the interacting site conditions, likely affecting the water dynamics in a catchment [6,16,33,37,38].

The objective of this study is to review the available literature on information for the SHP of forested sites in temperate zones. Our aim is to evaluate and assess the impact of tree species and age on soil hydraulic properties based on the reported values. We expect to find differences in SHP in terms of broad-leaved and needle-leaved tree species as reported by Archer et al. [14], Julich et al. [16], and Wahren et al. [17]. We hypothesize that in needle-leaved stands (at least for comparable soil geological settings), the overall hydraulic conductivity is lower, together with higher bulk density, compared to broad-leaved stands.

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

#### *2.1. Nature and Classification of the Studies*

For the review, a systematic search was conducted based on the PRISMA Statement [39]. The academic databases Scopus and Web of Science were used to search for information, with the last search performed on 24 March 2022. This review targeted peerreviewed publications in English until the year 2021. Publications in other languages were excluded as a thorough analysis requires an adequate understanding of the study site descriptions.

Figure 1 shows the search steps we followed to create our database. Firstly, we performed a broad search ("Search 1", Figure 1), considering the terms "soil hydraulic properties" and "forests", obtaining a total of 462 articles. The results were classified according to "temperate forests" and "tropical forests (+others)", with the latter group presenting the highest number of studies regarding soil properties. In this step the articles that were not related to the study of soil properties (i.e., soil fauna, soil fungi, root development), that were not specifically carried out in forested areas, that did not specify the type of climate, or articles where it was not clear, were excluded. Then, the focus of our study was specifically directed to temperate forests as our first research revealed that the studies conducted on these forests were deficient.

**Figure 1.** The flowchart on the left displays the steps performed in a first broad search, which shows the current state of the literature in relation to the variation in soil hydraulic properties in different climatic zones. The flow chart on the right illustrates the steps followed to filter the studies related to soil hydraulic properties in temperate forests (the focus) and how the final database based on current literature was compiled. The number of papers found (n) and the studies that were discarded at each step along the literature review are included. (\*) Some articles are repeated between the subdivisions.

After the first literature review, we carried out a complementary search using the same academic databases to aim for studies addressing SHP in temperate forests, including the words "hydraulic conductivity", "water retention", and "water content" to make certain of covering relevant soil hydraulic properties ("Search 2", Figure 1). For each of the beforementioned words, different combinations of keywords such as "temperate forest", "temperate climate", "coniferous", "deciduous", "stand", "species", "mixed forest(s)", and "forest conversion" were used.

Although we collected a total of 539 articles up to 2021, several of these studies referred to the term "soil properties" in a general context. Therefore, 267 studies were discarded before the assessment as they focused specifically on soil biogeochemical properties, traits, and root water uptake, rather than on soil hydraulic properties. Moreover, some of those studies were discarded since they concerned other land uses. After neglecting the studies that were not within the scope of this review, the compilation consisted of 272 studies. Since our analysis aimed at soil hydraulic properties obtained between 0 and 30 cm depth, 16 additional studies not reporting depth of measurement or performed at greater depths were excluded from the database. In addition, studies that used units deviating from the standard were discarded to avoid confusion and misinterpretation over converted values. Once all these steps were completed, the final database was subdivided into those articles containing studies related to soil bulk density, hydraulic conductivity, and soil water content.

To assess the effect of the different forest ecosystems on the soil properties, we defined a classification for this study according to needle-leaved, broad-leaved, and mixed forests. Within each forest group, we attempted to classify the studies in terms of site parameters such as soil texture, organic matter concentration, tree age, and tree mixing rate (in the case of mixed forests). Still, due to the limited information, we were able to classify the articles only in terms of soil texture and tree age. The classification of soil texture was performed based on fine for clay-dominated soils, medium for mostly loamy soils, and coarse for mostly sandy soils. Given the range in ages of the trees in the reviewed studies, and for the purpose of our study, we decided to group the values collected according to the following criteria: young forest for trees < 50 years old, mature stand for those containing trees of 50–100 years old, and old forest for trees older than 100 years.

When evaluating the articles that assess soil water content, we observed that most of the studies were concerned with variations in water content over a given period of time. Such results are of interest for assessing temporal fluctuations in soil moisture content; however, they are relative to seasonal and climatic factors. Therefore, we focused on measurements of water holding capacity, as this property is relevant to assess the effect of forest type. From the 92 articles collected that addressed water content, only 15 could be used for our research.

In the case of saturated hydraulic conductivity, 3 out of 22 articles performed their measurements under laboratory conditions, which were excluded since such measurements likely neglect the site characteristics. Regarding the experiments performed in the field, one article did not specify the measurement device, whereas 10 studies used permeameters and nine used infiltrometers. Although both methods work under different techniques, we considered all the 20 articles in our analysis. From these studies, we compiled a total of 145 values for analysis. Regarding soil bulk density, from a total of 56 studies, we collected 145 values to be further assessed.

#### *2.2. Statistical Analysis*

As most of the results came from diverse studies and were performed under specific conditions, the collected data could not be assumed to be uniformly distributed. Even when some of the reviewed papers presented comparable site parameters, they frequently differed in other aspects, such as measurement approaches, devices, and sampling designs. Therefore, the Kruskal–Wallis test was used to assess whether the data presented significant differences between the forest types under the assumption that all data points were obtained independently under the same experimental conditions.

#### **3. Results**

The compilation of studies was screened for values of SHP in terms of water holding capacity, saturated hydraulic conductivity, and soil porosity. However, the available information on the aforementioned parameters was limited. Thus, it was not possible to add most of these hydraulic parameters to the analysis. Only for soil hydraulic conductivity could sufficient information be extracted for further assessment. Since soil bulk density was present in almost all the reviewed articles, it was included in the analysis, as it provides relevant information regarding soil structure.

In order to appropriately compare each soil hydraulic property in terms of comparable site parameters, we classified the data according to forest type (i.e., tree species and age of the trees) and the soil texture where the parameters were measured. We found that not all the collected articles reported these relevant parameters, hence several studies were unsuitable for further analysis. As a result, we created two datasets containing all values found in the pool of collected articles. Database A contains all results obtained for a soil depth of 0–30 cm (independent of the site description) and database B contains the filtered information of database A, comprising only those articles reporting forest type and soil texture. Furthermore, all values were converted to the same units for further comparison.

Table 1 presents the number of studies included in databases A and B, the latter containing the results that are suitable to assess the collected bulk density and hydraulic conductivity values. To visualize the information, the distribution was performed by forest type, and each field was additionally sorted by soil texture and tree age. It can be seen that most studies comprised information regarding bulk density and fewer studies measured hydraulic conductivity. It is also evident that the majority of the studies were conducted in needle-leaved and broad-leaved forests.

**Table 1.** Classification of values found in the reviewed literature filtered by forest type as main field and further divisions by soil texture and tree age, which were defined by this study. The left side of the table (a) corresponds to the hydraulic conductivity (K) and right side (b) refers to the soil bulk density (BD). All values are measurements between 0 and 30 cm soil depth.


#### *3.1. Soil Water Retention*

Table 2 summarises the information regarding soil water retention characteristics found in the screened articles. From 15 studies a total of 30 values were found, yet not all studies referred to the same measurements of soil water holding capacity, resulting in limited information for the individual soil water types. Field capacity (FC) ranged from 0.12 m<sup>3</sup> m−<sup>3</sup> for sandy soils up to 0.44 m3 m−<sup>3</sup> for silty clay soils. Regarding plant available water (PAW), only a little amount of information was available with values ranging from 0.02 to 0.33 m<sup>3</sup> m−3. For water content at the permanent wilting point (PWP), values between 0.05 and 0.20 m3 m−<sup>3</sup> were reported.

#### *3.2. Soil Bulk Density*

Figure 2 compares the values for soil bulk density of the topsoil by using dataset A, which contains all values found in the literature screening, and dataset B containing those that were filtered regarding the available information on tree age and soil texture. Dataset A comprised a total of 163 bulk density values, whereas dataset B was reduced to 79, excluding about 52% of the information. For dataset A, the highest and lowest means were observed in the broad-leaved and needle-leaved forests, respectively, whereas dataset B showed comparable mean values between forest types (Figure 2).


**Table 2.** Detailed information on soil water holding capacity from the reviewed literature. In the division Forest Type, the abbreviations Nee, Broa, and Mix stand for needle-leaved, broad-leaved, and mixed forest type, respectively.

<sup>1</sup> Soil water at field capacity. <sup>2</sup> Soil plant available water. <sup>3</sup> Soil water at permanent wilting point.

Table 3 shows the number of studies and descriptive statistics of bulk density for datasets A and B. The *p*-values suggest that in dataset A there is an effect of forest type on the soil bulk density (*p*-values < 0.05); however, the *p*-values in dataset B show no effect of forest type on this soil property. The comparison with mixed forest is not adequate in both cases since the number of studies on this forest type is not sufficient in contrast to the others (Table 3).

**Figure 2.** Descriptive statistics of top mineral soil bulk density by forest type considering the two data sets: dataset A containing all values found and dataset B filtered values based on available information on tree age and the soil particle size. All values refer to the topsoil (0–30 cm depth). Black dots represent the outliers.

**Table 3.** Summary of statistics for soil bulk density and saturated hydraulic conductivity sorted according to forest type. Dataset A contains all results found in the screened articles; dataset B considers only those values reported with additional information on soil particle size and tree age. All values refer to the topsoil (0–30 cm depth). G Mean refers to the geometric mean. (\*) indicates the group of forest with significantly different BD values.


Figure 3 and Table 4 provide the descriptive statistics of dataset B regarding the bulk density of the topsoil. Most of the reviewed studies with information regarding tree age were conducted on broad-leaved forests, whereas for mixed forests only a few studies were carried out. Values for bulk density do not follow any trend (in terms of forest type or tree age) and reveal that mature and young stands are the dominant stages in the reviewed studies.

**Figure 3.** Descriptive statistics of the top mineral soil bulk density in relation to the tree age and forest type for dataset B. Scattered points represent the individual values classified by the different forest stages. Values correspond to the 0–30 cm soil depth. Grey points indicate the outliers.

**Table 4.** Summary of reported values on topsoil bulk density and hydraulic conductivity for the different forest types, data are subdivided based on the information on tree age: young (Y) less than 50 years, mature (M) 50 to 100 years, and old (O) more than 100 years.


The topsoil bulk density for needle-leaved and broad-leaved forests presented the highest values in the mature stage; however, comparison with mixed forests was not possible since values for the young stage were not found. Most values for the old stages under needle-leaved and broad-leaved forests were lower in contrast to the other stages.

Figure 4 and Table 5 describe the reported values for bulk density in relation to soil texture and forest type. Most of the screened studies with information on soil texture were conducted in broad-leaved forests, whereas only a few were conducted in mixed stands. The distribution of values reveals no apparent relationship between bulk density and soil texture for the different forest types. Figure 4 also shows that very few reported values correspond to the Oh-horizon in contrast to those measured in the first 30 cm of soil.

**Figure 4.** Descriptive statistics of the top mineral soil bulk density for dataset B in terms of forest type and soil texture: fine (clay-dominated), medium (mostly loamy), and coarse (mostly sandy). Scattered points represent the individual values between 0 and 30 cm soil depth, including the Oh-horizon. Grey points indicate the outliers.

**Table 5.** Topsoil bulk density and hydraulic conductivity values in terms of forest type, filtered by soil texture. Letters C, M, and F stand for soils with coarse (mostly sandy), medium (mostly loamy), and fine (clay-dominated) soil texture, respectively, between 0 and 30 cm soil depth. In the case of bulk density, the letter Oh represents the values for the Oh-horizon.


For needle-leaved and broad-leaved forests, the highest mean values for bulk density were observed in forest stands with a loam-dominated soil texture, whereas the lowest were found in sites with a fine soil texture. Overall bulk density in the Oh-horizon were in the order: broad-leaved > mixed > needle-leaved, with a limited number of values reported in the screened articles (Table 5). Moreover, it was often difficult to identify whether the soil density values corresponded to the mineral layer or the organic layer. Therefore, it is possible that there is a bias in the values reported in our study due to a lack of detail in terms of the depth of measurements.

#### *3.3. Soil Saturated Hydraulic Conductivity*

Figure 5 and Tables 3–5 summarise all the information regarding the topsoil saturated hydraulic conductivity obtained from the reviewed articles. Analogous to the bulk density, dataset A comprises 59 values, while dataset B contains 46 values. Between the two datasets, no change in the mean values in terms of forest type was noticeable: mixed > broad-leaved > needle-leaved (Table 3, Figure 5). Furthermore, no statistical effect of forest type on the soil hydraulic conductivity of the topsoil was observed (*p*-values > 0.05, Table 3).

**Figure 5.** Descriptive statistics of topsoil hydraulic conductivity in terms of forest type considering two datasets: dataset A containing all values found and dataset B containing filtered values based on available information on tree age and soil particle size. All values refer to the topsoil (0–30 cm depth). Black dots represent the outliers.

For needle-leaved forests, most studies concentrated on the young stage, whereas the old stages of broad-leaved forests dominated the screened studies (Table 4, Figure 6). For mixed forests, there were no data on hydraulic conductivity at the young stage, while most studies occurred in the stands of the mature age (Figure 6). In the case of needleleaved forests, the order of the mean values of hydraulic conductivity recorded was: mature > old > young, differing from that under broad-leaved forests: young > old > mature. However, the number of studies in the mature and old stages of needle-leaved forests was extremely low and therefore did not allow any statistical comparison.

Figure 7 and Table 5 describe the statistical distribution for dataset B with regard to the saturated hydraulic conductivity in relation to forest type and soil texture. Most of the reviewed studies with information on soil texture were conducted in broad-leaved forests and study sites were dominated by a coarse soil texture.

Under needle-leaved and broad-leaved forests, the order of magnitude of hydraulic conductivity was comparable between the coarse and medium soil textures, whereas soil dominated by fine texture presented higher values (Table 5). It has to be noted that the number of studies carried out in all forest types with regard to hydraulic conductivity in stands with fine-textured soils was very low.

**Figure 6.** Descriptive statistics of top mineral soil hydraulic conductivity for dataset B in relation to tree age and forest type. Scattered points represent the individual values classified by the different forest stages considered in this study. Values correspond to the 0–30 cm soil depth. Grey points represent the outliers.

**Figure 7.** Descriptive statistics for dataset B top mineral soil hydraulic conductivity in terms of forest type and soil texture: fine (clay-dominated), medium (mostly loamy), and coarse (mostly sandy). Grey points represent the outliers.

#### **4. Discussion**

Our literature review revealed that with respect to water retention characteristics of soils under forests, only limited information is available, with fewer than 20 studies reporting values (Table 2). Thus, the reported information does not allow any conclusion on the influence of tree species and age on water retention characteristics (which are part of SHP) for different site conditions (i.e., soil texture, tree species). Moreover, the reported studies differ in the applied methodology and experimental design, which allows only limited comparability among the studies. Dedicated studies such as Wahren et al. [17] and Julich et al. [16] indicate that for sites influenced by stagnant soil water conditions, differences in SHP for stands of different tree types and age can be expected. However, the transferability of those findings needs to be tested by more experiments on sites with different site conditions.

With respect to bulk density as an important indicator of soil structure, more information in the screened literature is available (Tables 1 and 3–5). Bulk density varies together with the changes in forest structure. Bakhshandeh-Navroud et al. [51] observed a soil bulk density of 1.33 g cm−<sup>3</sup> in a broad-leaved monoculture, whereas a higher value was observed when eight different broad-leaved tree species were mixed (1.48 g cm−3). Jost et al. [37] found a higher soil bulk density in a broad-leaved forest in contrast to a needle-leaved forest (1.13 vs. 0.97 g cm−3). Bło ´nska et al. [33] observed higher values in a broad-leaved forest (range: 0.35–0.81 g cm−3) compared to needle-leaved stands (range: 0.21–0.28 g cm−3), and in a mixed forest (0.36 g cm<sup>−</sup>3); however, it has to be considered that in this study, the organic layer instead of the top mineral layer was sampled. This emphasizes the need for a clear description of the sampled depth and horizon. Moragues-Saitua et al. [52] reported a higher bulk density under needle-leaved forests (1.06 g cm<sup>−</sup>3) and in broad-leaved forests (0.67 g cm−3) in contrast to the beforementioned studies. The results of the presented literature review suggest an effect of forest type (tree species) on soil bulk density (Table 3a, *p*-value < 0.05). Moreover, it is not only soil texture that has been described in the literature as a factor that might explain variations in soil bulk density but tree age as well. Polláková et al. [53] observed that bulk density was higher (1.15 and 1.01 g cm<sup>−</sup>3) in mature needle-leaved forests with comparable characteristics (110 years old, sandy soil), in contrast to old forests. Furthermore, they noticed that bulk density was higher in a silty soil young forest (50 years old, 1.31 g cm<sup>−</sup>3). Rosenkranz et al. [54] reported that a 100-year-old mature needle-leaved forest presented a soil bulk density of 0.82 g cm−3, but when mixed with a three-year-old broad-leaved tree species, it was reduced to 0.64 g cm<sup>−</sup>3. Baldocchi et al. [43] found that soil bulk density was 1.05 g cm−<sup>3</sup> in a 250-year-old needle-leaved forest with sandy soil, while Link et al. [55] reported lower values (0.8 g cm−3) in an older forest (450 years old) containing needle-leaved tree species. Frêne et al. [29] observed a higher bulk density in a 110-year-old needle-leaved forest in contrast to an old (350/400 years old) mixed forest (0.77 vs. 0.44 g cm<sup>−</sup>3, respectively); however, information on soil texture was not provided.

Regarding soil hydraulic conductivity, our literature review suggests that there is no apparent effect of forest type on soil hydraulic conductivity (Table 3). However, there are studies concluding that this parameter not only varies in terms of tree species, but also that site characteristics play a relevant role. Bauwe et al. [22] reported lower infiltration capacity under needle-leaved stands (349 cm d<sup>−</sup>1) in contrast to broad-leaved forests (both with mature tree stands). Bens et al. [21] observed the highest hydraulic conductivity in mixed forests (76-year-old pine with 34-year-old beech) with 1045 cm d−1, whereas needle-leaved and broad-leaved forests presented 57% and 38% lower infiltration rates, respectively. Similarly, Wahl et al. [7,8] reported the highest infiltration rates in mixed forests and the lowest in needle-leaved forests, both presenting comparable sandy soil. Julich et al. [16] reported a higher amount of well-connected macropores for a 170-year-old broad-leaved forest, and higher infiltration capacities, in contrast to a young 25-yearold needle-leaved forest, probably influenced by their deep-developed rooting system. Moreover, Frey et al. [56] found higher hydraulic conductivity in a mixed forest with loamy soil (320 cm d<sup>−</sup>1) in contrast to a mixed forest with loamy soil but containing a much higher clay content (63 cm d<sup>−</sup>1).

As reported by these authors, soil hydraulic conductivity is sensitive to other soil properties; still, many authors describe this parameter under varying conditions. For instance, Archer et al. [20] considered a wide range of parameters to assess the water infiltration rates such as the soil particle size, tree ages, and measurements at different depths. However, results were not able to be derived on the influence of such parameters since all forests differed in age and soil texture. They considered two mature stages (160/180 years old) and one old (500 years old) stage with broad-leaved stands and a young stage (~50 years old) with needle-leaved tree species. The sites differed in soil texture. In young needle-leaved forests the range was 82–101 cm d−<sup>1</sup> (mean 90.9 cm d−1), in mature broad-leaved sites the ranges were 12–17 cm d−<sup>1</sup> (mean 14.3 cm d<sup>−</sup>1) and 245–312 cm d−<sup>1</sup> (mean 280.2 cm d<sup>−</sup>1), and the range in the old broad-leaved site was 235–288 cm d−<sup>1</sup> (mean: 258.3 cm d−1). Similarly, Bens et al. [21] and Wahl et al. [7,8] reported higher values (648, 458, and 751 cm d−1, respectively, by author) in broad-leaved forests in contrast to needleleaved stands (451, 251, and 538 cm d−1, respectively, by author). However, the study considered a young (<40 years old) needle-leaved and a mature (91 years old) broad-leaved forest, which also differed in terms of soil texture.

#### **5. Conclusions**

The various soils under forests present an extensive assortment of site conditions: i.e., soil texture, porosity, biotic activity (in particular that of the meso- and macrofauna), amount and vertical distribution of soil organic matter, and tree species and ages, among many others. Nevertheless, in our literature review we found considerable limitations in assessing the effect of forest type on soil hydraulic properties, which led us to the following conclusions:


We suggest that the effect of forest type on soil water dynamics should be assessed based on studies with comparable site parameters to avoid overlooking the high sensitivity of SHP. For this reason, a much higher number of systematic studies following standardised methods is clearly needed. Such studies should include a minimum of information on tree age, mixture rate, depth of measurement, and soil texture in order to make future studies comparable and to create an understanding of the influence of forest type and forest structure on SHP under different site conditions.

**Author Contributions:** Conceptualisation S.J. and V.V.-R.; literature review, data analysis and writingoriginal draft preparation V.V.-R.; writing-review and editing S.J. and K.-H.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the scholarship programme of the Graduate Academy for the Promotion of Early-Career Female Scientists at the TU Dresden.

**Data Availability Statement:** The author, upon request, can provide information regarding the reviewed papers.

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

#### **References**

