# Oak stands along an elevation gradient have different molecular strategies for regulating bud phenology

Gregoire Le Provost, Céline Lalanne, Isabelle Lesur, Jean-Marc Louvet, Sylvain Delzon, Antoine Kremer, Karine Labadie, Jean-Marc Aury, Corinne da Silva, Thomas Moritz, et al.

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

## To cite this version:

Gregoire Le Provost, Céline Lalanne, Isabelle Lesur, Jean-Marc Louvet, Sylvain Delzon, et al. Oak stands along an elevation gradient have different molecular strategies for regulating bud phenology. BMC Plant Biology, 2023, 23 (108), 10.1186/s12870-023-04069-2. hal-04003618

HAL Id: hal-04003618

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

Submitted on 24 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH

Open Access

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

Check for
updates

# Oak stands along an elevation gradient have different molecular strategies for regulating bud phenology

Gregoire Le Provost<fnref n="1" />*, Céline Lalanne<fnref n="1" />, Isabelle Lesur<fnref n="1" />,<fnref n="2" />, Jean-Marc Louvet<fnref n="1" />, Sylvain Delzon<fnref n="1" />, Antoine Kremer<fnref n="1" />, Karine Labadie<fnref n="3" />, Jean-Marc Aury<fnref n="4" />, Corinne Da Silva<fnref n="4" />, Thomas Moritz<fnref n="5" /> and Christophe Plomion<fnref n="1" />

**Abstract**

**Background** Global warming raises serious concerns about the persistence of species and populations locally adapted to their environment, simply because of the shift it produces in their adaptive landscape. For instance, the phenological cycle of tree species may be strongly affected by higher winter temperatures and late frost in spring. Given the variety of ecosystem services they provide, the question of forest tree adaptation has received increasing attention in the scientific community and catalyzed research efforts in ecology, evolutionary biology and functional genomics to study their adaptive capacity to respond to such perturbations.

**Results** In the present study, we used an elevation gradient in the Pyrenees Mountains to explore the gene expression network underlying dormancy regulation in natural populations of sessile oak stands sampled along an elevation cline and potentially adapted to different climatic conditions mainly driven by temperature. By performing analyses of gene expression in terminal buds we identified genes displaying significant dormancy, elevation or dormancy-by-elevation interaction effects. Our Results highlighted that low- and high-altitude populations have evolved different molecular strategies for minimizing late frost damage and maximizing the growth period, thereby increasing potentially their respective fitness in these contrasting environmental conditions. More particularly, population from high elevation overexpressed genes involved in the inhibition of cell elongation and delaying flowering time while genes involved in cell division and flowering, enabling buds to flush earlier were identified in population from low elevation.

**Conclusion** Our study made it possible to identify key dormancy-by-elevation responsive genes revealing that the stands analyzed in this study have evolved distinct molecular strategies to adapt their bud phenology in response to temperature.

**Keywords** Response to temperature, Bud phenology, Elevation cline, Gene expression, Hormone quantification, Sessile oak

<footnote>*Correspondence:
Gregoire Le Provost
gregoire.le-provost@inrae.fr</footnote>

<footnote n="1">¹ INRAE, Univ. Bordeaux, BIOGECO, F-33610 Cestas, France</footnote>

<footnote n="2">² Helix Venture, F-33700 Mérignac, France</footnote>

<footnote n="3">³ Genoscope, Institut François Jacob, CEA, Université Paris-Saclay, Evry, France</footnote>

<footnote n="4">⁴ Génomique Métabolique, Genoscope, Institut François Jacob, CEA, CNRS, Univ Evry, Université Paris-Saclay, 91057 Evry, France</footnote>

<footnote n="5">⁵ Department of Forest Genetics and Plant Physiology, Umeå Plant Science Centre, Swedish University of Agricultural Sciences, 901 87 Umeå, Sweden</footnote>
## Background

Despite extensive gene flow between populations [1], tree species display clear adaptation to local environmental conditions [2–4]. Ample evidence of the efficiency of natural selection can be found in the large body of literature on forest tree provenance tests (reviewed by [5, 6]). In almost all tree species for which provenance tests have been performed, significant genetic variation between populations (i.e. differentiation) has been observed for fitness-related traits, including bud phenology [7]. These common garden experiments have also underpinned clinal patterns of phenotypic variation along temperature and/or photoperiod gradients. European white oaks illustrate how diversifying selection generated clinal variation of bud phenology variation and the rapidity with which forest tree populations became locally adapted over a short time span following the end of the last ice age [8, 9]. Indeed, extant sessile (Quercus petraea (Matt.) Liebl.) and pedunculate (Quercus robur L.) oak populations originating from the same or different (refugal) sources of glacial origin, but currently growing at different latitudes display strong phenotypic differentiation for the timing of leaf unfolding [10], but no differentiation for neutral genetic markers [11, 12]. The same holds true for populations distributed along elevation clines [10, 13]. Divergent selection along temperature gradients drives adaptation for early-late flushing dates along low–high elevation (or latitude) to avoid late frost damages [13, 14].

Environmental gradients (with continuous variations of environmental conditions) impose spatially variable selective pressures on populations within the range of the species, inducing genetic differentiation. The distribution of genetic diversity across environmental gradients provides information about the effects of past selection on the genetic composition of populations and facilitates the identification of parts of the genome that have undergone selection [15, 16]. Although elevation gradients can present significant climatic variations (i.e. precipitation, cloud cover, irradiance, light quality), such gradients offer ideal opportunities for the discovery of genes of importance for adaptation to temperature: (i) phenotypic divergence between populations is generated by a predominant environmental driver — temperature — avoiding confounding effects with other potential drivers of adaptation, (ii) extensive gene flow between populations leaves the genome with narrow genomic windows in which spatially varying selection is maintained and detectable, and (iii) the populations studied share the same evolutionary history, preventing confounding contributions between demographic dynamics and natural selection. Elevation gradients have, therefore, been

widely used to identify polymorphisms subject to natural selection [17], in diverse species, including forest trees [10, 18].

In the wild, Gene expression data have also proved powerful for the characterization of functional genetic variation and phenotypic plasticity between populations [19]. Comparisons of differences in expression profiles between locally adapted populations have been performed in many biological systems and have led to the discovery of genes of adaptive significance (reviewed by [20]). Moreover, the combination of allelic variation and gene expression data has proved very powerful for identifying signatures of adaptation in natural populations [21, 22]. However, it should be kept in mind that gene expression analysis in locally adapted populations to identify genes that matter for adaptation may have some limitations. Indeed, there is a shortcut between their differential expression and a clear sign of local adaption due to the possibility of epigenetic effects or possible population acclimation on gene expression. In Norway Spruce, [23] reported strong epigenetic effects on the expression of genes involved in bud phenology, cold acclimation and embryogenesis. Thus, genes identified by RNaseq must be validated by other approaches such as those described by [16].

As indicated above, bud phenology is a key trait contributing to local adaptation of forest trees. In temperate regions, the phenological cycle of the primary meristem can be split into two principal phases: a growing period (i.e. from bud burst to bud set), during which environmental conditions are favorable, and a resting period, also known as dormancy, from bud set to bud burst, when environmental conditions are unfavorable for stem elongation [24]. Dormancy begins when growth stops in late summer, in response to decreases in photoperiod and temperature [25], a phase known as paradormancy or summer dormancy [26]. This initial phase is followed by endodormancy, which is initiated by further decreases in temperature and photoperiod in the fall [27]. Endodormancy is generally associated with a cold acclimation process [28] and its induction is progressive, culminating in a complete lack of bud growth response under favorable conditions [29]. Endodormancy is the deepest stage of dormancy and is not released until chilling requirements have been fulfilled [30], generally around late January for oak populations in the Pyrenees [31]. In temperate regions, endodormancy generally starts in the early fall and peaks in early winter, with a timing dependent on the species considered and climatic conditions [32]. Once chilling requirements have been fulfilled, there is a transition from endodormancy to ecodormancy in the buds [24]. The development of ecodormant buds can be triggered by forced increases in air temperature in
late winter and early spring. The timing of the phenologi-
cal events observed in trees is correlated principally with
photoperiod (the main factor triggering entry into endo-
dormancy) and temperature (with chilling required to
break endodormancy and increased temperature trigger-
ing ecodormancy release) [33]. These two environmental
cues proved to be the most accurate in predictive models
of leaf unfolding [34].

In the ongoing context of global warming (which is
particularly pronounced in mountainous areas [35, 36]),
serious concerns have been raised about the long term
persistence of locally adapted populations at their cur-
rent locations. Indeed, with increasing temperatures,
the phenological cycle of tree species may be strongly
affected, because higher winter temperatures may pre-
vent endodormancy release if chilling requirements are
not met during winter [31], and late frosts in spring may
result in damage to trees if their buds flush too early [37].
Many molecular studies have been performed in trees,
to decipher the molecular mechanisms involved in dor-
mancy induction and release. These studies have revealed
that dormancy is regulated principally via molecular
pathways relating to phytohormones [38, 39], carbohy-
drates [40], temperature [41], photoperiod [42], oxygen
species [43], water [44] and cold acclimation. Most of
the studies aiming at deciphering the molecular players
underlying the dormancy phases have focused on poplar
[29], oak [45], spruce [46], chestnut [47] and beech [42].
However, these studies involved the use of a small num-
ber of genotypes harvested only during endodormancy
or ecodormancy, and therefore provide a limited view of
the internal and external drivers of dormancy induction
and release. The most detailed study performed to date
was published in poplar [48]. These authors identified key
molecular pathways involved in bud formation and endo-
dormancy induction and release. They also identified a
large set of genes commonly expressed during growth-to-
dormancy transitions, in poplar apical buds, cambium,
or *Arabidopsis thaliana* seeds, suggesting that similar
molecular mechanisms may underlie these processes in
different plant organs. This study was however conducted
under controlled temperature and photoperiod condi-
tions. Moreover, none of the studies cited above used
population genetics approaches to account for natural
variation in the ability of the trees to adapt their pheno-
logical cycles to fluctuating external cues.

Here, we used sessile oak populations sampled along
an elevation gradient in the Pyrenees to explore the
gene expression network involved in dormancy regula-
tion in this species. We addressed three main questions:
Does gene expression in vegetative buds reflect divergent
selection forces shaping population adaptation to eleva-
tion? Have stands sampled along the elevation gradient

evolved different molecular mechanisms for regulating
the endodormancy and ecodormancy stages? What gene
functions (i.e. molecular mechanism) have served as tar-
gets for adaptation at the level of gene expression? We
harvested terminal buds from trees growing at different
elevations, during endo- and ecodormancy. The avail-
ability of long term records of flushing dates along the
gradient allowed to account for differences in bud burst
date between populations. We quantified gene expres-
sion by RNAseq and used a likelihood ratio test to iden-
tify genes displaying significant dormancy, elevation or
dormancy-by-elevation interaction effects. The analysis
of dormancy-by-elevation interaction effect is of main
interest because this effect highlight different molecular
strategies in the oak stands sampled along the gradient
for adapting their bud phenology in response to environ-
mental variation (i.e. temperature, in this study).

