# Molecular plasticity to soil water deficit differs between sessile oak (Quercus Petraea (Matt.) Liebl.) high- and low-water use efficiency genotypes.

Gregoire Le Provost, Theo Gerardin, Christophe Plomion, Oliver Brendel

► To cite this version:

Gregoire Le Provost, Theo Gerardin, Christophe Plomion, Oliver Brendel. Molecular plasticity to soil water deficit differs between sessile oak (Quercus Petraea (Matt.) Liebl.) high- and low-water use efficiency genotypes.. Tree Physiology, 2022, 10.1093/treephys/tpac087/6648420. [https://hal.inrae.fr/hal-03739690](https://hal.inrae.fr/hal-03739690)

HAL Id: hal-03739690

https://hal.inrae.fr/hal-03739690

Submitted on 28 Jul 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Molecular plasticity to soil water deficit differs between sessile oak (Quercus
Petraea (Matt.) Liebl.) high- and low-water use efficiency genotypes.

Gregoire Le Provost¹, Theo Gerardin², Christophe Plomion¹ and Oliver Brendel².

1: Univ. Bordeaux, INRAE, BIOGECO, F-33610 Cestas, France.

2: Université de Lorraine, AgroParisTech, INRAE, UMR Silva, Nancy, France

† Corresponding author: Grégoire Le Provost

Univ. Bordeaux, INRAE, BIOGECO, F-33610 Cestas, France.

Tel: +33(0)535385334,

Fax: +33(0)557122881

E-mail: gregoire.le-provost@inrae.fr

Running Head: Molecular mechanisms of WUE in Sessile oaks.

Key words: WUE, GxE, gene expression, drought adaptation, sessile oak.

© The Author(s) 2022. Published by Oxford University Press. All rights reserved. For
permissions, please e-mail: journals permissions@oup.com
Abstract

Water use efficiency (WUE) is an important adaptive trait for soil water deficit. The molecular and physiological bases of WUE regulation in crops have been studied in detail in the context of plant breeding. Knowledge for most forest tree species lags behind, despite the need to identify populations or genotypes able to cope with the longer, more intense drought periods likely to result from climate warming. We aimed to bridge this gap in knowledge for sessile oak (Quercus petraeae (Matt.) L.), one of the most ecologically and economically important tree species in Europe, using a factorial design including trees with contrasted phenotypic values (low and high WUE) and two watering regimes (control and drought). By monitoring the ecophysiological response, we first qualified genotypes for their WUE (by using instantaneous and long-term measures). We then performed RNA-seq to quantify gene expression for the three most extreme genotypes exposed to the two watering regimes. By analyzing the interaction term, we were able to capture the molecular strategy of each group of plants for coping with drought. We identified putative candidate genes potentially involved in the regulation of transpiration rate in high WUE phenotypes. Regardless of water availability, trees from the high WUE phenotypic class overexpressed genes associated with drought responses, and in the control of stomatal density and distribution, and displayed a downregulation of genes associated with early stomatal closure and high transpiration rate.

Fine physiological screening of sessile oaks with contrasting WUE, and their molecular characterization i) highlighted subtle differences in transcription between low and high WUE genotypes, identifying key molecular players in the genetic control of this trait, and ii) revealed the genes underlying the molecular strategy that evolved in each group to potentially cope with water deficit, providing new insight into the within species diversity in drought adaptation strategies.
Introduction

Maintenance of the net primary production of trees, and, therefore, of the capacity of planted
and natural forests to buffer the effects of global warming, is dependent on various biotic and
abiotic factors, including their ability to adapt to water shortages during the growing period
(Keenan 2015).

It has been suggested that water use efficiency (WUE), defined as biomass production divided
by the water used through transpiration (Aranda et al. 2007, Brendel 2021), plays a functional
role in acclimation to water scarcity. Indeed, this trait, whether assessed instantaneously or
through integrated measurements, is determined environmentally (as a function of, for
example, soil water deficit, vapor pressure deficit or temperature) and genetically (Brendel
2021). Maritime pine, a tree species scattered throughout the Mediterranean basin, constitutes
a unique case study illustrating these findings. The phenotypic plasticity of WUE has been
shown to depend, among others, on soil water deficit and evaporative demand (Marguerit et
al. 2014, Plomion et al. 2016), and ample genetic variation has been found both within
breeding populations (Brendel et al. 2002, Marguerit et al. 2014) and between natural
populations adapted to local climatic or edaphic selection pressures (Guehl et al. 1995,
Aranda et al. 2010, Lamy et al. 2011). One study (Tognetti et al. 2011) also reported a
significant GxE interaction, suggesting that maritime pine gene pools have adopted different
WUE strategies in response to drought.

Sessile oak (*Quercus petraea* (Matt.) Liebl.) is a dominant species in Europe. It has a very
wide ecological niche, tolerating pH values from 3.5 to 9; it can tolerate all degrees of soil
wetness, from xeric to moderatly wet soils, and it is present in a zone extending from southern
Spain to Scandinavia, with a large climatic amplitude (Ducousos and Bordacs 2004). It is
relatively tolerant to drought and poor soils, but is sensitive to waterlogging conditions. It has
been reported to have a variable WUE in response to drought directly linked to stomatal conductance (Ponton et al. 2002). Both environmental (plasticity to drought) and genetic variation (within and between populations) have been shown to drive this variation in *Quercus robur* L., a close, sympatric species. (Ponton et al. 2002, Brendel et al. 2008, Roussel et al. 2009a).

Studies in oaks have contributed to our understanding of the physiological mechanisms underlying WUE diversity in trees. Studies in pedunculate oak genotypes presenting extreme phenotypic values for WUE (Roussel et al. 2009a, Roussel et al. 2009b) clearly demonstrated (in well-watered conditions) the genetic effect of stomatal conductance on leaf and whole plant WUE, particularly in terms of stomatal density and the diurnal dynamics of stomata.

Several studies have investigated the molecular mechanisms underlying WUE variation in annual plants. Most of these studies were conducted in the model plants *Arabidopsis thaliana* and *Oryza sativa* (reviewed in (Ruggiero et al. 2017)). Changes in stomatal aperture have been shown to be linked to both drought stress tolerance and WUE (Chaerle et al. 2005, Kim et al. 2010). Based on these findings, several genes involved in the stomatal development pathway have emerged as promising candidate genes for the control of WUE (reviewed in Yoo et al. 2010). In particular, the transcription factor GTL1 has been shown to regulate WUE by modifying stomatal density through repression of the SDD1 gene (Yoo et al. 2010).

The HARDY gene, encoding an AP2/ERF-like transcription factor, has been reported to affect WUE directly in rice, by enhancing photosynthetic assimilation and efficiency (Karaba et al. 2007). In *Arabidopsis*, the *ERECTA* gene, encoding a putative leucine-rich receptor-like protein kinase (LRR-RLK) regulates WUE through its action on leaf morphology (i.e. mesophyll proliferation) and stomatal density (Masle et al. 2005). Other authors have reported that the AtEDT1/HDG11 gene, encoding a Homeodomain START (HD-START) transcription factor, also regulates stomatal density and WUE through its interaction with the
ERECTA and E2Fa in both *Arabidopsis* and tobacco (Chaerle et al. 2005). In maize, the Asr1 transcription factor regulates WUE by modulating C4 (C4-PEPC) activity (Jeanneau et al. 2002). A possible role for carbonic anhydrase, which catalyzes CO₂ hydration, has also recently been suggested. Indeed, the WUE of double mutants for anhydrase carbonic genes is affected in drought stress conditions, suggesting an important role for these genes in maize (Kolbe et al. 2018). The molecular mechanisms involved in WUE regulation in forest trees have been little studied, and we are aware of only two studies, one in oak (Roussel et al. 2009a) and the other in poplar (Weng et al. 2012), confirming the regulation of WUE by known candidate genes (ERECTA and the GTL1), which had been identified previously in studies on annual plants.

In this study, we aimed to bridge this gap in our knowledge, by investigating molecular signatures in the leaves of sessile oak genotypes with extreme phenotypic values for WUE subjected to drought. We used a full factorial design in which oak seedlings with contrasting WUE (G) were raised in drought and control conditions (E) and characterized phenotypically and molecularly. The selection of seedlings with contrasting WUE were based on instantaneous as well as long-term estimators of WUE, thereby selecting genotypes were instantaneous and long-term integrated measures were coherent, Transcripts "presenting a E or G" effect can be used to identify important gene pathways responding to drought and underlying WUE variation, respectively. Moreover, using this experimental design, we were able to assess adaptive function from the analysis of genes presenting GxE interactions, revealing molecular strategies based on high-WUE genotypes in an environment in which water supply is limited. This study provides important clues to the adaptive value of WUE in sessile oak and the underlying molecular mechanisms.
# Material and Methods

## Plant material and experimental design

The plants analyzed here were part of a larger experiment, carried out at INRAE (Champenoux, France, 48°45′8″N, 6°20′28″E, elevation: 259 m asl) on 60 *Q. petraea* seedlings originating from six different French provenances (Table 1). The seeds were selected from stands growing in conditions ranging from moderately dry to moderately humid (annual precipitations ranging from 688mm to 872mm) to increase the within-species variation. All these populations were attributed to the same chloroplast DNA lineage (Petit et al. 2002), thereby minimizing the variation of the genetic background. The acorns were collected in autumn 2015 and sown during spring 2016 in 6 L pots filled with a 5/3/2 (V/V/V) mixture of sand, peat and silty-clay forest soil. This soil mixture was tested to a field capacity (FC) of 33% soil water content (SWC). The plants were grown in a greenhouse equipped with a robotic system for the automatic weighing and watering of plants (Bogeat-Triboulot et al. 2019). All the plants were initially subjected to the same non-limiting growth conditions: natural light, with fertilization and irrigation at 85% of relative extractable water (REW, with field capacity at 33% and a wilting point at 4% soil volumetric humidity). The volumetric soil water content was measured by regular time domain reflectrometry (TRIME-TDR; IMKO GmbH, Ettlingen, DE), with one to two measurements per week.

Plants were randomly assigned to the control and drought groups at day 193, with the REW decreased to 31% REW for the water stress group (reached at day 202). The final harvest of the leaves used in this experiment took place during day 258, whereas the final harvest of the whole plants for biomass analysis took place during day 266. For more details see Gerardin (2019). This strategy was specifically chosen to identify the molecular mechanisms involved in a long term constitutive response, rather than stress perception. The objective of the
experiment was to decipher the molecular mechanisms involved in differences in WUE,
where two of the used estimators of WUE are long term and time integrated. We are aware
that such sampling strategy allow us to identify only the molecular mechanism involved in the
long term constitutive response rather than that implicated in the perception of the stress. We
made the choice of long-term investigation, because our main goal was to decipher the
molecular mechanism involved in WUE upon well-established stress condition.

Monitoring of leaf gas exchange

