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

HAL
open science

# Genome-wide evolutionary response of European oaks during the Anthropocene

Dounia Saleh, Jun Chen, Jean-Charles Leplé, Thibault Leroy, Laura Truffaut, Benjamin Dencausse, Céline Lalanne, Karine Labadie, Isabelle Lesur, Didier Bert, et al.

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

## To cite this version:

Dounia Saleh, Jun Chen, Jean-Charles Leplé, Thibault Leroy, Laura Truffaut, et al., Genome-wide evolutionary response of European oaks during the Anthropocene. Evolution Letters, 2022, 6 (1), pp.4-20. 10.1002/evl3.269. hal-03553889

HAL Id: hal-03553889

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

Submitted on 3 Feb 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
doi:10.1002/evl3.269

# Genome-wide evolutionary response of European oaks during the Anthropocene

Dounia Saleh,<fnref n="1" />,* Jun Chen,<fnref n="2" />,* ID Jean-Charles Leplé,<fnref n="1" /> ID Thibault Leroy,<fnref n="3" /> ID Laura Truffaut,<fnref n="1" />

Benjamin Dencausse,<fnref n="1" /> Céline Lalanne,<fnref n="1" /> Karine Labadie,<fnref n="4" /> ID Isabelle Lesur,<fnref n="5" /> Didier Bert,<fnref n="1" /> ID

Frédéric Lagane,<fnref n="1" /> François Morneau,<fnref n="6" />,<fnref n="7" /> ID Jean-Marc Aury,<fnref n="8" /> ID Christophe Plomion,<fnref n="1" /> ID Martin Lascoux,<fnref n="9" /> ID and Antoine Kremer<fnref n="1" />,<fnref n="10" /> ID

Received May 25, 2021

Accepted December 2, 2021

The pace of tree microevolution during Anthropocene warming is largely unknown. We used a retrospective approach to monitor genomic changes in oak trees since the Little Ice Age (LIA). Allelic frequency changes were assessed from whole-genome pooled sequences for four age-structured cohorts of sessile oak (*Quercus petraea*) dating back to 1680, in each of three different oak forests in France. The genetic covariances of allelic frequency changes increased between successive time periods, highlighting genome-wide effects of linked selection. We found imprints of parallel linked selection in the three forests during the late LIA, and a shift of selection during more recent time periods of the Anthropocene. The changes in allelic covariances within and between forests mirrored the documented changes in the occurrence of extreme events (droughts and frosts) over the last 300 years. The genomic regions with the highest covariances were enriched in genes involved in plant responses to pathogens and abiotic stresses (temperature and drought). These responses are consistent with the reported sequence of frost (or drought) and disease damage ultimately leading to the oak dieback after extreme events. They provide support for adaptive evolution of long-lived species during recent climatic changes. Although we acknowledge that other sources (e.g., gene flow, generation overlap) may have contributed to temporal covariances of allelic frequency changes, the consistent and correlated response across the three forests lends support to the existence of a systematic driving force such as natural selection.

**KEY WORDS:** Anthropocene, evolution, linked selection, Little Ice Age, *Quercus petraea*.

<footnote n="1">¹UMR BIOGECO, INRAE, Université de Bordeaux, Cestas 33612, France</footnote>

<footnote n="2">²College of Life Sciences, Zhejiang University, Hangzhou 310058, China</footnote>

<footnote n="3">³Department of Botany and Biodiversity Research, University of Vienna, Vienna 1010, Austria</footnote>

<footnote n="4">⁴Genoscope, Institut de Biologie François Jacob, Commissariat à l'énergie atomique (CEA), Université de Paris-Saclay, Evry 91057, France</footnote>

<footnote n="5">⁵Helix Venture, Merignac 33700, France</footnote>

<footnote n="6">⁶Département Recherche Développement Innovation, Office National des Forêts, Boigny-Sur-Bionne 45760, France</footnote>

<footnote n="7">⁷Current Address: Service de l'Information Statistique Forestière et Environnementale, Institut National de l'Information géographique et Forestière, Nogent-sur-Vernisson 45290, France</footnote>

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

<footnote n="9">⁹Department of Ecology and Genetics, Evolutionary Biology Centre, Uppsala University, Uppsala SE-75236, Sweden</footnote>

<footnote n="10">¹⁰E-mail: antoine.kremer@inrae.fr</footnote>
## Impact Summary

A highly debated issue today in ecology is how fast evolutionary changes can occur, especially for long-lived organisms, and if this pace is sufficient to cope with ongoing climatic changes. Climate warming started after the Little Ice Age (LIA), at the mid-19th century; thus, multicentennial trees offer a unique opportunity to explore microevolution during documented climatic transitions. We implemented a retrospective approach in European white oak populations using whole genome sequences to track traces of selection. We found genome-wide evidence of parallel linked selection during the late LIA, and a shift of selection during more recent time periods. Our results provide the first evidence that selection has operated on long-lived trees over short time periods, and will likely continue to do so. Beyond oaks, other woody species sharing similar life history traits and attributes (large standing genetic, polygenicity of fitness, high fecundity, and severe selection screening at the juvenile stage) may as well be prone to linked selection and allow retrospective tracking of evolutionary and adaptive pathways. Such retrospective approaches may improve our understanding of future responses to ongoing climatic changes.

Our ability to forecast the response of species to climate change is limited by our lack of knowledge about the pace of adaptive evolution, particularly in long-lived species, such as trees. Multicentennial trees in the Northern Hemisphere that are old enough to have experienced the transition from the cold Little Ice Age (LIA; 1450–1850) to the warm Anthropocene (1850–today) (Luterbacher et al. 2004; Corona et al. 2010; Luterbacher et al. 2016; Anchukaitis et al. 2017) offer a unique opportunity to assess the extent of evolutionary responses to this extreme and well-described environmental challenge. The LIA (Tkachuck 1983) was a cold period characterized by climatic extremes such as long and harsh winters, but also severe droughts that led to plagues, famines, and ultimately revolutions (Pfister 1984; Fagan 2002; Parker 2013). The consequences of the LIA for plants are best illustrated by the recurrent poor crop harvests (Le Roy Ladurie 2004, 2006). Evidence of the impact of the LIA on forest trees is provided by comparisons of tree ring sizes and historical temperature records (Carrer and Urbinati 2006; Edouard et al. 2009). The decrease in temperature during the LIA resulted in a retreat of the forest tree line at high latitudes (Kullman 2005; MacDonald et al. 2008; Linderholm et al. 2014; Kullman 2015; Helama et al. 2020) or altitudes (Camarero et al. 2015)

and changes in species composition (Campbell and McAndrews 1993).

Here, we investigated whether climatic trends during and after the LIA and occurrences of extreme events had evolutionary consequences for forest tree populations. We addressed two major questions in this work. First, we investigated whether extreme events occurring during the late LIA left a genomic signature. Second, we investigated whether the shift in climate after the LIA also left a genomic imprint. We chose sessile oak (*Quercus petraea* (Matt.) Liebl.) as the model species for this study, as the oldest forests in Europe contain sessile oak stands that came into existence in the middle of the LIA (near 1650), a time at which the French statesman Colbert implemented even-aged management in French forests (Gallon 1752). Four age-structured cohorts of roughly 340, 170, 60, and 12 years old were sampled within three forests, to explore changes in allele frequencies over time. *Quercus petraea* is known to display considerable genetic diversity (Mariette et al. 2002; Kremer and Hipp 2020; Leroy et al. 2020) and to have a high genetic variance for fitness (Alexandre et al. 2020). The selection for viability or the demographic dynamics generated by extreme weather events would be expected to result in changes in allelic frequency.

Climate change-driven evolution over the course of a few generations has mostly been reported in invasive species (Chown et al. 2015) or in controlled experiments (Ravenscroft et al. 2015). By contrast, in this study, we explored the ability of a native species with high levels of standing genetic variation to respond to recent documented climatic changes. Our approach mirrors experiments monitoring the change in allelic frequencies over successive discrete generations in natura (Malaspinas 2016) or in controlled selection experiments (Schlotterer et al. 2015), except that our study made use of contemporary age-structured cohorts. Earlier genome-wide investigations performed with a synchronous approach in common gardens highlighted the multifaceted and unrepeatable signatures of natural selection, characterized by heterogeneous polygenicity (trait architecture determined by a large number of genes with small effects) and allelic heterogeneity (Alberto et al. 2011; Plomion et al. 2016; Rellstab et al. 2016).

In two recent publications, Buffalo and Coop (2019, 2020) showed that allele frequency trajectories can be shifted due to linked selection between selected and neutral loci and generate covariances between allelic frequency changes at successive time periods. Theoretically, the magnitude of covariances depends on the genetic variance of fitness, recombination, and linkage disequilibrium between selected and neutral loci and the strength of selection, whereas the sign of the covariances depends on the maintenance or fluctuation of selection pressures over time. We hypothesized that this approach would be suitable for detecting genomic footprints