Results

Differences in phytohormone concentrations
between elevations and dormancy stages

The endogenous concentrations of three main phytohor-
mones (indole-3-acetic acid (IAA), abscisic acid (ABA)
and cytokinin (CTK)) were determined on three biologi-
cal replicates for each dormancy stage, population and
valley. We first performed an ANOVA on the whole dataset,
to determine whether the two valleys could be con-
sidered as biological replicates. We used the following
fixed-effects model: Y<sub>ijk</sub> = μ + D<sub>i</sub> + E<sub>j</sub> + V<sub>k</sub> + ε<sub>ikj</sub>, where
D<sub>i</sub> is dormancy stage (i = "endodormancy" or "ecodor-
mancy"), E<sub>j</sub> is elevation (j = "low", "medium" and "high"
elevation) and V is the valley factor (k = "Ossau" or "Luz").
The results are shown in Additional File 1 Fig. S1 panel A.
No significant valley effect was detected, as for the date
of leaf unfolding (data not shown). In the second analysis,
we therefore used the whole dataset, with the following
fixed-effects model Y<sub>ijk</sub> = μ + D<sub>i</sub> + E<sub>j</sub> + D<sub>i</sub>*E<sub>j</sub> + ε<sub>ij</sub>, where
D<sub>i</sub> is the dormancy stage, and E<sub>j</sub> the elevation effect, to
identify significant dormancy, elevation and dormancy-
by-elevation interaction effects. The results are summa-
rized in Additional File 1 Fig. S1, Panel B for the three
phytohormones analyzed.

IAA were higher in EcoD samples (i.e. ecodormant
buds) at all elevations. CTK content were higher in EcoD
samples only at low and high elevations, respectively.
Finally, ABA concentrations were lower in EcoD samples
at all elevations (Additional File 1 Fig. S1, Panel B). These
results are consistent with those obtained by [49] and
[50] for ecodormant poplar and endodormant pear buds,
respectively, where similar hormone dynamics results
were obtained over the winter period.

A dormancy-by-elevation interaction effect was iden-
tified for both IAA and ABA (Additional File 1 Fig. S1,
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Evaluation of RNA-seq data quality using PCA (A) and clustering analysis (B) (using the “complete” method of the expander software) among the replicates of the different dormancy stages and elevations. For PCA blue and yellow dots represent EndoD and EcoD samples, respectively</figcaption></figure>

Panel B). IAA concentration increased with elevation in EndoD samples (i.e. endodormant buds), but the opposite pattern was observed in EcoD samples. ABA concentration increased with elevation in EndoD samples, but remained stable across the cline in EcoD samples.

RNA-seq analysis and identification of DEGs

The main goal of our study was to identify genes involved in dormancy regulation from the analysis of oak stands sampled at different elevations and potentially locally adapted to temperature (see methods section). We generated 24 cDNA libraries for this purpose. A general overview of the libraries generated is shown in Additional File 2 Table S1. More than 597 million reads in total were mapped onto the *Q. robur* reference genome. We made the choice to map our RNAseq reads on a closely related species (pedunculate oak) because currently, no reference genome is available for sessile oak. Mapping rates ranged from 66 to 82%. In total, 17,676 genes out of the 25,808 gene models available in the *Q. robur* reference genome were retained after quality-based filtering.

We used PCA (Principal Component Analysis) to explore transcriptome-wide changes in our dataset. This exploratory analysis showed close clustering within biological replicates, supporting the quality of the RNAseq data (Fig. 1A). It also showed clustering of the samples into two distinct groups corresponding to each dormancy stage (Blue and Yellow dots spread along the first PCA axis in Fig. 1A). Finally, it should be noticed that elevations separated also along the PCA 2 axis within a dormancy stage (except for two replicates sampled at 800 m of elevation) suggesting a stronger influence of temperature rather than photoperiod on dormancy regulation in our experimental design.

The adjusted data were also evaluated by a clustering approach, with the complete method available in Expander software (Fig. 1B). This analysis confirmed the high reproducibility of the RNAseq data and revealed that, for each treatment analyzed in this study, the biological replicates and dormancy stages were grouped into the same cluster (Fig. 1B). Both the PCA (second PCA axis) and clustering analyses exhibited interactions between dormancy and elevation.
A preliminary analysis was performed to determine whether the two valleys could be considered as biological replicates for RNAseq. For each dormancy stage, we compared gene expression levels between pairs of population samples from the same elevation (i.e. O-01 vs. L-01, O-08 vs. L-08, O-16 vs. L-16), resulting in six different gene sets (2 dormancy stages * 3 elevations). The results are summarized in Additional File 2 Table S2. Regardless of the dormancy stage considered, we identified very few genes differentially regulated between the two valleys, the maximum proportion of DEGs reaching 1.2% of the 17,676 genes between O-08 and L-08 in EndoD samples. These results confirmed that the two valleys could be considered as biological replicates. Increasing the number of biological replicates (from 2 to 4) decreased the false-positive rate, facilitating the distinction between noise and true biological signals in our dataset.

We then used likelihood ratio tests to identify genes displaying significant dormancy, elevation and interaction effects. Genes displaying differential expression (at least a two-fold difference in expression) with a corrected *P*-value < 0.01, and with adjustment for a false discovery rate (q-value) < 0.05 were considered to display significant differential expression. The results of the differential expression analysis are shown in Fig. 2, panel A. The overlap between the three effects is shown in Fig. 2, panel B. It should be noted that a very small number of genes (25) displayed all three effects simultaneously. Overall 2,084, 1,089 and 635 genes displayed significant dormancy, elevation and dormancy-by-elevation effects, respectively. We found that 397 genes were regulated by both dormancy and elevation, and that 203 genes displayed significant dormancy and dormancy-by elevation interaction effects (Fig. 2, panel B). Finally, 80 genes presented both significant elevation and dormancy-by-elevation interaction effects.

### Genes regulated during dormancy induction and release (Geneset#1)

By comparing models M1 and M2, using the threshold presented above, we identified 2,709 genes differentially regulated between EndoD and EcoD samples (Fig. 2, panel A, B and C). For this first gene set, most of the variation identified was quantitative (i.e. all the genes displaying differential expression were covered by reads from both the EndoD and EcoD libraries). We identified a single gene (Qrob_P0090160.2, encoding a riboflavin synthase) covered only by reads from EcoD samples. The results are presented in detail in Supplementary File 1. Gene set enrichment analysis was performed with the TopGo Package for this first gene set. Our analysis revealed different GO terms between EndoD and EcoD regulated genes. The results are shown only for the first

100 GO (ranked according to their *p*-value) in Supplementary File 1. A graphical representation (using bubble plots) of the first 20 ontologies is also available for a clearer view of the ontologies regulated.

We also performed subnetwork enrichment analysis independently for the genes regulated during EndoD and EcoD. We identified biological processes relating to cold tolerance, response to dehydration, freezing tolerance and drought tolerance in EndoD samples, whereas very different biological processes relating to meristem identity, flower development and floral organ identity were identified in EcoD samples. These results are consistent with those for dormancy induction and release reported in a previous study in oaks [45], validating the use of our experimental design and statistical model for finding dormancy-regulating genes and opening up promising new avenues for the identification of genes with expression profiles affected by elevation and interactions (see below).

### Elevation-regulated genes (Geneset#2)

By comparing models M1 and M3, we identified 1,941 genes regulated by elevation (Fig. 2, panel A, B and C and Supplementary File 2). We hypothesized that most of these genes may be involved in response to temperature. We clustered these genes according to their expression profiles. With the number of clusters set to 5, we obtained an overall homogeneity close to 0.98. Cluster#1 included 321 genes progressively upregulated by elevation; conversely cluster#2 contained 476 genes following the opposite pattern, whereas cluster#3 comprised 95 genes upregulated at intermediate elevations. Clusters #4 and #5 were strongly differentiated, with contrasting patterns of gene expression. Cluster#4 contained 715 genes that were overexpressed at the highest elevation, whereas cluster#5 contained 334 genes downregulated at highest elevation (Fig. 2, panel D and supplementary File 2).

We then performed gene set enrichment analysis on each cluster. GO terms enriched for each cluster in the three main ontologies are available in Supplementary File 2 both in tabular and graphical (bubble plot) representations. Subnetwork enrichment analysis was performed for each of the five clusters (Supplementary File 2). We identified hubs relating to leaf size, root growth and pollen development for the genes of cluster #1, senescence, cell death and cell division for cluster #2, root development, root growth and flower development for cluster#3, internode patterning, fatty acid omega oxidation and pedicel development for cluster #4, and heat tolerance for cluster #5.
## Genes displaying a significant dormancy-by-elevation interaction effect (Genest#3)

Using the likelihood ratio test to compare models M4 and M1, we identified 943 genes displaying a significant dormancy-by-elevation interaction effect (Fig. 2 panel A, B and D and Supplementary File 3). The majority of these genes displayed this effect only (635 of 943 genes, Fig. 2 panel B). Only 25 genes displayed all three main effects simultaneously. The remaining genes (i.e. 283 genes) displayed one of two significant effects: significant effects of dormancy for 203 genes and of elevation for 80 genes. We also applied the Kmeans method to this gene set to cluster genes according to their expression profiles. A stable homogeneity value of approximately 0.91 was obtained for a *k* value of five, giving five different clusters (Fig. 2 panel D and Supplementary File 3). Cluster #3 (139 genes) and cluster #2 (186 genes) had similar expression profiles characterized by a singular pattern in which certain genes were upregulated in EcoD samples at high elevation and in EndoD samples at low elevation. Cluster #1 (438 genes), cluster #5 (57 genes) and cluster #4 (118 genes) encompassed genes upregulated in EcoD samples at low elevation, but with stable or downregulated expression at the same developmental stage at high elevation.

GO terms enrichment analysis for each cluster is available in supplementary File 3 (both in tabular and graphical representations).

The subnetwork enrichment analysis is also shown for each cluster in Supplementary File 3. For cluster #2 and cluster #3, which were characterized by genes upregulated in EcoD samples in populations growing at high elevation, we identified significant enrichment for biological processes relating to the stress response: genes relating to “salinity”, “plant defense” and “plant development” in cluster #2, and genes involved in “salinity response”, “detoxification” and “cell elongation” in cluster #3. For clusters #1, #4 and #5, characterized by genes upregulated in EcoD samples in populations growing at low elevation, most of the biological processes identified were related to meristem functioning. Important hubs relating to “plant growth”, “mRNA splicing” and “transcription activation” were identified in cluster #1. For the other two clusters (#4 and #5), we