Regular *in situ* gas exchange monitoring was performed with a portable photosynthesis
system (Li-Cor 6200; Li-Cor, Lincoln, NE, USA). Net CO2 assimilation rate (A insitu), and
stomatal conductance to water vapor (gs insitu) were measured four times during the drought
period, on the same third-flush leaf. Leaf intrinsic water use efficiency (Wi) was estimated as
A/g, and an overall mean of all measurements during the drought period was then calculated
(Wi insitu).

Leaf gas exchange was also measured in a more standardized way using a Li-Cor 6400XT
(Li-Cor, Lincoln, NE, USA) equipped with a 2 cm² leaf chamber. All measurements were
performed between 10:00 and 19:00 h (Central European summertime), with the following
conditions in the leaf chamber: CO₂ concentration of 400 µmol mol⁻¹, temperature of 25°C,
~70% humidity and a PAR of 1200 µmol m⁻² s⁻¹. These measurements were performed three
times during the drought period for each plant. Wi was calculated as A/g and an overall mean
was then calculated for all three measure (A<sub>SS</sub>, g<sub>SS</sub>, W<sub>SS</sub>).

Final harvest
At the end of the experiment on the 266th day of the year, all the plants were harvested. Plants were oven-dried at 60 °C until they reached a constant dry mass. Leaves were weighed separately (LBM) to estimate the total leaf area (LA) according to the following allometry established in a previous study of (Gerardin, 2019), based on *Q, petraea* as well as *Q, robur* seedlings of similar size as used in this study and grown in the same greenhouse, : LA[cm²] = 146.2*LBm[g] (p<0.0001R² = 0.985, n = 59 plants).

The leaves used for standardized gas exchange were harvested and cut in half. One half of the leaf was used for stable carbon isotope analysis, and the other half was used for stomatal density measurements. The leaves for bulk leaf carbon isotope composition (δ¹³C) were ground in a ball mill and analyzed with an elemental analyzer (vario ISOTOPE cube, Elementar, Hanau, Germany) coupled to an isotope ratio mass spectrometer (IsoPrime 100, Isoprime Ltd, Cheadle, UK). We calculated δ¹³C according to the international standard (Vienna Pee Dee Belemnite, VPDB).

## Measurement of stomatal density

Oak trees have stomata only on the abaxial surface of the leaves. We collected 1 cm² segments of leaf and took nail varnish imprints of the abaxial surface with an adhesive film, which was then applied to a microscope slide for analysis. ImageJ2 software was used to count the stomata on the images obtained, excluding stomata overlapping the margins of the image.

## Estimation of transpiration efficiency

For the estimation of transpiration efficiency (TE), the cumulative transpiration (CT) of the plants was measured by summing the amount of water from all the watering cycles of each
plant over a four months period from mid-May to the final harvest in September, corrected for
the water evaporation from the soil surface from pots without plants (CT). This calculation
was also performed for the drought period only (CT<sub>d</sub>). Plant water consumption is also
expressed on a per unit leaf area basis (CT<sub>LA</sub> and CT<sub>d,LA</sub>), based on the leaf area at the end of
the experiment. Initial dry mass (BMi) was estimated using an allometric equation, based on a
large number of individuals from Q, petraea as well as Q, robur, sampled during several
experiments of (Gerardin, 2019). Overall 186 complete trees (including roots, leaves and
stems) had been harvested from the youngest seedling ages to plants over 2m in height
(H[mm]) and nearly 2cm in diameter (D[mm]). BMi was estimated using diameter and height
measurements as BMi[g] = 0,0598*[(D/2)<sup>2</sup>*pi*H*1/3]**0,7264 with a residual standard error
of &plusmn;0,44g. We also used this equation to estimate the dry biomass at the start of the drought
period. The increase in biomass was calculated for the whole period (BM) and for the drought
period only (BM<sub>d</sub>). TE was calculated by dividing accumulated biomass by water
consumption for the whole period (TE) and for the drought period only (TE<sub>d</sub>).

Selection of plants for RNAseq

At the final harvest, just before the plants were oven-dried, third-flush leaves close to those
used for gas exchange were harvested from all 60 plants and frozen in liquid nitrogen for
RNA extraction. We selected 12 plants based on Wi insitu, Wi ss, TE, and TEd. Individual trees
were chosen so as to constitute high- and low-WUE phenotypic groups consistent for all four
measurements of WUE, thereby selecting trees where the instantaneous measurements were
coherent with the resulting long term TE. First, for each trait and separately for control and
drought plants, the upper and lower quantiles were chosen. The trait values were transformed
into ranks and the mean rank over the four traits was calculated for each plant. The six plants
with the most extreme ranks in each group were selected. Finally, three plants per group were
chosen, by visual selection, to constitute groups homogeneous for all traits separately.

Statistical analysis of physiological traits

All statistical analyses were performed with R 3.4.3 (R Core Team (2015). The significance
of differences between treatments were analyzed by one-way analysis of variance (t-test, n =
60). Differences were considered significant if *P* < 0.05.

A principal components analysis (PCA) was performed for the four traits used for selection
(TE, TE<sub>d</sub>, Wi<sub>insitu</sub>, Wi<sub>SS</sub>) and also for all traits that were measured (N=18).

RNA extraction, library construction and sequencing

We first carefully removed all leaf veins, to focus the analysis exclusively on the leaf blade.
RNA was then extracted as previously described (Le Provost et al. 2007). Total RNA was
recovered from the 12 selected genotypes (as described above) with a CTAB-based extraction
buffer. Residual genomic DNA was then removed from each sample with 10 U of RNAse-
free DNaseI (Promega, Madison, WI, USA) according to the manufacturer's instructions.
After DNase digestion, the RNA was purified with the cleanup protocol from the RNeasy
plant minikit (Qiagen, Valencia, CA, USA). RNA quantity and quality was assessed for each
sample by both optical density measurements and electrophoresis on an Agilent 2100
Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). RNA integrity number
(RIN) values ranged from 7 to 8.5 for the 12 samples. We then generated cDNA libraries
from 2 µg of RNA according to the Illumina Kit (TruSeq Stranded mRNA Sample Prep Kit)
in which the Poly-A containing mRNA molecules are purified using poly-T oligo attached
magnetic beads We sequenced each library with 150 base reads, in a paired-end flow cell, on an Illumina HiSeq2000 (Illumina, San Diego, CA, USA). Over 50 million useable reads were generated for each library (Supplementary Table 1). The sequences were produced by Genewiz (Leipzig, Germany).

Trimming, mapping and identification of differentially expressed genes (DEGs)

Bioinformatic analyses were performed in a Galaxy web environment according to a published procedure (Le Provost et al. 2016). The quality of each batch of sequences was first assessed with FAsTQC, and low-quality reads (i.e. quality value<20) were removed. High-quality reads were then mapped onto the 25k pedunculate oak gene models (Plomion et al. 2018) with BWA V0.6.1 (Li and Durbin 2009). We finally selected only gene models covered by at least 60 reads (i.e. at least 5 reads per sample) for differential gene expression analysis. DEGs were identified with DESeq2 software (Love et al. 2014), with a p-value<0.05 after adjustment for multiple testing with a false discovery rate (FDR) of 5%. The effects of WUE (i.e. low vs. high WUE), treatment (i.e. control vs. drought) and their interaction were assessed in likelihood ratio tests implemented in the DESeq2 package. The WUE and treatment effects were investigated by comparing a model without interaction (M1) with two reduced models for the treatment (M2) and WUE (M3) effects. For estimation of the interaction effect, we compared M4 to M1.

$$M1: Y_{ijk} = \mu_1 + W_i + T_j + \varepsilon_{ijk}$$

where $W_i$ is the WUE effect (i= “high WUE” vs. “Low WUE”), $T_j$ is the treatment effect (j= “Control” vs. “drought”).

$$M2: Y_{jk} = \mu_2 + T_j + \varepsilon'_{jk}$$
$$M3: Y_{ik} = \mu_3 + W_i + \varepsilon''_{ik}$$

$$M4: Y_{ijk} = \mu + T_i + W_j + (T * W)_{ij} + \varepsilon''_{ijk}$$

Only three of the seven theoretical categories of genes displaying a single or combined effect (WUE, treatment or WUE x treatment interaction) included DEGs. Their annotations were recovered from the corresponding oak gene models (Plomion et al. 2018).

### Gene set and subnetwork enrichment analysis

GO term enrichment analysis was performed for the DEGs within each of the three categories, with topGO in R (Alexa 2010). We corrected the *P*-value for false discovery rate and considered ontology terms with a corrected *p*-value below 0.05 to be significantly enriched. Results were displayed for the first 10 ontologies only.

We then performed pathway enrichment analysis with Pathway Studio software (Pathway Studio®, Elsevier 2017), as previously described (Le Provost et al. 2016). We considered a *P*-value<0.05 to indicate significant enrichment in a biological process (BP) or a molecular function (MF) in our DEG datasets. Gene networks were then generated with the gene sets corresponding to the main and interaction effects.

### RT-qPCR validation

We quantified the expression of 21 DEGs. All the primer pairs used for the RT-qPCR were designed with Primer3 software (Rozen and Skaletsky 2000). All the genes analyzed are listed, together with their associated effects, in Supplementary Table 2. For each RT-qPCR reaction, we normalized the fluorescence data against two control genes (Qrob_P0530610 and
Qrob_P0426000) selected from the non-DEG and characterized by highly repeatable levels of expression (i.e. similar numbers of reads across the different samples analyzed). RT-qPCR was performed on a LightCycler480 Real-Time PCR detection system (Roche Applied Science, Penzberg, Germany) with standard PCR parameters, as previously described (Le Provost et al. 2016). The fluorescence data obtained were analyzed with StatQPCR software, using methods derived from the Genorm and qBASE programs. StatQPCR is available from the following URL: http://satqpcr.sophia.inra.fr/cgi/tool.cgi.

# Results

## *Ecophysiological variation*

A principal components analysis (PCA) was performed for four traits (TE, TE<sub>d</sub>, W<i>i</i><sub>insitu</sub>, W<i>i</i><sub>SS</sub>) used to select genotypes with extreme WUE phenotypic values and also for all traits (N=18) that were measured. For the four-trait PCA, the first axis explained 83.5%, and the second axis 12.4% of the observed variation. The first axis separated clearly the chosen genotypes in terms of WUE (Figure 1A), whereas the first and the second axis (upwards sloping diagonal) separated the drought groups. For the 18-trait PCA (Figure 1B), the first axis represented 51.9% and the second axis 15.9% of the observed variation. The first axis was strongly related to measurements of WUE (including δ<sup>13</sup>C), but also stomatal conductance measurements. The separation of drought groups (upwards sloping diagonal) corresponded mainly to whole plant traits such as CT and BM.