<footnote>*These two authors contributed equally to the work.</footnote>
of selection in the past, during the late LIA, and shifts in selection pressures due to warming after the end of the LIA. We also explored a more qualitative approach addressing the underlying functions of the multiple genes contributing to the covariances of allelic frequency changes. We thus had three objectives: (1) to retrace trends for the covariances of allelic frequencies between age-structured cohorts spanning the last three centuries in oak stands, (2) to determine whether these trends were repeatable over replicated observations in three different forests, and (3) to explore the gene networks involved in these responses to environmental change.

# Methods

## SAMPLING FORESTS AND AGE STRUCTURED COHORTS

We sampled three oak forests located in the central and western part of France (Bercé, Réno-Valdieu, and Tronçais; Fig. 1A). These forests include stands of up to 349 years of age when the study started in 2014, and are managed under even-aged silvicultural regimes (Supporting Information S1). The upper canopy consisted principally of *Quercus petraea*, and historical records and genetic evidence (based on chloroplast DNA haplotypes) indicated that the three forests were of natural origin (Petit et al. 2002). In each forest, we sampled individuals belonging to four age-class cohorts corresponding to ages of 340, 170, 60, and 12 years (born approximatively in 1680, 1850, 1960, and 2008) (Fig. 1; Table S1). These cohorts are referred to as cohorts 4, 3, 2, and 1, respectively, below. The regeneration period, during which mature trees mate and the stand is renewed by natural seeding, takes today about 10–20 years, but extended over longer periods in the past (up to 30 years). The age of the trees in a given cohort may, therefore, vary by up to 10–30 years. Cohorts were dated on the basis of management records, together with dendrochronological recordings for a few felled trees within each cohort (Fig. 1B). Knowledge of the demographic dynamics of even-aged forests is required to identify the periods during which selection was at its strongest. When a stand is renewed by natural seeding, a very dense cohort of seedlings develops (more than 100,000 seedlings/ha), the number of plants gradually decreasing to about 4000/ha by the age of 10 years, as a result of natural selection and chance events. Crucially, subsequent silvicultural thinning is applied to only the remaining 4% of the trees. Hence, an oak stand that is about 340 years old today probably underwent its strongest bout of natural selection in the late 17th century, when it was at the seedling stage. Such reasoning provided the rationale for sampling age-structured cohorts in even-aged stands for retrospective monitoring of the selective impact of past climatic changes. The different cohorts within a given forest are derived from the same founding population established more than

10,000 years ago (Giesecke 2016; Giesecke and Brewer 2018), but there is no direct traceable generation-to-generation link between the cohorts. The three forests as well share a common evolutionary history as they all derived from the iberic glacial refugial during the Last Glacial Maximum (Petit et al. 2002). Therefore, from an evolutionary stand point, the three forests cannot be considered as independent replicates.

High-resolution regional and temporal temperature reconstructions based on a combination of instrumental data (for the most recent period), documentary records, ice core, and tree ring proxy data revealed clear temperature trends over the last 340 years in the three forests (Luterbacher et al. 2004, 2016) (Supporting Information S2 and Figs. S1–S3).

## DNA EXTRACTION AND SEQUENCING

Leaf samples were collected from cohort 1 in the spring and summer 2014, whereas cambium samples were harvested at the base of the trunk with a punch for the three older cohorts. The sample size for each cohort was between 40 and 68 trees (Fig. 1B). DNA was extracted from the 639 samples with a Qiagen extraction kit. DNA yields were measured with a NanoDrop 1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and DNA samples were pooled in equimolar amounts within each cohort from each forest. The 12 pools (3 forests × 4 cohorts) were sequenced on an Illumina HiSeq4000 sequencer generating 150 bp paired-end reads.

## SNP DISCOVERY

The various steps in the SNP calling pipeline were as described in Altmann et al. (2012) and Pfeifer (2017). The adapters and primers were removed and reads were trimmed to remove nucleotides with a quality value below 20 from the two ends. The sequences between the second unknown nucleotide and the end of the read were removed, and reads of less than 30 nucleotides in length were discarded. Finally, the read pairs from the low-concentration spike-in Illumina PhiX Control library were removed. We obtained a mean of 348,863,070 reads per pool (Table 1). The reads were mapped onto the version 2.3 *Q. robur* genome assembly (Plomion et al. 2018) with bwa-mem, with a seed size of 39 (Li and Durbin 2010). Incorrectly paired reads and reads giving multiple alignments were removed with samtools (Li et al. 2009). Duplications were removed with Picard tools (no publication, Broad Institute). Base Alignment Quality (BAQ) was calculated with samtools (Li et al. 2009). Pileup files were generated for each scaffold over all forests, with samtools. These files were converted into synchronized pileups with a minimum alignment quality of 10, and allele frequencies were calculated for SNPs with a minimum count of two for minor alleles, a minimum coverage of 40×, and a maximum coverage of 10% of total coverage within each pool, with Popoolation2
A Geographic distribution of age structured cohorts in three even-aged managed forests

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

B Age, area, and size of age structured cohorts

<table><thead><tr><th>Forest</th><th>Cohort</th><th>Dates of birth according to management documents</th><th>Dates of birth according to dendrochronological data *</th><th>Surface of compartment (ha)</th><th>Sample size</th></tr></thead><tbody><tr><td rowspan="4">Bercé</td><td>B1</td><td>2003-2013</td><td></td><td>20</td><td>53</td></tr><tr><td>B2</td><td>1957</td><td>1954-1963 (3)</td><td>28</td><td>53</td></tr><tr><td>B3</td><td>1855</td><td>1856-1858 (3)</td><td>28</td><td>44</td></tr><tr><td>B4</td><td>1687</td><td>1675-1713 (18)</td><td>21</td><td>40</td></tr><tr><td rowspan="4">Réno-Valdieu</td><td>R1</td><td>2001-2010</td><td></td><td>10</td><td>55</td></tr><tr><td>R2</td><td>1959</td><td>1959-1973 (26)</td><td>10</td><td>68</td></tr><tr><td>R3</td><td>1854</td><td>1820-1857 (3)</td><td>12</td><td>54</td></tr><tr><td>R4</td><td>1704</td><td>1665-1731 (10)</td><td>13</td><td>56</td></tr><tr><td rowspan="4">Tronçais</td><td>T1</td><td>2002-2010</td><td></td><td>35</td><td>55</td></tr><tr><td>T2</td><td>1961</td><td>1953-1962 (10)</td><td>19</td><td>55</td></tr><tr><td>T3</td><td>1833</td><td>1841-1852 (3)</td><td>16</td><td>53</td></tr><tr><td>T4</td><td>1671</td><td>1654-1677 (4)</td><td>28</td><td>53</td></tr></tbody></table>

**Figure 1.** Sampling of forests and age structured cohorts of sessile oak. (A) Distribution of age-structured cohorts of sessile oak in three even-aged managed national forests in France. Each forest is subdivided in compartments (about 20 ha in size) limited by the black lines on the forest maps. Age class compartments are evenly distributed in the forests. Densities are extremely high at the seedling stage (>100,0000/ha) and decrease very rapidly due to natural selection during the early stage (≈4000 at age 10). (B) Age, area, and size of age-structured cohorts. Dendrochronological data of tree rings on felled trees in each cohort allowed to confirm documentary records of tree ages. About 50 trees were randomly sampled in each cohort for whole genome sequencing.

<footnote>*Minimum and maximum dates assessed on the rings of the stump. Numbers within brackets indicate the number of trees felled for the assessment.</footnote>
<table><caption>**Table 1.** SNP diversity statistics of the age structured cohorts.</caption><thead><tr><th>Forest</th><th>Cohort</th><th>Number of reads after postprocessing</th><th>Number of SNPs</th><th>π ± SD</th></tr></thead><tbody><tr><td rowspan="4">Bercé</td><td>B1</td><td>427,365,137</td><td>13,277,388</td><td>0.01202 ± 5 × 10<sup>-5</sup></td></tr><tr><td>B2</td><td>444,637,716</td><td>13,533,867</td><td>0.01208 ± 5 × 10<sup>-5</sup></td></tr><tr><td>B3</td><td>450,468,871</td><td>13,586,680</td><td>0.01334 ± 8 × 10<sup>-5</sup></td></tr><tr><td>B4</td><td>439,841,819</td><td>13,334,778</td><td>0.01447 ± 1 × 10<sup>-6</sup></td></tr><tr><td rowspan="4">Réno-Valdieu</td><td>R1</td><td>421,851,747</td><td>13,344,026</td><td>0.01211 ± 5 × 10<sup>-5</sup></td></tr><tr><td>R2</td><td>443,414,110</td><td>13,678,790</td><td>0.01256 ± 5 × 10<sup>-5</sup></td></tr><tr><td>R3</td><td>435,807,782</td><td>13,518,651</td><td>0.01203 ± 5 × 10<sup>-5</sup></td></tr><tr><td>R4</td><td>440,543,051</td><td>13,605,834</td><td>0.01206 ± 5 × 10<sup>-5</sup></td></tr><tr><td rowspan="4">Tronçais</td><td>T1</td><td>543,901,014</td><td>15,592,854</td><td>0.01351 ± 7 × 10<sup>-5</sup></td></tr><tr><td>T2</td><td>440,923,574</td><td>13,457,837</td><td>0.01208 ± 5 × 10<sup>-5</sup></td></tr><tr><td>T3</td><td>437,816,360</td><td>13,704,258</td><td>0.01215 ± 5 × 10<sup>-5</sup></td></tr><tr><td>T4</td><td>432,251,274</td><td>13,514,365</td><td>0.01208 ± 5 × 10<sup>-5</sup></td></tr></tbody></table>