identified hubs relating to “plant growth”, “flower development” and “seed germination” for cluster #4, and “cell elongation”, “developmental process” and “flower development” for cluster #5.

## RT-qPCR validation

In total, 12 genes were selected for validation of their expression profiles by real time reverse transcription-quantitative PCR. Two of these genes displayed either a multi-banding pattern on agarose gel electrophoresis or unsuccessful PCR amplification and were removed from the analysis. For the other 10 genes, PCR efficiency ranged from 95 to 110%, in accordance with Taq polymerase activity (Additional File 2 Table S3). The patterns of expression of the genes tested by RT-qPCR were similar to those obtained by RNAseq (Additional File 1 Fig. S2), validating our RNA-seq data.

## Discussion

The main goal of this study was to provide new insights into the key molecular pathways involved in dormancy regulation. We used sessile oak populations sampled along an elevation gradient in the Pyrenees to explore the gene expression network involved in dormancy regulation in this species. We considered the two main phases of dormancy — endodormancy and ecodormancy — and performed an analysis of differential gene expression on terminal vegetative buds by an RNAseq approach. We identified genes displaying significant changes in expression associated with dormancy, elevation and dormancy-by-elevation interaction. Our results highlighted a number of molecular processes at work in sessile oak populations along an elevationcline. Our findings provide also a set of candidate genes involved in the response of sessile oaks to temperature at different elevations. Investigating genomic signatures of local adaptation in and around these loci could help identify the genetic bases of adaptation along this elevation gradient. Further, inform predictive models aiming to infer the risk of maladaptation in assisted migration strategies in the framework of global warming [16].

The variation of dormancy-related traits has been well investigated along latitudinal clines, in which both photoperiod and temperature vary [51, 52], but elevation clines over a restricted area have major advantages for

(See figure on next page.)