## *Drought stress effect*
TE and Wi were higher in drought conditions, whereas BM (dry biomass), CT (cumulative transpiration), A (net CO₂ assimilation), and gs (stomatal conductance to water vapor) were lower (Table 2). The effect of drought was not significant for all whole-plant traits measured (see Supplementary File 1): TE, underlying BM and CT were estimated for the whole growing period and also for the drought period alone. However, the non-destructive estimation of initial biomass as used for BM and TE was more precise (because the plants were smaller and more homogeneous) than estimates of biomass when the desired level of drought stress was reached (as used for BM<sub>d</sub> and TE<sub>d</sub>). However, TE, BM and CT, included a period during which the plants were small and not subjected to drought, whereas TE<sub>d</sub>, BM<sub>d</sub> and CT<sub>d</sub> covered only the drought period, but may potentially contain more noise due to a less precise initial biomass estimation. This may explain why TE<sub>d</sub> displayed no significant drought effect, whereas TE displayed significant drought effects, as did BM<sub>d</sub>, CT and CT<sub>d</sub>, and BM displayed a tendency at *p*<0.1. A weakly significant interaction was detected for TE<sub>d</sub> (*p* = 0.046), which was not confirmed by Tukey analysis (Supplementary File 1). We also found no significant effect of drought for BM standardized per unit leaf area (for the period of drought only and for the whole period) or CT standardized per unit leaf area, but a significant effect was found for standardized CT for the drought period only (Ct<sub>d,LA</sub>). Wi displayed a significant drought effect both *in situ* and in steady state, as did the corresponding stomatal conductances. A<sub>ss</sub> was also significantly lower under drought conditions, whereas only a tendency was observed for A<sub>insitu</sub>.

## WUE effect

The selection of sessile oak genotypes based on Wi (Wi<sub>in situ</sub> and Wi<sub>SS</sub>) and TE (TE and TE<sub>d</sub>) resulted in highly significant differences between the WUE phenotypic groups (Table 2;
Figure 1A). Furthermore, the δ¹³C of bulk leaf material, an independent estimator of integrated Wi, also differed significantly between the high and low WUE phenotypic groups. We investigated the traits underlying TE and Wi. Biomass accumulation (BM, for the whole period or for the period of drought only, and after standardization by leaf area) did not differ significantly between the high and low WUE phenotypic groups. Similarly, no differences were detected for net CO₂ assimilation rate (A). By contrast, cumulative transpiration (CT) and stomatal conductance (gs) displayed significant differences or strong tendencies in the same direction (CT<sub>d</sub>, gs<sub>ss</sub>), with plants with a high WUE characterized by a lower CT and gs.

### Overview of the cDNA libraries generated

We generated 12 cDNA libraries corresponding to three genotypes for each of the four groups obtained (2 WUE phenotypes x 2 watering regimes).

An overview of the libraries generated is available in Supplementary Table 1. The number of reads ranged from 55 million to 68 million, corresponding to 27 and 34 million paired-sequences, respectively. All the libraries had quality scores above 35. Quality criteria were applied (see methods section), and 26 to 32 million paired reads were finally mapped onto the oak gene models (Supplementary Table 1).

The overall structure of the studied genotypes was investigated by both PCA and hierarchical clustering analysis (Figure 2). Unsurprisingly, there was as much variability within as there was between the high and low WUE phenotypic groups. This first observation indicates that only subtle transcriptional changes drove the variation of WUE between the selected genotypes. This was confirmed by the very small number of DEGs displaying a significant WUE phenotypic effect, consistent with the oligogenic control of WUE (see discussion section). It has been also shown that a major driver of genetic diversity in WUE can be related
to stomatal density and distribution (Chaerle et al. 2005, Kim et al. 2010), and RNAseq
analysis was performed on the whole leaf blade, which consists of several cell types. Thus, if
the differences in gene expression between genotypes with different WUE values is due both
to regulation in specific leaf structures and differences in physiological functioning, our
sampling procedure may have simply diluted the signal, homogenizing the transcription
profiles pools of the two groups (high vs. low WUE) of three contrasting genotypes. Despite
having used three different genotypes in the WUE groups and bulk leaf material,
homogenizing the transcription profiles, the few DEGs identified were also validated by RT-
qPCR, confirming the power of high-throughput sequencing and the robustness of RNA-seq
for "finding needles in a haystack". More striking was the lack of a clear pattern regarding the
second structuring factor considered in the factorial design: water availability. The applied,
controlled, drought level had been chosen to result in a significant 50% reduction of stomatal
conductance (based on gs<sub>insitu</sub>, Table 2), inducing a significant physiological response, but
allowing continued growth. After reaching this drought level, it was maintained for 56 days
before the leaves for expression analysis were harvested. Therefore we expected to only
detect changes in gene expression which are related to long term, constitutive responses. The
extremely small number of DEGs associated with watering regime corroborated the absence
of pattern. The different genotypes used in each group might have responded with slightly
different strategies, inducing over time more within group variability. These findings suggest
that the transcriptome of the tissue studied was affected only very mildly by the intensity of
the long-term water deficit imposed.

Identification of differentially expressed genes (DEGs)
We used likelihood ratio tests to identify DEGs displaying significant WUE (high vs. low WUE phenotypes), treatment (i.e. watering regime), and interaction effects. This analysis was performed on 15,576 gene models covered by at least five reads in each biological replicate (i.e. 60 reads in total for one gene model). Overall, 92 genes were found to be differentially expressed, with an adjusted *p*-value<0.05. Each DEG presented a unique signature with respect to the main or interaction effects (Figure 3).

