# Response of Downy Oak (Quercus pubescens Willd.) to Climate Change: Transcriptome Assembly, Differential Gene Analysis and Targeted Metabolomics

J.P. Mevy, Béatrice B. Loriod, Xi Liu, Erwan Corre, Magali Torres, Michael

Büttner, Anne Haguenauer, Ilja M. Reiter, Catherine Fernandez, Thierry

Gauquelin

<figure><img src="image_1.png" /></figure>

## To cite this version:

J.P. Mevy, Béatrice B. Loriod, Xi Liu, Erwan Corre, Magali Torres, et al., Response of Downy Oak (Quercus pubescens Willd.) to Climate Change: Transcriptome Assembly, Differential Gene Analysis and Targeted Metabolomics. Plants, 2020, 9 (9), pp.1149. 10.3390/plants9091149. [http://dx.doi.org/10.3390/plants9091149](http://dx.doi.org/10.3390/plants9091149)

HAL Id: hal-02945417

https://amu.hal.science/hal-02945417

Submitted on 22 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Article

# Response of Downy Oak (Quercus pubescens Willd.) to Climate Change: Transcriptome Assembly, Differential Gene Analysis and Targeted Metabolomics

Jean-Philippe Mevy ¹,*, Beatrice Loriod ², Xi Liu ³, Erwan Corre ³, Magali Torres ², Michael Büttner ⁴, Anne Haguenauer ¹, Ilja Marco Reiter ⁵, Catherine Fernandez ¹ and Thierry Gauquelin ¹

¹ CNRS, Aix-Marseille University, Avignon University, IRD, IMBE, 13331 Marseille, France; anne.haguenauer@criobe.pf (A.H.); catherine.fernandez@imbe.fr (C.F.); thierry.gauquelin@imbe.fr (T.G.)

² TGML-TAGC—Inserm UMR1090 Aix-Marseille Université 163 avenue de Luminy, 13288 Marseille, France; beatrice.loriod@inserm.fr (B.L.); magali.torres@univ-amu.fr (M.T.)

³ CNRS, Sorbonne Université, FR2424, ABiMS platform, Station Biologique, 29680 Roscoff, France; xi.liu@sb-ro recoff.fr (X.L.); erwan.corre@sb-ro recoff.fr (E.C.)

⁴ Metabolomics Core Technology Platform Ruprecht-Karls-University Heidelberg Centre for Organismal Studies (COS) Im Neuenheimer Feld 360, 69120 Heidelberg, Germany; michael.buettner@cos.uni-heidelberg.de

⁵ Research Federation ECCOREV FR3098, CNRS, 13545 Aix-en-Provence, France; ilja.reiter@web.de

* Correspondence: jean-philippe.mevy@univ-amu.fr; Tel.: +33-0413550766

Received: 20 July 2020; Accepted: 1 September 2020; Published: 4 September 2020

<figure><img src="image_4.png" /></figure>

check for
updates

**Abstract:** Global change scenarios in the Mediterranean basin predict a precipitation reduction within the coming hundred years. Therefore, increased drought will affect forests both in terms of adaptive ecology and ecosystemic services. However, how vegetation might adapt to drought is poorly understood. In this report, four years of climate change was simulated by excluding 35% of precipitation above a downy oak forest. RNASeq data allowed us to assemble a genome-guided transcriptome. This led to the identification of differentially expressed features, which was supported by the characterization of target metabolites using a metabolomics approach. We provided 2.5 Tb of RNASeq data and the assembly of the first genome guided transcriptome of *Quercus pubescens*. Up to 5724 differentially expressed transcripts were obtained; 42 involved in plant response to drought. Transcript set enrichment analysis showed that drought induces an increase in oxidative pressure that is mitigated by the upregulation of ubiquitin-like protein protease, ferrochelatase, oxaloacetate decarboxylase and oxo-acid-lyase activities. Furthermore, the downregulation of auxin biosynthesis and transport, carbohydrate storage metabolism were observed as well as the concomitant accumulation of metabolites, such as oxalic acid, malate and isocitrate. Our data suggest that early metabolic changes in the resistance of *Q. pubescens* to drought involve a tricarboxylic acid (TCA) cycle shunt through the glyoxylate pathway, galactose metabolism by reducing carbohydrate storage and increased proteolytic activity.

**Keywords:** transcriptome analysis; metabolism; drought; *Quercus pubescens*; adaptation; RNASeq

## 1. Introduction

Global change refers to a wide range of global phenomena such as warming, increased sea levels and extreme weather events. This includes significant changes in either the mean state of the climate or in its variability, persisting for an extended period as defined by the United Nations Intergovernmental
Panel on Climate Change [1]. Different biomes are impacted in different ways, but the consequences in terms of ecosystem functioning and ecosystem services are far from being well understood.

The Mediterranean basin stands at the crossroads of the influences of the arid climate of northern Africa and the temperate climate of central Europe, and this raises the question of the vulnerability of this region to climate change variables. The region accounts for about 25,000 autochthon plant species, among which more than half are endemic. The Mediterranean basin was thus recognized as one of the first 25 global biodiversity hotspots [2]. At the same time, annual rain fall is projected to decline by 20–30% by the end of this century [3,4] which confirms the urgency of setting up conservation strategies in relation to the vulnerability of these ecosystems to drought.

The limits of the distribution areas of the Mediterranean basin will be modified by climate change depending on the drought tolerance range of species. The most sensitive species are therefore those with specific tolerance thresholds for precipitation and other factors. However, little is known about how these species will cope with the predicted climate change.

The genus *Quercus* (Fagaceae, family) is one of the most important taxa of woody angiosperms in the northern hemisphere. In the Mediterranean basin, *Quercus pubescens*, or downy oak, is often found in coastal plains to about 1300 m altitude. It is a heliophilous and thermophilous species and is well adapted to moderate summer drought stress, although it avoids continental locations subjected to more frequent drought events.

Several studies have been carried out based on the tolerance of *Q. pubescens* to drought ranging from morphoanatomy to physiology. For instance, in drought-stressed individuals, a significant earlier senescence and increased amounts of soluble sugars in the leaves were observed [5]. This high sugar content is in accordance with the capacity of *Q. pubescens* to maintain much higher photosynthetic potential when compared to other species, such as *Quercus ilex*, *Pinus halepensis* and *Arbutus unedo* [6]. The capacity of *Q. pubescens* to control the antioxidant system response to drought was also demonstrated [7]. In general, the response of early plants to water deficit involves complex mechanisms related to the reduction in growth and stomatal conductance with the concomitant decrease in photosynthesis activity [8,9] in order to prevent excessive water loss and tissue oxidation. Subsequently, depending on the duration and intensity of water stress, specific metabolic processes are activated on the basis of transcription of several categories of genes. Overall, it is considered that the genes expressed in response to drought can be divided into two broad categories: those expressing stress proteins and those encoding regulatory ones [10]. Despite these advances, the genetic and metabolic basis of the adaptation of plants to water stress still remain unclear because of the complexity of the metabolic networks and the specificity of ecotypes.

The onset of NGS (next-generation sequencing) technologies as well as metabolomic tools offer new opportunities to explore the response of the whole genome to environmental constraints. The de novo transcriptome assembly of *Q. pubescens* plantlets cultivated in greenhouse conditions led to the identification of several genes associated with drought avoidance [11]. However, in controlled environments, plant growth often does not totally reflect the conditions experienced in their natural setting, especially with regards to the prevailing light conditions, biotic interactions and soil conditions. In addition, transcriptomic analyses should be associated with ecophysiological and metabolic data for a complete description of the response of plants to environmental factors. Indeed, few studies have attempted to analyze transcriptome from plants grown in situ and even fewer have underpinned these studies with thorough in situ physiological characterization [12]. To address these requirements, the concept of ecological genomics has been proposed as an integrative discipline aimed at understanding the genetic mechanisms in relation to the responses of organisms to the conditions of their natural environment [13]. Ecological genomics thus associate the study of ecophysiological components to those of the transcripts, proteins and metabolites as to understand gene function and coordination between cellular processes and biochemical pathways.

The aim of this in situ study was first to test the hypothesis of plant–climate coevolution through the transcription profile for drought response. Second, a multi-omic approach linking transcriptomics
and metabolomics was carried out to bridge the gap between gene expression, metabolic and ecophysiological responses, and to provide new insights into plant water-stress-dependent metabolic networks. For these purposes, a robust transcriptome was assembled, improving also the scarcity of the genomic resource of *Q. pubescens* by generating a clean 50 Gb RNA-Seq dataset.

## 2. Materials and Methods

### 2.1. Climate Change Simulation in Oak Forest

Leaf samples were collected at the Oak Observatory of the Observatoire de Haute Provence (O<sub>3</sub>HP) in the vicinity of Manosque in Southern France (43°56′06.74″ N, 5°42′38.57″ E). The O<sub>3</sub>HP site is situated in a 100-year old forest dominated by Downy oak (*Quercus pubescens*) mixed with Montpellier maple (*Acer monspessulanum*, 25% of leaf biomass). Since 2012, a rain exclusion device consisting of a retractable roof was built over a plot of 300 m<sup>2</sup>, thus allowing us to apply a precipitation scenario close to a reduction in annual rainfall of about 30–40%. The control plot (300 m<sup>2</sup>) is close to the exclusion plot. The upper and lower canopy is accessed by a scaffolding (Figure 1). Five mature leaves per tree were harvested at the upper canopy level at solar noon, then pooled and immediately frozen in liquid nitrogen and stored at –80° until used. Ten trees were selected for each plot and leaf collection was carried out in spring (26 May 2016) and summer (8 September 2016).

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Oak Observatory of the Observatoire de Haute-Provence (O<sub>3</sub>HP). 43°56′06.74″ N, 5°42′38.57″ E. Rain exclusion device and the scaffolding for canopy measurements.</figcaption></figure>

### 2.2. RNA Extraction

Total RNA was isolated using the RNeasy Plant Mini Kit (Qiagen, France SAS, Villebon-sur-Yvette, France) following the manufacturer's instructions with minor modifications. The extracted RNA samples were then qualitatively and quantitatively checked using the Agilent 2100 Bioanalyzer with an RNA Nano chip and stored at -80 °C until further analyses. Ten replicates per plot and date were realized, resulting in forty samples extraction.

### 2.3. Library Preparation and Sequencing

The cDNA library for RNA-seq was constructed using TruSeq™ RNA Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions, and was then sequenced using Next Seq (Illumina Inc., San Diego, CA, USA) to obtain ~150 bp sequences from both ends of each cDNA. Poly-A selection was applied and two runs were carried out, which generated 50 million reads per sample in paired-end sequencing.
## 2.4. Transcriptome Assembly and Completeness

After quality assessment, adaptors were trimmed and low-quality reads (average quality control QC < 30) were filtered with Trimmomatic v0.38 [14] raw reads, which were then used for several transcriptome assembly tests using Trinity [15], CLC-assembly-cell-4.1.0, BinPacker [16]. Reads were also remapped (HISAT2-2.0.5) [17] against the genome of the close species *Quercus robur* available in the European Nucleotide Archive (ENA) (PRJEB7855). Mapping results were used to build a genome-guided de novo transcriptome assembly using the Trinity assembler (v2.5.1). Finally, only transcripts with at least a Transcripts Per Million (TPM) value above 1 and isoforms corresponding to more than 1% of the total gene count were kept. Transcriptome assembly completeness was evaluated by pseudo-mapping the sequencing reads back to assembly using Kallisto via trinityrnaseq-2.5.1 [18] and BUSCO-3.0.2 analysis [19].

## 2.5. Functional Annotation

Peptide prediction was performed using TransDecoder-3.0.1 [20]. A similarity search (Diamond of the Transdecoder predicted peptides) was performed against the UniProt/SwissProt, the Uniref90 and the nr databases (release: 5 December 2017). Peptide signal prediction was performed using SignalP v4. [21]. Transmembrane peptide detection was performed using TMHMM v2.0c [22]. Protein domain search was performed using hmmscan from the hmmer v.3.1b2 suite against the Pfam-A database [23] version 29.0. Finally, transcriptome functional annotation was performed using the RNAmmer and Trinotate-3.0.2 pipeline [20].

## 2.6. Identification of Differentially Expressed Genes Involved in Drought Tolerance and Gene Ontology Enrichment

Differential gene expression between drought treatments was identified using R software [24] and the DESeq2, edgeR and pcaMethods packages. The RNA transcripts with adjusted (Bonferroni correction) *p*-values < 0.001 and absolute fold changes ≥ 4 were considered as differentially expressed. A GO (gene ontology) term enrichment analysis was conducted jointly using the goseq R package. All these analyses were conducted using companion scripts of the Trinity package. To visualize the interactions between the differentially expressed gene datasets, Venn diagrams were realized in R with the Vennerable package.

## 2.7. Targeted Metabolomics by Gas Chromatography/Mass Spectrometry (GC/MS) Analysis

**Extraction:** Sample material (50 mg) was extracted in 360 µL of 100% MeOH for 15 min at 70 °C with vigorous shaking. After the addition of 200 µL of chloroform, samples were shaken at 37 °C for 5 min. To separate polar and organic phases, 400 µL of water was added, and samples were centrifuged for 10 min at 11,000× g. As internal standard, 20 µL of ribitol (0.2 mg/mL) was added to each sample. For the derivatization, 300 µL of the polar (upper) phase was transferred to a fresh tube and dried in a vacuum concentrator (speed-vac, Thermo Scientific, Savant, Germany) without heating.

**Derivatization (Methoximation and Silylation):** Pellets were re-dissolved in 20 µL of methoximation reagent containing 20 mg/mL of methoxyamine hydrochloride (Sigma 226904) in pyridine (Sigma 270970) and incubated for 2 h at 37 °C with shaking. For silylation, 32.2 µL of N-methyl-N-(trimethylsilyl)trifluoroacetamide (MSTFA; Sigma M7891) and 2.8 µL of alkane standard mixture (50mg/mL C₁₀–C₄₀; Fluka 68281) were added to each sample. After incubation for 30 min at 37 °C, samples were transferred to glass vials for GC/MS analysis.

**Gas Chromatography/Mass Spectrometry (GC/MS) Analysis:** A GC/MS-QP2010 Plus (Shimadzu®, Shimadzu Corporation, Kyoto, Japan) fitted with a Zebron ZB 5MS column (Phenomenex®; Torrance, CA, USA, 30 m × 0.25 mm × 0.25 µm) was used for GC/MS analysis. A sample of 1 µL was injected at 230 °C in split mode (1:10). The program started with a 1 min hold at 70 °C, followed by a 6 °C/min ramp to 310 °C, a 20 °C/min ramp to 330 °C and a bake-out for 5 min at 330 °C using helium as carrier
gas with constant linear velocity. The MS was operated with an ion source, an interface temperatures of 250 °C, a solvent cut time of 7 min and a scan range (m/z) of 40–700 with an event time of 0.2 s.

The "GCMS solution" software (Shimadzu®) and the NIST11 Mass Spectral library were used for data processing. Five replicates were carried out for each plot (control and exclusion) in spring and summer. The differential metabolites data were analyzed with MetaboAnalyst 4.0 [25].

## 2.8. Data Accessibility

All sequencing data have been deposited in the NCBI Sequence Read Archive under the accession code PRJNA381515.

# 3. Results

## 3.1. *High Throughput Sequencing Output and De Novo Transcriptome Assembly*

Leaves from twenty downy oak trees were used to generate 1.97 Md of paired-end raw reads 150 bp long. After cleaning in order to remove adaptors, low quality and short sequences, this results in 1.79 Md reads representing about 2.5 Tb of data (Table 1). The clean reads were mapped over the previously published de novo transcriptome of *Q. pubescens* [11], but only a covering of 40% was obtained. The studies were thus conducted for de novo transcriptome assembly with different tools, such as BinPacker, CLC and Trinity. The results showed that the read covering rate was also low, ranging from 50 to 75%. This led to a genome-guided transcriptome assembly based on Trinity algorithms using *Quercus robur* contigs as a reference [26]. A high-quality mapping of the clean reads was obtained ranging from 80 to 90% depending on the tested samples (Figure 2). This genome-guided transcriptome assembled more sequences compared to the other assemblers, precisely 530,080 transcripts with a mean length of 595 bp (Table 1). N50 and E90N50 were 780 and 1544, respectively. The filtered version of the assembly (transcripts with TPM > 1) contains 156,986 transcripts. Ninety percent of the global expression is present in 34,572 transcripts. The completeness of the assembly in terms of expected genes was checked through the BUSCO algorithm. As shown in Table 2, among the 303 and 1444 orthologous gene sets of Eukaryota and Embryophyta respectively searched, 92.4% and 83.8% were complete BUSCO single copies or duplicated for Eukaryota and Embryophyta, respectively.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Comparison of pseudoalignment rates between different assemblies. Kallisto analysis with the FR (Forward/Reverse) strand-specific paired-end pseudoalignment option. *Quercus pubescens* leaf samples from control (C) and exclusion (E) plots harvested in spring (p) and summer (e).</figcaption></figure>
<table><caption>**Table 1.** Summary of RNA-Seq and genome-guided sequence assemblies (raw and filtered) of *Q. pubescens*. TPM (Transcripts Per Kilobase Million).</caption><tr><td>Total raw paired-end reads</td><td>1,971,431,570 × 2</td></tr><tr><td>Total clean paired-end reads</td><td>1,794,615,817 × 2</td></tr><tr><td>Number of transcripts/trinity "genes"—raw assembly</td><td>530,080/395,969</td></tr><tr><td>Mean contig length (bp)/N50—raw assembly</td><td>595/780</td></tr><tr><td>E90N50/E90 number of transcripts</td><td>1544/34,572</td></tr><tr><td>Number of transcripts/trinity "genes"—filtered transcriptome (TPM &lt; 1)</td><td>156,986/126,106</td></tr><tr><td>Mean contig length (bp)/N50—filtered transcriptome (TPM &lt; 1)</td><td>720/1133</td></tr></table>

<table><caption>**Table 2.** BUSCO run for genome-guided de novo transcriptome assembly of *Q. pubescens* and comparison with *Q. robur* de novo assembly (OCV4_assembly_final.fas) [27].</caption><thead><tr><th rowspan="2">Transcriptomes<br>Assemblies</th><th colspan="2">Eukaryota</th><th colspan="2">Embryophyta</th></tr><tr><th>Q. pubescens<br>Filtered</th><th>Q. pubescens<br>Raw</th><th>Q. robur</th><th>Q. pubescens<br>Filtered</th><th>Q. pubescens<br>Raw</th></tr></thead><tbody><tr><td><strong>BUSCO sets</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Complete BUSCOs (C)</td><td>280</td><td>288</td><td>140</td><td>1206</td><td>1266</td></tr><tr><td>Complete BUSCOs%</td><td>92.40</td><td>95.10</td><td>46.30</td><td>83.80</td><td>87.90</td></tr><tr><td>Complete and single-copy BUSCOs (S)</td><td>112</td><td>103</td><td>115</td><td>645</td><td>601</td></tr><tr><td>Complete and duplicated BUSCOs (D)</td><td>168</td><td>185</td><td>25</td><td>561</td><td>665</td></tr><tr><td>Fragmented BUSCOs (F)</td><td>14</td><td>12</td><td>20</td><td>105</td><td>113</td></tr><tr><td>Missing BUSCOs (M)</td><td>9</td><td>3</td><td>143</td><td>129</td><td>61</td></tr><tr><td>Total BUSCO groups searched</td><td>303</td><td>303</td><td>303</td><td>1440</td><td>1440</td></tr></tbody></table>

## 3.2. Functional Annotation of the Assembled Contigs

Two annotation strategies were carried out: similarity search based on sequence- or functional domain-alignments. As shown in Table 3, there were in the raw transcriptome 200,726 unigenes that exhibited homology to known protein sequences in the Uniref90 BLASTX database (38%) followed by SwissProt BLASTX (27%) and Prot ID (26%) databases. A total of 75,052 contigs (14%) were mapped to the PFAM protein domain/family. For the annotation queries, the default filtering E value was set at <1 × 10⁻⁵. Sequence homology searches were carried out on the Viridiplantae database, which allows us to identify the species closest to *Q. pubescens*. In order of decreasing homology: *Juglans regia* (88,912 hits) followed by *Vitis vinifera* (10,707 hits), *Ziziphus jujube* (5239 hits), *Hevea brasiliensis* (4492 hits) and *Theobroma cacao* (4129 hits) were the closest species (Figure 3). These results highlight the genome specificity of *Q. pubescens* since, except for the genus *Quercus*, it shares only 17% of the genes expressed with the nearest plant species, namely, walnut.
<table><caption>**Table 3.** Summary of annotations of raw and filtered *Q. pubescens* genome-guided sequence assembly. KEGG (Kyoto Encyclopedia of Genes and Genomes); eggNOG database (evolutionary genealogy of genes: Non-supervised Orthologous Groups); PFAM (Protein family); TMHMM Transmembrane Helices Hidden Markov Models; RNAmmer (consistent and rapid annotation of ribosomal RNA genes).</caption><thead><tr><th rowspan="2">Database</th><th colspan="2">Q. pubescens Raw</th><th colspan="2">Q. pubescens Filtered</th></tr><tr><th>Number of transcripts</th><th>Percentage of transcripts</th><th>Number of transcripts</th><th>Percentage of transcripts</th></tr></thead><tbody><tr><td>Number of transcripts</td><td>530,080</td><td></td><td>156,986</td><td></td></tr><tr><td>Uniref90 DiamondX</td><td>214,856</td><td>40.54</td><td>102,445</td><td>65.26</td></tr><tr><td>SwissProt DiamondX</td><td>155,684</td><td>29.37</td><td>79,367</td><td>50.56</td></tr><tr><td>GO Extract from Sequence Similarity Search</td><td>146,689</td><td></td><td>75,852</td><td></td></tr><tr><td>Number of Protein Predicted</td><td>135,957</td><td></td><td>71,166</td><td></td></tr><tr><td>KEGG</td><td>124,285</td><td>91.14</td><td>67,348</td><td>94.63</td></tr><tr><td>Uniref90 DiamondP</td><td>120,789</td><td>88.85</td><td>66,374</td><td>93.27</td></tr><tr><td>EggNOG Function</td><td>118,026</td><td>86.81</td><td>63,724</td><td>89.54</td></tr><tr><td>SwissProt DiamondP</td><td>95,098</td><td>69.94</td><td>54,967</td><td>77.24</td></tr><tr><td>PFAM</td><td>86,021</td><td>63.27</td><td>49,462</td><td>69.5</td></tr><tr><td>GO extract from PFAM Analysis</td><td>52,107</td><td></td><td>31,994</td><td></td></tr><tr><td>TmHMM</td><td>22,131</td><td></td><td>11,986</td><td></td></tr><tr><td>RNAMMER</td><td>110</td><td></td><td>68</td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Plant species distribution of DiamondX searches on NCBI nr database of *Q. pubescens* transcripts. Percentage of total top ten homologies. Filtering at E value threshold < 10 × 10⁻¹⁰.</figcaption></figure>
### 3.3. *Differentially Expressed Genes (DEG) Analysis*

In order to address the question of the identification of the genes involved in Downy oak's response to the reduced precipitation treatment, two strategies were used: the principal component analysis (PCA) and an analysis based on a model using the negative binomial distribution for the estimation of the differential expression of genes (R, packages DESeq2 and edgeR). The PCA representation of the individuals on the first factorial plane 1–2 is displayed in Figure 4. The first two components represent 90% of the total variance with 89 and 1.5% for axes 1 and 2, respectively. Axis 2 clearly exhibited a season effect in terms of gene expression profile, while a slight trend was shown regarding to rain exclusion impact when considering axis 1. Furthermore, the Pearson correlation matrix showed that the correlation coefficient between samples was very high (up to 0.85); most of the expression data are similar and no differences were found. Some samples in summer clearly exhibited a rain exclusion effect. The use of R packages results in the identification of up to 5724 genes. Among these, 42 (31 and 11) were associated with the exclusion treatment, while 4261 were related to the season based on edgeR algorithms (Table 4). Further analysis showed that DESeq2 and edgeR share 3828 features. The Venn diagram, based on EdgeR differential analysis data, identified transcripts specific or common to different modalities (Figure 5). In particular, four of the 4153 differentially expressed genes are specific to the summer drought (eC vs. eE). The effect of seasonal variation (eC vs. pC) is determined by 236 transcripts, while the combined effect of treatment and season (eE vs. pE) is regulated by 108 features. We will focus on the 42 differentially expressed genes due to the aridification of the climate which exhibited 4 common features between spring and summer (Figure S2). When available, the corresponding TAIR (The *Arabidopsis* Information Resource) gene ID of *Q. pubescens* differential genes or contigs was provided (Table S1).

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Principal component analysis (PCA) of foliar transcript profile expression of *Q. pubescens* trees; e and p denote sampling in summer and spring, respectively; C and E denote control and rain exclusion treatments, respectively. Circles denote a possible exclusion effect in spring.</figcaption></figure>
<table><caption>**Table 4.** Differentially expressed genes (DEG) of *Q. pubescens* leaves in response to the sampling season and the rain exclusion treatment. The thresholds of statistical tests were fixed at $p_{adj} < 10^{-3}$ and $log2FC ≥ 2$. 5724, and 4153 DEG were detected by DESeq2 (bold) and edgeR, respectively; e and p denote sampling in summer and spring, respectively; C and E denote control and rain exclusion treatments, respectively.</caption><thead><tr><th>Samples</th><th>eC</th><th>eE</th><th>pC</th><th>pE</th></tr></thead><tbody><tr><td>eC</td><td>/</td><td><strong>15</strong></td><td><strong>3122</strong></td><td><strong>3529</strong></td></tr><tr><td>eE</td><td>31</td><td>/</td><td><strong>3460</strong></td><td><strong>3709</strong></td></tr><tr><td>pC</td><td>2312</td><td>2613</td><td>/</td><td><strong>2</strong></td></tr><tr><td>pE</td><td>2569</td><td>1949</td><td>11</td><td>/</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Venn diagram of edegR data showing the intersections between the differentially expressed genes from *Q. pubescens*: eC vs. eE; eE vs. pE; eC vs. pC; pC vs. pE; eE vs. pC and eC vs. pC. Control plot in summer, eC; rain exclusion plot in summer, eE; control plot in spring pC and rain exclusion plot in spring, pE.</figcaption></figure>

The GO (Gene Ontology) terms of the 42 differential expressed transcripts obtained from the edgeR algorithm are presented in Figures 6 and 7 which reflect the effect of climate aridification according to the seasons.
**Spring DEG.** Soil moisture levels at the sampling date were about 25% and 20% in controls and rain exclusion, respectively (Figure S1). Differential transcripts were nevertheless expressed. Upregulated features concern GO categories related to molecular function and biological processes (Figure 6A). The main molecular function that is exacerbated because of aridity is the activity of cation channels, precisely, that of potassium leak channels. This results in stomatal closure mediated by the jasmonic acid signal transduction pathway.

Downregulated genes are mainly related to biological processes, molecular function and cellular components of GO categories (Figure 6B). The most downregulated features involve the biosynthesis of raffinose oligosaccharides (RFOs) through the activity of the inositol 3-α-galactosyltransferase.

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** GO annotation of differential transcripts of *Q. pubescens* in spring as a function of the enrichment *p* value *p* < 0.01. (A) Upregulated features and (B) downregulated features. MF, molecular function; BP, biological process; CC, cellular component.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** GO annotation of differential transcripts of *Q. pubescens* in summer as a function of the enrichment *p* value *p* < 0.01. (A) Upregulated features and (B) downregulated features. MF, molecular function; BP, biological process.</figcaption></figure>

**Summer DEG.** Soil moisture levels in controls and rain exclusion were 13% and 11%, respectively (Figure S1). Upregulated genes concern GO categories that are related to molecular function and biological process (Figure 7A). Among the molecular functions identified, lyase activity was upregulated mainly by that of 4-hydroxy-4-methyl-2-oxoglutarate aldolase. The enzyme catalyzes the breakdown of 4-hydroxy-4-methyl-2-oxoglutarate (HMG) into two molecules of pyruvate. Furthermore, it contains a secondary oxaloacetate (OAA) decarboxylase activity because of the common pyruvate enolate intermediate formed following the breakdown of the C-C bond in the retro-aldolization and decarboxylation reactions. HMG aldolase also exhibited putative ribonuclease inhibitor activity. Among other molecular functions, ferrochelatase was identified. Additionally, the thiol-dependent ubiquitin-specific protease activity was enhanced as in spring. Ubiquitin-specific proteases (UBPs) are a family of hydrolases that remove polypeptides covalently linked to the C-terminal glycine of ubiquitin. UBPs regulate the ubiquitin/26S proteasome proteolytic pathway through deubiquitination, thus recycling ubiquitins during the breakdown of ubiquitin-protein conjugates. In turn, the main biological processes deal with the response to sucrose and monosaccharides and also with heme metabolic pathways including porphyrin-containing compounds biosynthesis.
Of the 12 genes downregulated in summer, only two belong to the category of molecular functions—the rest being involved in biological processes (Figure 7B). Peptide methionine (S)-S-oxide reductase A (PMR-A) catalyzes the reduction of oxidized methionine (S) in proteins/peptides. Another protein whose activity was downregulated is a phosphate ion transmembrane transporter. The biological processes affected by drought in summer concern hormonal, amino acid, mineral, DNA and oxidative metabolism. For plant survival, the most important molybdenum enzyme is the cytosolic nitrate reductase, which catalyzes the first step of nitrate assimilation. A weak activity of this enzyme was shown upon drought as well as a modification of chromatin structure.

### 3.4. Targeted Metabolomic

Twenty-three metabolites were identified in leaf samples (Table 5). As in RNA profiling, several metabolites were drought- and/or season-specific. Spring rain exclusion promoted more than a 1.5-fold increase in lactic acid, pyroglutamic acid, xylulose, isocitric acid and gallic acid contents in leaves of the exclusion plot as compared to control. However, a decrease in the accumulation and/or biosynthesis of myo-inositol, quinic acid and shikimic acid was observed. In summer, very few metabolites were differentially accumulated in the drought condition, namely, oxalic acid (2-fold), pyroglutamic acid (2-fold) and malic acid (1.6-fold). It should be noted that both seasons share the same differential increase in pyroglutamic acid and also the decrease in catechin (0.5-fold) contents in terms of adaptation to water stress. When considering the effect of the season, xylulose and sorbitol were the main sugar and sugar alcohol markers. Through metabolic pathway enrichment, it appears that the response of *Q. pubescens* to drought involves principally galactose, glyoxylate, a tricarboxylic acid (TCA) cycle and amino acid metabolisms in spring (Figure 8A).

<table><caption>**Table 5.** Targeted metabolomics of the trimethylsilyl/ethoxime (TMS) derivatized organic acids and sugars from leaf extracts of *Quercus pubescens*. * According to the NIST (National Institute of Standards and Technology) database. Data of the ratios are means of five independent extractions. E, rain exclusion; C, control; n.d., no data. No change is considered for a fold change of one.</caption><thead><tr><th rowspan="3">Name</th><th colspan="4">Fold Change</th></tr><tr><th>Summer</th><th>Spring</th><th>Spring/Summer</th><th>E</th></tr><tr><th>E/C</th><th>E/C</th><th>C</th><th>E</th></tr></thead><tbody><tr><td>Lactic acid</td><td>1.0</td><td>1.6</td><td>0.7</td><td>0.4</td></tr><tr><td>Oxalic acid</td><td>1.8</td><td>1.2</td><td>0.4</td><td>0.6</td></tr><tr><td>Phosphoric acid</td><td>1.2</td><td>n.d.</td><td>n.d.</td><td>0.8</td></tr><tr><td>Glycerol</td><td>1.0</td><td>1.3</td><td>0.9</td><td>0.7</td></tr><tr><td>Succinic acid</td><td>1.0</td><td>1.2</td><td>0.4</td><td>0.4</td></tr><tr><td>Aspartic acid</td><td>n.d.</td><td>1.1</td><td>n.d.</td><td>0.6</td></tr><tr><td>Pyroglutamic acid</td><td>1.8</td><td>2.1</td><td>1.3</td><td>1.1</td></tr><tr><td>Xylulose</td><td>0.9</td><td>1.8</td><td>9.8</td><td>4.9</td></tr><tr><td>Isocitric acid</td><td>1.0</td><td>1.6</td><td>0.9</td><td>0.6</td></tr><tr><td>Sorbitol</td><td>1.4</td><td>1.1</td><td>2.1</td><td>2.5</td></tr><tr><td>Viburnitol *</td><td>0.9</td><td>1.2</td><td>0.9</td><td>0.7</td></tr><tr><td>Inositol, myo-</td><td>1.2</td><td>0.6</td><td>0.7</td><td>1.2</td></tr><tr><td>Sucrose</td><td>1.1</td><td>1.1</td><td>0.6</td><td>0.6</td></tr><tr><td>Quinic acid</td><td>0.8</td><td>0.6</td><td>0.6</td><td>0.8</td></tr><tr><td>D-Ketohexose sugars</td><td>1.2</td><td>0.9</td><td>0.8</td><td>1.2</td></tr><tr><td>D-Aldohexoses</td><td>1.0</td><td>0.9</td><td>1.0</td><td>1.1</td></tr><tr><td>Erythronic acid</td><td>0.8</td><td>0.7</td><td>0.9</td><td>1.0</td></tr><tr><td>Hexadecanoic acid</td><td>0.8</td><td>1.1</td><td>1.5</td><td>1.2</td></tr><tr><td>Octadecanoic acid</td><td>0.8</td><td>1.0</td><td>1.0</td><td>0.8</td></tr><tr><td>Malic acid</td><td>1.6</td><td>1.3</td><td>1.2</td><td>1.5</td></tr><tr><td>Shikimic acid</td><td>n.d.</td><td>0.5</td><td>n.d.</td><td>n.d.</td></tr><tr><td>Gallic acid</td><td>1.0</td><td>1.6</td><td>0.6</td><td>0.4</td></tr><tr><td>Catechin</td><td>0.5</td><td>0.5</td><td>0.3</td><td>0.3</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Map of enriched metabolic pathways from differential metabolites of *Quercus pubescens* leaves; rain exclusion versus control during spring (A) and summer (B). Topological analysis of impact factors and the *p* value by pathway enrichment analysis. The color intensity (white to red) indicates the increasing *p* value, while the circle diameter covaries with pathway impact.</figcaption></figure>

**4. Discussion**

Our study provided the first genome-guided transcriptome of *Q. pubescens* leaves. Samples originate from individuals that were grown in situ on a site that has been a forest since the 18th century. As such, this dataset may be considered as a fundamental reference. Although high-throughput sequencing now allows one to assemble transcriptomes for non-model species, our comparative study has shown that assembly errors such as the formation of partial or chimeric transcripts can deeply affect the quality of the assemblies, and we demonstrate the utility of using closely related species to
perform genome-guided de novo transcriptome assembly [28]. A sequence homology search revealed that *Q. pubescens* exhibited the highest similarity with *J. regia* (Figure 3).

Differential analysis allowed us to identify four key genes common to both spring and summer, of which three are annotated: upregulated ubiquitin carboxyl-terminal hydrolase 12-like; downregulated putative calcium-binding protein cml26; and nifU-like protein 1, chloroplastic (Figure S2).

In a previous study carried out on the site of O₃HP, a two-fold increase in stomata closure was found in spring from *Q. pubescens* of the rain exclusion treatment as compared to the control plot [29].

Among spring DEG, an ionic channel gene was upregulated (Figure 6A). Several reports have shown that the increase in K⁺ efflux in stomata cell guards results in stomata closure [30–32]. This closure is mediated by both abscisic acid (ABA)-dependent and -independent signaling pathways [33]. Drought also induces the occurrence of ROS (Reactive Oxygen Species), which in turn mediate abscisic acid (ABA)- and jasmonate-induced stomatal closures [34,35]. It appears that drought induces *Q. pubescens* stomata closure in the signaling network of both ABA and jasmonic acid.

RFO transcripts were downregulated (Figure 6B); they are described as the most abundant galactose containing oligosaccharides in the plant kingdom. They are α-galactosyl extensions of sucrose, whose predominant functions are related to transport and storage [36]. While RFOs were identified as differential genes from plants submitted to abiotic and biotic stresses, their role remains unclear [36]. In our study, the downregulation of the biosynthesis of RFOs and more generally galactose metabolism may be interpreted as the need of carbon resource as a consequence of stomatal closure during spring.

Concerning summer DEG, despite the slight difference in soil moisture levels in control compared to rain exclusion (Figure S1), differential transcripts were expressed. Because the exclusion device was set up since 2010, it is likely that the response obtained is rather the consequence of long-term drought application. One may consider that the metabolic activity of *Q. pubescens* in response to the summer drought is pyruvate rate-limiting. Since the substrate of HMG aldolase is also an oxaloacetate, an intermediate of the tricarboxylic acid cycle, it is likely that drought stimulates its accumulation through the pyruvate carboxylation by the anaplerotic pathway. Another key upregulated transcript is ferrochelatase, whose sequence exhibited a high similarity with *Quercus suber* chloroplastic ferrochelatase-2 (FC2). The enzyme has been functionally characterized as involved in the biosynthesis of the heme of photosynthetic cytochromes and stress response proteins [37]. Our data are in accordance with previous studies that showed that stress induced by salinity and water stress result in upregulated genes of ferrochelatase [38–40].

UBPs helps confer *Arabidopsis thaliana* resistance to abiotic stress [41], and the upregulated level of the enzyme found in our study may be interpreted as the need for metabolic control by increased protein turnover during water stress. Peptide-methionine (S)-S-oxide reductase A (PMR-A) was shown to enhance cell resistance to ROS (reactive oxygen species) and to be located in tobacco chloroplasts [42]. An increase in resistance to drought and plant-pathogen attacks was also proposed upon identification of the overexpression of the PMR gene [43–45]. Since methionine is known as the most sensitive amino acid to oxidation in proteins or peptides, the downregulation of PMR-A in *Q. pubescens* may be interpreted as an accelerated proteolytic activity in view of the fact that specific amino acids demand satisfaction.

Summer stomatal conductance was significantly reduced to 56% in control as compared to the rain exclusion plot with a concomitant reduction of net photosynthesis [29]. Hence, it is probable that the decrease in phosphate ion transmembrane transporter activity reflects a balance of the redox state between photosynthesis light reactions and the Calvin cycle. Similarly, a decrease in the respiratory intensity could be considered, thus limiting the risk of generating ROS due to the importance of the electronic transport flow. The plastic growth and adaptability of plants to the external environmental conditions are mainly mediated by hormones of which auxin plays a fundamental role. Its biological function is determined by the strict coordination of three complex processes: auxin metabolism,
auxin transport, and auxin signaling [46]. The downregulation of tryptophan, the precursor of auxin biosynthesis, as well as that of auxin efflux, suggest that drought impairs auxin biosynthesis and transport. Our results seem to be in contradiction with other works that show that the application of auxin improves plant resistance to water stress [47,48]. Indeed, auxin efflux is mediated by specific PIN-FORMED transporters regulating the asymmetric spatiotemporal distribution of auxin in the whole plant. The downregulation of auxin efflux may result in a high accumulation of auxin more than that would be expected from the tryptophan biosynthetic pathway. The weak activity of molybdenum incorporation into the molybdenum–molybdopterin complex may be explained by a reduction in ATP (Adenosine Triphosphate) input from photosynthetic light reactions at the expense of nitrification during drought. Light absorbed in excess could induce a light-dependent cell death, a form of programmed cell death (PCD) that is unique to plants. The requirement of PCD has often been associated with the production of reactive oxygen species during photosynthesis [49]. Hence, the mitigated singlet oxygen-mediated PCD likely upstream of jasmonate would reduce the deleterious effects of ROS in downy oak during water stress. In many plant species, drought was accompanied by DNA methylation [50]. The downregulation of methylation-dependent chromatin silencing is in accordance with the need for the recruitment of specific transcripts in response to aridification.

Differential gene analysis in more severe drought conditions reveals that several features were involved in the tricarboxylic acids cycle and carbohydrate and amino acid metabolisms. This motivated us to conduct a directed metabolome analysis accordingly.

The large accumulation of malate and isocitrate (Table 5) may be interpreted as a tricarboxylic acid (TCA) shunt in favor of the glyoxylate cycle. Indeed, the glyoxylate cycle is an alternative to the production of glucose from AcetylCoA given the low availability of CO2. Water stress experiments in rice demonstrated a high expression of genes encoding isocitrate lyase and malate synthase, two key enzymes of the glyoxylate cycle [51]. In addition, there is an increase in oxalic acid concentrations, an intermediate product of the glyoxylate cycle but also recognized as being an environmental stress marker [52]. The metabolic profiling thus produced shows similarities with the gene expression pattern. Indeed, during the summer, the increased expression of many genes confirms a switch in the energetic metabolism towards the cyclic pathway of glyoxylate and that of gluconeogenesis to the detriment of the TCA cycle. In particular, the upregulation of the genes encoding oxaloacetate decarboxylase causes a high accumulation of pyruvate, and that of the oxo-acid lyase genes is correlated with the increased production of isocitrate through isocitrate lyase activity. This latter activity is supported by the strong accumulation of malate from the reaction between AcetylCoA and glyoxylate catalyzed by malate synthase. Furthermore, the downregulation of the inositol 3-α-galactosyltransferase activity is correlated with the decrease in myo-inositol content in response to drought. Throughout the seasons, drought induces an increasing accumulation of pyroglutamic acid (Table 5). This result is in accordance with the finding that the application of pyroglutamic acid highly promotes water stress tolerance effects in deficit irrigated lettuce, probably as a precursor of proline biosynthesis [53]. However, a cyclization of free glutamine and glutamate to free pyroglutamic acid in electrospray ionization-based mass spectrometry is not excluded [54].

Of the differentially accumulated metabolites and transcripts, we thus proposed a hypothetical scheme of the early response of *Q. pubescens* to in situ water withholding conditions (Figure 9). One of the consequences of stomata closure is the imbalance in the redox state between light and dark reactions of photosynthesis. To promote the availability of carbohydrates, we hypothesized a shunt of the TCA cycle through the glyoxylate cycle in relation to the formation of oxaloacetate. The latter compound is decarboxylated by an oxaloacetate decarboxylase, leading to the formation of two molecules of pyruvate that enter either into the neoglucogenesis pathway or as a precursor of AcetylCoA formation. Simultaneously, the carbohydrate storage rate is reduced since the activity of the inositol 3-α-galactosyltransferase is downregulated. Furthermore, the upregulation of ubiquitin-like-protein-specific protease activity suggests that drought favors intense proteolytic activities coupled with the availability of free ubiquitin.
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Hypothetical scheme of the early metabolic pathways involved in the adaptation of *Quercus pubescens* to drought. Complex proteins and metabolic pathways (green); proteins of upregulated transcripts (blue circle); upregulated metabolites (blue); downregulated metabolite (red); protein of downregulated transcript (red circle). Oxalocetate (OAA); oxaloacetate decarboxylase (OAAD); ubiquitin-like-protein-specific protease activity (UBPA); inositol 3-α-galactosyltransferase activity (IGAT); phosphoenolpyruvate (PEP); uridin diphosphate (UDP); tricarboxylic acid (TCA); isocitrate lyase (ICL); coenzymeA (CoA).</figcaption></figure>

In conclusion, the present work provides an important dataset on downy oak genomics. We assembled for the first time the genome-guided transcriptome of *Q. pubescens*. This work generated 1.7 Md of paired-end long reads of 150 nucleotides, therefore contributing to feed the genomic database of the genus *Quercus* hitherto limited to a few species. Potential functions of about 102,000 filtered transcripts from *Q. pubescens* were described in this study. The application of long-term water stress (four years) at O3HP site triggers adaptive processes, which confirms the expected genomic plasticity in regard to global changes. The differential analysis of transcripts coupled with the targeted metabolomics allowed us to demonstrate the overlap of changes in gene expression with the occurrence of some metabolites, particularly in relation to the TCA cycle. That said, the Krebs cycle might be involved in the early metabolic changes in plant response to drought together with the glyoxylate cycle for de novo synthesis of glucose. The highly upregulated gene of ubiquitin carboxyl-terminal hydrolase 12-like plays an important role in ubiquitin recycling during protein degradation. Other candidate genes were identified, mainly those of ferrochelatase 2 and porphyrin metabolism involved in photosynthesis electron transport; galactose metabolism and carbohydrate storage; and antioxidative defense by downregulating singlet oxygen-mediated PCD. However, the initiating events that trigger the regulation of these genes are unclear, although an epigenetic mechanism related to the heterochromatic gene silencing was shown. Further studies of functional genomics will be necessary to confirm or invalidate the role of the different identified candidate genes. Specifically, this study, while aiming to elucidate
the metabolic pathways involved in plants tolerance to drought, can also serve as a resource for further
functional genomic and phylogenetic analyses on the group of species of the *Quercus* genus.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2223-7747/9/9/1149/s1, Figure S1: O3HP soil water status taking into account the different sampling dates. Figure S2: Differential transcripts of leaves of *Quercus pubescens* annotated from an NCBI Blastn search: (a) spring and (b) summer. Blue, upregulated features; red, downregulated features. Table S1: *Quercus pubescens* genes or contigs differentially expressed during rain exclusion treatment.

**Author Contributions:** J.-P.M., funding acquisition, planned and designed the research, annotation of differentially expressed genes, enriched metabolic pathways, wrote the manuscript; B.L., funding acquisition, designed the research, library construction, high throughput sequencing, revised the manuscript; X.L., bioinformatics data treatment, revised the manuscript; E.C., designed the research, bioinformatics data treatment, revised the manuscript; M.T., RNA extraction; M.B., metabolomics analysis; A.H., RNA extraction; I.M.R., produced O3HP site data and employed CC treatments, revised the manuscript; C.F., revised the manuscript; T.G., revised the manuscript. All authors have read and agreed to the published version of the manuscript.

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

**Acknowledgments:** The authors gratefully acknowledge ECCOREV for funding this research. High throughput sequencing was performed at the TGML Platform, supported by grants from Inserm, GIS IBiSA, Aix-Marseille Université, and ANR-10-INBS-0009-10. The authors thank the Observatoire de Haute-Provence (OHP) and the programs AnaEE and Cooperate for their support. Bruno Fady (INRA of Avignon) and Alex Bomel (IMBE) are gratefully acknowledged for their scientific advice. Christophe Plomion (INRA of Bordeaux) is acknowledged for providing *Quercus robur* genomic sequences. We are grateful to the Institut Français de Bioinformatique (ANR-11-INBS-0013) and the Roscoff Bioinformatics platform ABiMS (http://abims.sb-roscoff.fr) for providing computing and storage resources.

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

**References**

1. AR4 Climate Change 2007: The Physical Science Basis. Available online: https://www.ipcc.ch/report/ar4/wg1/ (accessed on 10 June 2020).

2. Myers, N.; Mittermeier, R.A.; Mittermeier, C.G.; Fonseca, G.A.B.; Kent, J. Biodiversity hotspots for conservation priorities. *Nature* **2000**, *403*, 853–858. [CrossRef]

3. Giorgi, F.; Lionello, P. Climate change projections for the Mediterranean region. *Glob. Planet. Chang.* **2008**, *63*, 90–104. [CrossRef]

4. Polade, S.D.; Pierce, D.W.; Cayan, D.R.; Gershunov, A.; Dettinger, M.D. The key role of dry days in changing regional climate and precipitation regimes. *Sci. Rep.* **2014**, *4*, 4364. [CrossRef]

5. Holland, V.; Koller, S.; Lukas, S.; Brüggemann, W. Drought- and frost-induced accumulation of soluble carbohydrates during accelerated senescence in *Quercus pubescens*. *Trees* **2016**, *30*, 215–226. [CrossRef]

6. Sperlich, D.; Chang, C.T.; Penuelas, J.; Gracia, C.; Sabate, S. Seasonal variability of foliar photosynthetic and morphological traits and drought impacts in a Mediterranean mixed forest. *Tree Physiol.* **2015**, *35*, 501–520. [CrossRef] [PubMed]

7. Contran, N.; Günthardt-Goerg, M.S.; Kuster, T.M.; Cerana, R.; Crosti, P.; Paoletti, E. Physiological and biochemical responses of *Quercus pubescens* to air warming and drought on acidic and calcareous soils. *Plant Biol.* **2013**, *15*, 157–168. [CrossRef] [PubMed]

8. Peñuelas, J.; Llusià, J. Linking photorespiration, monoterpenes and thermotolerance in *Quercus*. *New Phytol.* **2002**, *155*, 227–237. [CrossRef]

9. Prieto, P.; Peñuelas, J.; Llusià, J.; Asensio, D.; Estiarte, M. Effects of long-term experimental night-time warming and drought on photosynthesis, Fv/Fm and stomatal conductance in the dominant species of a Mediterranean shrubland. *Acta Physiol. Plant* **2009**, *31*, 729–739. [CrossRef]

10. Gugger, P.F.; Penaloza-Ramirez, J.M.; Wright, J.W.; Sork, V.L. Whole-transcriptome response to water stress in a California endemic oak, Quercus lobata. *Tree Physiol.* **2017**, *37*, 632–644. [CrossRef]

11. Torre, S.; Tattini, M.; Brunetti, C.; Fineschi, S.; Fini, A.; Ferrini, F.; Sebastiani, F. RNA-Seq Analysis of *Quercus pubescens* Leaves: De Novo Transcriptome Assembly, Annotation and Functional Markers Development. *PLoS ONE* **2014**, *9*, e112487. [CrossRef]
12. Yates, S.A.; Chernukhin, I.; Alvarez-Fernandez, R.; Bechtold, U.; Baeshen, M.; Baeshen, N.; Mutwakil, M.Z.; Sabir, J.; Lawson, T.; Mullineaux, P.M. The temporal foliar transcriptome of the perennial C3 desert plant Rhazya stricta in its natural environment. *BMC Plant Biol.* **2014**, *14*, 2. [CrossRef] [PubMed]

13. Ungerer, M.C.; Johnson, L.C.; Herman, M.A. Ecological genomics: Understanding gene and genome function in the natural environment. *Heredity* **2008**, *100*, 178–183. [CrossRef] [PubMed]

14. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. *Bioinform. Oxf. Engl.* **2014**, *30*, 2114–2120. [CrossRef] [PubMed]

15. Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q.; et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. *Nat. Biotechnol.* **2011**, *29*, 644–652. [CrossRef] [PubMed]

16. Liu, J.; Li, G.; Chang, Z.; Yu, T.; Liu, B.; McMullen, R.; Chen, P.; Huang, X. BinPacker: Packing-Based De Novo Transcriptome Assembly from RNA-seq Data. *PLoS Comput. Biol.* **2016**, *12*, e1004772. [CrossRef]

17. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. *Nat. Methods* **2015**, *12*, 357–360. [CrossRef]

18. Bray, N.L.; Pimentel, H.; Melsted, P.; Pachter, L. Near-optimal probabilistic RNA-seq quantification. *Nat. Biotechnol.* **2016**, *34*, 525–527. [CrossRef]

19. Waterhouse, R.M.; Seppey, M.; Simão, F.A.; Manni, M.; Ioannidis, P.; Klioutchnikov, G.; Kriventseva, E.V.; Zdobnov, E.M. BUSCO Applications from Quality Assessments to Gene Prediction and Phylogenomics. *Mol. Biol. Evol.* **2018**, *35*, 543–548. [CrossRef]

20. Haas, B.J.; Papanicolaou, A.; Yassour, M.; Grabherr, M.; Blood, P.D.; Bowden, J.; Couger, M.B.; Eccles, D.; Li, B.; Lieber, M.; et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. *Nat. Protoc.* **2013**, *8*, 1494–1512. [CrossRef]

21. Petersen, T.N.; Brunak, S.; von Heijne, G.; Nielsen, H. SignalP 4.0: Discriminating signal peptides from transmembrane regions. *Nat. Methods* **2011**, *8*, 785–786. [CrossRef]

22. Krogh, A.; Larsson, B.; von Heijne, G.; Sonnhammer, E.L. Predicting transmembrane protein topology with a hidden Markov model: Application to complete genomes. *J. Mol. Biol.* **2001**, *305*, 567–580. [CrossRef] [PubMed]

23. Finn, R.D.; Bateman, A.; Clements, J.; Coggill, P.; Eberhardt, R.Y.; Eddy, S.R.; Heger, A.; Hetherington, K.; Holm, L.; Mistry, J.; et al. Pfam: The protein families database. *Nucleic Acids Res.* **2014**, *42*, D222–D230. [CrossRef] [PubMed]

24. R: The R Project for Statistical Computing. Available online: https://www.r-project.org/ (accessed on 10 July 2020).

25. Chong, J.; Soufan, O.; Li, C.; Caraus, I.; Li, S.; Bourque, G.; Wishart, D.S.; Xia, J. MetaboAnalyst 4.0: Towards more transparent and integrative metabolomics analysis. *Nucleic Acids Res.* **2018**, *46*, W486–W494. [CrossRef] [PubMed]

26. Plomion, C.; Aury, J.-M.; Amselem, J.; Alaeitabar, T.; Barbe, V.; Belser, C.; Bergès, H.; Bodénès, C.; Boudet, N.; Boury, C.; et al. Decoding the oak genome: Public release of sequence data, assembly, annotation and publication strategies. *Mol. Ecol. Resour.* **2016**, *16*, 254–265. [CrossRef]

27. OAK GENOME SEQUENCING: Sequencing of the Oak Genome and Identification of Genes that Matter for Forest Tree Adaptation. Available online: http://www.oakgenome.fr/ (accessed on 10 July 2020).

28. McGettigan, P.A. Transcriptomics in the RNA-seq era. *Curr. Opin. Chem. Biol.* **2013**, *17*, 4–11. [CrossRef] [PubMed]

29. Saunier, A.; Ormeño, E.; Havaux, M.; Wortham, H.; Ksas, B.; Temime-Roussel, B.; Blande, J.D.; Lecareux, C.; Mévy, J.-P.; Bousquet-Mélou, A.; et al. Resistance of native oak to recurrent drought conditions simulating predicted climatic changes in the Mediterranean region. *Plant Cell Environ.* **2018**, *41*, 2299–2312. [CrossRef] [PubMed]

30. Hosy, E.; Vavasseur, A.; Mouline, K.; Dreyer, I.; Gaymard, F.; Porée, F.; Boucherez, J.; Lebaudy, A.; Bouchez, D.; Véry, A.-A.; et al. The Arabidopsis outward K+ channel GORK is involved in regulation of stomatal movements and plant transpiration. *Proc. Natl. Acad. Sci. USA* **2003**, *100*, 5549–5554. [CrossRef]

31. Zhang, A.; Ren, H.-M.; Tan, Y.-Q.; Qi, G.-N.; Yao, F.-Y.; Wu, G.-L.; Yang, L.-W.; Hussain, J.; Sun, S.-J.; Wang, Y.-F. S-Type Anion Channels SLAC1 and SLAH3 Function as Essential Negative Regulators of Inward K+ Channels and Stomatal Opening in Arabidopsis. *Plant Cell* **2016**, *28*, 949–965. [CrossRef]
32. Jin, Z.; Wang, Z.; Ma, Q.; Sun, L.; Zhang, L.; Liu, Z.; Liu, D.; Hao, X.; Pei, Y. Hydrogen sulfide mediates ion fluxes inducing stomatal closure in response to drought stress in *Arabidopsis thaliana*. *Plant Soil* **2017**, *419*, 141–152. [CrossRef]

33. Luan, S. Signalling drought in guard cells. *Plant Cell Environ.* **2002**, *25*, 229–237. [CrossRef]

34. Wang, P.; Song, C.-P. Guard-cell signalling for hydrogen peroxide and abscisic acid. *New Phytol.* **2008**, *178*, 703–718. [CrossRef] [PubMed]

35. Yin, Y.; Adachi, Y.; Nakamura, Y.; Munemasa, S.; Mori, I.C.; Murata, Y. Involvement of OST1 Protein Kinase and PYR/PYL/RCAR Receptors in Methyl Jasmonate-Induced Stomatal Closure in *Arabidopsis* Guard Cells. *Plant Cell Physiol.* **2016**, *57*, 1779–1790. [CrossRef] [PubMed]

36. Sengupta, S.; Mukherjee, S.; Basak, P.; Majumder, A.L. Significance of galactinol and raffinose family oligosaccharide synthesis in plants. *Front. Plant Sci.* **2015**, *6*, 656. [CrossRef] [PubMed]

37. Scharfenberg, M.; Mittermayr, L.; Roepenack-Lahaye, E.V.; Schlicke, H.; Grimm, B.; Leister, D.; Kleine, T. Functional characterization of the two ferrochelatases in *Arabidopsis thaliana*. *Plant Cell Environ.* **2015**, *38*, 280–298. [CrossRef]

38. Wu, H.; Shabala, L.; Zhou, M.; Shabala, S. Chloroplast-generated ROS dominate NaCl-induced K+ efflux in wheat leaf mesophyll. *Plant Signal. Behav.* **2015**, *10*. [CrossRef]

39. Wu, Y.; Jin, X.; Liao, W.; Hu, L.; Dawuda, M.M.; Zhao, X.; Tang, Z.; Gong, T.; Yu, J. 5-Aminolevulinic Acid (ALA) Alleviated Salinity Stress in Cucumber Seedlings by Enhancing Chlorophyll Synthesis Pathway. *Front. Plant Sci.* **2018**, *9*. [CrossRef]

40. Liu, D.; Kong, D.D.; Fu, X.K.; Ali, B.; Xu, L.; Zhou, W.J. Influence of exogenous 5-aminolevulinic acid on chlorophyll synthesis and related gene expression in oilseed rape de-etiolated cotyledons under water-deficit stress. *Photosynthetica* **2016**, *54*, 468–474. [CrossRef]

41. Yan, N.; Doelling, J.H.; Falbel, T.G.; Durski, A.M.; Vierstra, R.D. The Ubiquitin-Specific Protease Family from *Arabidopsis*. *AtUBP1* and 2 Are Required for the Resistance to the Amino Acid Analog Canavanine. *Plant Physiol.* **2000**, *124*, 1828–1843. [CrossRef]

42. Liu, L.; Wang, M.H. Expression and Biological Properties of a Novel Methionine Sulfoxide Reductase A in Tobacco (*Nicotiana tabacum*). *Protein J.* **2013**, *32*, 266–274. [CrossRef]

43. Romero, H.M.; Berlett, B.S.; Jensen, P.J.; Pell, E.J.; Tien, M. Investigations into the Role of the Plastidial Peptide Methionine Sulfoxide Reductase in Response to Oxidative Stress in *Arabidopsis*. *Plant Physiol.* **2004**, *136*, 3784–3794. [CrossRef]

44. Oh, S.-K.; Baek, K.-H.; Seong, E.S.; Joung, Y.H.; Choi, G.-J.; Park, J.M.; Cho, H.S.; Kim, E.A.; Lee, S.; Choi, D. CaMsrB2, Pepper Methionine Sulfoxide Reductase B2, Is a Novel Defense Regulator against Oxidative Stress and Pathogen Attack. *Plant Physiol.* **2010**, *154*, 245–261. [CrossRef] [PubMed]

45. Gao, Y.; Zhu, J.; Ding, P.; Chen, F.; Xia, G. The *Brachypodium distachyon* methionine sulfoxide reductase gene family. *Genes Genom.* **2017**, *39*, 975–985. [CrossRef]

46. Forestan, C.; Varotto, S. The Role of PIN Auxin Efflux Carriers in Polar Auxin Transport and Accumulation and Their Effect on Shaping Maize Development. *Mol. Plant* **2012**, *5*, 787–798. [CrossRef] [PubMed]

47. Kim, J.I.; Baek, D.; Park, H.C.; Chun, H.J.; Oh, D.-H.; Lee, M.K.; Cha, J.-Y.; Kim, W.-Y.; Kim, M.C.; Chung, W.S.; et al. Overexpression of Arabidopsis YUCCA6 in Potato Results in High-Auxin Developmental Phenotypes and Enhanced Resistance to Water Deficit. *Mol. Plant* **2013**, *6*, 337–349. [CrossRef] [PubMed]

48. Lecube, M.L.; Noriega, G.O.; Cruz, D.M.S.; Tomaro, M.L.; Batlle, A.; Balestrasse, K.B. Indole acetic acid is responsible for protection against oxidative stress caused by drought in soybean plants: The role of heme oxygenase induction. *Redox Rep.* **2014**, *19*, 242–250. [CrossRef] [PubMed]

49. Danon, A.; Coll, N.S.; Apel, K. Cryptochrome-1-dependent execution of programmed cell death induced by singlet oxygen in *Arabidopsis thaliana*. *Proc. Natl. Acad. Sci. USA* **2006**, *103*, 17036–17041. [CrossRef]

50. Deng, J.; Kou, S.; Zhang, C.; Zou, Q.; Li, P.; Yuan, P. DNA Methylation and Plant Stress Responses. *J. Plant Physiol. Pathol.* **2018**, *2018*. [CrossRef]

51. Maruyama, K.; Urano, K.; Yoshiwara, K.; Morishita, Y.; Sakurai, N.; Suzuki, H.; Kojima, M.; Sakakibara, H.; Shibata, D.; Saito, K.; et al. Integrated Analysis of the Effects of Cold and Dehydration on Rice Metabolites, Phytohormones, and Gene Transcripts. *Plant Physiol.* **2014**, *164*, 1759–1771. [CrossRef]

52. Yu, L.; Jiang, J.; Zhang, C.; Jiang, L.; Ye, N.; Lu, Y.; Yang, G.; Liu, E.; Peng, C.; He, Z.; et al. Glyoxylate rather than ascorbate is an efficient precursor for oxalate biosynthesis in rice. *J. Exp. Bot.* **2010**, *61*, 1625–1634. [CrossRef]
53. Jiménez-Arias, D.; García-Machado, F.J.; Morales-Sierra, S.; Luis, J.C.; Suarez, E.; Hernández, M.; Valdés, F.; Borges, A.A. Lettuce plants treated with L-pyroglutamic acid increase yield under water deficit stress. *Environ. Exp. Bot.* **2019**, *158*, 215-222. [CrossRef]

54. Purwaha, P.; Silva, L.P.; Hawke, D.H.; Weinstein, J.N.; Lorenzi, P.L. An artifact in LC-MS/MS measurement of glutamine and glutamic acid: In-source cyclization to pyroglutamic acid. *Anal. Chem.* **2014**, *86*, 5633-5637. [CrossRef] [PubMed]

<figure><img src="image_16.png" /></figure>

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