**Fig. 2** Illustration of the main results obtained from the DEGs analysis. Panel (A): Summary of the differential expression analysis. Results are shown for an adjusted *P*-value < 0.01 and a fold-change ratio > 2. The percentage indicated in the last column is given considering the 25,808 gene models found in the pedunculate oak genome. Panel (B): Venn-diagram showing the number of genes displaying main (Dormancy, Elevation) and dormancy-by-elevation-interaction effects or a combination of them. The Bar Chart illustrate the number of genes identified per effect and their associated clusters. Panel (C) Cluster identified for the elevation responsive genes and Panel (D) clusters identified for the dormancy-by-elevation responsive genes. For Panel (C) and (D) the X-axis represents the elevation (low: 100 m of elevation, mean: 800 m of elevation and high: 1,600 m of elevation. The Y-axis represents the normalized expression value for the genes belonging to the cluster (Mean and variance fixed to 0 and 1, respectively)
<table><thead><tr><th></th><th>Effect</th><th>Models comparison</th><th>Number of genes</th><th>Percentage</th></tr></thead><tbody><tr><td rowspan="3">(A)</td><td>Dormancy</td><td>M1 vs. M2</td><td>2,709</td><td>10%</td></tr><tr><td>Elevation</td><td>M1 vs. M3</td><td>1,941</td><td>7.5%</td></tr><tr><td>interaction</td><td>M4 vs. M1</td><td>943</td><td>3.5%</td></tr></tbody></table>

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

<figure><img src="image_6.png" /><figcaption>**Fig. 2** (See legend on previous page.)</figcaption></figure>
disentangling these two environmental factors and pro-
vide a unique opportunity to analyze the effect of tem-
perature independently from that of photoperiod [9].
Moreover, to our best knowledge, elevation clines have
never been used to analyze the molecular mechanisms
involved in bud dormancy regulation in a forest tree spe-
cies. The predicted temperature increase in the coming
years can be expected to affect the phenological cycle of
perennial species [28, 37]. We therefore mainly focused
the discussion on genes regulated by elevation (i.e. tem-
perature-responsive genes) and genes displaying a signifi-
cant dormancy-by-elevation interaction effect, because
these two sets of genes may include genes of importance
for forest tree adaptation. However, we have to bear in
mind that the molecular mechanisms identified in this
study were obtained using a single sampling campaign.
Dormancy induction and release is strongly influenced
by temperature and photoperiod. Thus, interannual cli-
mate variability can influence gene expression, leading
potentially to different gene networks identified if multi-
year sampling campaigns had been used.

Molecular plasticity during development underlying
the shift between endo- and ecodormancy

We first quantified three phytohormones (IAA, ABA and cytokinines) to characterize the samples phenotypically. As reported in fruit trees [53], ABA content was higher in EndoD samples. Several authors have also reported a role for ABA in the photoperiodic control of growth and, potentially, in the cessation of growth induced by short days [29]. Conversely, IAA and cytokinin concentrations were higher in EcoD samples. IAA concentration is known to increase during the reactivation of cambium activity [50] and cytokinins are activated during budburst [54]. These findings suggest that the higher concentrations of these two phytohormones in the buds during ecodormancy are correlated with the reinitiation of mitotic activity in the buds of the sampled trees.

Transcriptomic studies identified two biological net-
works corresponding to gene set#1 related to either
EndoD or EcoD samples (Supplementary File 1). As
expected, important hubs related to cold tolerance, freez-
ing tolerance, drought tolerance and defense response
were identified in the EndoD network, whereas biologi-
cal hubs related to meristem functioning (meristem iden-
tity, flower development, floral organ identity, etc.) were
found in the EcoD network. Similar results have already
been published for oak [45] and beech [42], but the
dataset reported here is much broader than those previ-
ously published, extending to genes with lower levels of
expression (e.g. genes encoding transcription factors).
The upregulation of molecular mechanisms involved in
dehydration and cold tolerance during endodormancy

is consistent with the greater tolerance to ice formation
of cold-tolerant dehydrated tissues, making it easier for
the plant to cope with cold winter temperatures [55].
The central hubs (i.e. biological processes) relating to cell
activity and meristem functioning identified in the eco-
dormancy network reflect the reinitiation of the cellular
metabolism required for bud burst when environmental
conditions become favorable [56, 57].

Environmental molecular plasticity in response
to temperature

In mainland France, regardless of the climatic scenario modeled, mean temperature is predicted to increase by 0.6 °C (in northern France) to 1.3 °C (in southern France) by 2050, with even greater temperature contrasts predicted for mountainous areas [58]. This increase in temperature may influence the distribution of forest ecosystems [59]. Thus, the elevationcline used here, with a temperature gradient extending over 6.9 °C, provided a opportunity to analyze the role of gene expression in responseto temperature.

(a) Populations from low and medium elevations

Two clusters encompassed genes upregulated at low
and/or medium elevations (Cluster #2 and #5, respec-
tively, Supplementary File 2, Fig. 2 panel C). Populations
of these two clusters exhibit an extended growing sea-
son characterized by a very early budburst date (Fig. 3).
This phenology may be explained by the main molecular
mechanisms identified in our subnetwork enrichment
analysis.

First, we identified a large number of genes in Clus-
ter #2 involved in cell death inhibition, suggesting an
extension of cell life expectancy, which would enable
the tree to extend its growing period. Indeed, half the
genes relating to the cell death biological process are
involved in cell death inhibition (CAS1, BON3, GLIP1,
TOPII, T3G21.18 and F5O11.34). For example, we iden-
tified a CAS1 gene encoding a cycloartenol synthase
1 (4.2 fold-change ratio -FC- between low and high
elevation populations). Kim et al. (2010) [60] showed
that *Arabidopsis* transgenic lines in which this gene
was downregulated were characterized by a cell death
phenotype, suggesting a key role of the CAS1 gene in
cell death inhibition. Similarly, *BON3* (3.8 FC), a gene
encoding a copine-like protein, has been reported to
be involved in inhibiting cell death [61]. A TOPII gene
(2.24 FC) encoding a topoisomerase-like protein, was
found to be downregulated in rice, this downregula-
tion being associated with delayed cell death [62]. We
also identified 15 genes encoding proteins involved in
the biological process “senescence”, including five genes
(MYB62 (2.45 FC), AGL15 (4.35 FC), LOX2 (2.15 FC),
ABI5 (2.81 FC) and SAG13 (2.23 FC)) involved in the inhibition of senescence, in agreement with our observations indicating that leaf fall occur later in populations growing at lower elevations. For example, the MYB62 gene, which encodes a R2R3 MYB transcription factor, is known to be involved in apical dominance, delayed flowering time and late senescence in sunflower [63]; and the AGL15 gene, encoding an agamous-like protein, has been clearly implicated in delayed senescence in *Arabidopsis thaliana* [64].

Finally, nine genes involved in cell division and 23 genes involved in root growth were identified, suggesting that populations growing at lower elevations may have higher levels of meristematic activity, enabling the trees to flush earlier, thereby extending their growing period. One of these genes, *ERF3*, encodes a protein similar to an ethylene-responsive factor that stimulates cell division in the root tip in rice [65]. Another, *MYB12*, encodes a protein known to promote cell division and elongation in *Arabidopsis* [66].

Cluster #5 was enriched in a single biological process related to “heat tolerance” (Supplementary File 2). All the genes in this hub are involved in the acquisition of heat tolerance. This mechanism enables higher plants to cope with heat stress and to maintain their productivity under unfavorable environmental conditions (reviewed by [67]). We hypothesized that higher heat tolerance acquisition in populations growing at lower elevations would be an effective strategy for coping with stressful conditions

more frequently encountered in such populations, which have an extended growing period. Two of the genes identified were considered to be of particular importance: *AGO1* (FC=2.26) and *TOR* (FC=2.55). *AGO1* encodes an RNA slice protein known to be essential for the maintenance of acquired heat tolerance [68], and the down-regulation of *TOR* (encoding a phosphatidylinositol 3-kinase family-like protein, TOR, AT1G50030) in the shoot apical meristem of *Arabidopsis thaliana* is known to weaken heat tolerance [69].

### (b) Populations from high elevations

Two other clusters (#1 and #4) contained genes upregulated at higher elevations (Supplementary File 2, Fig. 2 panel C). Oak populations from high elevations flush later and their vegetative phase is much shorter, occurring when environmental conditions are favorable [9]. Moreover, leaf size is smaller in high-elevation than in low-elevation populations [70]. The molecular mechanisms we identified in the subnetwork enrichment analysis may explain these observations. Indeed, biological processes relating to “leaf size”, “pollen development”, “root growth”, “internode patterning” and “pedicel development” displayed enrichment in clusters #1 and #4.

Ten of the genes in cluster #1 are involved in the biological process “leaf size”. Three of these genes (POL, T1M15.220 and GAI) have been implicated in the reduction of leaf size, the others being essential for normal leaf development. POL (FC=2.47 between

<figure><img src="image_7.png" /><figcaption>**Fig. 3** Mean date and standard error of leaf unfolding (y-axis, Date of the Year) for the three selected sessile oak populations in the Luz (L) and Ossau (O) valley at 100 (O-01, L-01), 800 (O-08, L-08) and 1,600 m (O-16, L-16) a.s.l</figcaption></figure>
high- and low-elevation populations) encodes a protein phosphatase 2C. It is known to be involved in reducing leaf size in *Arabidopsis* [71]. T1M15.220 (encoding a protein resembling SAUR, FC=2.21) is also known to be associated with smaller leaf size in plants (reviewed by [72]). Finally, the overexpression *GAI* (encoding a GRAS family transcription factor, FC=2.63) has been shown to lead to smaller leaf size in petunia [73].

Cluster #1 also included 16 genes relating to root growth and seven genes relating to meristem size. Most of the genes identified in this cluster are essential for meristem functioning, again illustrating the ability of high-elevation populations to produce leaves over a short growing period. For example, the BIG gene (encoding an auxin transporter protein, FC=2.6) is known to be essential for cell proliferation in response to sucrose and glucose in *Arabidopsis* roots and shoot meristems [74]. RPK2 (encoding a protein resembling a receptor-like kinase RPK2, FC=2.36) is known to be involved in regulating meristem size by controlling cell proliferation [75].

Eleven genes relating to the biological process “pollen development” were also identified, suggesting that the cellular machinery of high-elevation populations is activated to ensure rapid flushing and flowering when environmental conditions become favorable. These 11 genes included three GSL genes (GSL1 (FC=2.32), GSL8 (FC=2.49) and GSL10 (FC=2.62)) encoding callose synthases. Callose synthase has been reported to play a key role in pollen development in *Arabidopsis thaliana* [76]. We also identified two RGL genes (RGL2 (FC=2.77) and RGL1 (FC=3.47)) encoding the DELLA protein. Fleet & Sun, (2005) [77] reported that the DELLA gene repressed petal and pollen development, together with anther elongation, in *Arabidopsis*, suggesting a possible role of these genes in the delay of flowering time observed in high elevation populations.

Finally, only one biological process related to “internode patterning” and corresponding to three genes (KNAT1 (FC=2.63), STM (FC=4) and BLH8 (FC=2.43)) was identified in cluster#4 (Supplementary File 2). Two of these genes (STM and BLH8) control meristem formation and/or maintenance, organ morphogenesis, organ position, and several aspects of the reproductive phase in *Arabidopsis thaliana*, suggesting a possible key role in meristem function [78].

To summarize, our study indicates that the main molecular mechanisms associated with genes differentially expressed at low and high elevations may account for the typical phenology and features of adaptation to temperature observed in these oak stands. Indeed, at low elevation, populations display a higher basal level of expression for genes involved in the inhibition of cell

death and senescence, the acquisition of heat tolerance and the maintenance of cell division activity, enabling the trees to flush earlier when environmental conditions are favorable (Fig. 4). Conversely, populations from high elevations are characterized by an overexpression of genes involved in leaf size reduction, meristem function and delayed flowering, enabling the trees to cope with the shorter growing season at high elevations (Fig. 4). It should be also noticed that our results are also in accordance with those of [79]. Indeed, these authors reported that at low elevation natural selection favors vegetative growth than reproduction thus explaining the early flushing date observed at low elevation.

## Dormancy-by-elevation interaction-responsive genes reveal different molecular strategies for coping with temperature variation in the oak stands

The five clusters identified with the Kmeans approach were merged into two main groups: a group of genes upregulated in ecodormant buds at high elevation (clusters #2 and #3, Fig. 2 panel D) and a group of genes upregulated in ecodormant buds at low elevation (clusters #1, #4 and #5, Fig. 2 panel D).

The high-elevation gene pool corresponds to late flushing populations [80]. These populations have to cope with severe environmental conditions. A 6.9 °C temperature gradient has been reported along the studied elevationcline. Once chilling recruitment has been fulfilled, the buds in the high-elevation populations enter the ecodormancy phase. However, they still have to cope with damage due to the late frosts known to occur during early spring at this elevation. The main molecular mechanisms identified in clusters #2 and #3 (i.e. flowering time, plant development, salinity response and detoxification (i.e. cold tolerance)) may explain the typical bud phenology observed in these high-elevation populations.

By contrast, populations at low elevations consist of early-flushing trees with a longer growing period [80]. We thus assume that the main biological processes found in clusters #1, #4 and #5 (cell activity, meristem functioning and flower development) prepare the buds to flush as soon as the heat requirement is fulfilled, accounting for the longer growing period observed in these populations.

An analysis of the genes displaying an elevation x dormancy phase interaction effect therefore provides elements of the answer to the following additional questions: how do populations at high elevation cope with late frost damage during ecodormancy? Which molecular mechanisms account for the early flushing date observed in populations at low elevation?
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Overview of the results obtained in this study. Main molecular functions associated to elevation and dormancy-by-elevation interaction effect are indicated in color. Red color was used for population from low elevation, Blue color for population from high elevation. The number of genes displaying a significant (i) dormancy, (ii) elevation or (iii) dormancy-by-elevation interaction effect is indicated in grey</figcaption></figure>

(a) How do populations from high elevations cope with adverse environmental conditions during ecodormancy to avoid late frost damage (cluster#2 and #3)

The network enrichment analysis performed for clusters #2 and #3 (Supplementary File 3) identified several genes involved in the inhibition of cell elongation and/ or delayed flowering time, but also genes involved in the response to abiotic stresses. For the biological process "flowering time", we identified ABF4 (a b-zip transcription factor) known to be involved in delaying flowering time in *Arabidopsis* [81]. The WRKY71 transcription factor, also found in cluster#2, directly affects the flowering time of plants by regulating the Constance and Flowering locus genes [82]. The CRY1 gene (cryptochrome 1 apoprotein) encodes a protein involved in the photoperiodic regulation of flowering and the inhibition of cell elongation [83]. MFB16.6 (encoding the SPL13 protein) was also upregulated during ecodormancy in high-elevation populations. Gao et al. (2018) [84] reported a potential role of this gene in delaying flowering time in *Medicago*. We also identified many genes involved in the regulation of cell elongation, generally through its inhibition, suggesting a potential role for these genes in slowing meristem development. Two of these genes, SWN and At3g51890,

encode a swinger and the CLC3 protein, respectively. Footitt et al. (2015) [85] reported that peak SWN gene expression coincided with peak dormancy in *Arabidopsis*, suggesting a key role of this gene in dormancy release. Wang et al. (2013b) [86] identified the CLC3 gene as a key molecular player involved both in basipetal transport and the sensitivity and distribution of auxin. However, the principal role of CLC3 is in auxin transport, and it was difficult to relate its expression level to auxin content in our biological samples.

Moreover, cluster#3 contained genes potentially involved in the acceleration of cell elongation or in the activation of flowering, potentially accounting for the faster growth and more rapid reproduction of these populations. BZR1 (encoding a positive regulator of brassinosteroid signaling) is known to induce cell elongation at high temperatures, under the combined effect of auxin and brassinosteroids [87]. *ROT3*, encoding a cytochrome P-450, is known to be involved in determining the rate of leaf expansion in response to brassinosteroid [88]. We also identified a CIB1 gene encoding a protein similar to the cryptochrome-interacting-basic-helix-loop-helix protein. Liu et al. (2018) [89] showed that cryptochromes are blue light receptors that mediate light responses
in plants and animals and activate flowering once environmental conditions are favorable (temperature and photoperiod).

Finally, we identified several genes relating to “salinity response”, “plant defense” and ‘detoxication’. This is not surprising, because these processes may enable the bud to remain dehydrated during ecodormancy, and, therefore, to cope more effectively with the low temperatures frequently observed during early spring at high elevations. Several studies have also reported that detoxication processes are essential for the acquisition of cold tolerance in plants [90, 91]. For example, we identified an APX3 gene encoding an ascorbate peroxidase involved in ROS detoxication. Wang & Li (2006) [92] reported that APX genes are overexpressed during cold acclimation. We also found an OSSA1 gene encoding a protein similar to a cytosolic O-acetylserine(thiol)lyase known to be involved in ROS detoxication in plants [93]. Two GSTU genes (GSTU17 and GSTU9) encoding proteins resembling glutathione transferase were also identified. Fujino & Matsuda, (2009) [94] showed that the GSUT gene contributed to cold tolerance in rice. The PR4 gene encodes a pathogenesis-related protein, and Cabello et al. (2012) [95] reported that the induction of this gene led to cold tolerance in *Arabidopsis*.

**(b) Which molecular mechanisms underlie the early flushing observed in populations from low elevations?**

Our subnetwork analysis on clusters #1, #4 and #5 (Supplementary File 3) highlighted biological processes relating to “mitotic cell activity”, “flower” and “plant development”. These molecular mechanisms are potentially involved in the early bud break observed in low-elevation populations and may account for the longer growing period of these populations. Genes relating to “plant growth” were mostly found in cluster#1 (32 genes). Seven genes with a potential role in the early initiation of meristem activity were identified: (i) a GR-RBP2 and two emb1138 genes encoding a glycine- and lysine-rich RNA-binding protein and DEAD box RNA helicases, respectively. Members of these two gene families are known to be involved both in cold stress tolerance and accelerating seed germination and seedling growth at low temperature [96, 97], (ii) an FBL17 gene encoding a protein similar to an ubiquitin protein ligase FBL17 was also detected. Noir et al. (2015) [98] showed that a loss of FBL17 function strongly decreased cell division rates in the apical meristem. Similarly, an SHR gene encoding a protein resembling GRAS family transcription factors, which known to be key regulators of cell division and meristem activity, was also identified [99], (iii) an SE gene encoding a C2H2 zinc-finger protein, Serrate (SE), which is known to accelerate leaf

production and to decrease time-to-flowering, was also identified [100], and finally, (iv) we identified two genes (ARIA and F21P24.13) involved in ABA sensitivity, heat tolerance and seedling growth, suggesting a potential role of these genes in reinitiating meristem activity [101].

We also identified four and six genes relating to “flower development”. Notably, cluster#4 included a LOX3 gene and a GA20X2 gene involved in hormone signaling. The LOX3 gene encode a lipoxygenase implicated in the biosynthesis of jasmonic acid and flower development [102]. The GA2OX 2 protein resembles a gibberellin 2-beta-dioxygenase known to be involved in the gibberellin biosynthesis pathway. Yamauchi et al., (2007) [103] reported a role for this gene in partly suppressing germination during exposure to darkness after phytochrome inactivation. An F15K9.22 gene encoding a protein similar to the Fantastic Four (FAF) protein was also found in this cluster. This gene influences the activity of the shoot apical meristem and is upregulated once flowering is initiated [104].

Cluster#5 included three genes (RAP2.7, GAI and LP1) with potential roles in both the early budburst and early flowering observed in populations from low elevations: (i) RAP2.7, encoding a protein resembling Apetala2 protein. Okamuro et al. (1997) [105] showed that this gene is active in the meristem during both reproductive and vegetative development, (ii) a GAI gene encoding a protein resembling a transcription factor involved in gibberellin signaling. In *Arabidopsis thaliana*, the upregulation of this gene accelerates flowering [106], and (iii) an LP1 gene encoding a non-specific lipid transfer protein. Nieuwland et al. (2005) [107] showed, by RNAi in *Arabidopsis*, that the LP1 gene played an important role in meristem function. Its downregulation resulted in dwarfing and a disruption of flower development.

**Conclusion**

In summary, this study provides a comprehensive overview of the sessile oak bud transcriptome along an elevation gradient. Our monitoring of gene expression along altitude and during eco and endodormancy highlights + plastic responses to temperature, but also a potential genetic control of this plastic response. Indeed low- and high-altitude populations have evolved different molecular strategies for minimizing late frost damage and maximizing the growth period, thereby increasing their respective fitness in these contrasting environmental conditions. As summarized in Fig. 4, late budburst of populations at high elevation was found to be related to the overexpression of genes involved in the inhibition of cell elongation and delaying flowering time. Genes associated with cold tolerance were also found, enabling the buds to avoid late frost damage during early spring.
Conversely, in populations at low elevation, we found a higher expression of genes involved in accelerating cell division and flowering, enabling buds to flush earlier once environmental conditions become favorable.

A molecular understanding of how oak stands have fine-tuned their leaf unfolding process to a temperature gradient throughout the dormancy period, over micro-evolutionary time scales, provides a valuable basis for further investigation in light of the increasing temperatures due to global warming. The candidate genes that we identified could be validated in the future by studying their gene expression level in reciprocal transplantations. Such experimental designs are available [70] for the analyzed stands and will be helpful in the coming years to quantify the respective contributions of phenotypic plasticity and genetic variation.

## Methods

### Study area and populations

We studied natural populations of sessile oaks (*Quercus petraea* (Matt.) Liebl.) growing in two neighboring valleys (Ossau and Luz, denoted O and L, respectively) on the northern side of the Pyrenees in France (from 43°15′N, 00°44′W to 42°53′N, 00°06′E). In each valley, populations grew at elevations varying between 100 and 1,600 m a.s.l., corresponding to a temperature gradient of 6.9 °C between the populations at the lowest and highest elevations [108]. Functional traits related to phenology, morphology and physiology of these populations have been extensively monitored in situ and ex situ (i.e. in common gardens) over the last 15 years, to quantify the respective contributions of phenotypic plasticity and genetic variation to their within and between population variation [13, 31, 70, 109]. In addition, reciprocal transplant experiments conducted at a smaller scale confirmed the co-gradient variation [6] of the timing of bud burst along elevation, thus enhancing local adaptation [80]. A general overview of the populations sampled in this study is provided in Additional File 2 Table S4 and Fig. 5A.

### Bud phenology monitoring and sampling

The timing of leaf unfolding was recorded over a period of eight years (between 2002 and 2014) in six populations, consistently taking measurements from the same 20–34 adult trees per population each year, over the two elevation gradients. Leaf unfolding was monitored at 10-day intervals in each population, from March to June, by the same two observers with binoculars (magnifying power × 10, Fig. 5D), standing about 15 m away from the tree. At each visit, the developmental state of each tree was assessed and the date of leaf unfolding relative to January 1st (i.e. day of the year, DOY) was determined as

the date at which the leaves of 50% of the buds were fully unfolded, as described in a previous study [9].

We selected three populations per valley, based on the timing of leaf unfolding, as representative of the intraspecific phenological variability along these two elevation clines (Fig. 3). These populations were growing 100, 800 and 1600 m a.s.l. and differed considerably in terms of the timing of leaf unfolding (about 50 days' difference between the populations at the lowest and highest elevations).

Buds were sampled as follows: (i) a first sampling campaign was performed during the endodormancy period in the fall of 2013. Endodormancy is triggered by both light and temperature, but with photoperiod playing the dominant role [110], as demonstrated by short-day experiments [111]. We therefore harvested buds from all the trees during the same week. The samples collected were endodormant buds, referred to hereafter as EndoD buds, (ii) a second sampling campaign was performed in spring 2014, just before bud burst. As release from ecodormancy is driven primarily by increases in temperature (i.e. trees from populations growing at high elevations flush later than those from lower elevations), we used the leaf unfolding time series to determine the most relevant sampling dates at each elevation. Sampling thus spanned a period of four weeks, from early March to early April (Additional File 2 Table S4). The samples collected were ecodormant buds, referred to hereafter as EcoD buds.

Buds were sampled from the same 20 trees in each population for both developmental stages. A shotgun was used to detach the terminal branches for sampling. The terminal buds were then harvested and immediately frozen in liquid nitrogen, for storage at -80 °C until subsequent analyses (Fig. 5B and 5C). Taxonomic attribution of the trees used in this study was confirmed using diagnostic SNP markers described by [112].

A formal permission both from the National Park of the Pyrenees and to the different town halls was obtained in order to sample buds in the different populations used in this study.

### Quantification of hormones

For hormone quantification, we used three pools of six individuals each, corresponding to three biological replicates for each population and dormancy stage (i.e. 36 samples in total). Buds were independently ground to a fine powder in liquid nitrogen and freeze dried, and then mixed in equimolar ratios to constitute pools. Phytohormones were quantified as previously described [113].

### RNA extraction and sequencing

RNA was extracted independently for each individual and dormancy stage, according to the procedure described by
<figure><img src="image_9.png" /><figcaption>**Fig. 5** Sampling buds from sessile oak trees in the Pyrenees. (A): Location of the three selected populations along two elevation clines (Ossau and Luz valleys). Population ID follows that of Additional file 2 Table S4. (B): Shotgun to harvest terminal branches on three canopee. (C): Bud storage in liquid nitrogen. (D): Leaf unfolding observation with binoculars</figcaption></figure>

[114]. Residual genomic DNA was removed before purification with DNase RQ1 (Promega, Madison, WI, USA), according to the manufacturer's instructions. The quantity and quality of each extract were then determined with an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). We generated 24 libraries, corresponding to 2 dormancy phases × 3 populations × 2 valleys × 2 biological replicates. For each biological replicate, we pooled equimolar amounts of RNA extracted from 10 individual trees. Libraries were generated with

the Illumina protocol (TruSeq Stranded mRNA Sample Prep Kit, Illumina, San Diego, CA, USA). Briefly, we selected mRNA from a 2 µg sample of total RNA. mRNA was chemically fragmented and reversed transcribed with random hexamer primers. The second-strand cDNA was generated to create double-stranded cDNA. The cDNA was 3'-adenylated, and Illumina adapters were added. We amplified DNA fragments (with adapters) by PCR with Illumina adapter-specific primers. We quantified the libraries with a Qubit Fluorometer (Life
Technologies, NY, USA). We estimated library size with Agilent 2100 Bioanalyzer technology (Agilent). We then sequenced each library by 101 base-length read chemistry, in a paired-end flow cell, on an Illumina HiSeq2000 sequencer (Illumina).

## Trimming, mapping and identification of differentially expressed genes (DEG)

Cleaning and mapping were performed as described by [115]. Briefly, we first removed low-quality reads and retained sequences with a mean quality value above 20. We mapped reads onto the 25,808 oak gene models published with the reference oak genome [116], using BWA (V.0.6.1) with the default parameters. We then identified differentially expressed genes (DEG) with the DESeq2 package, using a P-value < 0.01 after adjustment for multiple testing with a false discovery rate (FDR) of 5%. We also considered two-fold changes in expression ratio as a threshold for identifying the genes with the highest degree of differential expression. We investigated the effects of dormancy stage, elevation and their interaction in likelihood ratio tests. The dormancy and elevation effects were assessed by comparing a statistical model without interaction (M1) with two reduced models for the dormancy (M2) and elevation (M3) effects, respectively.

$$M1: Y_{ijk} = \mu + D_i + E_j + \varepsilon_{ijk}$$

where Dᵢ is the dormancy stage (i = "Endodormancy" or "Ecodormancy"), and Eⱼ the elevation effect (j = "low-100 m", "medium-800 m" or "high-1,600 m" elevation).

$$M2: Y_{jk} = \mu + E_j + \varepsilon_{jk}$$

$$M3: Y_{jk} = \mu + D_i + \varepsilon_{jk}$$

For estimation of the interaction effect, we compared the following a complete model:

$$Y_{ijk} = \mu + D_i + E_j + (D * E)_{ij} + \varepsilon_{ijk} \quad (M4) \text{ to M1}$$

Three gene sets were thus generated: (i) DEGs (geneset#1) corresponding to dormancy regulation (regardless of elevation), (ii) DEGs (geneset#2) corresponding to differences in elevation (regardless of dormancy stage), and (iii) DEG (geneset#3) displaying a significant dormancy-by-elevation interaction. Annotations for each DEG were recovered from the published pedunculate oak genome sequence [116]. Below, we focus particularly on geneset#2 and #3, which should include the key molecular players involved in response to temperature (and potentially to local adaptation) to temperature for geneset#2, and should reveal differences in the strategies of oak stands analysed across bud

phenological stages to cope with temperature variation for geneset#3. The elevation term used here does not allow to disentangle the effect of phenotypic plasticity from those of genetic differentiation. Indeed, [117] reported that the adaptive response of populations to divergent selection pressures along the elevation systematically results from a combination of non-optimal phenotypic plasticity and genetic differentiation. Thus, other experimental design such as reciprocal transplantations are needed to separate these two effects in the stands analyzed.

DEGs from geneset#2 and #3 were analyzed with EXPANDER software [118], which clusters genes according to their expression profile, using a Kmeans algorithm [118]. For both gene sets, we set the number of clusters to 5 (k=5) to maximize the homogeneity of each cluster. The genes from each cluster were then used for independent gene set and subnetwork enrichment analysis (see below).

## Gene set and subnetwork enrichment analysis

GO term enrichment analyses were first performed with the TopGO R package [119, 120]. We considered corrected P-values below 0.05 (after adjustment for multiple testing with a false discovery rate (FDR of 5%)) to indicate significant enrichment in a particular GO term in our dataset. We then performed subnetwork enrichment analysis with Pathway Studio software (Pathway Studio®, Elsevier 2017), as described by [115].

## RT-qPCR validation

For the elevation and dormancy-by-elevation effects, we selected six genes each for RT-qPCR validation. All 12 genes analyzed are listed, together with their associated effect, in Additional File2 Table S3. RT-qPCR quantification was performed with 1 µg total RNA, as described by [121], on a Lightcycler® LC 384 (Roche Lifescience, Germany, EU) with standard parameters, as described by [115]. Fluorescence data were normalized against two control genes (Qrob_P0530610 and Qrob_P0426000) and analyzed by STATQPCR (http://satqpcr.sophia.inra.fr/cgi/home.cgi). All the primer pairs used here were designed with Primer3 software [122].

### Abbreviations

<table><tr><td>EndoD</td><td>Endodormant buds</td></tr><tr><td>EcoD</td><td>Ecodormant buds</td></tr><tr><td>L</td><td>Luz valley</td></tr><tr><td>O</td><td>Ossau valley</td></tr><tr><td>O (L)-01</td><td>Ossau (or Luz) valley 100 m a.s.l</td></tr><tr><td>O (L)-08</td><td>Ossau (or Luz) valley 800 m a.s.l</td></tr><tr><td>O (L)-16</td><td>Ossau (or Luz) valley 1,600 m a.s.l</td></tr></table>
RT-qPCR Quantitative Real Time Polymerase Chain Reaction

GO Gene Ontology and

PCA Principal Component Analysis

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12870-023-04069-2.

**Additional file 1: Figure S1.** Evolution phytohormone content in each population according to the Dormancy stage. Panel A: Table for ANOVA results for phytohormone analysis. *P* value is indicated in each cell. *P* value<0.05, **P* value<0.01 and ***P* value<0.0001. NS stands for not significant. Panel B: Graphical representation of their accumulation over the Dormancy period. We used blue and orange color for EndoD (i.e. Endormancy) and EcoD (i.e. Ecodormancy) samples, respectively. Standard deviations were obtained from the 3 measurements performed in each population (Low=O-01+L-01, Mean=O-08+L-08 and High=O-18+L-18). Effects identified in the linear model where also indicated. Abbreviations correspond to: D: Dormancy effect, E: Elevation effect, D\*E: interaction effect. *P* value<0.05, **P* value<0.001 and ***P* value<0.0001. **Figure S2.** qPCR validation. Panel (A) genes displaying a significant elevation effect. Panel (B) genes displaying a significant Dormancy-byelevation effect. Abbreviations correspond to EndoD: Endodormant buds, EcoD: Ecodormant buds, Low: 100 mts (i.e. O-01+L-01), Mean: 800 mts (i.e. O-08+L-08) and High: 1,600 mts (i.e. O-16+L-16). Standard deviations were obtained from the four biological replicates.

**Additional file 2: Table S1.** Overview of the cDNA libraries generated in this study. **Table S2.** Gene expression level comparison between the two valleys for a specific dormancy stage. The comparison was performed for the population harvested at a same elevation in the two valleys considered. In each cell, we indicated the number of differentially expressed genes; the percentage is indicated in parenthesis. **Table S3.** List of the primer pairs used for qPCR analysis. Abbreviations: Tm: annealing temperature, E: Elevation effect, D\*E: Dormancy-byelevation effect. For each effect, the cluster ID is indicated in parenthesis in the second column. Qrob_IDs were retrieved from the oak genome available in Plomion et al. (2018). **Table S4.** Overview of the sessile oak populations used in this study.

**Additional file 3: Supplementary Files 1, 2, 3.** are available online at the INRAE dataverse portal: G. Le Provost, 2021. Oak stands along an elevation gradient have different molecular strategies for regulating bud phenology. https://doi.org/10.15454/XMEKFX - Portail Data INRAE, v3.0. These files includes normalized values for RNAseq data, Fold change Ratio, Gene set enrichment analysis and subnetwork enrichment analysis for genes displaying a significant dormancy, elevation and dormancy-by-elevation interaction effect respectively.

## Acknowledgements

We thank the Genotoul Bioinformatics Platform Toulouse Occitanie (Bioinfo Genotoul, https://doi.org/10.15454/1.5572369328961167E12) for providing computing resources.

## Authors' contributions

GLP and CP designed the study. GLP, IL, AK and CP wrote the manuscript. GLP, CP and JML sampled the populations. GLP performed bioinformatics analysis with IL. GLP and CL were involved the RT-qPCR experiments. cDNA libraries were constructed and sequenced by KL, CDS and JMA. SD, AK and JML were involved in monitoring bud burst date and CP selected the populations for this study. TM performed hormone quantification. All the authors have read and approved the manuscript.

## Funding

This work was supported by France Génomique (ANR-10-INBS-09-08, Oakadapt project). IL received funding from the European Research Council under the European Union's Seventh Framework Program (FP/2007-2013) / ERC Treepeace Grant Agreement n. 339728. Illumina sequencing was performed at the Genoscope sequencing facility (www.genoscope.cns.fr). The current publication benefitted also fundings from the ANR (EPITREE ANR-17-CE32-0009-01).

## Availability of data and materials

All the sequences generated for this publication have been deposited in the short-read archive of NCBI (https://www.ncbi.nlm.nih.gov/sra) under accession PRJEB17876. All data generated or analysed during this study are included in this published article and its supplementary online information files (Additional File 1 and 2 and Supplementary File 1,2 and 3).

## Declarations

### Ethics approval and consent to participate

Our research complies local, national or international guidelines and legislation on plant. According to Nagoya protocol on access and benefit-sharing, a formal request both to the National Park of the Pyrenees and to the different town halls were made in order to sample buds in the different populations analyzed.

All of the materials employed in this study are preserved in the -80 °C freezer at Biogeco research unit, Cestas, France.

### Consent for publication

Peoples visible in the Fig. 5B were involved in this study. The informed consent to publish was obtained from those participants.

### Competing interests

The authors declare that they have no competing interests.

Received: 29 September 2022 Accepted: 16 January 2023

Published online: 23 February 2023

## References

1. Kremer A, Ronce O, Robledo-Arnuncio JJ, Guillaume F, Bohrer G, Nathan R, et al. Long-distance gene flow and adaptation of forest trees to rapid climate change. Ecol Lett. 2012;15:378-92.

2. Kawecki TJ, Ebert D. Conceptual issues in local adaptation. Ecol Lett. 2004;7:1225-41.

3. Savolainen O, Pyhjärvi T, Knürr T. Gene flow and local adaptation in trees. Annu Rev Ecol Evol Syst. 2007;38:595-619.

4. Aitken SN, Yeaman S, Holliday JA, Wang T, Curtis-Mclane S. Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol Appl. 2008;1:95-111.

5. Alberto FJ, Derory J, Boury C, Frigerio J-M, Zimmermann NE, Kremer A. Imprints of natural selection along environmental gradients in phenology-related genes of Quercus petraea. Genetics. 2013;195:495-512.

6. Kremer A, Potts BM, Delzon S. Genetic divergence in forest trees: understanding the consequences of climate change. Funct Ecol. 2014;28:22-36.

7. White TL, Adams WT, Neale DB. Genetic markers - morphological, biochemical and molecular markers. Forest genetics. 2007. p. 53-74. https://doi.org/10.1079/9781845932855.00.

8. Ducousso A, Guyon JP, Krémer A. Latitudinal and altitudinal variation of bud burst in western populations of sessile oak (Quercus petraea (Matt) Liebl). Ann For Sci. 1996;53:775-82.

9. Vitasse Y, Delzon S, Dufrêne E, Pontailler J-Y, Louvet J-M, Kremer A, et al. Leaf phenology sensitivity to temperature in European trees: do within-species populations exhibit similar responses? Agric For Meteorol. 2009;149:735-44.

10. Leroy T, Louvet J-M, Lalanne C, Provost GL, Labadie K, Aury J-M, et al. Adaptive introgression as a driver of local adaptation to climate in European white oaks. New Phytol. 2020;226:1171-82.

11. Zanetto A, Kremer A. Geographical structure of gene diversity in Quercus petraea (Matt.) Liebl. I. Monolocus patterns of variation. Heredity. 1995;75:506-17.

12. Kremer A, Dupouey JL, Deans JD, Cottrell J, Csaikl U, Finkeldey R, et al. Leaf morphological differentiation between Quercus robur and Quercus petraea is stable across western European mixed oak stands. Ann For Sci. 2002;59:777-87.
13. Firmat C, Delzon S, Louvet J-M, Parmentier J, Kremer A. Evolutionary dynamics of the leaf phenological cycle in an oak metapopulation along an elevation gradient. J Evol Biol. 2017;30:2116–31.

14. Soularue J-P, Kremer A. Assortative mating and gene flow generate clinal phenological variation in trees. BMC Evol Biol. 2012;12:79.

15. Sork VL, Squire K, Gugger PF, Steele SE, Levy ED, Eckert AJ. Landscape genomic analysis of candidate genes for climate adaptation in a California endemic oak. Quercus lobata Am J Bot. 2016;103:33–46.

16. Rellstab C, Zoller S, Walthert L, Lesur I, Pluess AR, Graf R, et al. Signatures of local adaptation in candidate genes of oaks (Quercus spp.) with respect to present and future climatic conditions. Mol Ecol. 2016;25:5907–24.

17. Lobréaux S, Miquel C. Identification of Arabis alpina genomic regions associated with climatic variables along an elevation gradient through whole genome scan. Genomics. 2020;112:729–35.

18. Capblancq T, Morin X, Gueguen M, Renaud J, Lobreaux S, Bazin E. Climate-associated genetic variation in Fagus sylvatica and potential responses to climate change in the French Alps. J Evol Biol. 2020;33:783–96.

19. Mead A, Ramirez JP, Bartlett MK, Wright JW, Sack L, Sork VL. Seedling response to water stress in valley oak (Quercus lobata) is shaped by different gene networks across populations. Mol Ecol. 2019;28:5248–64.

20. Sork VL. Genomic studies of local adaptation in natural plant populations. J Hered. 2018;109:3–15.

21. Lasky JR, Des Marais DL, Lowry DB, Povolotskaya I, McKay JK, Richards JH, et al. Natural variation in abiotic stress responsive gene expression and local adaptation to climate in arabidopsis thaliana. Mol Biol Evol. 2014;31:2283–96.

22. Cooper RD, Shaffer HB. Allele-specific expression and gene regulation help explain transgressive thermal tolerance in non-native hybrids of the endangered California tiger salamander (Ambystoma californiense). Mol Ecol. 2021;30:987–1004.

23. Johnsen Ø, Fossdal CG, Nagy N, Mølmann J, DæHLEN OG, Skrøppa T. Climatic adaptation in Picea abies progenies is affected by the temperature during zygotic embryogenesis and seed maturation. Plant, Cell Environ. 2005;28:1090–102.

24. Rohde A, Bhalerao RP. Plant dormancy in the perennial context. Trends Plant Sci. 2007;12:217–23.

25. Balandier P, Bonhomme M, Rageau R, Capitan F, Parisot E. Leaf bud endodormancy release in peach trees: evaluation of temperature models in temperate and tropical climates. Agric For Meteorol. 1993;67:95–113.

26. Cline MG, Deppong DO. The role of apical dominance in parador-mancy of temperate woody plants: a reappraisal. J Plant Physiol. 1999;155:350–6.

27. Heide OM. Interaction of photoperiod and temperature in the control of growth and dormancy of Prunus species. Sci Hortic. 2008;115:309–14.

28. Tanino KK, Kalcits L, Silim S, Kendall E, Gray GR. Temperature-driven plasticity in growth cessation and dormancy development in deciduous woody plants: a working hypothesis suggesting how molecular and cellular function is affected by temperature during dormancy induction. Plant Mol Biol. 2010;73:49–65.

29. Maurya JP, Bhalerao RP. Photoperiod- and temperature-mediated control of growth cessation and dormancy in trees: a molecular perspective. Ann Bot. 2017;120:351–60.

30. Luedeling E, Girvetz EH, Semenov MA, Brown PH. Climate change affects winter chill for temperate fruit and nut trees. PLoS ONE. 2011;6:e20155.

31. Dantec CF, Vitasse Y, Bonhomme M, Louvet J-M, Kremer A, Delzon S. Chilling and heat requirements for leaf unfolding in European beech and sessile oak populations at the southern limit of their distribution range. Int J Biometeorol. 2014;58:1853–64.

32. Naor A, Faishman M, Stern R, Moshe A, Erez A. Temperature effects on dormancy completion of vegetative buds in apple. J Amer Soc Hort Sci. 2003;128:636–41.

33. Doi H, Katano I. Phenological timings of leaf budburst with climate change in Japan. Agric For Meteorol. 2008;148:512–6.

34. Schaber J, Badeck F-W. Physiology-based phenology models for forest tree species in Germany. Int J Biometeorol. 2003;47:193–201.

35. Lenoir J, Gégout JC, Marquet PA, de Ruffay P, Brisse H. A significant upward shift in plant species optimum elevation during the 20th century. Science. 2008;320:1768–71.

36. Lenoir J, Svenning J-C. Climate-related range shifts – a global multidimensional synthesis and new research directions. Ecography. 2015;38:15–28.

37. Olsen JE, Lee Y, Junttila O. Effect of alternating day and night temperature on short day-induced bud set and subsequent bud burst in long days in Norway spruce. Front Plant Sci. 2014;5:691.

38. Khalil-Ur-Rehman M, Sun L, Li C-X, Faheem M, Wang W, Tao J-M. Comparative RNA-seq based transcriptomic analysis of bud dormancy in grape. BMC Plant Biol. 2017;17:18.

39. Chao WS, Doğramacı M, Horvath DP, Anderson JV, Foley ME. Comparison of phytohormone levels and transcript profiles during seasonal dormancy transitions in underground adventitious buds of leafy spurge. Plant Mol Biol. 2017;94:281–302.

40. Min Z, Zhao X, Li R, Yang B, Liu M, Fang Y. Comparative transcriptome analysis provides insight into differentially expressed genes related to bud dormancy in grapevine (Vitis vinifera). Sci Hortic. 2017;225:213–20.

41. Paul A, Jha A, Bhardwaj S, Singh S, Shankar R, Kumar S. RNA-seq-mediated transcriptome analysis of actively growing and winter dormant shoots identifies non-deciduous habit of evergreen tree tea during winters. Sci Rep. 2014;4:5932.

42. Lesur I, Bechade A, Klopp C, Noirot C, Leplé JC, Kremer A, et al. A unigene set for European beech Fagus sylvatica L and its use to decipher the molecular mechanisms involved in dormancy regulation. Mol Ecol Res. 2015;15:1192. https://doi.org/10.1111/1755-0998.12373.

43. Takemura Y, Kuroki K, Shida Y, Araki S, Takeuchi Y, Tanaka K, et al. Comparative transcriptome analysis of the less-dormant Taiwanese pear and the dormant Japanese pear during winter season. PLoS ONE. 2015;10:e0139595.

44. Lesur I, Provost GL, Bento P, Silva CD, Leplé J-C, Murat F, et al. The oak gene expression atlas: insights into Fagaceae genome evolution and the discovery of genes regulated during bud dormancy release. BMC Genomics. 2015;16:112.

45. Ueno S, Klopp C, Leplé JC, Derory J, Noirot C, Léger V, et al. Transcriptional profiling of bud dormancy induction and release in oak by next-generation sequencing. BMC Genomics. 2013;14:236.

46. Karlgren A, Gyllenstrand N, Clapham D, Lagercrantz U. FLOWERING LOCUS T/TERMINAL FLOWER1-Like Genes Affect Growth Rhythm and Bud Set in Norway Spruce. Plant Physiol. 2013;163:792–803.

47. Santamaría ME, Rodríguez R, Cañal MJ, Toorop PE. Transcriptome analysis of chestnut (Castanea sativa) tree buds suggests a putative role for epigenetic control of bud dormancy. Ann Bot. 2011;108:485–98.

48. Rutink T, Arend M, Morreel K, Storme V, Rombauts S, Fromm J, et al. A molecular timetable for apical bud formation and dormancy induction in poplar. Plant Cell. 2007;19:2370–90.

49. Yordanov VS, Ma C, Strauss SH, Busov VB. EARLY BUD-BREAK 1 (EBB1) is a regulator of release from seasonal dormancy in poplar trees. PNAS. 2014;111:10001. https://doi.org/10.1073/pnas.1405621111.

50. Li W-F, Ding Q, Cui K-M, He X-Q. Cambium reactivation independent of bud unfolding involves de novo IAA biosynthesis in cambium regions in Populus tomentosa Carr. Acta Physiol Plant. 2013;35:1827–36.

51. Ingarsson PK, García MV, Hall D, Luquez V, Jansson S. Clinical variation in phyB2, a candidate gene for day-length-induced growth cessation and bud set, across a latitudinal gradient in European aspen (Populus tremula). Genetics. 2006;172:1845–53.

52. Hall D, Luquez V, Garcia VM, St Onge KR, Jansson S, Ingarsson PK. Adaptive population differentiation in phenology across a latitudinal gradient in European aspen (Populus tremula, L.): a comparison of neutral markers, candidate genes and phenotypic traits. Evolution. 2007;61:2849–60.

53. Li J, Xu Y, Niu Q, He L, Teng Y, Bai S. Abscisic Acid (ABA) promotes the induction and maintenance of pear (Pyrus pyrifolia White Pear Group) flower bud endodormancy. Int J Mol Sci. 2018;19:310.

54. Dierck R, Keyser ED, Riek JD, Dhooghe E, Huylenbroeck JV, Prinsen E, et al. Change in auxin and cytokinin levels coincides with altered expression of branching genes during axillary bud outgrowth in chrysanthemum. PLoS ONE. 2016;11:e0161732.
55. Horvath D. Bud Dormancy and Growth. Plant Developmental Biology - Biotechnological Perspectives. Springer Berlin Heidelberg; 2010. p. 53-70.

56. Zhong W, Gao Z, Zhuang W, Shi T, Zhang Z, Ni Z. Genome-wide expression profiles of seasonal bud dormancy at four critical stages in Japanese apricot. Plant Mol Biol. 2013;83:247-64.

57. Shim D, Ko J-H, Kim W-C, Wang Q, Keathley DE, Han K-H. A molecular framework for seasonal growth-dormancy regulation in perennial plants. Hortic Res. 2014;1:14059.

58. Allamano P, Claps P, Laio F. Global warming increases flood risk in mountainous areas. Geophys Res Lett. 2009;36(24):1-5.

59. Tylianakis JM, Didham RK, Bascompte J, Wardle DA. Global change and species interactions in terrestrial ecosystems. Ecol Lett. 2008;11:1351-63.

60. Kim HB, Lee H, Oh CJ, Lee H-Y, Eum HL, Kim H-S, et al. Postembryonic seedling lethality in the sterol-deficient arabidopsis cyp51A2 mutant is partially mediated by the composite action of ethylene and reactive oxygen species. Plant Physiol. 2010;152:192-205.

61. Kim SY, Shang Y, Joo S-H, Kim S-K, Nam KH. Overexpression of BAK1 causes salicylic acid accumulation and deregulation of cell death control genes. Biochem Biophys Res Commun. 2017;484:781-6.

62. Yang X, Yu Y, Jiang L, Lin X, Zhang C, Ou X, et al. Changes in DNA methylation and transgenerational mobilization of a transposable element (mPing) by the Topoisomerase II inhibitor, Etoposide, in rice. BMC Plant Biol. 2012;12:48.

63. Moschen S, BengoaLuoni S, Paniego NB, Hopp HE, Dosio GAA, Fernandez P, et al. Identification of candidate genes associated with leaf senescence in cultivated sunflower (Helianthus annuus L.). PLoS One. 2014;9:e104379.

64. Fang S-C, Fernandez DE. Effect of regulated overexpression of the MADS domain factor AGL15 on flower senescence and fruit maturation. Plant Physiol. 2002;130:78-89.

65. Zhao Y, Cheng S, Song Y, Huang Y, Zhou S, Liu X, et al. The interaction between rice ERF3 and WOX11 promotes crown root development by regulating gene expression involved in cytokinin signaling. Plant Cell. 2015;27:2469-83.

66. Tan H, Man C, Xie Y, Yan J, Chu J, Huang J. A crucial role of GA-regulated flavonol biosynthesis in root growth of arabidopsis. Mol Plant. 2019;12:521-37.

67. Burke JJ. Identification of genetic diversity and mutations in higher plant acquired thermotolerance. Physiol Plant. 2001;112:167-70.

68. Stief A, Altmann S, Hoffmann K, Pant BD, Scheible W-R, Bäurle I. Arabidopsis miR156 regulates tolerance to recurring environmental stress through SPL transcription factors. Plant Cell. 2014;26:1792-807.

69. Sharma M, Banday ZZ, Shukla BN, Laxmi A. Glucose-regulated HLP1 acts as a key molecule in governing thermomemory. Plant Physiol. 2019;180:1081-100.

70. Bresson CC, Vitasse Y, Kremer A, Delzon S. To what extent is altitudinal variation of functional traits driven by genetic adaptation in European oak and beech? Tree Physiol. 2011;31:1164-74.

71. DeYoung BJ, Bickle KL, Schrage KJ, Muskett P, Patel K, Clark SE. The CLAVATA1-related BAM1, BAM2 and BAM3 receptor kinase-like proteins are required for meristem function in Arabidopsis. Plant J. 2006;45:1-16.

72. Ren H, Gray WM. SAUR Proteins as effectors of hormonal and environmental signals in plant growth. Mol Plant. 2015;8:1153-64.

73. Gargul JM, Mibus H, Serek M. Characterization of transgenic Kalanchoë and petunia with organ-specific expression of GUS or GA2ox genes led by the deletion BOX-I version (dBI) of the PAL1 promoter. J Plant Growth Regul. 2017;36:424-35.

74. Ayala-Rodríguez JA, Barrera-Ortiz S, Ruiz-Herrera LF, López-Bucio J. Folic acid orchestrates root development linking cell elongation with auxin response and acts independently of the TARGET OF RAPAMYCIN signaling in arabidopsis thaliana. Plant Sci. 2017;264:168-78.

75. Racolta A, Nodine MD, Davies K, Lee C, Rowe S, Velazco Y, et al. A common pathway of root growth control and response to CLE peptides through two receptor kinases in arabidopsis. Genetics. 2018;208:687-704.

76. Zhu Y, Shi Z, Li S, Liu H, Liu F, Niu Q, et al. Fine mapping of the novel male-sterile mutant gene ms39 in maize originated from outer space flight. Mol Breeding. 2018;38:125.

77. Fleet CM, Sun T. A DELLAcate balance: the role of gibberellin in plant morphogenesis. Curr Opin Plant Biol. 2005;8:77-85.

78. Hamant O, Pautot V. Plant development: a TALE story. CR Biol. 2010;333:371-81.

79. Caignard T, Kremer A, Bouteiller XP, Parmentier J, Louvet J-M, Venner S, et al. Counter-gradient variation of reproductive effort in a widely distributed temperate oak. Funct Ecol. 2021;35:1745-55.

80. Vitasse Y, Delzon S, Bresson CC, Michalet R, Kremer A. Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden. Can J Forest Res. 2009;39:1259.

81. Wang Y, Li L, Ye T, Lu Y, Chen X, Wu Y. The inhibitory effect of ABA on floral transition is mediated by ABI5 in arabidopsis. J Exp Bot. 2013;64:675-84.

82. Jue D, Sang X, Liu L, Shu B, Wang Y, Liu C, et al. Identification of WRKY gene family from dimocarpus longan and its expression analysis during flower induction and abiotic stress responses. Int J Mol Sci. 2018;19:2169.

83. He S-B, Wang W-X, Zhang J-Y, Xu F, Lian H-L, Li L, et al. The CNT1 domain of arabidopsis CRY1 alone is sufficient to mediate blue light inhibition of hypocotyl elongation. Mol Plant. 2015;8:822-5.

84. Gao R, Gruber MY, Amyot L, Hannoufa A. SPL13 regulates shoot branching and flowering time in Medicago sativa. Plant Mol Biol. 2018;96:119-33.

85. Footitt S, Müller K, Kermode AR, Finch-Savage WE. Seed dormancy cycling in Arabidopsis: chromatin remodelling and regulation of DOG1 in response to seasonal environmental signals. Plant J. 2015;81:413-25.

86. Wang C, Yan X, Chen Q, Jiang N, Fu W, Ma B, et al. Clathrin light chains regulate clathrin-mediated trafficking, auxin signaling, and development in arabidopsis. Plant Cell. 2013;25:499-516.

87. Bouré N, Kumar SV, Arnaud N. The BAP module: a multisignal integrator orchestrating growth. Trends Plant Sci. 2019;24:602-10.

88. Bancos S, Nomura T, Sato T, Molnár G, Bishop GJ, Koncz C, et al. Regulation of transcript levels of the Arabidopsis cytochrome p450 genes involved in brassinosteroid biosynthesis. Plant Physiol. 2002;130:504-13.

89. Liu Y, Li X, Ma D, Chen Z, Wang J-W, Liu H. CIB1 and CO interact to mediate CRY2-dependent regulation of flowering. EMBO Rep. 2018;19:e45762.

90. Koehler G, Wilson RC, Goodpaster JV, Sønsteby A, Lai X, Witzmann FA, et al. Proteomic study of low-temperature responses in strawberry cultivars (Fragaria × ananassa) that differ in cold tolerance. Plant Physiol. 2012;159:1787-805.

91. Yang Q-S, Wu J-H, Li C-Y, Wei Y-R, Sheng O, Hu C-H, et al. Quantitative proteomic analysis reveals that antioxidant mechanisms contribute to cold tolerance in plantain (Musa paradisiaca L.; ABB group) seedlings. Mol Cell Proteomics. 2012;11:1853-69.

92. Wang L-J, Li S-H. Salicylic acid-induced heat or cold tolerance in relation to Ca2+ homeostasis and antioxidant systems in young grape plants. Plant Sci. 2006;170:685-94.

93. Shirzadian-Khorramabad R, Jing HC, Everts GE, Schippers JHM, Hille J, Dijkwel PP. A mutation in the cytosolic O acetylserine (thiol) lyase induces a genome-dependent early leaf death phenotype in Arabidopsis. BMC Plant Biol. 2010;10:80.

94. Fujino K, Matsuda Y. Genome-wide analysis of genes targeted by qLTG3-1 controlling low-temperature germinability in rice. Plant Mol Physiol. 2009;72:137.

95. Cabello JV, Arce AL, Chan RL. The homologous HD-Zip I transcription factors HaHB1 and AtHB13 confer cold tolerance via the induction of pathogenesis-related and glucanase proteins. Plant J. 2012;69:141-53.

96. Kim JY, Park SJ, Jang B, Jung C-H, Ahn SJ, Goh C-H, et al. Functional characterization of a glycine-rich RNA-binding protein 2 in Arabidopsis thaliana under abiotic stress conditions. Plant J. 2007;50:439-51.

97. Gu L, Xu T, Lee K, Lee KH, Kang H. A chloroplast-localized DEAD-box RNA helicase AtRH3 is essential for intron splicing and plays an important role in the growth and stress response in Arabidopsis thaliana. Plant Physiol Biochem. 2014;82:309-18.

98. Noir S, Marrocco K, Masoud K, Thomann A, Gusti A, Bitrian M, et al. The control of arabidopsis thaliana growth by cell proliferation and Endoreplication requires the F-box protein FBL17. Plant Cell. 2015;27:1461-76.
99. Wang J, Andersson-Gunnerås S, Gaboreanu I, Hertzberg M, Tucker MR, Zheng B, et al. Reduced expression of the SHORT-ROOT gene increases the rates of growth and development in hybrid poplar and arabidopsis. PLoS One. 2011;6:e28878.

100. Wilson MD, Wang D, Wagner R, Breysens H, Gertsenstein M, Lobe C, et al. ARS2 is a conserved eukaryotic gene essential for early mammalian development. Mol Cell Biol. 2008;28:1503–14.

101. Lee S, Cho D-I, Kang J, Kim M-D, Kim SY. AtNEK6 interacts with ARIA and is involved in ABA response during seed germination. Mol Cells. 2010;29:559–66.

102. Cao S, Chen H, Zhang C, Tang Y, Liu J, Qi H. Heterologous expression and biochemical characterization of two Lipoxygenases in oriental melon, Cucumis melo var. Makuwa Makino. PLoS One. 2016;11:e0153801.

103. Yamauchi Y, Takeda-Kamiya N, Hanada A, Ogawa M, Kuwahara A, Seo M, et al. Contribution of gibberellin deactivation by AtGA2ox2 to the suppression of germination of dark-imbibed arabidopsis thaliana seeds. Plant Cell Physiol. 2007;48:555–61.

104. Wahl V, Brand LH, Guo Y-L, Schmid M. The FANTASTIC FOUR proteins influence shoot meristem size in arabidopsis thaliana. BMC Plant Biol. 2010;10:1–12.

105. Okamuro JK, Caster B, Villarroel R, Montagu MV, Jofuku KD. The AP2 domain of APETALA2 defines a large new family of DNA binding proteins in arabidopsis. PNAS. 1997;94:7076–81.

106. Petty LM, Harberd NP, Carré IA, Thomas B, Jackson SD. Expression of the Arabidopsis gai gene under its own promoter causes a reduction in plant height in chrysanthemum by attenuation of the gibberellin response. Plant Sci. 2003;164:175–82.

107. Nieuwland J, Feron R, Huisman BAH, Fasolino A, Hilbers CW, Derksen J, et al. Lipid transfer proteins enhance cell wall extension in tobacco. Plant Cell. 2005;17:2009–19.

108. Vitasse Y, Porté AJ, Kremer A, Michalet R, Delzon S. Responses of canopy duration to temperature changes in four temperate tree species: relative contributions of spring and autumn leaf phenology. Oecologia. 2009;161:187–98.

109. Alberto F, Bouffier L, Louvet J-M, Lamy J-B, Delzon S, Kremer A. Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient. J Evol Biol. 2011;24:1442–54.

110. Horvath DP, Anderson JV, Chao WS, Foley ME. Knowing when to grow: signals regulating bud dormancy. Trends Plant Sci. 2003;8:534–40.

111. Yang Q, Gao Y, Wu X, Moriguchi T, Bai S, Teng Y. Bud endodormancy in deciduous fruit trees: advances and prospects. Hortic Res. 2021;8:1–11.

112. Guichoux E, Lagache L, Wagner S, Léger P, Petit RJ. Two highly validated multiplexes 12-plex and 8-plex for species delimitation and parentage analysis in oaks (Quercus spp.). Mol Ecol Resour. 2011;11:578–85.

113. Björklund S, Antti H, Uddestrand I, Moritz T, Sundberg B. Cross-talk between gibberellin and auxin in development of Populus wood: gibberellin stimulates polar auxin transport and has a common transcriptome with auxin. Plant J. 2007;52:499–511.

114. Le Provost G, Herrera R, Paiva JA, Chaumeil P, Salin F, Plomion C. A micromethod for high throughput RNA extraction in forest trees. Biol Res. 2007;40:291–7.

115. Le Provost G, Lesur J, Lalanne C, Da Silva C, Labadie K, Aury JM, et al. Implication of the suberin pathway in adaptation to waterlogging and hypertrophied lenticels formation in pedunculate oak Quercus robur L. Tree Physiol. 2016;36:1330. https://doi.org/10.1093/treephys/tpw056.

116. Plomion C, Aury J-M, Amselem J, Leroy T, Murat F, Duplessis S, et al. Oak genome reveals facets of long lifespan. Nature Plants. 2018;4:440.

117. Soularue J-P, Firmat C, Caignard T, Thöni A, Arnoux L, Delzon S, et al. Antagonistic effects of assortative mating on the evolution of phenotypic plasticity along environmental gradients. 2022. https://hal.archives-ouvertes.fr/hal-03813187.

118. Sharan R, Maron-Katz A, Shamir R. CLICK and EXPANDER: a system for clustering and visualizing gene expression data. Bioinformatics. 2003;19:1787–99.

119. Tavazoie S, Hughes JD, Campbell MJ, Cho RJ, Church GM. Systematic determination of genetic network architecture. Nat Genet. 1999;22:281–5.

120. Alexa A. Gene set enrichment analysis with topGO. R packageversion 2.8.0. Springer Berlin Heidelberg.

121. Le Provost G, Sulmon C, Frigerio JM, Bodénès C, Kremer A, Plomion C. Role of waterlogging-responsive genes in shaping interspecific differentiation between two sympatric oak species. Tree Physiol. 2012;32:119–34.

122. Rozen S, Skaletsky H. Primer3 on the WWW for general users and for biologist programmers. Methods Mol Biol. 2000;132:365–86.

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

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