(Kofler et al. 2011). Further filtering was applied to the three types of pileup files, to select biallelic SNPs with a minimum minor allele frequency of 0.02.

## DIVERSITY WITHIN AND BETWEEN COHORTS

Genetic diversity was estimated on nonoverlapping genomic windows of 10 kb spanning the whole oak genome. For each pool, we therefore generated new pileup files, including monomorphic sites. To take into account the variance in pool size and coverage among populations, we used the *subsample-pileup.pl* script from Popoolation2 (Kofler et al. 2011) to target a minimum coverage of 30 at each position and for each pool. These parameters were chosen after optimization to minimize the amount of missing data between pools. Tajima's π (Tajima 1989) was then calculated with the *variance-sliding.pl* script. As for the SNP sets, sites with a minimum alignment quality of 10 were retained, and were considered polymorphic if at least two copies of the minor allele were detected among all reads. Windows with a minimum covered fraction of 50% were retained for the calculation of Tajima's π. F<sub>ST</sub> values were calculated for each SNP, between each pair of cohorts, using Popoolation2 (Kofler et al. 2011), and were averaged over the whole genome.

## TEMPORAL COVARIANCES

Buffalo and Coop (2019, 2020) have developed a method for testing for genomic signals of selection on polygenic traits due to linked selection, based on temporal covariances of allelic frequency changes. Allele frequency changes were calculated for different time spans separating the cohorts in the three replicated forests. We used the CVTK software package available from http://github.com/vsbuffalo/cvtk to calculate the temporal co-

variances between allelic frequency changes (Buffalo and Coop 2020). Covariances were standardized by sample heterozygosity.

We filtered SNPs by removing sites with a depth below the number of alleles sampled for each cohort, and with a minor allele frequency below 0.02. Contigs were also filtered out when shorter than 200 kb long (excluding “N”s). Read depth and sample size were used to correct for bias in variance estimates (Buffalo and Coop 2020). Corrections were also done for bias caused by sampling noise common to adjacent time points. This bias stems from the allele frequency at a given time point common to the terms of the covariances between adjacent time periods (e.g., for Cov(Δ<sub>1680-1850</sub>, Δ<sub>1850-1960</sub>), allele frequency at 1850 is included in both terms).

Genome-wide temporal covariances were then computed between pairs of nonoverlapping time spans within each forest. Covariances were calculated based on all SNPs, and for all three forests separately. The mean and 95% confidence intervals of the covariance were obtained by bootstrapping, with 5000 iterations.

The change in the variance of allelic frequency changes over the whole period, from 1680 to 2008, becomes:

$$ 
\begin{aligned}
V \Delta_{1680-2008} = & V \Delta_{1680-1850} + V \Delta_{1850-1960} + V \Delta_{1960-2008} \\ 
& + 2 \text{Cov}(\Delta_{1680-1850}, \Delta_{1850-1960}) \\ 
& + 2 \text{Cov}(\Delta_{1680-1850}, \Delta_{1960-2008}) \\ 
& + 2 \text{Cov}(\Delta_{1850-1960}, \Delta_{1960-2008}).
\end{aligned} 
\quad (1) $$

The contribution of the temporal covariances to VΔ<sub>1680-2008</sub> due to linked selection thus becomes (G equation [1] of Buffalo and Coop (2020)):
$$G(\Delta_{1680-2008}) = \frac{2\mathrm{Cov}(\Delta_{1680-1850}, \Delta_{1850-1960}) + 2\mathrm{Cov}(\Delta_{1680-1850}, \Delta_{1960-2008}) + 2\mathrm{Cov}(\Delta_{1850-1960}, \Delta_{1960-2008})}{V\Delta_{1680-2008}} \quad (2)$$

G estimates the contribution of linked selection to the variance of the allelic frequency changes from the starting point to the most recent time period. Buffalo and Coop (2020) suggested that G can also be understood as the decrease in neutral diversity due to linked selection.

We extended this approach to the calculation of covariances between pairs of forests, (Cov $\Delta X_{i-j}$, $\Delta Y_{i-j}$) over the same period $(i-j)$ or at different nonoverlapping time periods. The maintenance of a positive covariance would indicate similar responses and directions of selection (parallel selection) in the different forests, whereas negative values would provide evidence of differences in the direction of selection. The extension to different forests made it possible to estimate covariance over contemporary time periods, or disconnected time scales (adjacent or distant), thereby making it possible to determine whether parallel selection occurred over the same time period, or whether fluctuating selection occurred over different time periods.

## GENE ANNOTATION AND ONTOLOGY ENRICHMENT ANALYSIS

We divided the oak genome into tiles of 100,000 bp in length and calculated temporal variances and covariances for all tiles. Genomic regions with strong responses to linked selection were identified on the basis of covariance values between time intervals above 0.01 for all three replicated forests.