In total, 12 genes (Gene-set #1(WUE-set), Supplementary File 2) presented a significant “WUE” effect. These genes displayed only quantitative variation. Half were up-regulated in the high- or low-WUE phenotypic group. The six genes up-regulated in the low WUE phenotypic group were: NAD(P)-linked oxidoreductase super family protein (Qrob_P0056900.2), endoxyloglucan transferase 28 protein (Qrob_P0320680.2), RPL/BEL1-like homeodomain 9 (Qrob_P0256370.2), ALY4 protein (Qrob_P0532680.), nicotinate phosphoribosyltransferase protein (Qrob_P0162910.2) and MAH-like protein (Qrob_P0125520.2). The six genes up-regulated in the high WUE phenotypic group were: Qrob_P0389480, which is similar to an ankyrin protein, Qrob_P0162930.2, which encodes a protein kinase, Qrob_P0412020.2, which is similar to a DNA-binding protein, a homeodomain-like superfamily protein, Qrob_P0244690.2, which encodes receptor-like protein 15, Qrob_P0489730.2, which is similar to a G-type lectin S-receptor-like serine/threonine-protein kinase, and Qrob_P0021870.2, which encodes a leucine-rich repeat receptor-like protein kinase family protein.

We also identified 27 genes (Gene-set #2 (Drought-set), Supplementary File 3) responding to the watering regime. In total, 17 genes were up-regulated in drought conditions, three of which (Qrob_P0268900.2, Qrob_P0104190.2 and Qrob_P0625010.2) were specifically expressed in these conditions. The other 10 genes were up-regulated in the control condition.
Finally, 53 genes presented an interaction effect (Gene-set #3 (interaction-set), Supplementary File 4) with clear and opposite patterns of expression between the low- and high-WUE phenotypic groups as a function of watering regime. The Kmeans algorithm in Expander software (Sharan et al. 2003) identified two opposite clusters: in cluster #1, 16 genes were up-regulated in low-WUE genotypes and down-regulated in high-WUE genotypes) with irrigation (i.e. control condition), whereas, in cluster #2, 37 genes displayed the opposite pattern (i.e. up-regulated in high-WUE genotypes and down-regulated in low-WUE genotypes with irrigation (i.e. control condition)).

Functional annotations for the three gene-sets were retrieved from the TAIR database.

## Analysis of GO term enrichment

We performed GO term enrichment analysis for the three gene-sets described above, corresponding to the two main and one interaction effects (available in Supplementary Files 2, 3 and 4). The highest degree of molecular function (MF) enrichment was observed for terms: (i) for Gene-set#1 (WUE-set) relating to cytokinin receptor activity, alpha-amylase activity and enone reductase activity, (ii) for Gene-set#2 (drought-set), relating to drug binding, organic cyclic compound binding and heterocyclic compound binding, and (iii) for Gene-set#3 (interaction-set), relating to quercetin 3-O-glucosyltransferase activity, UDP-glucosyltransferase activity and quercetin 7-O-glucosyltransferase activity. For biological process (BP) ontology, we identified terms relating to cellular polysaccharide catabolic processes and ion transport for Gene-set#1, response to external stimulus and cellular processes for Gene-set#2, and defense response and regulation of raffinose metabolic processes for Gene-set#3.
GO term enrichment analysis was also performed independently for the two contrasting clusters obtained from Gene-set#3. We identified significant enrichment only for cluster#1 (Supplementary File 4) for terms relating to quercetin 3-O-glucosyltransferase activity, UDP-glucosyltransferase activity and glucosyltransferase activity. Enrichment was observed for MF terms relating to the negative regulation of protein catabolic process, and the regulation of vesicle fusion, and for BP terms relating to the regulation of SNARE complex assembly.

## Subnetwork enrichment analysis

We performed subnetwork enrichment analysis independently for the three gene-sets, with Pathway Studio software. Results are also available in Supplementary Files 2, 3 and 4.

For Gene-set#1 (WUE-set effect), we identified three genes interacting with at least one of the four cellular processes identified (flower development, flowering time, plant development, cell differentiation): a RPL gene (Qrob_P0244690.2) interacting with all four of these processes, an ALY4 gene (Qrob_P0532680.2) and a XTH28 (Qrob_P0320680.2) gene interacting with three of the processes (flower development, plant development and flowering time).

For genes displaying a significant effect of treatment (Gene-set#2, drought-set), we identified six DEGs known to be involved in five BP. Four of these DEGs (FER (Qrob_P0672440.2), FRK1 (Qrob_P0625010.2), FRO6 (Qrob_P0222260.2) and Akinbeta1 (Qrob_P0268900.2)) were up-regulated in the control treatment and regulated three of the five BP identified in this analysis (sugar response, sugar metabolism and chlorophyll content). Two genes (AVP1, (Qrob_P0221580.2) and WOL, (Qrob_P0247500.2)) were over-expressed in the drought treatment and interacted with three BP (turgor, organ formation and chlorophyll content).
We focused in particular on the DEGs displaying a significant WUE-by-treatment interaction effect (Gene-set #3, interaction-set), revealing different molecular strategies for dealing with water constraints in low- and high-WUE genotypes. Subnetwork enrichment analysis was performed for the whole dataset, and separately for clusters #1 and #2. Unfortunately, no significant enrichment was identified for cluster #2. We therefore discuss below only the network obtained for the 53 genes of cluster #1 displaying a significant WUE-by-treatment interaction. We identified nine DEGs (RLK1 (Qrob_P0430150.2), CPR30 (Qrob_P0255490.2), SAG101 (Qrob_P0003420.2), ACA1 (Qrob_P0090430.2), ERH3 (Qrob_P0303640.2), T13D8.29 (Qrob_P0542920.2), VPS41 (Qrob_P0111100.2), UGT74D1 (Qrob_P0574050.2) and CNCG1 (Qrob_P0085080.2)) known to be involved in six BP (regulation of cell shape, plant immunity, cell morphogenesis, photosynthesis, pollen fertility and positive gravitropism). Four of the genes identified (SAG101, 1CA1, CPR30 and ERH3) are involved in regulation of the biological process photosynthesis, and three (RLK1, CPR30 and SAG101) are involved in plant immunity. All the other biological processes identified (pollen fertility, regulation of cell shape, positive gravitropism and cell morphogenesis) are known to interact with two of our DEGs.

*Gene expression – phenotypic trait correlation analysis*

A correlation analysis was done between all the ecophysiological traits available in Table 2 and the levels of expression of all genes in the three gene-sets, WUE-set (gene-set#1), Drought-set (Gene-set#2) and Interaction-set (Gene-set#3). Given the small number of data points in each regression (N=12) and the large number of correlations tested (18 traits and 92 genes, i.e. 1,656 correlations) a stringent type I error of 0.0005 was applied. Eight significant correlations were detected, for six different genes (see Supplementary File 5), all with a
correlation coefficient above 0.85: Qrob_P0056900.2 was negatively correlated with Wi<sub>SS</sub> and positively correlated with gs<sub>SS</sub>, Qrob_P0532680.2 was also negatively correlated with Wi<sub>insitu</sub> and gs<sub>insitu</sub>. Qrob_P0162910.2 and Qrob_P0320680.2 were negatively correlated with TE, whereas Qrob_P0412020.2 was positively correlated with TE<sub>d</sub>. All of these genes belong to Gene-set#1 (WUE-set). Qrob_P0226850.2 was the only gene from Gene-set#2 (Drought-set) for which a negative correlation was detected with Wi<sub>insitu</sub> and no significant correlation was detected for Gene-set#3 (Interaction-set). All of these genes displayed weaker significant correlations (p<0.05, Supplementary File 5) with at least six other traits from Table 2, with the exception of traits related to photosynthesis, for which no significant correlations were detected.

## Real-time quantitative PCR validation

For confirmation of the accuracy of RNA-seq results, we performed RT-qPCR analyses for 21 genes: six displaying significant WUE effects and 15 displaying WUE-by-treatment effects. Six of the genes displaying significant phenotype-by-treatment interaction effects belonged to Cluster #1 and nine belonged to Cluster #2. We excluded genes displaying either a multiple-band amplification pattern on agarose gels or inconsistent PCR efficiency. We then quantified the expression of the remaining 16 genes by RT-qPCR (5 genes for the phenotype effect, 3 for cluster #1 and 7 for cluster #2, Supplementary Table 2). The patterns of expression of the genes tested by RT-qPCR were similar to those obtained by RNAseq (Supplementary Figure 1), indicating that both the RNA-seq data and the bioinformatics procedure were reliable and that the pattern of DEG expression revealed by RNA-seq could be used for further biological interpretation.
Discussion

The molecular mechanisms involved in regulating WUE are well documented for crops (Ruggiero et al. 2017).. However, the elucidation of the molecular mechanisms regulating WUE remains in its infancy for forest trees. In this study, an analysis of differential gene expression enabled us to identify three sets of genes. Below, we focus mostly on the genes of Gene-set#1 (i.e.showing a WUE effect, G, Supplementary File2) and Gene-set #3(i.e displaying a treatment by WUE interaction effect, G*E, Supplementary File4), which encode proteins characterizing the molecular strategies of low- and high-WUE genotypes for coping with drought.

Our differential gene expression analysis was performed on the whole leaf blade and on different genotypes (i.e. the four WUE by watering treatment combinations are represented by three independent genotypes, each, with a specific genetic background). We have taken care that seedlings analysed in this study belong to the same recolinization lineage ( i.e. minimizing the variation of the genetic background) and that no provenance effect was identified for the selected populations (data not shown). Finally, it should be noticed that Wi is a variable trait over time and it is therefore difficult to characterize a consistent phenotype (i.e. High WUE efficiency for example) which could be included in both water regimes (well watered vs. non irrigate). For this reason, we have made the choice to select the individuals according to the ecophysiological data, combining instant and time integrated measurements of WUE, resulting in twelve different genotypes to be analysed in this study. This strategy may account for the relatively small number of DEGs identified, however it enlarges the validity towards the natural variability in this species.. One previous study (Brendel et al. 2008) reported an oligogenic genetic determinism of WUE in a pedunculate oak family, mainly driven by variation in stomatal conductance (Roussel et al. 2009a), thus suggesting that variation of WUE in oak species is influenced by a small number of genes, each with a
large effect, consistent with our observations for sessile oak (see the high fold-change ratios below). The subnetwork enrichment analysis identified key genes of potential importance for sessile oak WUE (see below), demonstrating that the RNaseq approach is a method of choice for accelerating gene discovery for fine ecophysiological traits (reviewed in Nguyen et al. (2019)), even in the presence of confounding effects.

### (i) Genes regulated by soil water deficiency

In total, 27 genes in Gene-set#2 (Drought-set) were found to correspond to the effect of soil water deficit, and the subnetwork enrichment analysis identified BPs (turgor, organ formation, chlorophyll content and sugar metabolism) known to be regulated by drought stress (Hsiao, 1990). The overexpression of two genes (WOL, FC=0.22 and AVP1, FC=0.28) regulating turgor in plants is consistent with several studies that have already reported that cell turgor regulation during water depletion is essential to cope with stressful conditions (Kim et al. 2010). We also observed a decrease in transcript accumulation for four genes involved in sugar metabolism, sugar response and chlorophyll content, probably due to the decrease in primary metabolism classically observed during the abiotic stress response (Du et al. 2020). We found a highly significant correlation between Qrob_P0226850.2 (FC=2.24), which encodes an abscisic acid-responsive gene (AT1G02130), and intrinsic water use efficiency, highlighting the long-term coordination of the drought response by ABA at leaf level.

### (ii) Genes regulated in high and low WUE Genotypes

Of considerable relevance to our research question, the second set included 12 genes constitutively expressed at higher levels in either low- or high-WUE phenotypic group..
In total, six genes were significantly correlated with TE, Wi or stomatal conductance (Supplementary File 5), which were all also detected for constitutive expression differences in Gene-set#1. These genes were not correlated with assimilation rate, suggesting that the genetic diversity of WUE observed in these sessile oak genotypes is driven principally by variation in stomatal conductance (Roussel et al. 2009b, Roussel et al., 2009a).

Six genes up-regulated in high WUE individuals were considered of particular interest: (i) the receptor-like protein 15 gene (RPL15) which presented a 4.9-fold change in expression. RPL genes encode cell-surface receptors involved primarily in defense against pathogens. Some RPL genes are also known to regulate stomatal density and distribution by disrupting their patterning. For example, the too many mouth (TMM) gene is highly similar to RLP16 (Wang et al. 2008). Other authors have identified TMM as a key gene for improving the WUE of crops through molecular genetics approaches (Ruggiero et al. 2017); (ii) we also identified a leucine-rich repeat receptor-like protein kinase gene (validated by RT-qPCR) displaying a 5.2-fold change in expression. Leucine-rich repeat receptor-like protein kinases are receptor kinases located on the surface of the plant cell and involved in the perception of signals from the environment. A study in *Arabidopsis* by (Does et al. 2017) reported a key role for this gene in both regulating several aspects of growth and the response to abiotic stress; (iii) we identified a G-type lectin S-receptor-like serine/threonine protein kinase gene (validated by RT-qPCR) with a 4.8-fold change in expression. The overexpression of this gene in *Arabidopsis* (Zhao et al. 2013) enhances salinity and drought stress tolerance; (iv) a homeodomain-like superfamily protein gene, encoding a DNA-binding protein, was also identified, with a 2.02-fold change in expression, for which expression was significantly positively correlated with TE<sub>d</sub> (and, to a lesser extent, with other estimators of WUE and whole-plant transpiration). In *Quercus lobata*, it has been suggested that such genes may underlie the differences in WUE between oak populations (Gugger et al. 2017); and (v)
finally, we identified an ankyrin protein gene with a 4.7-fold change in expression. Ankyrin repeat proteins play essential roles in cell growth, development and the response to hormones and environmental stresses (Garcion et al. 2006). It has recently been reported (Zhao et al. 2020) that some ankyrin genes confer resistance to drought and salinity in *Arabidopsis*.

Five of the six genes up-regulated in low-WUE genotypes (or down-regulated in high-WUE genotypes) were also considered of considerable interest based on published findings or the results of our subnetwork analysis:

(i) an MAH-like protein gene (FC=-0.83) potentially involved in wax and cutin biosynthesis (Sakai et al. 2018). It has been suggested that changes in cutin thickness may be associated with WUE variation in pine, through effects on transpiration rate (Bert et al. 2021). However, no clear relationship is available in the literature on the influence of cuticle thickness on WUE; (ii) the NAPTR1 gene (FC=-0.32), which encodes a nicotinate phosphoribosyltransferase known to play a key role in both the biosynthesis and homeostasis of NAD (Hashida et al. 2009). NAD homeostasis has recently been shown to be a key component in stomatal movement (Hashida et al. 2009). The authors also showed, with transgenic plants, that mutations of this gene led to the stomata partially losing their ability to close, resulting in higher drought stress sensitivity (Hashida et al. 2010). NAPTR1 expression was inversely correlated with TE (and, to a lesser extent, with other estimators of WUE, transpiration and stomatal conductance), suggesting a major role in stomatal conductance, driving differences in WUE. Similarly, the correlation analysis also showed an inverse correlation between the expression of Qrob_P0056900.2, which is similar to an NAD(P)-linked oxidoreductase superfamily protein (AT1G59960) and WUE, and positively correlation with gs<sub>SS</sub>. The three last genes were associated with the functional network (Supplementary File 2): (iii) ALY genes (FC=-0.26) encode key proteins involved in the transport of mRNA from the nucleus to the cytosol and are essential for correct plant growth
and development (Pfaff et al. 2018). It has also been shown (Teng et al. 2014) that the ALY4 gene is involved in stomatal closure in *Arabidopsis*. We also observed two highly significant and strong correlations between ALY4 expression and $W_i_{insitu}$ (inverse correlation) and to $g_{s_{insitu}}$ (positive correlation), suggesting that higher levels of expression for this gene resulted in greater stomatal conductance in low-WUE genotypes; (iv) the XTH28 gene ($F_c$=-0.71) encodes a xyloglucan endotransglycosylase (XTH) protein. In plants, the members of the XTH multigene family are involved in the metabolism of xyloglucan, an important cell-wall compound (Hara et al. 2014). In *Solanum lycopersicum*, a role for the product of the XTH3 gene in regulating cell wall remodeling in guard cells, and, potentially, in regulating transpiration rates during abiotic stress has been reported (Choi et al. 2011), suggesting a potential role of this gene in WUE. The significant correlations observed between the expression of this gene and traits related to stomatal conductance and cumulative water consumption suggest that a stronger expression of this gene will decrease WUE by increasing stomatal conductance and thereby the cumulative water consumption; (v) finally, we identified an RPL gene encoding a Bell-like homeodomain ($F_c$=-2.14). In rice, Bel-1 has been involved in the regulation both of organ size and water loss, through the regulation of stomatal density. In studies of mutations of Bel-1 genes (Gong et al. 2018), a key role was identified for the product of the RPL gene in cell proliferation, expansion and response to environmental clues. A higher stomatal density, leading to senescence associated with underdeveloped chloroplasts and low levels of photosynthesis, was observed in this mutant, consistent with a greater sensitivity to drought stress. Many of the genes down-regulated for high WUE individuals suggest strongly that long term, constitutive differences in WUE among sessile oak genotypes were linked to either differences in the molecular functioning of stomatal guard cells or the distribution of stomata on the leaf. This had already been suggested by (Roussel et al. 2009b) for pedunculate oak.
The identification of these DEGs highlights a series of genes for which up- or down-
regulation is associated with high WUE in sessile oaks regardless of the watering regime. An
overview of the genes identified is available in Figure 4, with the corresponding BPs, which
include stomatal density and distribution, cuticle biosynthesis, regulation of transpiration or
stomatal movement. This study thus validates current scientific knowledge regarding the
importance of water loss regulation in oak genotypes with a high water-use efficiency.

(iii) *Gene regulation triggering high WUE under drought conditions: toward adaptive molecular mechanisms combining drought tolerance and WUE*

The third set of genes, displaying GxE (i.e. corresponding to transcripts for which plasticity in expression was dependent on the level of WUE of the genotype studied), is of even greater interest in the context of this study, as these molecular players represent key candidate genes of importance in the long-term adaptation of oak to drought.

Two clusters were identified by the Kmeans approach. Cluster #1 contained 16 genes, which were down-regulated in the low-WUE genotypes and up-regulated in high-WUE genotype in the drought treatment. and Cluster #2 comprised (37 genes) with the opposite pattern of expression (i.e. up- and down-regulated in the low- and high-WUE genotypes, respectively, in the drought treatment). These two clusters of genes correspond to molecular mechanisms enabling oak trees to maintain a high WUE in conditions of water deficit.These genes show also that high and low WUE individuals of a same species can have opposite molecular drought response patterns, suggesting a within species diversity of drought responses (Brendel and Epron 2022). This gene pool may, therefore, correspond to genes of importance for the adaptation of forests to a drier environment, therefore worth for conserving their genetic
diversity. A description of the biological processes and associated genes identified is available
in Figure 5.

In cluster #1, we identified two genes validated by RT-qPCR (CES101, and TPS02) and known to be regulated by drought stress. The first, CES101, is highly similar to a lectin S-receptor-like serine/threonine protein kinase, the expression of which is modulated during salt stress in *Arabidopsis* (Sun et al. 2013). The second, TPS 02, encodes a terpene synthase. Monoterpenes are important aromatic molecules widespread in nature. The levels of monoterpene and terpene synthases increase massively in plants exposed to drought stress (Radwan et al. 2017). The up-regulation of a TPS gene in high-WUE genotypes may result in better drought tolerance. We identified eight other genes (AP2, CPR30, ERH3, ERF12, EN121, and three UGT genes (UGT73B5, UGT71C4, UGT 74D1) known to be involved in the regulation of stomatal movement or stomatal aperture by regulating IAA levels in guard cells (Jin et al. 2013). For example, AP2 is similar to aspartic protease 2, which is specifically expressed in guard cells. In *Arabidopsis*, overexpression of the AP2 gene resulted in a smaller mean stomatal aperture, suggesting a possible role for the AP2 gene in the fine regulation of transpiration rate (Guo et al. 2015) and, indirectly, in WUE. EN121, another of the genes identified, encodes a basic helix-loop-helix (bHLH) DNA-binding protein. In wheat, bHLH proteins are transcription factors known to be regulated by both ABA and drought stress. A transgenic line overexpressing the bHLH gene was found to display early stomatal closure in response to drought stress, enabling the plant to regulate its transpiration more effectively (Yang et al. 2016). We identified an ERF12 gene displaying similarity to ethylene-responsive factor 12 gene. A transgenic line of cotton overexpressing the ERF38 gene displayed both a larger guard cell stomatal aperture and a lack of sensitivity to ABA in terms of stomatal closure under drought stress (Ma et al. 2017), suggesting a key role of ERFs in regulating transpiration rate. The last two genes (CPR30 and ERH3) were identified by
functional network analysis (Supplementary File 3). The ERH3 gene is similar to a katanin P60 subunit protein. In *Arabidopsis*, an ERH3 loss-of-function mutant had small stomata, suggesting a possible role for the katanin protein in both cell differentiation and stomatal development (Bouquin et al. 2003). We hypothesized that the differential expression of genes involved in stomatal size would affect WUE by regulating both transpiration and assimilation rate. However, no significant differences in stomatal size were observed in our study. Finally, we identified a CP30 gene highly similar to an F box protein as differentially expressed. The CPR 30 gene was also found to be up-regulated in response to drought stress in another European white oak, *Quercus pubescens* (Madritsch et al. 2019). CPR30 is known to interact with ASK proteins, which promote stomatal opening, suggesting a possible functional role of CPR30 in the regulation of WUE (Li et al. 2012), by increasing the transpiration rate under drought stress. Last, we identified three UGT genes (UGT73B5, UGT71C4 and UGT74D1). UGT74D1 is a key gene involved in auxin glycosylation, a key mechanism regulating the levels of auxin in cells (Jin et al. 2013). Indole-3-acetic acid (IAA) is a key auxin hormone involved in the regulation of several biological processes. It has been reported to stimulate stomatal opening (ŘicÁnek and VicherkovÁ 1992). The regulation of IAA levels by genes in cluster #1 may explain the lower WUE, through the modulation of stomatal movements. However, further studies would be required to confirm this, as no IAA determinations were performed in this study.

Fourteen of the 37 genes in cluster #2 were considered of potential interest based on the available literature or on our network analysis. We first identified genes (4) known to improve resistance to drought stress: two CNGC1 genes (Qrob_P0085080.2 validated by RT-qPCR) and an SD1-13 (validated by RT-qPCR) gene known to be up-regulated during drought stress in *Arabidopsis* (Yang et al. 2010, Jha et al. 2016). We also identified a DUF
247 gene (validated by RT-qPCR) known to be differentially regulated in response to drought stress in *Quercus lobata* (Gugger et al. 2017).

The other 10 genes have been involved in WUE and/or in the regulation of stomatal movement. We identified an ACA1 gene encoding a carbonic anhydrase. ACA genes encode key enzymes involved in the transfer of CO2 to the catalytic site of ribulose 1.5 biphosphate carboxylase/oxygenase in the cell, which may act as important regulators of WUE in plants (Fabre et al. 2007). We identified two late embryogenesis abundant (LEA) proteins. Barley transgenic lines with LEA mutations perform better in terms of both biomass productivity and WUE under drought stress (Sivamani et al. 2000). We also identified two genes encoding RLK proteins (CrRLK and GbRLK, validated by RT-qPCR). In rice, CrRLK proteins are known to be involved in the regulation of the circadian clock and, potentially, in WUE (Nguyen et al. 2015). GbRLK proteins are key cell-surface receptors that may regulate stomatal movements in response to environmental cues in cotton (Jun et al. 2015). Indeed, several authors reported that loss of function mutants for some GbRLK genes were defective for stomatal closure (Hua et al. 2012). We identified an SHT protein in this cluster, which was up-regulated (validated by RT-qPCR). SHT proteins are involved in the polyamine biosynthesis pathway. Several authors have reported that polyamines are key components of the plant response to environmental variations. Indeed, they may be involved in stomatal movement through regulation of the voltage-dependent inward K+ channel in the plasma membrane of the guard cell (Liu et al. 2000). We identified an ABC2 gene (validated by RT-qPCR) encoding an ATP-binding cassette. In *Arabidopsis*, loss-of-function mutations affected WUE under drought stress, suggesting a key role for this gene in water uptake and loss (Lee et al. 2008, Ruggiero et al. 2017). We also identified a WAK2 gene encoding a serine/threonine kinase. In *Arabidopsis*, loss-of-function mutations of this serine/threonine kinase gene strongly decreased WUE, suggesting a possible role for this gene in the
regulation of stomatal movement (Fujii and Zhu 2009). Although the role of this gene under drought conditions is not very clear for high WUE genotypes, its upregulation in control conditions may explain the higher WUE observed. Finally, we identified two LST8-1 genes. LST8-1 is a partner of the target of rampamycin kinase involved in several biological processes in plants. The LST8-1 gene was recently reported to be strongly expressed in *Arabidopsis* guard cells, suggesting a possible role in stomatal regulation (Moreau et al. 2012).

## Conclusion

We first identified genes differentially regulated between individuals with contrasting WUE (low vs. high WUE, Gene-set#1). Low-WUE genotypes displayed a regulation of genes involved in cuticle remodeling or the rapidity of stomatal closure in response to drought stress, whereas high-WUE genotypes were characterized by a regulation of genes involved in stomatal density and distribution, or genes encoding surface receptors potentially involved in the regulation of transpiration. These genes may contribute to fine regulation of the tree's transpiration rate, improving WUE in these genotypes.

We also identified genes displaying a significant WUE-by-treatment interaction, highlighting different molecular strategies between low- and high-WUE genotypes for coping with a long-term soil water deficit. Brendel and Epron (2022) have suggested that a within species diversity of WUE might be linked to different drought response strategies. The expression profile of low-WUE genotypes responding to drought stress highlighted an important role for genes potentially involved in early stomatal closure or in increasing transpiration rate, whereas that obtained for high-WUE genotypes responding to drought highlighted an
important role for genes involved in the CO₂ flux to the carboxylation sites, in the regulation
of stomatal movements or in the water loss under drought stress.

This comprehensive investigation constitutes a first step toward understanding complex
molecular processes which play a role in trees leaf-level responses to a long term drought and
highlight a within species diversity of adaptive strategies. Although this study has already
focused on the natural diversity found within the species, using a detailed ecophysiological
characterization of the genotypes, the number of genotypes involved was too small for a
genetic analysis, Additional insights are emerging from ongoing investigations focusing on
the association between WUE and gene polymorphisms in a pedigree population and
association panels of unrelated sessile oak genotypes. Our future investigations planned in the
coming years will be to preform Whole Genome Sequencing Analysis in hundreds of sessile
oak genotypes harvested in a common garden experiment and to analyse the association
between the polymorphism identified within the candidate genes in this studt and the WUE
related traits.
**Declarations**

Availability of data and materials: The datasets generated and/or analyzed in this study are available from the SRA-NCBI repository under accession number [PRJNA763825](https://www.ncbi.nlm.nih.gov/SRA/?term=PRJNA763825). All the scripts used for bioinformatics analysis and the identification of differentially expressed genes are available from the corresponding author upon request.

**Conflict of interest:** none to declare

**Funding:** This work was supported by the ANR H2oak project (2014 14-CE02-0013-02, "Diversity for adaptive traits related to water use in two temperate European white oaks"). Illumina sequencing was performed by the Genewiz Company.

**Acknowledgments**

We thank the Genotoul bioinformatics facility in Toulouse for providing computing resources. RT-qPCR analysis was performed at the sequencing and genotyping facility of Bordeaux, supported by Conseil Régional d'Aquitaine grants nos. 20030304002FA and 20040305003FA, European Union FEDER grant no. 2003227 and Investissements d'Avenir ANR-10-EQPX-16-01. We also thank Cyril Buré for help with the preparation and the running of the experiment and the management of the robotic irrigation system. The author(s) would like to thank SILVATECH (Silvatech, INRAE, 2018. Structural and functional analysis of tree and wood Facility, doi: 10.15454/1.5572400113627854E12) from UMR 1434
SILVA, 1136 IAM, 1138 BEF and 4370 EA LERMAB from the research center INRAE
Grand-Est Nancy for the stable carbon isotope measurements. SILVATECH facility is
supported by the French National Research Agency through the Laboratory of Excellence
ARBRE (ANR-11-LABX-0002-01

**Authors' contributions:** GLP, TG and OB designed the study. GLP and CP wrote the manuscript with OB. GLP performed bioinformatics analysis. GLP was involved in the RT-qPCR experiments. OB and TG were responsible for running the experiment and the ecophysiological measurements and their statistical analysis. All the authors have read and approved the manuscript.
## References

* Alexa A (2010) Gene set enrichment analysis with topGO. http://www.mpi-sb.mpg.de/~alexa

* Aranda I, Alía R, Ortega U, Dantas ÂK, Majada J (2010) Intra-specific variability in biomass partitioning and carbon isotopic discrimination under moderate drought stress in seedlings from four Pinus pinaster populations. Tree Genetics & Genomes 6:169–178.

* Aranda I, Pardos M, Puértolas J, Jiménez MD, Pardos JA (2007) Water-use efficiency in cork oak (Quercus suber) is modified by the interaction of water and light availabilities. Tree Physiol 27:671–677.

* Bert D, Le Provost G, Delzon S, Plomion C, Gion J-M (2021) Higher needle anatomic plasticity is related to better water-use efficiency and higher resistance to embolism in fast-growing Pinus pinaster families under water scarcity. Trees 35:287–306.

* Bogeat-Triboulot MB, Buré C, Gerardin T, Chuste PA, Le Thiec D, Hummel I, Durand M, Wildhagen H, Douthe C, Molins A, Galmés J, Smith HK, Flexas J, Polle A, Taylor G, Brendel O (2019) Additive effects of high growth rate and low transpiration rate drive differences in whole plant transpiration efficiency among black poplar genotypes. Environmental and Experimental Botany 166:103784.

* Bouquin T, Mattsson O, Næsted H, Foster R, Mundy J (2003) The Arabidopsis luel mutant defines a katanin p60 ortholog involved in hormonal control of microtubule orientation during cell growth. Journal of Cell Science 116:791–801.

* Brendel O (2021) The relationship between plant growth and water consumption: a history from the classical four elements to modern stable isotopes. Annals of Forest Science 78:1–16. https://link.springer.com/article/10.1007/s13595-021-01063-2.

* Brendel O, Epron D (2022) Are differences among forest tree populations in carbon isotope composition an indication of adaptation to drought? Tree Physiology 42:26–31.

* Brendel O, Le Thiec D, Scotti-Saintagne C, Bodénès C, Kremer A, Guehl J-M (2008) Quantitative trait loci controlling water use efficiency and related traits in Quercus robur L. Tree Genetics & Genomes 4:263–278.

* Brendel O, Pot D, Plomion C, Rozenberg P, Guehl J-M (2002) Genetic parameters and QTL analysis of δ13C and ring width in maritime pine. Plant, Cell & Environment 25:945–953.

* Chaerle L, Saibo N, Van Der Straeten D (2005) Tuning the pores: towards engineering plants for improved water use efficiency. Trends in Biotechnology 23:308–315.

* Choi JY, Seo YS, Kim SJ, Kim WT, Shin JS (2011) Constitutive expression of CaXTH3, a hot pepper xyloglucan endotransglucosylase/hydrolase, enhanced tolerance to salt and drought stresses without phenotypic defects in tomato plants (Solanum lycopersicum cv. Dotaerang). Plant Cell Rep 30:867–877.
Does DV der, Boutrot F, Engelsdorf T, Rhodes J, McKenna JF, Vernhettes S, Koevoets I, Tintor N, Veerabagu M, Miedes E, Segonzac C, Roux M, Breda AS, Hardtke CS, Molina A, Rep M, Testerink C, Mouille G, Höfte H, Hamann T, Zipfel C (2017) The Arabidopsis leucine-rich repeat receptor kinase MIK2/LRR-KISS connects cell wall integrity sensing, root growth and response to abiotic and biotic stresses. PLOS Genetics 13:e1006832.

Du Y, Zhao Q, Chen L, Yao X, Zhang W, Zhang B, Xie F (2020) Effect of drought stress on sugar metabolism in leaves and roots of soybean seedlings. Plant Physiology and Biochemistry 146:1-12.

Ducousos A, Bordacs S (2004) EUFORGEN Technical Guidelines for genetic conservation and use for Pedunculate and sessile oaks (Quercus robur) and (Quercus petraea). Bioversity International.

Fabre N, Reiter IM, Becuwe-Linka N, Genty B, Rumeau D (2007) Characterization and expression analysis of genes encoding α and β carbonic anhydrases in Arabidopsis. Plant, Cell & Environment 30:617-629.

Fujii H, Zhu J-K (2009) Arabidopsis mutant deficient in 3 abscisic acid-activated protein kinases reveals critical roles in growth, reproduction, and stress. PNAS 106:8380-8385.

Garcion C, Guilleminot J, Kroj T, Parcy F, Giraudat J, Devic M (2006) AKRP and EMB506 are two ankyrin repeat proteins essential for plastid differentiation and plant development in Arabidopsis. The Plant Journal 48:895-906.

Gerardin T (2019) Plasticité et diversité de l'efficience d'utilisation de l'eau chez deux espèces de chêne blanc d'Europe : les chênes pédonculé (Quercus robur L.) et sessile (Quercus petraea (Matt.) Liebl.) : approche descriptive de la dynamique de réponse stomatique aux changements environnementaux. These de doctorat, Université de Lorraine. http://www.theses.fr/2019LORR0120.

Gong X, Zhang Z, Yue J, Tang W, Tang X, Zeng Z, Niu X, Chen D, Sang X, Xiao F, He G, Liu Y (2018) Phenotypic Characterization, Fine Mapping, and Altered Expression Profiling of Roses1 Mutation That Affects Organ Size and Water Loss Through Regulating Stomatal Density in Rice. Crop Science 58:486-506.

Guehl J-M, Fort C, Ferhi A (1995) Differential response of leaf conductance, carbon isotope discrimination and water-use efficiency to nitrogen deficiency in maritime pine and pedunculate oak plants. New Phytologist 131:149-157.

Gugger PF, Peñaloza-Ramírez JM, Wright JW, Sork VL (2017) Whole-transcriptome response to water stress in a California endemic oak, Quercus lobata. Tree Physiol 37:632-644.

Guo R, Zhao J, Wang X, Guo C, Li Z, Wang Y, Wang X (2015) Constitutive expression of a grape aspartic protease gene in transgenic Arabidopsis confers osmotic stress tolerance. Plant Cell Tiss Organ Cult 121:275-287.

Hara Y, Yokoyama R, Osakabe K, Toki S, Nishitani K (2014) Function of xyloglucan endotransglucosylase/hydrolases in rice. Annals of Botany 114:1309-1318.
Hashida S, Itami T, Takahashi H, Takahara K, Nagano M, Kawai-Yamada M, Shoji K, Goto F, Yoshihara T, Uchimiya H (2010) Nicotinate/nicotinamide mononucleotide adenyltransferase-mediated regulation of NAD biosynthesis protects guard cells from reactive oxygen species in ABA-mediated stomatal movement in Arabidopsis. Journal of Experimental Botany 61:3813–3825.

Hashida S, Takahashi H, Uchimiya H (2009) The role of NAD biosynthesis in plant development and stress responses. Annals of Botany 103:819–824.

Hsiao T Measurement of plant water status. In: Irrigation of Agricultural Crops. Steward, B.A., Nielsen, D.R, pp 243–249.

Hua D, Wang C, He J, Liao H, Duan Y, Zhu Z, Guo Y, Chen Z, Gong Z (2012) A Plasma Membrane Receptor Kinase, GHR1, Mediates Abscisic Acid- and Hydrogen Peroxide-Regulated Stomatal Movement in Arabidopsis. The Plant Cell 24:2546–2561.

Jeanneau M, Gerentes D, Foueillassar X, Zivy M, Vidal J, Toppan A, Perez P (2002) Improvement of drought tolerance in maize: towards the functional validation of the Zm-Asr1 gene and increase of water use efficiency by over-expressing C4–PEPC. Biochimie 84:1127–1135.

Jin S-H, Ma X-M, Han P, Wang B, Sun Y-G, Zhang G-Z, Li Y-J, Hou B-K (2013) UGT74D1 Is a Novel Auxin Glycosyltransferase from Arabidopsis thaliana. PLOS ONE 8:e61705.

Jun Z, Zhang Z, Gao Y, Zhou L, Fang L, Chen X, Ning Z, Chen T, Guo W, Zhang T (2015) Overexpression of GbRLK, a putative receptor-like kinase gene, improved cotton tolerance to Verticillium wilt. Scientific Reports 5:15048.

Jha S, Sharma M, K. Pandey G (2016) Role of Cyclic Nucleotide Gated Channels in Stress Management in Plants. Current Genomics 17:315–329.

Karaba A, Dixit S, Greco R, Aharoni A, Trijatmiko KR, Marsch-Martinez N, Krishnan A, Nataraja KN, Udayakumar M, Pereira A (2007) Improvement of water use efficiency in rice by expression of HARDY, an Arabidopsis drought and salt tolerance gene. Proceedings of the National Academy of Sciences 104:15270–15275.

Keenan RJ (2015) Climate change impacts and adaptation in forest management: a review. Annals of Forest Science 72:145–167.

Kim T-H, Böhmer M, Hu H, Nishimura N, Schroeder JI (2010) Guard Cell Signal Transduction Network: Advances in Understanding Abscisic Acid, CO2, and Ca2+ Signaling. Annual Review of Plant Biology 61:561–591.

Kolbe AR, Brutnell TP, Cousins AB, Studer AJ (2018) Carbonic Anhydrase Mutants in Zea mays Have Altered Stomatal Responses to Environmental Signals. Plant Physiology 177:980–989.

Lamy J-B, Bouffier L, Burlett R, Plomion C, Cochard H, Delzon S (2011) Uniform Selection as a Primary Force Reducing Population Genetic Differentiation of Cavitation Resistance across a Species Range. PLOS ONE 6:e23476.
Lata C, Prasad M (2011) Role of DREBs in regulation of abiotic stress responses in plants. Journal of Experimental Botany 62:4731–4748.

Le Provost G, Herrera R, Paiva JAP, Chaumeil P, Salin F, Plomion C (2007) A micromethod for high throughput RNA extraction in forest trees. Biol Res 40:291–297.

Le Provost G, Lesur I, Lalanne C, Da Silva C, Labadie K, Aury JM, Leple JC, Plomion C (2016) Implication of the suberin pathway in adaptation to waterlogging and hypertrophied lenticels formation in pedunculate oak (Quercus robur L.). Tree Physiology:tpw056.

Lee M, Choi Y, Burla B, Kim Y-Y, Jeon B, Maeshima M, Yoo J-Y, Martinoia E, Lee Y (2008) The ABC transporter AtABCB14 is a malate importer and modulates stomatal response to CO2. Nat Cell Biol 10:1217–1223.

Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754–1760.

Li C, Liu Z, Zhang Q, Wang R, Xiao L, Ma H, Chong K, Xu Y (2012) SKP1 is involved in abscisic acid signalling to regulate seed germination, stomatal opening and root growth in Arabidopsis thaliana. Plant, Cell & Environment 35:952–965.

Liu K, Fu H, Bei Q, Luan S (2000) Inward Potassium Channel in Guard Cells As a Target for Polyamine Regulation of Stomatal Movements. Plant Physiology 124:1315–1326.

Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15. http://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8.

Ma L, Hu L, Fan J, Amombo E, Khaldun ABM, Zheng Y, Chen L (2017) Cotton GhERF38 gene is involved in plant response to salt/drought and ABA. Ecotoxicology 26:841–854.

Madritsch S, Wischnitzki E, Kotrade P, Ashoub A, Burg A, Fluch S, Brüggemann W, Sehr EM (2019) Elucidating Drought Stress Tolerance in European Oaks Through Cross-Species Transcriptomics. G3 Genes|Genomes|Genetics 9:3181–3199.

Marguerit E, Bouffier L, Chancerel E, Costa P, Lagane F, Guehl J-M, Plomion C, Brendel O (2014) The genetics of water-use efficiency and its relation to growth in maritime pine. Journal of Experimental Botany 65:4757–4768.

Masle J, Gilmore SR, Farquhar GD (2005) The ERECTA gene regulates plant transpiration efficiency in Arabidopsis. Nature 436:866–870.

Moreau M, Azzopardi M, Clément G, Dobrenel T, Marchive C, Renne C, Martin-Magniette M-L, Taconnat L, Renou J-P, Robaglia C, Meyer C (2012) Mutations in the Arabidopsis Homolog of LST8/GβL, a Partner of the Target of Rapamycin Kinase, Impair Plant Growth, Flowering, and Metabolic Adaptation to Long Days. The Plant Cell 24:463–481.

Nguyen KL, Grondin A, Courtois B, Gantet P (2019) Next-Generation Sequencing Accelerates Crop Gene Discovery. Trends in Plant Science 24:263–274.
Nguyen Q-N, Lee Y-S, Cho L-H, Jeong H-J, An G, Jung K-H (2015) Genome-wide identification and analysis of Catharanthus roseus RLK1-like kinases in rice. *Planta* 241:603-613.

Petit RJ, Latouche-Hallé C, Pemonge M-H, Kremer A (2002) Chloroplast DNA variation of oaks in France and the influence of forest fragmentation on genetic diversity. *Forest Ecology and Management* 156:115-129.

Pfaff C, Ehrnsberger HF, Flores-Tornero M, Sørensen BB, Schubert T, Längst G, Griesenbeck J, Sprunck S, Grasser M, Grasser KD (2018) ALY RNA-Binding Proteins Are Required for Nucleocytosolic mRNA Transport and Modulate Plant Growth and Development. *Plant Physiology* 177:226-240.

Plomion C, Aury J-M, Amselem J, Leroy T, Murat F, Duplessis S, Faye S, Francillonne N, Labadie K, Provost GL, Lesur I, Bartholomé J, Faivre-Rampant P, Kohler A, Leplé J-C, Chantret N, Chen J, Diévart A, Alaeitabar T, Barbe V, Belser C, Berges H, Bodénès C, Bogeat-Triboulot M-B, Bouffaud M-L, Brachi B, Chancerel E, Cohen D, Couloux A, Silva CD, Dossat C, Ehrenmann F, Gaspin C, Grima-Pettenati J, Guichoux E, Hecker A, Herrmann S, Hugueney P, Hummel I, Klopp C, Lalanne C, Lascoux M, Lasserre E, Lemainque A, Desprez-Loustau M-L, Luyten I, Madoui M-A, Mangenot S, Marchal C, Maumus F, Mercier J, Michotey C, Panaud O, Picault N, Rouhier N, Rué O, Rustenholz C, Salin F, Soler M, Tarkka M, Velt A, Zanne AE, Martin F, Wincker P, Quesneville H, Kremer A, Salse J (2018) Oak genome reveals facets of long lifespan. *Nature Plants* 4:440-452.

Plomion C, Bartholomé J, Bouffier L, Brendel O, Cochard H, De Miguel M, Delzon S, Gion JM, González-Martínez SC, Guehl JM, Lagraulet H, Le Provost G, Marguerit E, Porté A (2016) Understanding the genetic bases of adaptation to soil water deficit in trees through the examination of water use efficiency and cavitation resistance: maritime pine as a case study Cochard H (ed). *Journal of Plant Hydraulics* 3:e008.

Ponton S, Dupouey J-L, Bréda N, Dreyer E (2002) Comparison of water-use efficiency of seedlings from two sympatric oak species: genotype x environment interactions. *Tree Physiol* 22:413-422.

Radwan A, Kleinwächter M, Selmar D (2017) Impact of drought stress on specialised metabolism: Biosynthesis and the expression of monoterpene synthases in sage (Salvia officinalis). *Phytochemistry* 141:20-26.

ŘicÁnek M, VicherkovÁ M (1992) Stomatal responses to ABA and IAA in isolated epidermal strips of Vicia faba L. *Biol Plant* 34:259.

Roussel M, Dreyer E, Montpied P, Le-Provost G, Guehl J-M, Brendel O (2009a) The diversity of 13C isotope discrimination in a Quercus robur full-sib family is associated with differences in intrinsic water use efficiency, transpiration efficiency, and stomatal conductance. *Journal of Experimental Botany* 60:2419-2431.

Roussel M, Le Thiec D, Montpied P, Ningre N, Guehl J-M, Brendel O (2009b) Diversity of water use efficiency among Quercus robur genotypes: contribution of related leaf traits. *Ann For Sci* 66:408-408.
Rozen S, Skaletsky H (2000) Primer3 on the WWW for general users and for biologist programmers. Methods Mol Biol 132:365-386.

Ruggiero A, Punzo P, Landi S, Costa A, Van Oosten MJ, Grillo S (2017) Improving Plant Water Use Efficiency through Molecular Genetics. Horticulturae 3:31.

Sakai K, Taconnat L, Borrega N, Yansouni J, Brunaud V, Paysant-Le Roux C, Delannoy E, Martin Magniette M-L, Lepiniec L, Faure JD, Balzergue S, Dubreucq B (2018) Combining laser-assisted microdissection (LAM) and RNA-seq allows to perform a comprehensive transcriptomic analysis of epidermal cells of Arabidopsis embryo. Plant Methods 14:10.

Sharan R, Maron-Katz A, Shamir R (2003) CLICK and EXPANDER: a system for clustering and visualizing gene expression data. Bioinformatics 19:1787-1799.

Sivamani E, Bahieldin A, Wraith JM, Al-Niemi T, Dyer WE, Ho T-HD, Qu R (2000) Improved biomass productivity and water use efficiency under water deficit conditions in transgenic wheat constitutively expressing the barley HVA1 gene. Plant Science 155:1-9.

Sun X-L, Yu Q-Y, Tang L-L, Ji W, Bai X, Cai H, Liu X-F, Ding X-D, Zhu Y-M (2013) GsSRK, a G-type lectin S-receptor-like serine/threonine protein kinase, is a positive regulator of plant tolerance to salt stress. Journal of Plant Physiology 170:505-515.

Teng W, Zhang H, Wang W, Li D, Wang M, Liu J, Zhang H, Zheng X, Zhang Z (2014) ALY proteins participate in multifaceted Nep1Mo-triggered responses in Nicotiana benthamiana and Arabidopsis thaliana. J Exp Bot 65:2483-2494.

Tognetti R, Michelozzi M, Lauteri M, Brugnoli E, Giannini R (2011) Geographic variation in growth, carbon isotope discrimination, and monoterpene composition in Pinus pinaster Ait. provenances. Canadian Journal of Forest Research. https://cdnsciencepub.com/doi/abs/10.1139/x00-096 (24 February 2021, date last accessed).

Wang G, Ellendorff U, Kemp B, Mansfield JW, Forsyth A, Mitchell K, Bastas K, Liu C-M, Woods-Tör A, Zipfel C, Wit PJGM de, Jones JDG, Tör M, Thomma BPHJ (2008) A Genome-Wide Functional Investigation into the Roles of Receptor-Like Proteins in Arabidopsis. Plant Physiology 147:503-517.

Welander NT, Ottosson B (2000) The influence of low light, drought and fertilization on transpiration and growth in young seedlings of Quercus robur L. Forest Ecology and Management 127:139-151.

Weng H, Yoo CY, Gosney MJ, Hasegawa PM, Mickelbart MV (2012) Poplar GTL1 Is a Ca2+/Calmodulin-Binding Transcription Factor that Functions in Plant Water Use Efficiency and Drought Tolerance. PLoS One 7. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3292583/.

Xing HT, Guo P, Xia XL, Yin WL (2011) PdERECTA, a leucine-rich repeat receptor-like kinase of poplar, confers enhanced water use efficiency in Arabidopsis. Planta 234:229-241.
Yang L, Ji W, Zhu Y, Gao P, Li Y, Cai H, Bai X, Guo D (2010) GsCBRLK, a calcium/calmodulin-binding receptor-like kinase, is a positive regulator of plant tolerance to salt and ABA stress. Journal of Experimental Botany 61:2519–2533.

Yang T, Yao S, Hao L, Zhao Y, Lu W, Xiao K (2016) Wheat bHLH-type transcription factor gene TabHLH1 is crucial in mediating osmotic stresses tolerance through modulating largely the ABA-associated pathway. Plant Cell Rep 35:2309–2323.

Yoo CY, Pence HE, Jin JB, Miura K, Gosney MJ, Hasegawa PM, Mickelbart MV (2010) The *Arabidopsis* GTL1 Transcription Factor Regulates Water Use Efficiency and Drought Tolerance by Modulating Stomatal Density via Transrepression of *SDD1*. Plant Cell 22:4128–4141.

Zhang J, Marshall JD (1994) Population differences in water-use efficiency of well-watered and water-stressed western larch seedlings. Can J For Res 24:92–99.

Zhao J, Gao Y, Zhang Z, Chen T, Guo W, Zhang T (2013) A receptor-like kinase gene (GbRLK) from Gossypium barbadense enhances salinity and drought-stress tolerance in Arabidopsis. BMC Plant Biology 13:110.

Zhao J-Y, Lu Z-W, Sun Y, Fang Z-W, Chen J, Zhou Y-B, Chen M, Ma Y-Z, Xu Z-S, Min D-H (2020) The Ankyrin-Repeat Gene GmANK114 Confers Drought and Salt Tolerance in Arabidopsis and Soybean. Front Plant Sci 11. https://www.frontiersin.org/articles/10.3389/fpls.2020.584167/full.
Figures

**Figure 1:** Biplot (individuals and variable vectors) of the PCA analysis of a) the four traits used for genotype selection (TE, TE_d, Wi insitu, Wi ss) and b) all of the traits presented in Table 2. The four groups are represented as confidence ellipses around the group mean.

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

<figure><img src="image_3.png" /></figure>
**Figure 2:** PCA analysis (panel A) and clustering analysis (using the “complete” method of the expender software, panel B) among the 12 different libraries.

<figure><img src="image_4.png" /></figure>
**Figure 3:** Venn diagram showing the overlap between the three effects analyzed in this study. We used an adjusted p-value<0.05 to declare a gene differentially expressed. Abbreviation correspond to: G: WUE effect (i.e. phenotypic group), E: Treatment effect (i.e. water regime) and G\*E: WUE-by-treatment interaction effect.

<figure><img src="image_5.png" /></figure>
**Figure 4:** Illustration of the main biological processes and associated genes regulated by WUE in sessile oaks. Genes labelled in red and blue are respectively up- or down-regulated in high WUE genotypes. +: positive action of the gene, -: negative action.

<figure><img src="image_6.png" /></figure>
**Figure 5:** Illustration of the main biological processes and associated genes triggering high WUE under drought stress in sessile oak. Genes labelled in red belong cluster #1 (gene over-expressed in drought condition in high WUE genotypes), while those in blue are specific to cluster #2 (gene under-expressed in drought condition in high WUE genotypes). +: positive action, -: negative action.

<figure><img src="image_7.png" /></figure>
<table><thead><tr><th>Forest Name</th><th>Latitude</th><th>Longitude</th><th>T annual mean</th><th>T summer mean</th><th>P annual sum</th><th>P summer sum</th></tr></thead><tbody><tr><td>Cuverville</td><td>49°57'37" N</td><td>1°24'00" E</td><td>10.3</td><td>14.2</td><td>872.0</td><td>119.6</td></tr><tr><td>Les Hogues</td><td>49°25'42" N</td><td>1°23'49" E</td><td>10.0</td><td>14.9</td><td>769.4</td><td>124.2</td></tr><tr><td>Belleme</td><td>48°22'33" N</td><td>0°31'27" E</td><td>10.4</td><td>15.4</td><td>777.1</td><td>120.2</td></tr><tr><td>Boulogne</td><td>47°32'23" N</td><td>1°30'01" E</td><td>11.1</td><td>16.6</td><td>688.4</td><td>120.1</td></tr><tr><td>Loches</td><td>47°09'46" N</td><td>1°09'59" E</td><td>11.4</td><td>16.6</td><td>722.4</td><td>120.6</td></tr><tr><td>Chateauroux</td><td>46°44'19" N</td><td>1°43'38" E</td><td>11.5</td><td>16.8</td><td>803.3</td><td>136.6</td></tr></tbody></table>

**Table 1:** Forests from which the seeds were harvested and used in this study. Temperature (T) in °C and precipitation are shown as annual or summer (mai, june and july), means and sums, respectively. Data are from the SAFRAN meteorological analysis system<fnref n="1" />.

<footnote n="1">¹Quintana-Seguí, P., Le Moigne, P., Durand, Y., Martin, E., Habets, F., Baillon, M., Canellas, C., Franchisteguy, L., Morel, S., 2008. Analysis of near-surface atmospheric variables: Validation of the SAFRAN analysis over France. J. Appl. Meteorol. Climatol. 47, 92–107. https://doi.org/10.1175/2007JAMC1636.1</footnote>
<table><caption>**Table 2:** Results of the two-way ANOVA. Means and their standard errors (SEM) for the principal factors tested (Tmt) and their interaction. Significant differences between groups were identified in Tukey's highly significant differences test.</caption><thead><tr><th rowspan="2">Tmt</th><th rowspan="2">p</th><th rowspan="2">Phe</th><th rowspan="2">Int</th><th>Control</th><th>Drought</th><th>High</th><th>Low</th><th rowspan="2">Selection traits</th></tr><tr><th>mean ± SEM</th><th>mean ± SEM</th><th>WUE mean</th><th>WUE ± SEM</th><th>mean ± SEM</th></tr></thead><tbody><tr><td>TE</td><td>*</td><td>***</td><td></td><td>5.50±0.43</td><td>b</td><td>5.99±0.31</td><td>a</td><td>6.54±0.12</td><td>a</td><td>4.95±0.20</td><td>b</td><td>S</td></tr><tr><td>TE<sub>d</sub></td><td></td><td>***</td><td>*</td><td>4.91±0.45</td><td>a</td><td>5.11±0.28</td><td>a</td><td>5.77±0.15</td><td>a</td><td>4.24±0.16</td><td>b</td><td>S</td></tr><tr><td>BM</td><td>+</td><td></td><td></td><td>145.2±7.9</td><td>a</td><td>115.1±12.2</td><td>a</td><td>130.8±14.0</td><td>a</td><td>129.5±10.3</td><td>a</td><td></td></tr><tr><td>BM<sub>d</sub></td><td>*</td><td></td><td></td><td>116.5±8.5</td><td>a</td><td>79.4±7.2</td><td>b</td><td>100.7±12.3</td><td>a</td><td>95.1±10.3</td><td>a</td><td></td></tr><tr><td>CT</td><td>*</td><td>*</td><td></td><td>27.4±2.9</td><td>a</td><td>19.3±2.0</td><td>b</td><td>19.9±2.0</td><td>b</td><td>26.8±3.2</td><td>a</td><td></td></tr><tr><td>CT<sub>d</sub></td><td>**</td><td>+</td><td></td><td>24.5±2.5</td><td>a</td><td>15.7±1.6</td><td>b</td><td>17.4±2.0</td><td>a</td><td>22.9±3.1</td><td>a</td><td></td></tr><tr><td>BM<sub>LA</sub></td><td></td><td></td><td></td><td>0.0481±0.0025</td><td>a</td><td>0.0466±0.0046</td><td>a</td><td>0.0439±0.0029</td><td>a</td><td>0.0507±0.0038</td><td>a</td><td></td></tr><tr><td>BM<sub>d.LA</sub></td><td></td><td></td><td></td><td>0.0383±0.0018</td><td>a</td><td>0.0324±0.0031</td><td>a</td><td>0.0340±0.0033</td><td>a</td><td>0.0367±0.0021</td><td>a</td><td></td></tr><tr><td>CT<sub>LA</sub></td><td></td><td>**</td><td></td><td>9.05±0.90</td><td>a</td><td>7.93±0.99</td><td>a</td><td>6.70±0.38</td><td>b</td><td>10.28±0.70</td><td>a</td><td></td></tr><tr><td>CT<sub>d.LA</sub></td><td>*</td><td>**</td><td></td><td>8.10±0.74</td><td>a</td><td>6.49±0.85</td><td>b</td><td>5.85±0.45</td><td>b</td><td>8.74±0.69</td><td>a</td><td></td></tr><tr><td>d13C</td><td>*</td><td></td><td></td><td>-30.3±0.4</td><td>a</td><td>-30.2±0.4</td><td>a</td><td>-29.6±0.2</td><td>a</td><td>-30.8±0.4</td><td>b</td><td></td></tr><tr><td>Wi</td><td>***</td><td>*</td><td></td><td>58.8±9.0</td><td>b</td><td>108.1±9.9</td><td>a</td><td>102.4±12.8</td><td>a</td><td>64.5±10.7</td><td>b</td><td>S</td></tr><tr><td>A<sub>insitu</sub></td><td>+</td><td></td><td></td><td>14.8±1.3</td><td>a</td><td>11.4±1.3</td><td>a</td><td>12.5±1.7</td><td>a</td><td>13.7±1.1</td><td>a</td><td></td></tr><tr><td>gs<sub>insitu</sub></td><td>**</td><td>**</td><td></td><td>0.301±0.048</td><td>a</td><td>0.155±0.034</td><td>b</td><td>0.157±0.038</td><td>b</td><td>0.299±0.046</td><td>a</td><td></td></tr><tr><td>Wi<sub>SS</sub></td><td>**</td><td>***</td><td></td><td>76.5±13.7</td><td>b</td><td>110.9±9.1</td><td>a</td><td>116.2±6.6</td><td>a</td><td>71.2±11.9</td><td>b</td><td>S</td></tr><tr><td>A<sub>SS</sub></td><td>*</td><td></td><td></td><td>14.7±1.6</td><td>a</td><td>9.4±1.4</td><td>b</td><td>11.4±2.0</td><td>a</td><td>12.8±1.8</td><td>a</td><td></td></tr><tr><td>gs<sub>SS</sub></td><td>*</td><td>+</td><td></td><td>0.267±0.077</td><td>a</td><td>0.094±0.019</td><td>b</td><td>0.107±0.024</td><td>a</td><td>0.254±0.081</td><td>a</td><td></td></tr><tr><td>SD</td><td></td><td></td><td></td><td>567.7±25.8</td><td>a</td><td>606.8±28.0</td><td>a</td><td>608.7±34.8</td><td>a</td><td>565.8±14.2</td><td>a</td><td></td></tr></tbody></table>

TE: Transpiration efficiency [g/kg]; BM: dry biomass [g]; CT: cumulative transpiration [kg]; d13C: carbon isotope composition; Wi: intrinsic water use efficiency; A: net CO₂ assimilation rate; gs: stomatal conductance to water vapor; SD: stomatal density (mm⁻¹) indices: d: traits estimated for the drought period only; LA: traits standardized per unit plant leaf area; *in situ*: gas exchange measured under greenhouse conditions; SS: gas exchange performed under steady-state conditions; *p* values: + < 0.1.
* < 0.05. ** < 0.01. *** < 0.001; Tmt: watering treatment (control or water deficit); Phe: high- or low-WUE phenotype.