The genes located in the regions with the highest covariances were annotated with various bioinformatic tools. *Arabidopsis thaliana* was initially used as a reference. Orthologous genes between oak and *A. thaliana* were identified for each protein sequence (file "Qrob_PM1N_CDS_aa_20161004.fa.gz," available at https://urgi.versailles.inra.fr/download/oak/) by selecting the best blast hit (blastP, pvalue < 1e+05) to represent the oak gene model, and *Arabidopsis* GO terms from TAIR were used for annotation (https://www.arabidopsis.org/). For oak genes without orthologs in *Arabidopsis*, we performed eggNOG-mapper version 2 functional annotation based on fast orthology assignments, using precomputed eggNOG version 5.0 clusters and phylogenies (Huerta-Cepas et al. 2017, 2019). Gene ontology (GO) terms were inferred from eggNOG orthologous groups (OGs).

Gene sets associated with Biological Processes (BP), Molecular Functions (MF), and Cellular Components (CC) were tested for GO terms enrichment with the “topGO” R package and the “weight” algorithm associated with a Fisher’s exact test to select the most relevant terms (Alexa et al. 2006; Alexa and

Rahnenfuhrer 2020). A P-value < 0.05 was applied for the statistical test, and no FDR was calculated because the P-values returned by the “weight” method are interpreted by Alexa and Rahnenfuhrer (2020) as corrected or not affected by multiple testing. In addition to the annotation and enrichment analysis, we implemented a more integrated approach by exploring the biological relationships and interactions between genes and biological processes. This step made use of Pathway Studio™ Plant that assembles knowledgebase coming mostly from model plant species but using here oak gene *Arabidopsis* homologs of oak genes as gene entries.

## Results

### DIVERSITY AND DIFFERENTIATION WITHIN AND BETWEEN AGE-STRUCTURED COHORTS

On average, more than 13 million SNPs were called in each cohort (Table 1). Within age-structured cohorts, nucleotide diversity was high ($\pi$ ~ 0.01205). Levels of diversity were similar across cohorts and forests and consistent with previous estimates for oak (Plomion et al. 2018). The mean pairwise $F_{ST}$ between cohorts ranged from 0.010 to 0.015, with no detectable structure between cohorts and forests (Fig. 2).

### TEMPORAL COVARIANCES OF ALLELIC FREQUENCY CHANGES WITHIN FORESTS

We calculated the covariances of allelic frequency changes between proximal time periods (Cov($\Delta_{1680-1850}$, $\Delta_{1850-1960}$) and Cov($\Delta_{1850-1960}$, $\Delta_{1960-2008}$)) and distant time periods (Cov($\Delta_{1680-1850}$, $\Delta_{1960-2008}$)). All proximal and distant covariances were significant in all three forests (Fig. 3). The covariances between the earliest adjacent time periods Cov($\Delta_{1680-1850}$, $\Delta_{1850-1960}$) were positive and significant in each forest (Fig. 3A). The covariances between the more recent adjacent time periods Cov($\Delta_{1850-1960}$, $\Delta_{1960-2008}$) varied between the three forests, being positive in Bercé, slightly lower in Réno-Valdieu, and negative in Tronçais (Fig. 3B). Finally, the covariances between distant time periods were much lower than those between adjacent time periods in Bercé and Réno-Valdieu, and also partially in Tronçais (Fig. 3C). Overall, the patterns illustrated in Figure 3 show a shift of covariances between adjacent time periods from the earliest periods considered (1680–1850, 1850–1960) to the more recent time periods (1850–1960, 1960–2008).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** *F*<sub>ST</sub> values between age-structured cohorts in the three forests (B: Bercé; R: Réno-Valdieu; T: Tronçais). Subscripts to forest acronyms indicate the ages of the cohorts: 4, age ~340 or year of birth ~1680; 3, age ~170 or year of birth ~1850; 2, age ~60 year of birth ~1960; 1, age ~12 or year of birth ~2008.</figcaption></figure>

We also estimated the contribution of the temporal covariances to $V\Delta_{1680-2008}$ by calculating $G\Delta_{1680-2008}$ (eq. 2; Table 2). The proportion of the variance of allelic frequency change due to linked selection increased from 1680 to 2008, due to the positive covariances in Bercé (9–16%) and Réno Valdieu (18–22%). It decreased in Tronçais from 12% to –1% (Table 2), due to the overall decrease in covariances between distant time periods. The largest contribution was that of the covariances for the earliest periods $Cov(\Delta_{1680-1850}, \Delta_{1850-1960})$ that were positive in all three forests, resulting in $G\Delta_{1680-1960}$ values of 9%, 18%, and 12%, respectively. We explored the genomic distribution of the temporal covariances between the two earliest adjacent periods by calculating the covariances for tiles of 100 Kb and plotting their distribution as Manhattan plots along chromosomes (Fig. 4). Although the overall pattern depicted a rather random distribution, there were a few islands of higher values (especially on chromosomes 2, 3, and 7; Figs. 4 and S5) that led us to investigate the gene composition within these tiles (paragraph GENE ONTOLOGY AND ENRICHMENT ANALYSIS). However, the overall distribution of mean covariances over the three forests across all tiles showed a clear bell-shaped curve shifted toward positive values, supporting a genome-wide contribution to the covariances (blue shaded distribution on Fig. S4). These positive covariances exhibit a genome-wide distribution across all chromosomes (Figs. 4 and S5).

## TEMPORAL COVARIANCES OF ALLELIC FREQUENCY CHANGES BETWEEN FORESTS

Significant positive temporal covariances between forests would indicate parallel linked selection. Considering different time periods for calculating the covariances can provide an indication as to when parallel selection occurred. Does it occur at contemporary or distant time periods? We therefore calculated the covariances of allelic frequency changes between the three forests for three different time periods (Fig. 5):

* Contemporary time periods (Fig. 5A);

* Adjacent time periods (Fig. 5B);

* Distant time periods (Fig. 5C).

We did not consider overlapping time periods to avoid the inflation of covariances due to allelic frequency changes common to the covariance terms.

Overall, covariances between forests remained positive and significant for contemporary time periods, particularly between Bercé and Réno-Valdieu and between Réno-Valdieu and Tronçais (Fig. 5A). Covariances between Bercé and Tronçais were close to zero regardless of the time period considered. Leaving these latter covariances aside, there was a trend over time toward a decrease in contemporary covariances from more ancient to more recent time periods (Fig. 5A). The covariances between two forests, between adjacent or distant time periods, were strikingly different from contemporary covariances, as they were all close to zero (Fig. 5B, C).

## GENE ONTOLOGY AND ENRICHMENT ANALYSIS

We used the genomic distribution of covariances at the 100 kb tile level (Manhattan plot; Fig. 4) to identify putative regions (tiles) under linked selection. We restricted our analysis to the two earliest time intervals (yielding positive values in all three forests (Fig. 3A), for which strong covariances were also observed between forests (Fig. 5). Tiles exhibiting covariances between the two earliest time intervals $Cov(\Delta_{1680-1850}, \Delta_{1850-1960})$ larger than 0.01 in all three replicates were considered as outliers (open bar distribution on Fig. S4). In total, we identified 104 tiles exhibiting temporal covariances above the threshold, which corresponded to 1% of the tiles (Fig. S4). The outlier tiles were distributed across all chromosomes as shown by the Manhattan plots (Fig. 4). Interestingly, outlier tiles that exhibited the highest covariances between the oldest time periods ($Cov(\Delta_{1680-1850}, \Delta_{1850-1960})$) comprise SNPs that show temporal patterns of differentiation. Indeed, $F_{ST}$ values of SNPs between proximal older time periods are larger than between proximal recent time periods, within and between forests, and they constantly decrease from older to recent time periods (Fig. 6). The $F_{ST}$ values of those SNPs exhibited temporal patterns that parallel the
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Temporal covariances of allelic frequency changes between different time periods and occurrences of extreme climatic events since the Little Ice Age. Mean and 95% confidence intervals of the covariances were obtained by bootstrapping with 5000 iterations. (A) Temporal covariances of allelic frequency changes between 1680–1850 and 1850–1960 in the three forests. (B) Temporal covariances of allelic frequency changes between 1850–1960 and 1960–2008 in the three forests. (C) Temporal covariances of allelic frequency changes between 1680–1850 and 1960–2008 in the three forests. (D) Timeline subdivided in decades. On the right side of the timeline in blue bars, number of extreme winters per decade according to instrumental temperatures recorded at the Observatory of Paris between 1676 and 2010 (Rousseau 2012) (more details in Fig. S2). On the left side of the timeline in red bars, number of extreme summer droughts per decade according to Cook's data base of Old World megadroughts (Palmer 1965; van der Schrier et al. 2013; Cook et al. 2015) (for more details, see Fig. S3). Highlighted decades in yellow correspond to periods when the cohorts became installed after natural regeneration.</figcaption></figure>
<table><caption>**Table 2.** Contribution of the covariances between allelic frequency changes to the variance of allelic frequency changes between two time points.</caption><thead><tr><th></th><th>Bercé</th><th>Réno-Valdieu</th><th>Tronçais</th></tr></thead><tbody><tr><td>G Δ<sub>1680-1960</sub></td><td>0.085</td><td>0.179</td><td>0.119</td></tr><tr><td>G Δ<sub>1680-2008</sub></td><td>0.163</td><td>0.221</td><td>-0.009</td></tr></tbody></table>

Note: G Δ<sub>1680-1960</sub>: contribution of temporal covariances to the variance of allelic frequency changes from 1680 to 1960. G Δ<sub>1680-1850</sub> equals 0 as there is no covariance for the first time period.

covariance temporal pattern (Fig. 3). In a few cases (e.g., on chromosomes 2, 3, and 7), we found outlier tiles exhibiting systematically higher values in the three forests (Fig. S5). However, the comparisons of the covariances of the outlier tiles with their $F_{ST}$ values (between cohorts delimiting the two earliest time periods) indicate that covariances build up as a result of repeated numerous subtle allelic frequency changes over two successive time periods, rather than by a few recurrent extreme shifts (albeit with a few exceptions in Réno-Valdieu; Fig. S5).

We inventoried 280 protein coding genes in these regions. Functional annotations of the 280 genes putatively under linked selection are summarized in Table S2. Briefly, out of these 280 genes, 248 received a GO annotation from *A. thaliana*, and two also received eggNOG GO terms, leaving 30 genes without a GO annotation.

The 280 genes revealed significant enrichments in the different gene ontologies (Tables S3–S5). Enrichment analysis identified several Biological Processes (BP), most of which are related to the “plant-type hypersensitive response,” “defense response to fungus,” “wax and cutin biosynthetic processes,” and “anther de-

hiscence,” with higher connectivity between the two first terms that gather 15 and 13, genes respectively (Fig. 7; Table S3).

Additional enrichment analysis accounting for other gene ontologies as Molecular Function or Cellular Component and gene network analysis using Pathway Studio™ Plant knowledge-base highlighted as well genes and functions contributing to resistance to biotic or abiotic factors (Supporting Information 3, Fig. S6, and Tables S3–S6).

## Discussion

### TEMPORAL COVARIANCES AS IMPRINTS OF GENOME-WIDE LINKED SELECTION DURING AND AFTER THE LIA

Temporal allelic frequency changes for genes underlying adaptive polygenic traits have been studied in detail through theoretical approaches (Hollinger et al. 2019; Stephan and John 2020), simulations (Franssen et al. 2017), and empirical case studies (Long et al. 2015). Convergent predictions highlight small and subtle changes for genes targeted by directional selection, and have been supported by multigeneration selection experiments in model organisms (Burke et al. 2010; Barghi et al. 2019; Michalak et al. 2019). Our results suggest that genomic signatures of recent selection are also detected in natural populations of long-lived nonmodel species and over limited numbers of generations. Rapid evolution was reported at the genomic level in previous studies in nonmodel species in the context of rapid environmental change (Franks et al. 2018; Hamann et al. 2021), strong selection pressures (Dayan et al. 2019), or domestication (Guo et al. 2018), and this process was fueled by the presence of large amounts of standing genetic variation (Bitter et al. 2019; Dayan et al. 2019).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Manhattan plots of the temporal covariances of allelic frequency changes between the two oldest time periods (Cov(Δ<sub>1680-1850, Δ1850-1960</sub>)) calculated at the tile level, over the whole genome. Black dots correspond to the mean covariance over the three forests. Outliers tiles (red bars) are tiles for which covariances are larger than 0.01 in each of the three forests.</figcaption></figure>
The results we obtained in oaks meet theoretical expectations and show that tree populations are under genome-wide linked selection. Buffalo and Coop (2019) showed that under directional selection (based on an exponential fitness function), the covariances of neutral allele frequency changes over time depend on the genetic variance of fitness (V<sub>A</sub>), recombination distance, and linkage disequilibrium. We recently reported that *Q. petraea* forest stands display substantial heritable differences in reproductive success (Alexandre et al. 2020). The genetic variance of fitness in a forest located between two of the forests studied here, Bercé and Réno-Valdieu, was 0.468 (Alexandre et al. 2020), which is at the upper end of the range of values reported in a recent review of the literature (89% of the reported values were below 0.20) (Hendry et al. 2018). Furthermore, QTL studies of fitness-related traits in oaks (growth, phenology, and reproduction) have shown that all these traits depend on a large number of loci (Brendel et al. 2008; Derory et al. 2010; Caignard et al. 2019), suggesting the genetic architecture of fitness is highly polygenic. Sessile oak and most likely trees in general are therefore prone to linked selection, provided that environmental constraints are strong enough for selection to operate.

## SHIFTS OF TEMPORAL COVARIANCES MIRROR CLIMATIC TRANSITIONS BETWEEN THE LIA AND MODERN TIMES

If *Q. petraea* populations provide the necessary substrate for linked selection, what were the environmental constraints that triggered viability selection during the time periods considered? There is evidence that the LIA was a period with frequent extreme events in both summer and winter (Dobrovolny et al. 2010; Berchtold et al. 2012; Moreno-Chamarro et al. 2017) (Figs. S2 and S3). As illustrated in Supporting Information 2, extreme frosts and droughts were more frequent up to the late 19th century than during more recent periods. Frost damage, including tree decline and deaths, was reported after extreme winters in the late LIA. Hausendorff (1940) summarized historical writings about forest districts in Northern Germany mentioning recurrent episodes of adult oak tree mortality during the 1740–1748 period, following the extreme winter of 1739–1740. Many historical documents from Catholic parishes reported extensive tree losses following the winter of 1708–1709, which is known to be the most severe winter in European records (Avila and Avila 1987; Luterbacher et al. 2004). In Eastern France, a historical atlas dating from 1758 reported “The dead material was mainly oak. For example, forty-six oak high forests belonging to all categories were totally dead.” This high mortality was presumably attributed to the extreme winter of 1708–1709 (Schnitzler 2020). However, the climate improved after 1850. There were only three reported ex-

treme winters (1942, 1947, and 1963) and three extreme droughts (1921, 1934, and 1976) during the 20th century (Figs. S2 and S3), during which oak decline was also recorded (Delatour 1983). These extreme events and the shift in their occurrence over time (during the second half of the 19th century) may be responsible for generating the pattern of linked selection observed in this study.

In even-aged managed oak stands, selection intensity at the seedlings stage is extremely high and the size of the seedling population very large. Theory (Buffalo and Coop 2019) and empirical data from selection experiments (Buffalo and Coop 2020) suggest that, if the direction of selection is maintained over time, covariances tend to be positive, whereas shifts of directional selection are likely to decrease covariances, possibly resulting in negative values. In experimental selection, in which the direction of selection is maintained over successive generations, temporal covariances are predicted to decay as a result of recombination (Buffalo and Coop 2019). In natural populations, heritable fitness differences, the strength and shifts of direction of selection, and recombination are the drivers of the temporal covariances. The widespread and frequent occurrence of extreme frosts and drought events up to the late 19th century in Europe may therefore be responsible for the positive covariance observed between the two earliest time periods (1680–1850 and 1850–1960; Fig. 3A) within each forest and for the positive covariances at contemporary time scales between forests (Fig. 5A). Conversely, the lower and even negative covariances between the two more recent time periods (1850–1960 and 1960–2008; Fig. 3B) may be a genomic signature of the change in the frequency of extreme events between the mid-19th century and today and of the shift of directional selection. These observed patterns of temporal covariances were predicted by simulations mimicking different strengths of selection showing that stronger directional selection generates stronger temporal covariances (Buffalo and Coop 2020). The maintenance of positive (or negative) covariances between adjacent time periods is reinforced by the limited impact of recombination in our study case. The short time interval between cohorts would have prevented opportunities for recombination events to break the linked polymorphisms. The time trends for covariances between forests support the interpretation that the strength and shifts of selection since the late LIA have shaped signatures of linked selection. Between-forest covariances were higher and more consistent across forests during the late LIA, when climatic conditions were harsher than the milder climatic conditions observed in the late 20th century (Fig. 5). Such trends suggest fluctuating selection as also supported by the temporal trends of pairwise F<sub>ST</sub> values of SNPs located in the outlier tiles (Fig. 6).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Temporal covariances of allelic frequency changes between the different forests for different time periods. Mean and 95% confidence intervals of the covariances were obtained by bootstrapping with 5000 iterations. Colors of the arrows on the left diagram indicate the time periods considered in the graphs. B: Bercé, R: Rêno-Valdieu, T: Tronçais. (A) Temporal covariances of allelic frequency changes between forests for contemporary time periods. (B) Temporal covariances of allelic frequency changes between forests for adjacent time periods. (C) Temporal covariances of allelic frequency changes between forests for distant time periods.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Matrix of *F*<sub>ST</sub> values of SNPs located in the outlier tiles between all pairs of age-structured cohorts.</figcaption></figure>

BIOTIC INTERACTIONS DRIVEN BY ADAPTIVE RESPONSE TO CLIMATE

The functional analysis highlighted significant enrichments for genes involved mainly in plant defense responses to pathogens, or contributing to abiotic stress responses (temperature and drought) (Fig. 7). Many of the genes encoding resistance proteins (R proteins) belong to the NBS-LRR families. These R-genes are widely represented in the oak genome, in which they appear as expanded groups (Plomion et al. 2018). Such NBS-LRR receptors detect effectors used by pathogens to facilitate infection, and participate in the signaling cascade leading to various responses preventing further infection, such as the hypersensitive response, the production of reactive oxygen species, and cell wall modification (Boyes et al. 1998; Roux 2010; El Kasmi et al. 2017). Interestingly two homologs of genes RPM1 and EFR were both involved in plant defense against *Pseudomonas syringae* (Boyes et al. 1998), which can infect a wide range of herbaceous and woody plants that have suffered frost and freezing damage (Young et al. 1988; Luisetti et al. 1991; Balestra et al. 2009). In addition, five genes encoding beta-glucosidases that could confer freezing/cold tolerance were identified (Thorlby et al. 2004; Fourrier et al. 2008; Ambroise et al. 2020). These genes may also participate in responses to other biotic/abiotic stresses (Baba et al. 2017; Vassao et al. 2018). Our functional analysis is con-

sistent with the oak responses described in cases of oak decline ultimately leading to the death of the tree. Botanic and pathological descriptions of oak decline during the LIA are lacking, but reports of recurrent sparse oak dieback in more recent decades can be used to retrace the steps leading to oak death following severe winter or drought events. In a review of oak decline in Europe, Thomas et al. (2002) highlighted the combined effects of climatic extremes (drought or frost) and defoliating insects and pathogenic fungi (Thomas et al. 2002). The starting point is an extreme climatic event that weakens the trees (Vanoni et al. 2016), which is then followed by pest and insect attacks, which ultimately kill the tree. For example, episodes of oak decline in the first half of the 20th century in various parts of Europe were caused by a combination of winter frost, summer drought, insect defoliation (caused by oak leaf roller [*Tortrix viridana*] and oak processionary moth [*Thaumetopoea processionea*], and several pathogenic fungi (mainly powdery mildew [*Erysiphe alphioides*] and root pathogens (*Armillaria* and *Phytophthora*) (Delatour 1983; Donaubauer 1998; Thomas et al. 2002). In summary, extreme events, such as severe frost or drought, expose trees not only to abiotic stresses, but also to biotic selection pressures. These, in turn, trigger resistance responses, which were identified at genomic level in the functional analysis. We suspect that these processes also operated during the late LIA and that exposure
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Network plot of the 10 most significant Biological Process (BP terms from gene ontology enrichment analysis). The size of the nodes is proportional to their degree of connectivity. The labels correspond to the name of the *Quercus robur* gene model followed by the locus name of the best *Arabidopsis* homolog and the corresponding *Arabidopsis* gene name (from TAIR: https://www.arabidopsis.org/) between brackets when available. When no gene name exists, a short description (from TAIR or EggNOG databases) is added. If no description is available, the description is set to "unknown."</figcaption></figure>

to selection was more stringent, leading more widespread death (Hausendorff 1940; Schnitzler 2020), and, ultimately, to parallel linked selection across the three forests.

CAVEATS AND LIMITATIONS

The temporal allelic frequency changes assessed in our particular case study of age-structured cohorts are also potentially subject to other sources of variation that must be considered. Although the theory of temporal covariances generated by linked selection has been developed for discrete generations in isolated populations (Buffalo and Coop 2019), it is unknown how covariances may be shaped in natural settings (Dehasque et al. 2020). We see

four (nonexclusive) mechanisms that may have contributed co-
variances in our study in addition to climate triggered selection:
generation overlap, gene flow between cohorts (within the stud-
ied forests), gene flow from external populations (not included in
our settings), and development-related changes of fitness. Some
of these mechanisms are not independent: gene flow between
cohorts reinforces generation overlap and hence their joint effect
should be considered. Gene flow may potentially generate posi-
tive covariances if the same differentiated populations pollinate
repeatedly the studied cohort populations. A realistic scenario
to consider in this case is, for example, older cohorts pollinating
systematically younger cohorts. This will create a so-called
age-driven unidirectional gene flow within each forest. The
occurrence of age-driven unidirectional gene flow is minimized
in our case because successful matings leading to the recruitment
include also pollination from younger trees to the older seed trees
in even-aged managed forests (Supporting Information S1). The
large census size of the age classes corresponding to the cohorts,
their random spatial distribution (Fig. 1), and the extent of pollen
dispersal within forests (Gerber et al. 2014) all reduce the likeli-
hood of age-driven unidirectional gene flow. Generation overlap
may also be enhanced by peculiar demographic scenarios where
a limited number of large trees (corresponding in our case to the
oldest cohorts) would repeatedly contribute to recruitments to
the next generations. In even-aged managed forests, recruitment
is implemented at the level of each compartment, using as seed
trees those located in the compartment, thus excluding repeated
recruitment from the same seed trees. Repeated unidirectional
gene flow from external populations can similarly lead to a
covariance signal. For example, consistent pollination from a
differentiated population through time could increase covari-
ance. However, highly differentiated populations are unlikely
to be present in the study area. Genetic surveys conducted with
genome-wide SNP data show that population differentiation is
extremely low (FST values usually lower than 2 %) throughout
western lowland Europe (Leroy et al. 2020). Our own data in this
study show similar results (Fig. 2): cohorts within forests (and
between forests) exhibit extremely low differentiation. The low
differentiation is due to limited drift effects and extensive pollen
flow as a result of the continuous distribution of *Q. petraea*
forests and their large population sizes (also shown by the large
π values in Table 1) in this part of the natural distribution. It is
thus unlikely that the same “dissimilar” population may have
pollinated repeatedly the studied cohorts resulting in positive
covariances.

Lastly, the confounding effects of development-related
changes in fitness may have contributed to increases in tempo-
ral covariances too. If traits contributing to fitness change over
time in a long-lived species, then genetic covariances of fitness-
related traits will increase between traits assessed at similar ages.
Such temporal serial autocorrelations and their decay over time
have been reported for growth in pines (Kremer 1992). Such
development-related covariances would be entirely confounded
with temporal covariances in our case. However, we would argue
that developmental covariances occurring at the adult stage are
likely to have a limited impact on temporal covariances, as se-
lection is overwhelmingly more severe at the juvenile stage, with
more than 90% of oak seedlings in a stand eliminated before the
age of 10 years (Jarret 2004).

Theoretical investigations that have tentatively addressed
the combined effects of gene flow, generation overlap, and
development-related variation of fitness are missing. Although

we cannot exclude such consequences, we think that they were
minimized in our experimental settings. Despite these caveats,
and although we acknowledge that the three forests are not evolu-
tionary independent, the similar patterns of temporal covariances
observed across the three forests (Figs. 3 and 5) are a strong sig-
nal of systematic driving forces operating in the same direction
in the three forests. Among the different potential causes that we
just discussed, we believe that selection is the most likely one.

ACKNOWLEDGMENTS

This research was supported by the European Research Council through
an Advanced Grant (project TREEPEACE # FP7-339728), by France
Génomique (project EVOL-OAK, ANR-10-INBS-09-08), and by the
French Forest Service (ONF) (INRAE-ONF TREEPEACE contract).
JC was financed by the National Natural Science Foundation of China
(31972946). We thank the staff of the ONF Research and Development
Department (M. Legay, L. Arnaudet) and the ONF staff at the National
Forests of Bercé, Réno-Valdieu, and Tronçais (A. Jeanneau, V. Breton,
B. Laurendeau, L. Nicolas, A. Durin) for technical assistance during
the sampling of cohorts and trees. We acknowledge the contribution
of N. Zimmermann (WSL, Birmensdorf, Switzerland) for downscaling
the historical climatic data. We thank L. Arnoux, S. Delzon, P. Reynet,
and F. Routier (UMR BIOGECO) for collecting bud, leaf, and cambium
samples for the 12 cohorts. We thank A. Ducoussou (UMR BIOGECO)
for assistance with forest selection, and B. Brachi and S. Gonzalez-
Martinez (UMR BIOGECO) for fruitful discussions during data analy-
sis. We thank the Genotoul Bioinformatics Platform (Bioinfo Genotoul)
at Toulouse (France) for providing computing and storage resources, and
the Genoscope (CEA, Commissariat à l'Energie Atomique et aux En-
ergies Alternatives) at Evry (France) for sequencing the oak samples.
We thank V. Buffalo for his helpful comments on an earlier version of
this manuscript. We thank three anonymous reviewers, the associate, and
editor-in-chief for their constructive remarks and suggestions for improv-
ing the manuscript.

AUTHOR CONTRIBUTIONS

Conception and coordination of the study: AK, CP, and ML. Sampling
of forests and cohorts: AK, FM, and LT. Collection of samples and DNA
extraction: LT, BD, CL, AK, TL, and DB. Dendrochronological analy-
sis: DB and FL. Whole genome sequencing: JMA, KL, and CP. Bioin-
formatic analysis: DS, IL, JC, and TL. Enrichment analysis: JCL. Data
analysis: DS, JC, AK, and ML. Writing of the manuscript: AK, DS, JC,
JCL, ML, and TL. All authors reviewed the manuscript.

DATA ARCHIVING

The data that support the findings of this study are openly available in the
public accessible data repository of INRAE: https://data.inrae.fr/dataset.
xhtml?persistentId=doi:10.15454/JLA1A0.

REFERENCES

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

Alexa, A., and J. Rahnenfuhrer. 2020. topGO: enrichment analysis for gene ontology. R package version 2.40.0. Available via https://bioconductor.org/packages/release/bioc/html/topGO.html.
Alexa, A., J. Rahnenfuhrer, and T. Lengauer. 2006. Improved scoring of func-
tional groups from gene expression data by decorrelating GO graph
structure. Bioinformatics 22:1600–1607.

Alexandre, H., L. Truffaut, E. Klein, A. Ducoussou, E. Chancerel, I. Lesur, et al. 2020. How does contemporary selection shape oak phenotypes? Evol. Appl. 13:2772–2790.

Altmann, A., P. Weber, D. Bader, M. Preuss, E. B. Binder, and B. Muller-Myhsok. 2012. A beginners guide to SNP calling from high-throughput DNA-sequencing data. Hum. Genet. 131:1541–1554.

Ambroise, V., S. Legay, G. Guerriero, J. F. Hausman, A. Cuypers, and K. Sergeant. 2020. The roots of plant frost hardiness and tolerance. Plant Cell Physiol. 61:3–20.

Anchukaitis, K. J., R. Wilson, K. R. Briffa, U. Buntgen, E. R. Cook, R. D'Arrigo et al. 2017. Last millennium Northern Hemisphere summer temperatures from tree rings: part II, spatially resolved reconstructions. Quat. Sci. Rev. 163:1–22.

Avila, F., and M. Avila. 1987. Le froid de Janvier 1985 et les grands hivers passés dans le Sud-Ouest de la France. Rev. Géogr. Pyrén. Sud Ouest. 58:5–22.

Baba, S. A., R. A. Vishwakarma, and N. Ashraf. 2017. Functional character-
ization of CsBGlu12, a beta-glucosidase from Crocus sativus, provides
insights into its role in abiotic stress through accumulation of antioxi-
dant flavonols. J. Biol. Chem. 292:4700–4713.

Balestra, G. M., L. Perestrelo, A. Mazzaglia, and A. Rossetti. 2009. First report of blossom blight caused by *Pseudomonas syringae* on kiwi fruit in Portugal. J. Plant Pathol. 91:231–231.

Barghi, N., R. Tobler, V. Nolte, A. M. Jaksic, F. Mallard, K. A. Otte et al. 2019. Genetic redundancy fuels polygenic adaptation in *Drosophila*. PLoS Biol. 17:e3000128.

Berchtold, J., E. Le Roy Ladurie, J.-P. Sermain, and A. Vasak. 2012. Canicules et froids extrêmes. Hermann, Paris.

Bitter, M. C., L. Kapsenberg, J. P. Gattuso, and C. A. Pfister. 2019. Stand-
ing genetic variation fuels rapid adaptation to ocean acidification. Nat.
Commun. 10:5821.

Boyes, D. C., J. Nam, and J. L. Dangl. 1998. The Arabidopsis thaliana RPM1 disease resistance gene product is a peripheral plasma membrane protein that is degraded coincident with the hypersensitive response. Proc. Natl. Acad. Sci. USA. 95:15849–15854.

Brendel, O., D. Le Thiec, C. Scotti-Saintagne, C. Bodenes, A. Kremer, and J. M. Guehl. 2008. Quantitative trait loci controlling water use efficiency and related traits in *Quercus robur* L. Tree Genet. Genomes. 4:263–278.

Buffalo, V., and G. Coop. 2019. The linked selection signature of rapid adap-
tation in temporal genomic data. Genetics. 213:1007–1045.

—. 2020. Estimating the genome-wide contribution of selection to tem-
poral allele frequency change. Proc. Natl. Acad. Sci. USA. 117:20672–
20680.

Burke, M. K., J. P. Dunham, P. Shahrestani, K. R. Thornton, M. R. Rose, and A. D. Long 2010. Genome-wide analysis of a long-term evolution experiment with *Drosophila*. Nature. 467:587–U111.

Caignard, T., S. Delzon, C. Bodenes, B. Dencausse, and A. Kremer. 2019. Heritability and genetic architecture of reproduction-related traits in a temperate oak species. Tree Genet. Genomes 15:1.

Camarero, J. J., J. M. Garcia-Ruiz, G. Sangusa-Barreda, J. D. Galvan, A. Q. Alla, Y. Sanjuan et al. 2015. Recent and intense dynamics in a formerly static Pyrenean treeline. Arct. Antart. Alp. Res. 47:773–783.

Campbell, I. D., and J. H. McAndrews. 1993. Forest disequilibrium caused by rapid Little Ice Age cooling. Nature. 366:336–338.

Carrer, M., and C. Urbinati. 2006. Long-term change in the sensitivity of tree-
ring growth to climate forcing in *Larix decidua*. New Phytol. 170:861–
871.

Chown, S. L., K. A. Hodgins, P. C. Griffin, J. G. Oakeshott, M. Byrne, and A. A. Hoffmann. 2015. Biological invasions, climate change and ge-nomics. Evol. Appl. 8:23–46.

Cook, E. R., R. Seager, Y. Kushnir, K. R. Briffa, U. Buntgen, D. Frank et al. 2015. Old World megadroughts and pluvials during the Common Era. Sci. Adv. 1:e1500561.

Corona, C., J. Guiot, J. L. Edouard, F. Chalie, U. Buntgen, P. Nola et al. 2010. Millennium-long summer temperature variations in the European Alps as reconstructed from tree rings. Clim. Past 6:379–400.

Dayan, D. I., X. Du, T. Z. Baris, D. N. Wagner, D. L. Crawford, and M. F. Oleksiak. 2019. Population genomics of rapid evolution in natural populations: polygenic selection in response to power station thermal effluents. BMC Evol. Biol. 19:61.

Dehasque, M., M. C. Avila-Arcos, D. Diez-del-Molino, M. Fumagalli, K. Guschanski, E. D. Lorenzen et al. 2020. Inference of natural selection from ancient DNA. Evol. Lett. 4:94–108.

Delatour, C. 1983. Les dépérissements de chênes en Europe. Rev. Forest. Fr. 35:265–282.

Derory, J., C. Scotti-Saintagne, E. Bertocchi, L. Le Dantec, N. Graignic, A. Jauffres et al. 2010. Contrasting relations between diversity of candidate genes and variation of bud burst in natural and segregating populations of European oaks. Heredity. 105:401–411.

Dobrovolny, P., A. Moberg, R. Brazdil, C. Pfister, R. Glaser, R. Wilson et al. 2010. Monthly, seasonal and annual temperature reconstructions for Central Europe derived from documentary evidence and instrumental records since AD 1500. Clim. Change. 101:69–107.

Donaubauer, E. 1998. Die Bedeutung von Krankheitserregern beim gegen-
wärtigen Eichensterben in Europa-eine Literaturübersicht. Eur. J. Forest
Pathol. 28:91–98.

Edouard, J.-L., C. Corona, A. Thomas, F. Guibal, and N. Denelle. 2009. Le petit âge glaciaire dans les cernes de croissance des arbres des Alpes françaises. Approche dendrochronoligical. Archeol. Midi Médiév. 27:169–177.

El Kasmi, F., E. H. Chung, R. G. Anderson, J. Y. Li, L. Wan, T. K. Eitas et al. 2017. Signaling from the plasma-membrane localized plant immune re-ceptor RPM1 requires self-association of the full-length protein. Proc. Natl. Acad. Sci. USA. 114:E7385–E7394.

Fagan, B. 2002. The little ice age. How climate made history (1300-1850). Basic Books, New York.

Fourrier, N., J. Bedard, E. Lopez-Juez, A. Barbrook, J. Bowyer, P. Jarvis et al. 2008. A role for sensitive to freezing2 in protecting chloroplasts against freeze-induced damage in Arabidopsis. Plant J. 55:734–745.

Franks, S. J., E. Hamann, and A. E. Weis. 2018. Using the resurrection approach to understand contemporary evolution in changing environments. Evol. Appl. 11:17–28.

Franssen, S. U., R. Kofler, and C. Schlotterer. 2017. Uncovering the genetic signature of quantitative trait evolution with replicated time series data. Heredity 118:42–51.

Gallon, D. 1752. Conférence de l'ordonnance de Louis XIV du mois d'Aout 1969 sur le fait des Eaux et Forêts. Brunet, Paris.

Gerber, S., J. Chadœuf, F. Gugerli, M. Lascoux, J. Buiteveld, J. Cottrell et al. 2014. High rates of gene flow by pollen and seed in oak populations across Europe. PLoS ONE 9:e85130

Giesecke, T. 2016. Did thermophilous trees spread into central Europe during the Late Glacial? New Phytol. 212:15–18.

Giesecke, T., and S. Brewer. 2018. Notes on the postglacial spread of abundant European tree taxa. Veg. Hist. Archaeobot. 27:337–349.

Guo, J. Z., H. X. Tao, P. F. Li, L. Li, T. Zhong, L. J. Wang et al. 2018. Whole-genome sequencing reveals selection signatures associated with important traits in six goat breeds. Sci. Rep. 8:10405.
Hamann, E., C. S. Pauli, Z. Joly-Lopez, S. C. Groen, J. S. Rest, N. C. Kane et al. 2021. Rapid evolutionary changes in gene expression in response to climate fluctuations. Mol. Ecol. 30:193–206.

Hausendorff, E. 1940. Frostschäden an Eichen. Aus der Geschichte des Forstamtes Grimmitz/Uckermarck. Z. Forst Jagdw. 72:3–35.

Helama, S., M. Kuoppamaa, and R. Sutinen 2020. Subaerially preserved remains of pine stemwood as indicators of late Holocene timberline fluctuations in Fennoscandia, with comparisons of tree-ring and 14C dated depositional histories of subfossil trees from dry and wet sites. Rev. Palaeobot. Palynol. 278:104223

Hendry, A. P., D. J. Schoen, M. E. Wolak, and J. M. Reid 2018. The contemporary evolution of fitness. Annu. Rev. Ecol. Evol. Syst. 49:457–476.

Hollinger, I., P. S. Pennings, and J. Hermisson 2019. Polygenic adaptation: from sweeps to subtle frequency shifts. PLoS Genet. 15:e1008035.

Huerta-Cepas, J., K. Forslund, L. P. Coelho, D. Szklarczyk, L. J. Jensen, C. von Mering et al. 2017. Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol. Biol. Evol. 34:2115–2122.

Huerta-Cepas, J., D. Szklarczyk, D. Heller, A. Hernandez-Plaza, S. K. Forslund, H. Cook et al. 2019. eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 47:D309–D314.

Jarret, P. 2004. Chênaie atlantique. Lavoisier, Paris.

Kofler, R., P. Orozco-terWengel, N. De Maio, R. V. Pandey, V. Nolte, A. Futschik et al. 2011. PoPoolation: a toolbox for population genetic analysis of next generation sequencing data from pooled individuals. PLoS One. 6:e15925.

Kremer, A. 1992. Predictions of age-age correlations of total height based on serial correlations between increments in Maritime pine (Pinus pinaster Ait.). Theor. Appl. Genet. 85:152–158.

Kremer, A., and A. L. Hipp 2020. Oaks: an evolutionary success story. New Phytol. 226:987–2011.

Kullman, L. 2005. Pine (Pinus sylvestris) treeline dynamics during the past millennium - a population study in west-central Sweden. Ann. Bot. Fenn. 42:95–106.

———. 2015. Higher than present Medieval pine (Pinus sylvestris) treeline along the Swedish Scandes. Landsc. Online 42:1–14.

Le Roy Ladurie, E. 2004. Histoire humaine et comparée du climat. Canicules et glaciers XIIIième-XVIIIième siècles. Fayard, Paris.

———. 2006. Histoire humaine et comparée du climat. Disettes et révolutions 1740–1860. Fayard, Paris.

Leroy, T., J. M. Louvet, C. Lalanne, G. Le Provost, K. Labadie, J. M. Aury et al. 2020. Adaptive introgression as a driver of local adaptation to climate in European white oaks. New Phytol. 226:1171–1182.

Li, H., and R. Durbin. 2010. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 26:589–595.

Li, H., B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer et al. 2009. The sequence alignment/map format and SAMtools. Bioinformatics. 25:2078–2079.

Linderholm, H. W., P. Zhang, B. E., J., B. Gunnarson, E. Farahat, M. Fuentes et al. 2014. Growth dynamics of tree-line and lake-shore Scots Pine (Pinus sylvestris L.) in the Central Scandinavian Mountains during the Medieval Climate anomaly and the early Little Ice Age. Front. Ecol. Evol. 2:20.

Long, A., G. Liti, A. Luptak, and O. Tenaillon. 2015. Elucidating the molecular architecture of adaptation via evolve and resequence experiments. Nat. Rev. Genet. 16:567–582.

Luisetti, J., J. L. Gaignard, and M. Devaux. 1991. Pseudomonas syringae pv syringae as one of the factors affecting the ice nucleation of grapevine buds in controlled conditions. J. Phytopathol. 133:334–344.

Luterbacher, J., D. Dietrich, E. Xoplaki, M. Grosjean, and H. Wanner. 2004. European seasonal and annual temperature variability, trends, and extremes since 1500. Science 303:1499–1503.

Luterbacher, J., J. P. Werner, J. E. Smerdon, L. Fernandez-Donado, F. J. Gonzalez-Rouco, D. Barriopedro et al. 2016. European summer temperatures since Roman times. Environ. Res. Lett. 11:024001.

MacDonald, G. M., K. V. Kremenetski, and D. W. Beilman 2008. Climate change and the northern Russian treeline zone. Philos. Trans. R. Soc. B Biol. Sci. 363:2285–2299.

Malaspinas, A. S. 2016. Methods to characterize selective sweeps using time serial samples: an ancient DNA perspective. Mol. Ecol. 25:24–41.

Mariette, S., J. Cottrell, U. M. Csaikl, P. Goikoechea, A. Konig, A. J. Lowe et al. 2002. Comparison of levels of genetic diversity detected with AFLP and microsatellite markers within and among mixed Q-petraea (MATT.) LIEBL. and Q-robur L. stands. Silvae Genet. 51:72–79.

Michalak, P., L. Kang, M. F. Schou, H. R. Garner, and V. Loeschcke. 2019. Genomic signatures of experimental adaptive radiation in *Drosophila*. Mol. Ecol. 28: 600–614.

Moreno-Chamarro, E., D. Zanchettin, K. Lohmann, J. Luterbacher, and J. H. Jungclaus. 2017. Winter amplification of the European Little Ice Age cooling by the subpolar gyre. Sci. Rep. 7:9981.

Palmer, W. C. 1965. Meteorological drought. US Department of Commerce, Weather Bureau, Washington, D.C.

Parker, G. 2013. Global crisis. War, climate change and catastrophe in the seventeenth century. Yale Univ. Press, New Haven, CT.

Petit, R. J., U. M. Csaikl, S. Bordacs, K. Burg, E. Coart, J. Cottrell et al. 2002. Chloroplast DNA variation in European white oaks - phylogeny and patterns of diversity based on data from over 2600 populations. Forest Ecol. Manag. 156:5–26.

Pfeifer, S. P. 2017. From next-generation resequencing reads to a high-quality variant data set. Heredity 118:111–124.

Pfister, C. 1984. *Das Klima der Schweiz von 1525–1860 und seine Bedeutung in der Geschichte von Bevölkerung und Landwirtschaft. Volume 2. Bevölkerung, Klima und Agrarmodernisierung 1525–1860*. Paul Haupt, Bern, Switzerland.

Plomion, C., J. M. Aury, J. Amselem, T. Alaeitabar, V. Barbe, C. Belser et al. 2016. Decoding the oak genome: public release of sequence data, assembly, annotation and publication strategies. Mol. Ecol. Resour. 16:254–265.

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

Ravenscroft, C. H., R. Whitlock, and J. D. Fridley 2015. Rapid genetic divergence in response to 15 years of simulated climate change. Global Change Biol. 21:4165–4176.

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

Rousseau, D. 2012. Identification des grands hivers de 1676 à 2010 à l'aide de séries thermométriques de Paris. Pp. 345–360 in J. Berchtold, E. Le Roy Ladurie, J.-P. Sermain, and A. Vasak, eds. Canicules and froids extrêmes. Hermann, Paris.

Roux, M. 2010. Identification and characterization of EFR interacting proteins. The Univ. of East Anglia, Norwich, U.K.

Schlotterer, C., R. Kofler, E. Versace, R. Tobler, and S. U. Franssen. 2015. Combining experimental evolution with next-generation sequencing: a powerful tool to study adaptation from standing genetic variation. Heredity. 114:431–440.
Schnitzler, A. 2020. Spontaneous landscape dynamics in the Pays de Bitche, Lorraine (France), during the Little Ice Age. Landsc. History. 41:89–104.

Stephan, W., and S. John. 2020. Polygenic adaptation in a population of finite size. Entropy. 22:907.

Tajima, F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 123:585-595.

Thomas, F. M., R. Blank, and G. Hartmann. 2002. Abiotic and biotic factors and their interactions as causes of oak decline in Central Europe. Forest Pathol. 32:277-307.

Thorlby, G., N. Fourrier, and G. Warren. 2004. The sensitive to freezing2 gene, required for freezing tolerance in *Arabidopsis thaliana*, encodes a beta-Glucosidase. Plant Cell. 16:2192-2203.

Tkachuck, R. D. 1983. The Little Ice Age. Origins 10:51-65.

van der Schrier, G., J. Barichivich, K. R. Briffa, and P. D. Jones. 2013. A scPDSI-based global data set of dry and wet spells for 1901–2009. J. Geophys. Res. Atmos. 118:4025–4048.

Vanoni, M., H. Bugmann, M. Notzli, and C. Bigler. 2016. Drought and frost contribute to abrupt growth decreases before tree mortality in nine temperate tree species. Forest Ecol. Manag. 382:51–63.

Vassao, D. G., N. Wielsch, A. Gomes, S. Gebauer-Jung, Y. Hupfer, A. Svatos et al. 2018. Plant defensive beta-glucosidases resist digestion and sustain activity in the gut of a lepidopteran herbivore. Front. Plant Sci. 9:1389.

Young, J. M., G. J. Cheesmur, F. V. Welham, and W. R. Henshall 1988. Bacterial blight of kiwi fruit. Ann. Appl. Biol. 112:91-105.

Associate Editor: Z. Gompert

## Supporting Information

Additional supporting information may be found online in the Supporting Information section at the end of the article.

**Figure S1.** Mean yearly temperature trends in the three studied forests.

**Figure S2.** Number of extreme winters per decade.

**Figure S3.** Number of extreme summer droughts per decade near each studied forest according to Cook's data base of Old World megadroughts (Cook et al., 2015).

**Figure S4.** Distribution of the temporal covariances between allelic frequency changes between the two oldest time periods (Cov(Δ₁₆₈₀-₁₈₅₀, Δ₁₈₅₀-₁₉₆₀)), calculated for tiles of 100kb.

**Figure S5.** Manhattan plot of the temporal covariances of allelic frequency changes between the two earliest time periods (Cov(Δ₁₆₈₀-₁₈₅₀, Δ₁₈₅₀-₁₉₆₀)) and of Fst between the cohorts delimiting the time periods (Fst₁₆₈₀-₁₈₅₀ and Fst₁₈₅₀-₁₉₆₀).

**Figure S6.** Biological network showing the 15 selected Cell Processes and their 74 connected *Arabidopsis* genes (according to the pathway Studio™ Plant database).

**Table S1.** Geographic coordinates of sampled forests and age structured cohorts

**Table S2.** Functional annotation of the 280 oak genes located in outlier tiles

**Table S3.** Significance of Biological Process GO terms according to topGO "weight" method and Fisher's exact test.

**Table S4.** Significance of Molecular Function GO terms according to topGO "weight" method and Fisher's exact test.

**Table S5.** Significance of Cellular Component GO terms according to topGO "weight" method and Fisher's exact test.

**Table S6.** List of the 74 *Arabidopsis* homologs known to participate to at least one cell process with their name used in Pathway Studio™ Plant, their locus name as defined in TAIR (https://www.arabidopsis.org) and alias commonly used in scientific publications
Supporting Information S1, S2, S3