1. Genetics and Genomics
  2. Plant Biology
Download icon

Genomic basis for drought resistance in European beech forests threatened by climate change

  1. Markus Pfenninger  Is a corresponding author
  2. Friederike Reuss
  3. Angelika Kiebler
  4. Philipp Schönnenbeck
  5. Cosima Caliendo
  6. Susanne Gerber
  7. Berardino Cocchiararo
  8. Sabrina Reuter
  9. Nico Blüthgen
  10. Karsten Mody
  11. Bagdevi Mishra
  12. Miklós Bálint
  13. Marco Thines
  14. Barbara Feldmeyer
  1. Molecular Ecology, Senckenberg Biodiversity and Climate Research Centre, Germany
  2. Institute for Organismic and Molecular Evolution, Johannes Gutenberg University, Germany
  3. LOEWE Centre for Translational Biodiversity Genomics, Germany
  4. Institute of Human Genetics, University Medical Center, Johannes Gutenberg University, Germany
  5. Conservation Genetics Section, Senckenberg Research Institute and Natural History Museum Frankfurt, Germany
  6. Ecological Networks lab, Department of Biology, Technische Universität Darmstadt, Germany
  7. Department of Applied Ecology, Hochschule Geisenheim University, Germany
  8. Biological Archives, Senckenberg Biodiversity and Climate Research Centre, Germany
  9. Functional Environmental Genomics, Senckenberg Biodiversity and Climate Research Centre, Germany
  10. Agricultural Sciences, Nutritional Sciences, and Environmental Management, Universität Giessen, Germany
  11. Institute for Ecology, Evolution and Diversity, Johann Wolfgang Goethe-University, Germany
Research Article
  • Cited 0
  • Views 1,265
  • Annotations
Cite this article as: eLife 2021;10:e65532 doi: 10.7554/eLife.65532

Abstract

In the course of global climate change, Central Europe is experiencing more frequent and prolonged periods of drought. The drought years 2018 and 2019 affected European beeches (Fagus sylvatica L.) differently: even in the same stand, drought-damaged trees neighboured healthy trees, suggesting that the genotype rather than the environment was responsible for this conspicuous pattern. We used this natural experiment to study the genomic basis of drought resistance with Pool-GWAS. Contrasting the extreme phenotypes identified 106 significantly associated single-nucleotide polymorphisms (SNPs) throughout the genome. Most annotated genes with associated SNPs (>70%) were previously implicated in the drought reaction of plants. Non-synonymous substitutions led either to a functional amino acid exchange or premature termination. An SNP assay with 70 loci allowed predicting drought phenotype in 98.6% of a validation sample of 92 trees. Drought resistance in European beech is a moderately polygenic trait that should respond well to natural selection, selective management, and breeding.

eLife digest

Climate change is having a serious impact on many ecosystems. In the summer of 2018 and 2019, around two thirds of European beech trees were damaged or killed by extreme drought. It is critical to keep these beech woods healthy, as they are central to the survival of over 6,000 other species of animals and plants.

The level of damage caused by the drought varied between forests. However, not all the trees in each forest responded in the same way, with severely damaged trees often sitting next to fully healthy ones. This suggests that the genetic make-up of each tree determines how well it can adapt to drought rather than its local environment.

To investigate this further, Pfenninger et al. studied the genome of over 400 European beech trees from the Hesse region in Germany. The samples came from pairs of neighbouring trees that had responded differently to the droughts. The analysis found more than 80 parts of the genome that differed between healthy and damaged trees.

Pfenninger et al. then used this information to create a genetic test which can quickly and inexpensively predict how well an individual beech tree might survive in a drought. Applying this test to another 92 trees revealed that it can reliably detect which ones were healthy and which ones were damaged.

Beech forests are typically managed by private owners, agencies or breeders that could use this genetic test to select and reproduce trees that are better adapted to drought. The goal now is to develop the test so that it can be used more widely to manage European beech trees and potentially other species.

Introduction

Climate change comes in many different facets, amongst which are prolonged drought periods (Christensen et al., 2007). The Central European droughts in the years 2018 and 2019 caused severe water stress in many forest tree species, leading to the die-off of many trees (Schuldt et al., 2020). Among the suffering tree species was European beech, Fagus sylvatica L. As one of the most common deciduous tree species in Central Europe, F. sylvatica is of great ecological importance: beech forests are a habitat for more than 6000 different animal and plant species (Brunet et al., 2010; Dorow et al., 2010). The forestry use of beech in 2017 generated a turnover of more than €1 billion in Germany alone (Thünen_Institute, 2020), without taking the economic and societal value of the ecosystem services of woods into account (Elsasser et al., 2016). However, the drought years 2018 and 2019 severely impacted the beech trees in Germany (Paar and Dammann, 2019). Official reports on drought damage in beech recorded 62% of trees with rolled leaves and 20–30% of small leaves, mainly in the crown, resulting in 7% of badly damaged or dead trees. As shown before (Bressem, 2008), mainly medium- to old-aged trees were affected by drought stress (>60 years).

Under favourable conditions, beech is a competitive and shade-tolerant tree species, dominating mixed stands (Pretzsch et al., 2013). High genetic diversity within populations supports adaptation to local conditions (Kreyling et al., 2012). Significant differences between local populations in tolerance to various stress factors such as early frost (Czajkowski and Bolte, 2006), drought (Cocozza et al., 2016; Harter et al., 2015), or air pollution (Müller-Stark, 1985) are known. The distribution of F. sylvatica is mainly limited by water availability as the tree does not tolerate particularly wet or dry conditions (Sutmöller et al., 2008). Therefore, it is quite conceivable that the species could suffer even more under the predicted future climatic conditions than today (Sutmöller et al., 2008).

Despite the widespread, severe drought damage, a pattern observed in all beech forests was very noticeable (personal observations). Using crown deterioration as a significant indicator for drought damage (Choat et al., 2018), not all trees in a beech stand were equally damaged or healthy. The damage occurred rather in a mosaic-like pattern instead. Even though the extent of drought damage varied among sites, apparently completely healthy trees immediately neighboured severely damaged ones and vice versa. This observation gave rise to the hypothesis that not the local environmental conditions might be decisive for the observed drought damage, but rather the genetic make-up of the individual trees.

We decided to draw on this natural ‘experimental set-up’ to infer the genomic basis underlying the drought susceptibility in F. sylvatica. We identified more than 200 neighbouring pairs of trees with extreme phenotypes and used a Pool-GWAS approach (Bastide et al., 2013) to infer associated single-nucleotide polymorphism (SNP) loci by contrasting allele frequencies with replicated pools of drought-susceptible and -resistant individuals. In addition, we individually resequenced a subset of 51 pairs of susceptible and resistant trees. If the observed pattern indeed has a genetic basis, identifying the associated loci would enable the genomic prediction of drought resistance (Stocks et al., 2019). Constructing an SNP assay from the most highly phenotype-associated SNPs, we validated 70 identified loci by predicting the drought phenotype of an additional set of beech trees from their genotype at these loci using linear discriminant analysis (LDA) and a new Machine Learning approach (Horenko, 2020). These accurate genomic prediction tools, for example, the choice of drought-resistant seed-producing trees and selective logging, could help accelerate and monitor natural selection and thus harness beech forests against climate change (Waldvogel et al., 2020).

Results

Sampling, climate development, and phenotyping

Damaged and healthy beech tree pairs were sampled from woods in the lowland Rhein-Main plain, the adjacent low mountain ranges of Odenwald and Taunus, and mountain ranges from Central and Northern Hessen (Figure 1A). When summarising the climatic conditions from 1950 to 2019 for the sampling sites in a principal component analysis (PCA), the sites were divided into two groups by axis 1, a temperature gradient. The Taunus mountain sites grouped with those from the northern part of Hessen, while the Rhein-Main plain clustered with the Odenwald sites (Figure 1B). This grouping was also used to construct the GWAS pools (see below). Comparing the climate from the 1950s, when most of the trees sampled were already in place, with the decade from 2010 to 2019, showed that all local conditions changed substantially and similarly in the direction and extent of warmer and drier conditions (Figure 1B). The steepest temperature increase occurred in the 1980s, while precipitation patterns mainly changed in the last decade (Figure 1—figure supplement 1). A wide range of parameters, potentially relevant as selection pressures, changed drastically during this period: the mean January daily minimum temperature at the sampling sites increased by 1.49°C from −2.64°C (s.d. 1.68°C) in the 1950s to −1.15°C (s.d. 2.50°C) during the last decade. The mean August daily maximum temperatures increased even more by 2.37°C from 22.06°C (s.d. 1.95°C) to 24.43°C (s.d. 2.35°C). Simultaneously, mean annual precipitation decreased by 40.5 mm or 5.5% from 741.2 mm (s.d. 85.8 mm) to 700.7 mm (s.d. 70.9 mm). Most of the precipitation loss (84%) occurred during the main growth period between April and September, with a decrease of 33.9 mm from 410.4 mm (s.d. 36.1 mm) to 376.5 mm (s.d. 25.6 mm).

Figure 1 with 2 supplements see all
Location of sampling sites and climate history.

(A) Locations of sampling sites in Hessen, Germany. For abbreviations, see Supplementary file 1A. The sites where confirmation individuals were sampled are designated in grey. (B) Principal component analysis of monthly climate data 1950–2019. (C) Development of main growth period drought indicator from 1991 to 2019. Shown is the difference mean monthly evaporation potential in mm from April to September relative to the 1991 level. Climate and drought data obtained from https://opendata.dwd.de/climate_environment/CDC/grids_germany/monthly/.

Mean monthly evaporation potential, available from 1991 onwards, showed that, compared to the beginning of the 1990s, the main growth period of beech from April to September became increasingly drier, with up to 30 mm more evaporation per month. The drought dynamics suggested that the years 2018 and 2019 were not outliers, but rather part of a long-term, accelerating trend (Figure 1C), following the overall global pattern (Büntgen et al., 2021; Trenberth et al., 2014).

There was a strong negative correlation (r = 0.695) between the drought strength during the main growth period (Apr– Sept) and a proxy for (green) leaf cover (leaf area index [LAI]) for the sampled plots in the years 2015–2019 (Figure 1—figure supplement 2). This observation suggested that leaf loss and dried leaves are good indicators for drought stress.

The mean distance between paired trees was 5.1 m (s.d. 3.4 m, Figure 2—figure supplement 1). Phenotypic measurements generally confirmed the study design and selection of trees: healthy and damaged trees within each tree pair did not differ significantly in trunk circumference, tree height, canopy closure, and competition index (Figure 2A–D, Supplementary file 1B). Hence, these parameters were not considered in further analyses. As expected, and confirming the assignment of damage status, the quantity of dried leaves and leaf loss differed substantially between damaged and healthy ones (Figure 2E, F, Supplementary file 1B). A sample of photographs contrasting damaged and healthy paired trees can be found in Figure 2—figure supplement 2.

Figure 2 with 2 supplements see all
Comparison of sampled beech pairs.

(A) Trunk circumference, (B) canopy closure, (C) tree height, (D) competition index, (E) dried leaves, and (F) leaf loss. Boxplots with indicated means, the boxes represent one standard deviation, the whiskers are the 95% confidence intervals. Damaged trees in ochre, healthy trees in green. Except for (E) and (F), the difference of means among damaged and healthy trees is insignificant between the groups.

Linkage disequilibrium, population structure, and genome-wide association study

For a subsample of 300 out of the 402 sampled beech trees, we generated four DNA pools from two climatically distinct regions (North and South Hessen, Figure 1B), contrasting trees that were either healthy or highly drought damaged, respectively (Supplementary file 1A). The ‘South’ pools consisted of 100 individuals each, whereas the ‘North’ pools contained 50 individuals each. We created ~50 GB 150 bp-paired end reads with insert size 250–300 bp on an Illumina HiSeq 4000 system per pool. More than 96% of the reads mapped against the repeat-masked chromosome-level beech reference genome (accession no. PRJNA450822). After filtering the alignment for quality and a coverage between 15× and 70×, and removing indels, allele frequencies for 9.6 million SNPs were scored. All 100 individuals from the North population were additionally individually resequenced to ~20× coverage each (for more details, see Materials and methods). This data was used to (i) determine individual variability in allele frequencies and (ii) validate the information content of the candidate SNP set.

Using all individually resequenced individuals, we inferred the extent of genome-wide linkage disequilibrium (LD). The plot of LD r² against the distance from the focal SNP showed that LD fell to r² ~ 0.3 within <120 bp, which means that genome positions such a distance apart are on average effectively unlinked (Figure 3—figure supplement 1A). The PCA on SNP variation of the individually resequenced trees from the North population explained 12.3% of accumulated variation on the first two axes (Figure 3—figure supplement 1B). Trees from the same sampling site (within the North population) did not tend to cluster together (Figure 3—figure supplement 1B). FST estimates among pools for non-overlapping 1 kb windows were virtually identical among healthy/damaged pools within region as compared to between regions (Figure 3—figure supplement 2). Trees within a phenotypic class were genomically not more similar than between classes (Figure 3—figure supplement 3, ANOSIM R = −0.008, p=0.76, 9999 permutations).

Pool-GWAS analysis identified 106 SNPs significantly associated with the drought damage status using a Cochran–Mantel–Haenszel test on the two pairs of damaged and healthy pools after false discovery rate correction and a cutoff at 1 × 10−2 (Figure 3A, Figure 3—figure supplement 4). Some of the 106 SNPs were in close physical proximity (<120 bp) and thus probably linked. Taking this into account, 80 independent genomic regions were associated with the drought damage status. None of the significantly differentiated SNP loci was mutually fixed; the observed allele frequency differences between healthy and damaged trees at associated loci ranged between 0.12 and 0.51 (Figure 3B).

Figure 3 with 7 supplements see all
Single-nucleotide polymorphisms (SNPs) significantly associated with damaged phenotypes and combined results of SNP assay and discriminant analyses.

(A) Manhattan plot of false discovery rate (FDR)-corrected –log10 probability values from Cochran–Mantel–Haenszel test. The black horizontal line indicates the chosen significance threshold. SNPs on different chromosomes alternate in colour (black and blue). (B) Mean allele frequency difference at significantly associated SNP loci between healthy and damaged phenotypes. The loci are ordered according to amount of change. (C) The centre of the figure depicts the genotypogram of the SNP assay. Each column represents one of 70 loci, each row one of 92 beech individuals. The scored genotypes are colour-coded, with red squares = homozygous reference allele, light blue = homozygous alternate allele, white = heterozygous SNP, and grey squares = locus could not be scored in the respective individual. The left bar indicates the observed phenotype for each tree individual with ochre rectangles for damaged and green for healthy trees. Below the genotypogram, the relative contribution of each locus to the predictive model of the discriminant analysis is indicated, ordered from high to low. On the right side, first the genotype model scores for each individual are given, with the according predicted phenotype (ochre = damaged; green = healthy).

Associated genes and gene function

Of the 106 significant SNPs, 24 were found in 20 protein coding genes (Table 1). Forty-nine genes were the closest genes to the remaining 82 SNPs. For 61 of these genes, the best BLAST hit was with a tree, mainly from the Fagales genera Quercus and Castanea (Table 1, Supplementary file 1C). Among the 24 SNPs in genes, we observed 13 non-synonymous changes. In 11 of these changes, the alternate allele was associated with the damaged phenotype and only in two cases with the healthy phenotype. Three of the non-synonymous substitutions resulted in a stop codon. Of the remaining 10, 8 exchanges caused a major change in amino acid characteristics and thus probably in protein folding or function (Table 1). One gene, a PB1 domain-containing protein tyrosine kinase, contained four non-synonymous changes, suggesting that the allele version associated with the damaged phenotype lost its function (Table 1). From the 20 genes with significant SNPs, functional information could be obtained from the UniProt database for 14 (Supplementary file 1C). Of these, 10 genes were associated in previous studies with either environmental stress response (two) or specifically with drought stress response (eight; Supplementary file 1C). Of the 49 predicted genes closest to the remaining significant SNPs (Table 1), 16 could be reliably annotated (Supplementary file 1B). Twelve had been directly related to drought in previous studies, while three were previously associated with other environmental stress responses (Supplementary file 1C).

Table 1
Genes with significantly associated single-nucleotide polymorphisms (SNPs).

Given are the chromosome number (CHR), nucleotide position (position), the gene ID for Fagus sylvatica (gene), the UniProt ID of the closest match (UniProt ID), the name of the gene (name), the nucleotide base in the reference (ref DNA base), and the alternate base (alt DNA base), if applicable, the amino acid of the reference (ref AA) and the alternate base (non-synonymous change), functional change (effect), and the phenotype associated with the alternate base.

CHRPositionGeneUniProt IDNameRef DNA baseAlt DNA baseRef AANon-synoynmous changeEffectPhenotype associated with alternate base
1403747621 .g3851.t1NoneAGCRSH side chain > positive chargeHealthy
103229064510 .g3914.t1NoneCTF--
112047962811 .g2467.t1EXOS5_ORYSJExosome complex exonuclease RRP46 homologTCTAPolar > hydrophobicDamaged
112372230711 .g2832.t1PCN_ARATHWD repeat-containing protein PCNAGIVHydrophobic > hydrophobicDamaged
121390103412 .g1695.t1F4I5S1_ARATHPB1 domain-containing protein tyrosine kinaseCAPQHydrophobic > polarDamaged
13901063CTQStopTerminationDamaged
13901082ATHLPositive charge > hydrophobicDamaged
13901094TAINHydrophobic > polarDamaged
2433265712 .g4736.t1NoneGARCPositive charge > SH side chainDamaged
3312269403 .g3590.t1NoneCTQStopTerminationHealthy
4340770174 .g3980.t1CKX1_ARATHCytokinin dehydrogenase 1CAGCNo side chain > SH side chainDamaged
5163595875 .g1807.t1GDI2_ARATHGuanosine nucleotide diphosphate dissociation inhibitor 2TCP--
6198653116 .g2227.t1NDUS7_ARATHNADH dehydrogenase (ubiquinone) iron-sulfur protein 7, mitochondrialTCD--
6263831726 .g2921.t1NoneTCA--
714939047 .g177.t1TLP10_ARATHTubby-like F-box protein 10CTG--
7202420237 .g2350.t1PRK4_ARATHPollen receptor-like kinase 4AGP--
745047997 .g1655CTWStopTerminationDamaged
7314566947 .g3617.t1LSH4_ARATHProtein LIGHT-DEPENDENT SHORT HYPOCOTYLS 4GTR--
7331100007 .g3816.t1NoneGAL--
745048137 .g552.t1NoneGAG--
4504831CTL--
8292951398 .g3494.t1VATC_ARATHV-type proton ATPase subunit CGAG--
9255388279 .g3080.t1PPA14_ARATHProbable inactive purple acid phosphatase 14GCKNPositive charge > polarDamaged
9379557159 .g4504.t1TBL33_ARATHProtein trichome birefringence-like 33GCMIHydrophobic > hydrophobicDamaged

Genomic prediction

We furthermore set out to determine how many SNPs were needed to successfully predict the drought susceptibility of individual trees, that is, to develop a genotyping assay. All Pool-GWAS SNPs in addition to the top 20 individual resequencing SNPs were used to create an SNP combination to reach a genotyping success threshold of min 90%. After excluding loci due to technical reasons and filtering for genotyping success, 70 loci proved to be suitable for reliable genotyping with an SNP assay. We genotyped only individuals sampled in 2019 that were not used to identify the SNPs in the first place plus paired individuals sampled in August 2020. On average, each of the 95 individuals was successfully genotyped at 67.7 loci (96.7%). We coded the genotypes as 0 for homozygous reference allele, 1 for heterozygous, and 2 for the homozygous alternate allele, thus assuming a linear effect relationship. Figure 3B shows the genotypogram for the tested individuals.

Linear discriminant analysis (LDA) correctly predicted the observed phenotype from the genotype in 91 of 92 cases (98.9%). Prediction success decreased to 65% when successively removing loci from the analysis (Figure 3—figure supplement 5). Nevertheless, ordering the individuals according to the LDA score of axis 1 revealed no clear genotype pattern that distinguished healthy from damaged trees (Figure 3C). Observed heterozygosity at loci used in the SNP assay of individuals in the upper half of predictive values for a healthy phenotype was not significantly different from heterozygosity of the lower half (Figure 3—figure supplement 6). Ordering the loci according to their squared loadings showed that loci’s contribution to the genomic prediction differed substantially (Figure 3). As expected, the histogram of LDA scores showed two peaks, corresponding to the two phenotypes (Figure 3—figure supplement 7).

To validate the results of the LDA prediction and circumvent potential overfitting due to the small sample size, we also applied a non-parametric Machine Learning algorithm for feature selection and clustering that was especially designed for small sample sizes (Gerber et al., 2020; Horenko, 2020). The method identified the 20 most-significant SNPs allowing to make an almost 85% correct classification that distinguished healthy from damaged trees (Supplementary file 1G).

Discussion

Over the last two decades, increasing drought periods caused severe damage to European forests (Schuldt et al., 2020; Etzold et al., 2019; Pretzsch et al., 2013). Conifers seem to suffer the most, but deciduous trees were also strongly affected (Schuldt et al., 2020). Weather data from our study area from 1950 onwards suggested that the climatic conditions for beech trees in the area investigated changed dramatically during this period. Roughly estimating the tree age from their trunk circumference (Bošeľa et al., 2014), more than a third of the trees were already in place at the beginning of this period. About 60% were recruited prior to the acceleration of temperature change from the 1980s onwards. As a result, trees in the mountainous regions of the study area today experience climatic conditions comparable to those experienced by lowland trees in the 1950s, which in turn now experience a climate that used to be typical for regions much further South. Given the documented propensity of beech for local adaptation (Gárate‐Escamilla et al., 2019; Pluess et al., 2016; Aranda et al., 2015), including drought (Bolte et al., 2016), it is therefore conceivable that current conditions exceed the reaction norm of some previously locally well-adapted genotypes with detrimental consequences for their fitness. If the trend of an increasingly drier vegetation period persists, this will likely affect an even larger proportion of the currently growing beeches.

Evolutionary genomics will be indispensable to predict and manage the impact of global change on biodiversity (Waldvogel et al., 2020). As already shown for other partially managed (tree) species (Stocks et al., 2019), in particular pool-GWAS approaches (Endler et al., 2016) have proven to be useful in guiding conservation management.

Our strictly pairwise sampling design avoided many pitfalls of GWAS studies, arising, for example, from cryptic population structure and shared ancestry (Hoban et al., 2016; Wellenreuther and Hansson, 2016). Despite presented evidence from this and other studies (Schuldt et al., 2020) that the observed crown damages in large parts of Central Europe used for phenotyping here are directly or indirectly due to the severe drought years 2018 and 2019, we must acknowledge that we have no direct physiological proof that the trees surveyed here indeed suffered from drought stress. In addition, the observed diagnostic symptoms are not specific to drought stress. Nevertheless, an unknown independent stressor would have needed to accidentally co-occur spatially and temporally with the drought. The phenotypical drought response of individual trees may also be influenced by microspatial variation (Carrière et al., 2020). In the present study, however, the mean distance between sampled paired trees of about 5 m assured that their roots systems largely overlapped. Thus, environmental variation in soil quality, rooting depth, water availability, or other factors should have been minimal. Please note that any phenotypical misclassification due to such microspatial differences would have rather dissimulated the genotypic differences found in GWAS than enhanced them artificially.

As expected from previous studies (Rajendra et al., 2014), we found no population structure among the sampling sites. Applying relatively strict significance thresholds, we found systematic genomic differences between the healthy and damaged trees. In all cases, these differences were quantitative and not categorical, that is, we found allele frequency changes but no fixed SNPs between phenotypes. Significant SNPs were mostly not clustered – we found on average 1.4 selected SNPs in a particular genomic region. These findings were in line with the observed very short average LD in F. sylvatica, indicating that polymorphisms associated with the two phenotypes were likely old-standing genetic variation (Harris and Nielsen, 2013). Moreover, such SNPs are mostly detached from the background in which they arose and they are therefore often the actual causal variants or at least in very close proximity (Schaid et al., 2018). This observation is underlined by the high proportion of non-synonymous significant SNPs within genes, which in most cases caused substitution to an amino acid with different properties or even premature termination. Such deviant variants with likely substantial functional or conformational changes in the resulting proteins may be selectively neutral or nearly neutral under ancestral benign conditions, but may become selectively relevant under changing conditions (Paaby and Rockman, 2014). Interestingly, most of the allelic variants associated with a healthy phenotype were also the variants in the reference genome. This might be due to the choice of the F. sylvatica individual from which the reference genome was gained (Mishra et al., 2018). This more than 300-year-old individual is standing at a particularly dry site on a rocky outcrop on the rim of a scarp where precipitation swiftly runs off. Trees at such sites were likely selected for drought tolerance.

Even though the area sampled for this study was limited relative to the species distribution range, it comprised its core area. In addition, the climatic variation covered by the sampling sites for this study is representative for large parts of the species range (Baumbach et al., 2019). The relatively limited population structure over large parts of the species range (Magri et al., 2006), together with the propensity for long-range gene flow (Belmonte et al., 2008), suggested that the genomic variation responsible for drought tolerance identified here is widely distributed (Lander et al., 2021). Nevertheless, an assessment of the geographic distribution of the drought-related genomic variants over the entire distribution range would yield general insight into the species-wide architecture of this important trait.

None of the genes found here was involved in a transcriptomic study on drought response in beech saplings (Müller et al., 2017). However, most of the reliably annotated genes with or close to SNP loci significantly associated with drought phenotypes had putative homologs in other plant species previously shown to be involved in drought or different environmental stress response (for citations, see Supplementary file 1C, D). For the remaining, not annotated genes, it remained unclear whether they had really never been associated with drought before, or whether we were just unable to make this link due to the lack of (ecological) annotation and standardised reporting (Waldvogel et al., 2021). The involvement of in total 67 genes, together with the relatively flat effect size distribution, suggested that drought resistance in F. sylvatica is a moderately polygenic trait, which should respond well to artificial breeding attempts and natural selection. However, given the relatively strict threshold criteria, it is likely that more yet undetected loci contribute to the respective phenotypes. The low LD in beech predicts that an adaptation to drought will not compromise genome-wide genetic diversity and thus adaptation potential to other stressors. We achieved a high level of accuracy using genomic data to predict the drought phenotype from individuals not used to identify drought-associated SNP loci. However, due to the small sample size, LDA might have resulted in overfitting (Hawkins, 2004). We therefore also used a non-parametric Machine Learning algorithm that has been shown to produce more robust results, especially for small sample sizes (Horenko, 2020). Both analyses confirmed that we mainly identified alleles widespread throughout the sampled range and not locally specific. Besides, we confirmed a considerable level of genetic variation in the sampled regions. The observation that trees with the highest predictive values showed no loss of heterozygosity indicated that there is still adaptive potential for drought adaptation in the species (Gienapp et al., 2017). With the SNP assay, we therefore created a tool that can (i) support the choice of seed trees for reforestations, (ii) provide decision guidance for selective logging, and (iii) monitor whether natural selection on this quantitative trait is already acting in the species. The current study can also serve as a starting point for molecular and physiological research on how the identified loci or variants may, alone or in concert, confer resilience or tolerance to a range of drought stress symptoms.

Materials and methods

Sampling and phenotyping

Request a detailed protocol

In August/early September 2019, we sampled leaf tissue of 402 F. sylvatica trees from 32 locations in Hessen/Germany (set 1, Figure 1), of which 300 were used for the (pool)GWAS analysis. 43, plus additional 53 trees which were sampled in August 2020, additional 52 trees from four sites were sampled (set 2, Figure 1) which made up the confirmation set. The coordinates and characteristics of each site can be found in Supplementary file 1A. The sampling was performed in a strictly pairwise design. The pairs consisted of one tree with heavy damage of the crown (lost or rolled up, dried leaves) and one with an unaffected crown, respectively. This categorisation into least and most damaged trees was taken compared to the other trees in the respective forest patch. The pairs were a priori chosen such that the two trees were (i) mutually the closest neighbours with contrasting damage status (i.e. no other tree in the direct sight line), (ii) free from apparent mechanical damage, fungal infestations, or other signs of illness, similar (iii) in tree height, (iv) trunk circumference, (v) light availability, and (vi) canopy closure. In addition, each pair was situated at least 30 m from the closest forest edge. For each tree of the chosen pairs, we recorded the exact position, distance to the pair member and the estimated tree height (in 1 m increments), measured the trunk circumference at 150 cm height above the ground (in 10 cm increments), and estimated the leaf loss of the crown and the proportion of dried leaves (in 5% increments). We also recorded the estimated distance (in 1 m resolution) and the specific identity of the two closest neighbour trees for each pair member and calculated a competition index C as follows: C = S1/D1 +S2/D2, where S1and S2 are the trunk diameter at 150 cm and D1and D2 the distances of the nearest and second nearest neighbour tree of the same size or larger than the focal tree. Photographs from the crown and the trunk were taken from the trees sampled in 2019.

From each tree, we sampled 5–10 fully developed leaves from low branches. The leaves sampled from each tree were placed in paper bags. After returning from the field, they were dried at 50°C for 30–90 min and then kept on salt until they could be stored at −80°C.

Climate and remote sensing data

Request a detailed protocol

Monthly daily mean minimum and maximum temperature values and precipitation data were obtained for the 1 × 1 km grid cells harbouring the sampling sites for the period between 1950 and 2019. Data on the accumulated potential evapotranspiration during the growth season was obtained for the same grid cells. The data is publicly available from https://opendata.dwd.de/climate_environment/CDC/grids_germany/monthly/.

LAI data for the above grid cells was obtained from Copernicus (http://www.copernicus.eu) for the period 2014–2019, considering only the month of August. To see whether drought conditions influenced leaf coverage of the woods at the sampling sites, we calculated the relative annual deviation of LAI from the 2014 value. We correlated it to the relative deviation of the cumulated potential evatransporation over the growth season from 2014. The year 2014 was used as a baseline because of the significant drought increase since then (Büntgen et al., 2021). Please note that the absolute level of LAI depends on the wood coverage, vegetation density, and species composition of each plot. Changes in LAI are thus not exclusively due to drought damage in beech.

DNA extraction, construction of GWAS pools, and sequencing

Request a detailed protocol

DNA was extracted from 12.5 mm² of a single leaf from each tree following the NucleoMag Plant Kit (Macherey Nagel, Düren, Germany) protocol. We setup four DNA pools for Pool-GWAS by pooling equal amounts of DNA from each individual: damaged individuals from the Southern part (dSouth), healthy individuals from the South (hSouth), damaged North (dNorth), and healthy North (hNorth). The Southern pools consisted of 100 individuals each, the Northern pools of 50 individuals each. The pools were sent to Novogene (Cambridge, UK) for library construction and 150 bp paired end sequencing with 350 bp insert size with 25 Gb data for the Northern and 38 Gb data for the Southern samples. The 100 individuals used to construct the Northern pools were also individually resequenced. The exact composition of the genomic pools can be found in Supplementary file 1A. All sequence information can be found on the European Nucleotide Archive (ENA) under project accession number PRJEB24056.

Reference genome improvement

Request a detailed protocol

We used an improved version of the recently published reference genome for the European beech (Mishra et al., 2018). Contiguity was improved to chromosome level using Hi-C reads with the help of the allhic software after excluding the probable organelle backbones from the earlier assembly that was generated from the Illumina-corrected PacBio reads using Canu assembler (Mishra et al., 2021) accession no. PRJNA450822.

Mapping and variant calling

Request a detailed protocol

Reads of pools and individual resequencing were trimmed using the wrapper tool autotrim v0.6.1 (Waldvogel et al., 2018) that integrates trimmomatic (Bolger et al., 2014) for trimming and fastQC (Andrews, 2010) for quality control. The trimmed reads were then mapped on the latest chromosome-level build of the F. sylvatica genome using the BWA mem algorithm v.0.7.17 (Li and Durbin, 2009). Low-quality reads were subsequently filtered and SNPs were initially called using samtools v.1.10 (Li et al., 2009). A PCA was conducted on unlinked SNPs using the R package Factoextra v.1.0.7 (Kassambara and Mundt, 2017).

Pool GWAS and PLINK

Request a detailed protocol

The PoPoolation pipeline 2_2012 (Kofler et al., 2011a; Kofler et al., 2011b) was used to call SNPs and remove indels from the four pools. Allele frequencies for all SNPs with a coverage between 15× and 100× with a minimum allele count of 3 were estimated with the R library PoolSeq v. 0.35 (Taus et al., 2017).

The statistical test to detect significant allele frequency differences among damaged and healthy trees was the Cochran–Mantel–Haenszel test. With this test, a 2 × 2 table was created for each variable position and region with two phenotypes (healthy and damaged). The read counts of each allele for each phenotype were treated as the dependent variables. We controlled for false discovery rate using the Benjamini–Hochberg correction in R package p.adjust.

For the individual resequencing data, we followed the GATK-pipeline 4.1.3.0 (DePristo et al., 2011). In short, Picard tools v.2.20.8 was used to mark duplicates. GVCF files were created with HaplotypeCaller and genotyped with GEnotypeGVCFs. Since we did not have a standard SNP set, we hard-filtered SNPs with VariantFiltration QD < 2.0, MQ < 50.0, MQRankSum < 12.5, ReadPosRankSum < 8.0, FS > 80.0, SOR > 4.0, and QUAL < 10.0. This conservative SNP set was used for base recalibration before running the HaplotypeCaller pipeline a second round. Finally, the genotyped vcf-files were filtered using vcftools with --maf 0.03 --max-missing 0.9 --minQ 25 --min-meanDP 10 --max-meanDP 50 --minDP 10 --maxDP 50. The detailed pipeline can be found in Supplementary file 1G.

To conduct the GWAS association on the above-generated SNP set with phenotypes being either damaged or healthy and to generate a PCA on the SNP positions of the individually resequenced trees, we used PLINK 1.9 (Purcell et al., 2007). The detailed workflow can be found in Supplementary file 1G. We calculated a non-parametric ANOSIM on an inter-individual Euclidean distance matrix based on the first 10 principal components to infer whether the trees within phenotype groups are overall genetically more similar than within groups (9999 permutations; Hammer et al., 2001).

Inference of linkage disequilibrium

Request a detailed protocol

The expected length of segregating haplotypes in a species depends on the recombination rate and their age. The former can be approximated by an estimate of LD. To determine LD decay based on individually resequenced data, we used the software LDkit v 1.0.0 (Tang et al., 2020), in 1 kb and 100 kb windows.

Identification substitution type and gene function

Request a detailed protocol

We inferred whether significantly differentiated SNPs within genes lead to a (non-) synonymous amino acid substitution using tbg-tools v0.2 (https://github.com/Croxa/tbg-tools; Schoennenbeck et al., 2021). The protein sequences of the identified genes were used in a blastp search against all non-redundant GenBank CDS translations, PDB, SwissProt, PIR, PRF to infer potential gene functions. For each search, only the single BLAST- top hit was considered.

Selection of SNP loci for SNPtype assay design

Request a detailed protocol

For the design of SNPtype assays, we used the web-based D3 assay design tool (Fluidigm Corp.). We aimed in first preference for the most significant SNPs of each genomic region identified by Pool-GWAS (80 loci). If this was technically impossible and the region harboured more than a single significant SNP, we opted for the second most significant SNP and so forth. This resulted finally in 76 suitable loci. The remaining 20 loci were recruited from the 20 most significant SNPs of the PLINK analysis that were not scored in the Pool-GWAS.

SNP genotyping procedure

Request a detailed protocol

For validation of drought susceptibility-associated SNPs, we conducted SNP genotyping on 96.96 Dynamic Arrays (Fluidigm) with integrated fluidic circuits (Wang et al., 2009) (N = 96) to validate the effectiveness of the identified SNPs in discriminating healthy from damaged trees. Prior to genotyping PCR, DNA extracts were normalised to approximately 5–10 ng/µl. They underwent a pre-amplification PCR (Specific Target Amplification [STA]) according to the manufacturer's protocol to enrich target loci. PCR products were diluted 1:10 with DNA suspension buffer (TEKnova, PN T0221) before further use. Genotyping was performed according to the manufacturer's recommendations. Four additional PCR cycles were added to accommodate for samples of lower quality or including inhibitors (von Thaden et al., 2020). Fluorescent data were measured using the EP1 (Fluidigm) and analysed with the SNP Genotyping Analysis Software version 4.1.2 (Fluidigm). The automated scoring of the scatter plots was checked visually and, if applicable, manually corrected.

Genomic prediction

Request a detailed protocol

To predict drought susceptibility from genotype data, we used an LDA on 92 genotypes scored with the Fluidigm assay at 70 loci. Genotypes homozygous for the reference allele were scored as 0, heterozygous as 1 and homozygous alternate alleles as 2. We used the LDA option implemented in PAST v. 4.05. (Hammer et al., 2001) .

We also used a non-parametric entropy-based Scalable Probabilistic Analysis framework (eSPA). This method allows simultaneous solution of feature selection and clustering problems, meaning that does not rely on a particular choice of user-defined parameters and has been shown to produce more robust results, especially for small sample sizes (Gerber2020, Horenko2020). Following the suggestion of the user manual, eSPA was run 100 times with independent cross-validations of the area under the curve (AUC) on the validation data.

Data availability

Sequencing data have been deposited at ENA under project code PRJEB41889. The genome assembly including the annotation is available under the accession PRJNA450822.

The following data sets were generated
    1. Pfenninger M
    2. Feldmeyer B
    (2020) European Nucleotide Archive
    ID PRJEB41889. Genomic basis of drought resistance in Fagus sylvatica by PoolGWAS.

References

  1. Report
    1. Andrews S
    (2010)
    FastQC: A Quality Control Tool for High Throughput Sequence Data
    Cambridge, United Kingdom: Babraham Bioinformatics, Babraham Institute.
    1. Bressem U
    (2008)
    Komplexe erkrankungen an buche. Complex diseases in beech
    Ergebnisse Angewandter Forschung Zur Buche 3:87.
    1. Brunet J
    2. Ö F
    3. Richnau G
    (2010)
    Biodiversity in european beech forests-a review with recommendations for sustainable forest management
    Ecological Bulletins 53:77–94.
  2. Report
    1. Christensen JH
    2. Hewitson B
    3. Busuioc A
    (2007)
    Regional Climate Projections Climate Change, 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change.
    University Press.
    1. Czajkowski T
    2. Bolte A
    (2006)
    Frosttoleranz deutscher und polnischer herkünfte der buche (Fagus sylvatica L.) und ihre beeinflussung durch trockenheit
     Archiv Für Forstwesen Und Landschaftsökologie 40:119–126.
  3. Report
    1. Elsasser P
    2. Meyerhoff J
    3. Weller P
    (2016)
    An Updated Bibliography and Database on Forest Ecosystem Service Valuation Studies in Austria, Germany and Switzerland
    Thünen Working Paper.
    1. Hammer Ø
    2. Harper DA
    3. Ryan PD
    (2001)
    PAST: paleontological statistics software package for education and data analysis
    Palaeontologia Electronica 4:9.
    1. Hawkins DM
    (2004) The Problem of Overfitting
    Journal of Chemical Information and Computer Sciences 44:1–12.
    https://doi.org/10.1021/ci0342472
    1. Müller-Stark G
    (1985)
    Genetic differences between" tolerant" and" sensitive" beeches (Fagus sylvatica L.) in an environmentally stressed adult forest stand
    Silvae Genetica 34:241–246.
  4. Book
    1. Paar U
    2. Dammann I
    (2019)
    Waldzustandsbericht Hessen 2019
    Wiesbaden: Hessisches Ministerium für Umwelt K, Landwirtschft und Verbraucherschutz.
    1. Sutmöller J
    2. Spellmann H
    3. Fiebiger C
    4. Albert M
    (2008)
    Der klimawandel und seine auswirkungen auf die buchenwälder in deutschland the effects of climate change on beech forests in Germany
    Ergeb. Angew. Forsch. Buche 3:135.
  5. Book
    1. Thünen_Institute
    (2020)
    Clusterstatistik
    Forst & Holz.
  6. Conference
    1. Waldvogel A-M
    2. Schreiber D
    3. Pfenninger M
    4. Feldmeyer B
    (2021)
    Climate change genomics calls for standardized data reporting
    Coping with Climate Change: A Genomic Perspective on Thermal Adaptation.

Decision letter

  1. Meredith C Schuman
    Senior and Reviewing Editor; University of Zurich, Switzerland
  2. Lucienne de Witte
    Reviewer; Institute for Applied Plant Biology, Switzerland
  3. Meredith C Schuman
    Reviewer; University of Zurich, Switzerland

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Acceptance summary:

This study is a pioneering example of something like "personalized medicine" in the sense of diagnosis, applied to forests. The authors study a tree species that forms a large proportion of European forests, using real trees in medium- to old-growth forests. Using a paired sampling design combined with ecological and climatological observations, they are able to diagnose variants indicative of stress resistance or tolerance which vary among neighboring trees. In doing so the authors present new quantitative evidence that European beech forests are strongly affected by climate change, and that the genome of beech harbors resistant genes which may be broadly distributed in trees across its range.

Decision letter after peer review:

Thank you for submitting your article "Genomic basis of drought resistance in Fagus sylvatica" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Meredith C Schuman as the Reviewing and Senior Editor and Reviewer #3. The following individual involved in review of your submission has agreed to reveal their identity: Lucienne de Witte (Reviewer #1).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this evaluation to help you prepare a revised submission.

Summary:

The reviewers agree that the study is of broad interest and that the genetic work represents a substantial advance on existing tools in European beech, an economically and ecologically important forest species with an uncertain future under global change. The reviewers value the case-control design using trees in standing populations and taking advantage of existing environmental data. There are questions about many important details both regarding tree selection and phenotyping, and also the genetic analysis, which are missing from the manuscript and need to be addressed. The results provide a valuable basis for the development of hypotheses regarding the adaptive potential within beech populations facing global change, but many additional steps would be required to articulate and test specific hypotheses emerging from these results. Thus we would consider a revised version of this manuscript, addressing these essential revisions and considering the other points raised, for publication as a Tools and Resources article formatted as a Short Report, to accompany the sequencing dataset that the authors have already deposited at ENA, including an update of the beech genome (also at ENA/beechgenome.net) corresponding to the revisions made in the course of this study. The revision would need to include a full description (including specification of all data, software and parameter settings used) of how the assembly was improved, as well as details of the annotation (both genes and repeats).

Essential Revisions:

This is a list of revisions considered essential, referring to more detailed information from the individual reviews where relevant. The individual reviews are also appended in full for your information, including public reviews and a public summary designed to accompany your preprint, and recommendations for the authors from each reviewer.

1) The authors state that they selected trees which were the nearest neighbors showing distinct phenotypes (apparently unaffected canopies versus canopies with dead and curled, dried leaves) and not displaying any apparent indications of insect infestation or disease which would cause such differences. This stressed versus healthy classification used by the authors is not specific to drought. To the extent that their identified genetic targets could be annotated, most of these are also generally stress responsive and not specific to drought stress. Drought is a complex stressor and it is certainly plausible that drought could be the main cause of different stress status in this study, but the authors do not present evidence to support this other than the general knowledge that beech forests suffered under the European 2018 summer drought and other recent droughts.

– More direct indicators of water stress would provide better support for their classification, if any such data exist. Please refer to comments and specific suggestions from Reviewer #1 – Recommendations for authors point 1, and Reviewer #3 – Public review, "major weaknesses" section.

– Information on the relative health status of the chosen trees prior to the 2018 extreme drought could also be valuable, and perhaps the different climatic conditions of their study populations could be better leveraged to support the importance of drought in explaining different apparent stress status in the pairs of trees they sampled.

– To the extent that the authors cannot provide data on study trees which is more specific to drought stress, they must acknowledge that their stress classification cannot be clearly attributed to drought susceptibility of the trees, and may incorporate other (multiple) stress factors.

2) Beyond the concerns regarding interpretation of tree stress status, sufficient information must be provided regarding the trees studied.

– The authors should provide evidence to support their classification of trees as healthy versus stressed (such as images – even example photographs, but better systematic photographs or aerial imagery of pairs; and/or supporting measurements; field rankings of e.g. crown transparency or apparent stress based on specific criteria).

– The authors should provide at least a rough measurement of how far pairs of trees were allowed to be apart (e.g. minimum and maximum distance in m) and roughly how many other trees were allowed to stand in between (to what extent can these trees be considered neighbors). Please refer to Reviewer #3 – Public review, "major weaknesses" section.

– The authors must either provide information on known or expected within-site variation of water availability and other relevant factors such as soil quality, possible rooting depth, and differences in biotic interactions which could have affected drought status or responses, on the scale of distance between their selected pairs; or at a minimum, they should discuss possible environmental variation on this scale based on literature such as that suggested by Reviewer #1.

– There is conflicting information regarding numbers of individuals sampled, and included in each group, which must be resolved. Please refer to comments from Reviewer #2 under Recommendations for authors – Major suggestions – Number of individuals used.

3) There is also important information missing regarding sample and data analysis, and this must be provided, and concerns regarding non-standard methods addressed. Please see detailed comments from Reviewer #2 (Recommendations for authors – Major suggestions) regarding the following:

– Insufficient information is provided regarding the improvement of the reference genome.

– The authors used unevenly sized pools for poolGWAS and equal membership of pooled samples can also not be guaranteed based on the information provided in the Materials and methods.

– Insufficient information is provided regarding how GWAS was conducted.

– It is not clear which indicators of statistical significance are reported and how these are to be interpreted.

– The reported very low LD may be an artifact of interpreting LD from poolGWAS data, and should be assessed using the individual genomes which were sequenced at high coverage.

– It is not clear how group membership prediction was conducted and how the very high success rate is to be interpreted (GP analysis).

– The functional interpretation of SNPs is complicated by some apparent contradictions in genome locations and gene models, some of which may be resolved by more information from the improved genome assembly.

– Insufficient information is provided regarding software versions and parameters used, and also reasoning for choosing some non-standard software tools.

– Additionally, from Reviewer # 3, Recommendations for authors: The leaf handling protocol is unusual and it is not very precisely reported: collected fresh from field, 50 deg C in lab 30-90 min, then to salt until storage at -80. It seems like the DNA would be better preserved if leaves were harvested directly to silica gel and not heated as is commonly done (in which case storage at -80 deg C can be done but is generally not required).

4) One important hypothesis suggested by this study is that beech populations locally embody sufficient variation in drought-responsive genes to adapt to increased drought severity and frequency. However, this cannot be supported as a conclusion of the study. The authors should re-state this as a hypothesis and propose work required to test this hypothesis rigorously. As noted by Reviewer 1, the candidate genes revealed often do not overlap with those reported from other studies, including studies in beech, and as noted by Reviewer 3, the geographic limitations of the study preclude conclusions regarding to what extent the study populations embody the potential of members of the species to adapt to drought. Furthermore, there is insufficient phenotypic data to link the identified loci to specific response curves (e.g. change in ability to tolerate or recover from a given severity of water deficit and associated stresses such as salinity and heat). The ability to predict membership in a group (healthy versus stressed) by genotype demonstrates the consistency of associations identified within this dataset, but it does not indicate that the study populations are resilient to predicted future increases in drought frequency and severity. Again, rather, these results suggest hypotheses about how specific loci may, alone or in concert, confer resilience or tolerance to a range of drought stress, which could be better articulated based on knowledge in the beech system and gene annotation; and the authors could suggest follow-on studies, increasing the value of their own study as a resource.

Reviewer #1 (Recommendations for the authors):

1. The authors study the genomic basis of drought susceptibility and found systematic genomic differences. However, selection criteria of damaged trees are unclear because trait variables used for phenotype differentiation seem not consistent with protocols from internationally recognized monitoring networks.

In the first paragraph of the introduction (lines 55-58 and 68-74), it is stated that drought damage in beech includes rolled leaves and small leaves. But these two symptoms rather are resistance/stress symptoms, see also Wohlgemuth et al. 2020 (Schweiz Z Forestwes 171). Beech trees roll their leaves and produce smaller leaves to avoid to much evapotranspiration. The literature cited, Paar and Dammann 2019, state stress symptoms (Trockenstresssymptome) versus damage, the later defined as decreased growth, increased crown transparency and increased mortality. These definitions are according to standardized monitoring protocols used by international networks such as ICP forests (ICP forest manual). Why were these variables not included for tree selection?

In the Results chapter, lines 113-118, the damage status was assigned according to the quantity of dried leaves and leaf loss (Figure 2 E-F, Suppl. Table 2 is referred to but that table does not include any info on these variables, only the results from the Mann-Whitney U-test on differences between damaged and healthy trees). What do dried leaves and leaf loss mean? These are not usual terms or variables used to describe or name damage in beech trees in standardized monitoring. Do they translate to leaf discoloration and crown transparency?

In the Materials section (lines 237-254) you state again differently named traits: "heavy drought damage of the crown (lost or rolled up, dried leaves) (…) leaf loss of the crown and the proportion of dried leaves (in 5% steps)."

On lines 75-86 you use the terms „drought susceptibility" versus „resistant" trees which is good. The wording is more concise in this paragraph than first paragraph of the introduction.

Attention should be paid also to the damages recorded during monitoring: crown transparency, discoloration, decrease in growth and increase in mortality cannot be explained only by drought stress. Modelling analysis of long-term data for example shows that drought, next to other reasons, is a main reason for the actual damages observed in central European beech forests (Braun et al. 2020, Schweiz Z Forstwes 171). Observed damages in beech trees are not the consequence only of drought (temperature increase and precipitation decrease), but also of storm events (root damages!) and environmental pollution (including ozone and nitrogen deposition impact soil quality and growth of roots and associated microbiome), and because these factors also interact with drought. These factors should at least be considered and discussed (see also Pflug et al. 2018).

Reviewer #2 (Recommendations for the authors):

This is an interesting study on the genomic basis of drought resistance in an important forest tree species. However, I have concerns about some of the analyses and there are several aspects of the paper where there is a lack of detail or where further clarification is needed. I outline these points below.

Number of individuals used

I am confused about the number of individuals included in some of the analyses because there are conflicts between the details provided in different parts of the manuscript. For example, Lines 76-77 of the introduction suggest that more than 400 individuals were used for the pool-GWAS ("more than 200 neighbouring pairs of trees"), but Line 120 states that "300 beech trees" were used for the pool-seq. Lines 122-123 say that there were 100 individuals in each of the South pools, but, Table S1 indicates there are 135 trees in each South pool. On Lines 127-128 it states "All 100 individuals from the North population" were individually re-sequenced, but the Materials and Methods say 102 individuals were re-sequenced (Line 271) and Table S1 indicates that it was 95 individuals. On Lines 238-239 it says 53 trees were sampled for "set 2", but Table S1 indicates there are 80 individuals in set 2. There are various other similar examples, which altogether create a confusing picture. The Authors need to double-check all such details and correct them as necessary.

Reference genome

More details about the methods used for improving the reference genome are needed. For example, please provide more details of how the Hi-C reads were generated (i.e. source of DNA, extraction method, sequencing coverage) and add details of the version number and parameter settings used for the allhic software. Also, was the improved assembly reannotated? If so, how? If the original annotation was transferred to the new assembly, this should also be explained.

GWAS

The guidance for popoolation (https://sourceforge.net/p/popoolation2/wiki/Manual/) states that, for the experimental design, the most important thing is "the amount of DNA per individual in a single pool should be constant". In the Materials and Methods section that deals with DNA extraction (starting on Line 264), there is no mention of quantification of DNA being done to ensure that the amount of DNA contributed to a pool by each individual is equal. Please can the Authors provide details of how they quantified DNA from each individual prior to pooling to ensure they contributed equally? The guidance also states that, "the number of individuals per pool should be similar". However, the Southern pools are much larger than the Northern pools (more than double the size according to the information in Table S1). Could the Authors comment on why they used such unevenly sized pools and what impact this might have had on the results of the GWAS?

With regard to the potential impact of population stratification, which could create spurious associations in the GWAS, the Authors state on Lines 201-202 "As expected from previous studies, we found no population structure among the sampling sites". However, I cannot see anywhere an analysis of potential structure between all populations sampled for the main GWAS analysis (i.e. the pool-GWAS). Could the Authors comment on this further please?

Also, in relation to the control for FDR, the method used by the Authors (the R package qvalue) does not calculate "corrected" p-values. Normally for GWAS, false discovery rates are used to set an appropriate p-value significance threshold to control for the number of false discoveries, rather than to correct p-values as with multiple comparison correction. I.e. a threshold can be set under which the number of false positives is deemed to be acceptable. In Figure 4 the Authors present "corrected" -log10 p-values and indicate a p-value significance threshold. I think the values in this plot are actually the q-value estimates that are generated by the qvalue package. I suggest the Authors present the original p-values from the association analysis, but indicate a significance threshold depending on the number of false positives they are willing to allow within the set of "significant" SNPs. If this changes which SNPs are considered to be significant, then clearly any subsequent analyses and results that are presented should account for this fact.

Furthermore, for the GWAS analysis using data from the re-sequenced individuals (mentioned in Lines 303-308), please provide full details for what was actually done, including any additional filtering of SNPs done in PLINK (e.g. LD pruning, which is mentioned in the legend for Figure 3, but not in the Materials and Methods) and which statistical test for association was performed. Moreover, there don't seem to be any details for this GWAS analysis in the Results section, although it is mentioned on Lines 326-327 of the Materials and Methods that some of the significant SNPs were from this GWAS were included in the SNP array.

LD analysis

The Authors estimated LD in beech on the basis of a subset of the pool-seq data. This suggests a very short size of linkage block on average (c.120bp), which led the Authors to conclude that most of the significant SNPs they detected are "the actual causal variants" as they are unlikely to be linked to other SNPs. However, it has been suggested that the estimates of LD obtained from analysis of pool-seq data are not very accurate (Schlötterer et al., 2014; doi:10.1038/nrg3803). The finding that SNPs that are further than c. 120bp apart are likely to be unlinked is surprising to me given the length of linkage blocks that have been detected in some other tree species. As the Authors have whole genome sequence data for individuals from the Northern pools, I would recommend that they also use these data to estimate LD to see if this produces a result that is congruent with that based on the pool-seq data and supports their assertion that they have likely identified causative SNPs.

GP analysis

The details of the GP analysis need further clarification. In particular, I am confused by the reported level of accuracy of GP with the set of 70 variants. The Authors say "Linear discriminant analysis (LDA) correctly predicted the observed phenotype from the genotype in 91 of 92 cases (98.9%)". However, in the Materials and Methods it is reported that 80% of individuals were used in training and only 15 individuals were used in the test set. Can the Authors please make clear what they actually did to test the accuracy of the LDA analysis in discriminating trees with the different phenotypes? If only 15 trees were included in the test set, then accuracy should be reported based on the results of these. Also, was any replication done to see if results differ with different subsamples of 15 individuals (i.e. using multiple random samples of 15 individuals from the full set)? If not, then I think this would be worth considering. Alternatively, have the Authors considered also using some of the pool-seq data to train a GP model? This could be used to test the accuracy of prediction for the individually sequenced trees, which would provide a larger number of trees for testing. Further, the details of the number of individuals used for the GP analysis and numbers and percentages of individuals used for training and testing need to be double-checked (lines 342-346), as there seem to be some inaccuracies here.

Functional significance of SNPs

tbg-tools was used to check whether significant SNPs from the GWAS cause non-synonymous substitutions in non-coding regions. Could the Authors explain why they chose to use this tool (which is described on its GitHub page as "not in its finished version and there are surely still some bugs to be found") rather than more established software such as SnpEff (https://pcingola.github.io/SnpEff/)? Also, even if the software correctly predicts the impacts of the variants, these predictions are dependant on the accuracy of the genome annotation. As noted above, there is a lack of detail about the annotation for the improved reference genome (and also no indication of the availability of the annotation file). Could the Authors comment on whether they manually checked the gene models for those cases where a significant SNP fell within a gene (i.e. results in Table 1) to see if there was any evidence for gene model errors that could impact the results? For example, two variants are reported to be within gene model 7.g2350.t1, one of which causes a premature stop codon. However, the position of these variants suggests they are >15Mb apart; assuming the positions have been correctly reported, I would be surprised if this gene model is correct, especially as the best matching gene from A. thaliana is only c. 3kb. Also, in the case of 12.g1695.t1, four variants are reported with a c. 60bp region of the gene, including one inducing a premature stop codon in individuals with the alternate allele. This leads the Authors to suggest that the gene "lost its function" in individuals that carry the non-reference version of the allele. However, another possible cause of such a pattern would be incorrect specification of intron/exon boundaries in the gene model. Have the Authors double-checked this to confirm that the variants are actually within a protein-coding region? Even if no evidence of gene model errors is found, it cannot be concluded from the sequence data alone that the gene will be non-functional, as even with a premature stop codon a truncated product with some functionality could be formed. So, I suggest the Authors need to be a bit more cautious about their conclusions (e.g "is likely to have lost its function, or be only partially functional").

Details of software

For each piece of software used, please provide details of the exact version number used and any parameter settings; if default parameter settings were used, please state this. At the moment, these details are given for some of the software, but not all. For example, no version number is given for BWA, PoPoolation and qvalue; autotrim and samtools are missing details of the parameter settings, etc.

Reviewer #3 (Recommendations for the authors):

The genetic analysis suggests important hypotheses regarding the drought resistance or resilience of European beech. However, these hypotheses are neither very precisely articulated, nor tested in this study. Some of the information necessary to better articulate the hypotheses seems to be missing, and to articulate and test them would require substantial additional work. For this reason it seems most suitable as a Tools and Resources article accompanying publication of the dataset and improved reference genome, which could be a very valuable resource as described by the authors and the reviewers.

The leaf handling protocol is unusual and it is not very precisely reported: collected fresh from field, 50 deg C in lab 30-90 min, then to salt until storage at -80. It seems like the DNA would be better preserved if leaves were harvested directly to silica gel and not heated (in which case storage at -80 deg C can be done but is generally not required).

https://doi.org/10.7554/eLife.65532.sa1

Author response

Summary:

The reviewers agree that the study is of broad interest and that the genetic work represents a substantial advance on existing tools in European beech, an economically and ecologically important forest species with an uncertain future under global change. The reviewers value the case-control design using trees in standing populations and taking advantage of existing environmental data. There are questions about many important details both regarding tree selection and phenotyping, and also the genetic analysis, which are missing from the manuscript and need to be addressed. The results provide a valuable basis for the development of hypotheses regarding the adaptive potential within beech populations facing global change, but many additional steps would be required to articulate and test specific hypotheses emerging from these results. Thus we would consider a revised version of this manuscript, addressing these essential revisions and considering the other points raised, for publication as a Tools and Resources article formatted as a Short Report.

We followed the reviewers’ suggestions wherever possible, which, however, lead to a longer rather than shorter text of the manuscript

To accompany the sequencing dataset that the authors have already deposited at ENA, including an update of the beech genome (also at ENA/beechgenome.net) corresponding to the revisions made in the course of this study. The revision would need to include a full description (including specification of all data, software and parameter settings used) of how the assembly was improved, as well as details of the annotation (both genes and repeats).

The updated genome assembly is part of another project and is published separately. This paper in the publication process and already available as preprint. It contains all methodological details.

A chromosome-level genome assembly of the European Beech (Fagus sylvatica) reveals anomalies for organelle DNA integration, repeat content and distribution of SNPs

Bagdevi Mishra, Bartosz Ulaszewski, Joanna Meger, Markus Pfenninger, Deepak K Gupta, Stefan Wötzel, Sebastian Ploch, Jaroslaw Burczyk, Marco Thines

bioRxiv 2021.03.22.436437; doi: https://doi.org/10.1101/2021.03.22.436437

Essential Revisions:

This is a list of revisions considered essential, referring to more detailed information from the individual reviews where relevant. The individual reviews are also appended in full for your information, including public reviews and a public summary designed to accompany your preprint, and recommendations for the authors from each reviewer.

1) The authors state that they selected trees which were the nearest neighbors showing distinct phenotypes (apparently unaffected canopies versus canopies with dead and curled, dried leaves) and not displaying any apparent indications of insect infestation or disease which would cause such differences. This stressed versus healthy classification used by the authors is not specific to drought. To the extent that their identified genetic targets could be annotated, most of these are also generally stress responsive and not specific to drought stress. Drought is a complex stressor and it is certainly plausible that drought could be the main cause of different stress status in this study, but the authors do not present evidence to support this other than the general knowledge that beech forests suffered under the European 2018 summer drought and other recent droughts.

– More direct indicators of water stress would provide better support for their classification, if any such data exist. Please refer to comments and specific suggestions from Reviewer #1 – Recommendations for authors point 1, and Reviewer #3 – Public review, "major weaknesses" section.

Please note that we presented evidence that the all trees within the 1 x 1 km grid cell of all sampling sites experienced significantly increased drought stress in the years 2018 and 2019 compared to the years before (Figure 1C). Moreover we cite two publications showing that the droughts in the years 2018 and 2019 caused severe water stress in many forest tree species (Schuldt et al. 2020a) and even more specifically lead to drought damage in German beech trees (Paar and Dammann 2019). Therefore, there can be no doubt that all trees sampled experienced substantial and unusual drought stress.

– Information on the relative health status of the chosen trees prior to the 2018 extreme drought could also be valuable, and perhaps the different climatic conditions of their study populations could be better leveraged to support the importance of drought in explaining different apparent stress status in the pairs of trees they sampled.

The categorization of each individual tree to damaged versus healthy was done opportunistically and represents only a snapshot of the acute symptoms expressed at the time of the study; we lack data on the status of these tree individuals in preceding or following years, let alone their growth or mortality. We acknowledge that in addition to drought events there could be multiple reasons for damage of branches or entire trees, e.g. attacks by insects, pathogens or tree competition, either independently or even as a response to drought. Trees are particularly vulnerable to insect or pathogen infestation in the following growing season after drought events (see Schuldt et al. 2020), hence such infestations can be both cause or effect of stress symptoms, associated with reduced defence potential. We cannot exclude that affected trees may already have experienced stress or infections in previous years, and our understanding of legacy effects of multiple symptoms is just beginning (Pflug et al. 2018, Schuldt et al. 2020). We thus reported a list of additional symptoms for a broader description of the trees’ status and their environment (Figure 2; Suppl. Table 2); and to get a broader overview of their responses, although important traits such as root damage or changes over time could not be considered. We included an additional paragraph into the discussion addressing the above issues.

In addition, we provide now data on Leaf Area Index from remote sensing for the sampled areas (thus including the sampled pairs) and correlate these with the drought stress experienced (see also previous remark) which shows a clear negative correlation of “Leaf Area Index” and “cumulated growth season evatranporation potential (Suppl. Figure 2). This shows the close relation between drought experienced and leaf loss of the woods surveyed over a longer period.

– To the extent that the authors cannot provide data on study trees which is more specific to drought stress, they must acknowledge that their stress classification cannot be clearly attributed to drought susceptibility of the trees, and may incorporate other (multiple) stress factors.

We now acknowledge in the discussion that we do not fully causative evidence that the observed phenotypes are caused by drought alone. While it is of course possible that all the sampling sites likewise experienced simultaneously another stress and that this unknown stressor alone or in conjunction with drought led to the observed stress phenotypes, it is highly unlikely that yet another stressor other than drought. The nature of this stressor that would have needed to increase simultaneously yet independently of the observed drought, is, however, completely unclear. The discussion was modified accordingly.

2) Beyond the concerns regarding interpretation of tree stress status, sufficient information must be provided regarding the trees studied.

– The authors should provide evidence to support their classification of trees as healthy versus stressed (such as images – even example photographs, but better systematic photographs or aerial imagery of pairs; and/or supporting measurements; field rankings of e.g. crown transparency or apparent stress based on specific criteria).

Before going into detail, please note that any misclassification of trees due to unaccounted environmental differences would tend to dissimulate systematic genomic differences between the phenotypic classes – not enhance them. The reported results are therefore conservative in the sense that potential misclassifications may have prevented the discovery of additional associated loci, but cannot have produced spurious results. We now provide a sample of photographs trees from both phenotypic classifications, Suppl. Figure 4, which is also depicted in Figure 2 E-F.

– The authors should provide at least a rough measurement of how far pairs of trees were allowed to be apart (e.g. minimum and maximum distance in m) and roughly how many other trees were allowed to stand in between (to what extent can these trees be considered neighbors). Please refer to Reviewer #3 – Public review, "major weaknesses" section.

We provide now the mean distance (5.1 m) in the text and a distance distribution between all pairs in the Supplement. From this distance distribution, it becomes obvious that the trees were literally standing next to each other and that the root systems of most pairs overlapped and thus at least partially experienced the same soil quality, rooting depth etc. (see following comment).

– The authors must either provide information on known or expected within-site variation of water availability and other relevant factors such as soil quality, possible rooting depth, and differences in biotic interactions which could have affected drought status or responses, on the scale of distance between their selected pairs; or at a minimum, they should discuss possible environmental variation on this scale based on literature such as that suggested by Reviewer #1.

Suggested literature was included in the discussion, but see above – any such micro-heterogeneity in environmental conditions among paired trees cannot have contributed to artificially increase systematic genomic differences between them.

– There is conflicting information regarding numbers of individuals sampled, and included in each group, which must be resolved. Please refer to comments from Reviewer #2 under Recommendations for authors – Major suggestions – Number of individuals used.

We want to stress that there was no conflicting information regarding numbers. We admit that the (correct) information about the sizes of different sets might have been confusing. We tried to clarify these issues in the revision.

3) There is also important information missing regarding sample and data analysis, and this must be provided, and concerns regarding non-standard methods addressed. Please see detailed comments from Reviewer #2 (Recommendations for authors – Major suggestions) regarding the following:

– Insufficient information is provided regarding the improvement of the reference genome.

The improvement of the genome is now described in a separate article:

A chromosome-level genome assembly of the European Beech (Fagus sylvatica) reveals anomalies for organelle DNA integration, repeat content and distribution of SNPs

Bagdevi Mishra, Bartosz Ulaszewski, Joanna Meger, Markus Pfenninger, Deepak K Gupta, Stefan Wötzel, Sebastian Ploch, Jaroslaw Burczyk, Marco Thines

bioRxiv 2021.03.22.436437; doi: https://doi.org/10.1101/2021.03.22.436437

– The authors used unevenly sized pools for poolGWAS and equal membership of pooled samples can also not be guaranteed based on the information provided in the Materials and methods.

The pools that were actually contrasted in poolGWAS were of identical size; their common analysis with a CMH test (the multivariate version of a Fisher’s exact test) is not affected by this. See also comments to reviewers.

Regarding the individual contribution to the pools, Gaultier et al. 2013 showed that even for substantial unequal contributions of each individual to the final pool of sequence reads, the estimation of allele frequencies is at least of the same accuracy as individual‐based analyses. We provide information that we used equal amounts of DNA from each individual (which was mainly assured by the already included information that we used a similar amount of tissue from each tree).

– Insufficient information is provided regarding how GWAS was conducted.

More detailed information on how the individual GWAS was conducted is now included in the Materials and Methods section and the complete pipeline added to Suppl. Info 2.

– It is not clear which indicators of statistical significance are reported and how these are to be interpreted.

Issue resolved, see comment to the reviewer.

– The reported very low LD may be an artifact of interpreting LD from poolGWAS data, and should be assessed using the individual genomes which were sequenced at high coverage.

Please note that the cited reference (Schlötterer et al. NRG 2014) does NOT call the method or accuracy of genome-wide (short range) LD estimates from PoolSeq data into question. On the contrary: the method is listed in a table for (recommended) PoolSeq software. It simply states that analyses requiring (long range) LD information should not rely on PoolSeq data, if haplotypes are longer than the read length. We have nevertheless added an average LD estimate from the individual resequenced data. The outcome is qualitatively identical, quantitatively, the individual reseq data suggested an even steeper drop of correlation among SNPs (Figure 3a)

– It is not clear how group membership prediction was conducted and how the very high success rate is to be interpreted (GP analysis).

See below.

– The functional interpretation of SNPs is complicated by some apparent contradictions in genome locations and gene models, some of which may be resolved by more information from the improved genome assembly.

The single apparent contradiction was probably caused by the page-break in the Table, that somehow let a line vanish unnoticed. This mishap then propagated further, which we beg to excuse. The incriminated position belongs to a different gene model. The table and all consequential errors were corrected.

– Insufficient information is provided regarding software versions and parameters used, and also reasoning for choosing some non-standard software tools.

Parameters and versions were exhaustively added. What are (non)standard software tools?

– Additionally, from Reviewer # 3, Recommendations for authors: The leaf handling protocol is unusual and it is not very precisely reported: collected fresh from field, 50 deg C in lab 30-90 min, then to salt until storage at -80. It seems like the DNA would be better preserved if leaves were harvested directly to silica gel and not heated as is commonly done (in which case storage at -80 deg C can be done but is generally not required).

Please note that we simply and faithfully reported each step of the sampling and leaf handling procedure. As we encountered no problems in terms of quantity and quality of the DNA with our low-cost, field-proof method, we do not see the need to justify the reported proceeding. What is meant by “not very precisely reported”? Even the telegram-style summary above would allow exact replication.

4) One important hypothesis suggested by this study is that beech populations locally embody sufficient variation in drought-responsive genes to adapt to increased drought severity and frequency. However, this cannot be supported as a conclusion of the study.

This is based in our undisputed inference that “trees with the highest predictive values showed no loss of heterozygosity [which] indicated that there is still adaptive potential for drought adaptation in the species”. We therefore provided evidence that genetic variation for this quantitative trait is not depleted and strengthen this argument with a citation. We agree that this is no conclusive evidence and needs to be confirmed in future studies, but it is also more than a mere hypothesis. By the way, this claim occurred only in the one-sentence-twitter-style Impact Statement, not the scientific parts of the text. We down-toned it to: “European beech harbours substantial genetic variation at genomic loci associated to drought resistance and the loci identified in this study can help to accelerate and monitor this process.”

The authors should re-state this as a hypothesis and propose work required to test this hypothesis rigorously. As noted by Reviewer 1, the candidate genes revealed often do not overlap with those reported from other studies, including studies in beech,

We have refuted this critique in our direct response to Reviewer 1. The majority of all genes with known function were already implicated in drought response in previous studies. It remained unclear whether the remaining genes had not yet been associated with drought, or whether we were just unable to find them due to the lack of annotation. Respective sentence inserted.

and as noted by Reviewer 3, the geographic limitations of the study preclude conclusions regarding to what extent the study populations embody the potential of members of the species to adapt to drought.

We addressed this issue by adding a respective paragraph in the discussion.

Furthermore, there is insufficient phenotypic data to link the identified loci to specific response curves (e.g. change in ability to tolerate or recover from a given severity of water deficit and associated stresses such as salinity and heat). The ability to predict membership in a group (healthy versus stressed) by genotype demonstrates the consistency of associations identified within this dataset, but it does not indicate that the study populations are resilient to predicted future increases in drought frequency and severity.

Such functional assays to link genomic variation with phenotypic response curves of more or less wild forest trees over their life cycle would slightly exceed our current possibilities. Nevertheless, there is enough current literature, including a paper suggested by a reviewer that demonstrate this possibility from genomic data (e.g. Capblancq, T., Fitzpatrick, M. C., Bay, R. A., Exposito-Alonso, M., and Keller, S. R. (2020). Genomic prediction of (mal) adaptation across current and future climatic landscapes. Annual Review of Ecology, Evolution, and Systematics, 51, 245-269.)

Again, rather, these results suggest hypotheses about how specific loci may, alone or in concert, confer resilience or tolerance to a range of drought stress, which could be better articulated based on knowledge in the beech system and gene annotation; and the authors could suggest follow-on studies, increasing the value of their own study as a resource.

The main aim of this study was to test whether there is heritable genetic variation for this quantitative trait, i.e. if the species is likely to persist under climate change by natural selection and/or assisted evolutionary management. To study physiological or molecular details on whether and how specific loci may contribute to drought resilience or tolerance was not our aim. We now acknowledge that our study can be used as starting point for such research. We added a respective sentence.

Reviewer #2 (Recommendations for the authors):

This is an interesting study on the genomic basis of drought resistance in an important forest tree species. However, I have concerns about some of the analyses and there are several aspects of the paper where there is a lack of detail or where further clarification is needed. I outline these points below.

Number of individuals used

I am confused about the number of individuals included in some of the analyses because there are conflicts between the details provided in different parts of the manuscript. For example, Lines 76-77 of the introduction suggest that more than 400 individuals were used for the pool-GWAS ("more than 200 neighbouring pairs of trees"), but Line 120 states that "300 beech trees" were used for the pool-seq. Lines 122-123 say that there were 100 individuals in each of the South pools, but, Table S1 indicates there are 135 trees in each South pool. On Lines 127-128 it states "All 100 individuals from the North population" were individually re-sequenced, but the Materials and Methods say 102 individuals were re-sequenced (Line 271) and Table S1 indicates that it was 95 individuals. On Lines 238-239 it says 53 trees were sampled for "set 2", but Table S1 indicates there are 80 individuals in set 2. There are various other similar examples, which altogether create a confusing picture. The Authors need to double-check all such details and correct them as necessary.

See above.

Reference genome

More details about the methods used for improving the reference genome are needed. For example, please provide more details of how the Hi-C reads were generated (i.e. source of DNA, extraction method, sequencing coverage) and add details of the version number and parameter settings used for the allhic software. Also, was the improved assembly reannotated? If so, how? If the original annotation was transferred to the new assembly, this should also be explained.

See above.

GWAS

The guidance for popoolation (https://sourceforge.net/p/popoolation2/wiki/Manual/) states that, for the experimental design, the most important thing is "the amount of DNA per individual in a single pool should be constant". In the Materials and Methods section that deals with DNA extraction (starting on Line 264), there is no mention of quantification of DNA being done to ensure that the amount of DNA contributed to a pool by each individual is equal. Please can the Authors provide details of how they quantified DNA from each individual prior to pooling to ensure they contributed equally?

See comment above. We used the same amount of tissue from each individual by cutting an equally sized piece out of the leafs and measuring DNA content and quantified it. Details now added.

The guidance also states that, "the number of individuals per pool should be similar". However, the Southern pools are much larger than the Northern pools (more than double the size according to the information in Table S1). Could the Authors comment on why they used such unevenly sized pools and what impact this might have had on the results of the GWAS?

It is important that the pools, which are actually contrasted (i.e. with the same region) have similar sizes. This is given in our study design contrasting pools with healthy and damaged pools of 100 individuals from the South and 50 individuals each from the North. The multivariate CMH test, on which the identification of associated SNPs is based, is not affected.

With regard to the potential impact of population stratification, which could create spurious associations in the GWAS, the Authors state on Lines 201-202 "As expected from previous studies, we found no population structure among the sampling sites". However, I cannot see anywhere an analysis of potential structure between all populations sampled for the main GWAS analysis (i.e. the pool-GWAS). Could the Authors comment on this further please?

We now include FST distributions for genome-wide 1 kb windows among all pools, showing that i) the FST level among all pools is very low (standard deviation includes 0) and ii) not much different among regions compared to within regions.

Also, in relation to the control for FDR, the method used by the Authors (the R package qvalue) does not calculate "corrected" p-values. Normally for GWAS, false discovery rates are used to set an appropriate p-value significance threshold to control for the number of false discoveries, rather than to correct p-values as with multiple comparison correction. I.e. a threshold can be set under which the number of false positives is deemed to be acceptable. In Figure 4 the Authors present "corrected" -log10 p-values and indicate a p-value significance threshold. I think the values in this plot are actually the q-value estimates that are generated by the qvalue package. I suggest the Authors present the original p-values from the association analysis, but indicate a significance threshold depending on the number of false positives they are willing to allow within the set of "significant" SNPs. If this changes which SNPs are considered to be significant, then clearly any subsequent analyses and results that are presented should account for this fact.

The reviewer is right that qvalue does not calculate corrected p-values. However, we actually used p.adjust with Benjamini-Hochberg correction. The reported values are therefore correct. Please excuse this confusion.

Furthermore, for the GWAS analysis using data from the re-sequenced individuals (mentioned in Lines 303-308), please provide full details for what was actually done, including any additional filtering of SNPs done in PLINK (e.g. LD pruning, which is mentioned in the legend for Figure 3, but not in the Materials and Methods) and which statistical test for association was performed. Moreover, there don't seem to be any details for this GWAS analysis in the Results section, although it is mentioned on Lines 326-327 of the Materials and Methods that some of the significant SNPs were from this GWAS were included in the SNP array.

Details are now included.

LD analysis

The Authors estimated LD in beech on the basis of a subset of the pool-seq data. This suggests a very short size of linkage block on average (c.120bp), which led the Authors to conclude that most of the significant SNPs they detected are "the actual causal variants" as they are unlikely to be linked to other SNPs. However, it has been suggested that the estimates of LD obtained from analysis of pool-seq data are not very accurate (Schlötterer et al., 2014; doi:10.1038/nrg3803). The finding that SNPs that are further than c. 120bp apart are likely to be unlinked is surprising to me given the length of linkage blocks that have been detected in some other tree species. As the Authors have whole genome sequence data for individuals from the Northern pools, I would recommend that they also use these data to estimate LD to see if this produces a result that is congruent with that based on the pool-seq data and supports their assertion that they have likely identified causative SNPs.

See above.

GP analysis

The details of the GP analysis need further clarification. In particular, I am confused by the reported level of accuracy of GP with the set of 70 variants. The Authors say "Linear discriminant analysis (LDA) correctly predicted the observed phenotype from the genotype in 91 of 92 cases (98.9%)". However, in the Materials and Methods it is reported that 80% of individuals were used in training and only 15 individuals were used in the test set. Can the Authors please make clear what they actually did to test the accuracy of the LDA analysis in discriminating trees with the different phenotypes? If only 15 trees were included in the test set, then accuracy should be reported based on the results of these. Also, was any replication done to see if results differ with different subsamples of 15 individuals (i.e. using multiple random samples of 15 individuals from the full set)? If not, then I think this would be worth considering.

Due to some communication failure, the method was wrongly reported. We used LDA as implemented in PAST 4.05 with the entire data set and post-hoc prediction. Clarified in the manuscript.

Alternatively, have the Authors considered also using some of the pool-seq data to train a GP model? This could be used to test the accuracy of prediction for the individually sequenced trees, which would provide a larger number of trees for testing. Further, the details of the number of individuals used for the GP analysis and numbers and percentages of individuals used for training and testing need to be double-checked (lines 342-346), as there seem to be some inaccuracies here.

In addition we used a new machine learning method, eSPa (Horenko 2020) which avoids the problems arising from overfitting of a small data set. Here, we had 85% classification success. Method, results and discussion accordingly modified.

Functional significance of SNPs

tbg-tools was used to check whether significant SNPs from the GWAS cause non-synonymous substitutions in non-coding regions. Could the Authors explain why they chose to use this tool (which is described on its GitHub page as "not in its finished version and there are surely still some bugs to be found") rather than more established software such as SnpEff (https://pcingola.github.io/SnpEff/)?

The software is now in the publication process and the manuscript is available as preprint: tbg – a new file format for genomic data Philipp Schönnenbeck, Tilman SchellSusanne Gerber, Markus Pfenninger doi: https://doi.org/10.1101/2021.03.15.435393

Also, even if the software correctly predicts the impacts of the variants, these predictions are dependant on the accuracy of the genome annotation. As noted above, there is a lack of detail about the annotation for the improved reference genome (and also no indication of the availability of the annotation file). Could the Authors comment on whether they manually checked the gene models for those cases where a significant SNP fell within a gene (i.e. results in Table 1) to see if there was any evidence for gene model errors that could impact the results?

We manually checked the gene-models. All gene-models are correct and the intron/exon borders are supported by RNA-Seq data. Please refer to the new genome publication/preprint.

For example, two variants are reported to be within gene model 7.g2350.t1, one of which causes a premature stop codon. However, the position of these variants suggests they are >15Mb apart; assuming the positions have been correctly reported, I would be surprised if this gene model is correct, especially as the best matching gene from A. thaliana is only c. 3kb.

Thanks for remarking this, there a line simply vanished in manuscript preparation. The incriminated position belongs to a different gene model. Results accordingly corrected.

Also, in the case of 12.g1695.t1, four variants are reported with a c. 60bp region of the gene, including one inducing a premature stop codon in individuals with the alternate allele. This leads the Authors to suggest that the gene "lost its function" in individuals that carry the non-reference version of the allele. However, another possible cause of such a pattern would be incorrect specification of intron/exon boundaries in the gene model. Have the Authors double-checked this to confirm that the variants are actually within a protein-coding region?

See above.

Even if no evidence of gene model errors is found, it cannot be concluded from the sequence data alone that the gene will be non-functional, as even with a premature stop codon a truncated product with some functionality could be formed. So, I suggest the Authors need to be a bit more cautious about their conclusions (e.g "is likely to have lost its function, or be only partially functional").

The association of these variants with the damaged trees implies at least that the resulting variants are detrimental to fitness, even if some functionality is retained. Changed to “is possible…”.

Details of software

For each piece of software used, please provide details of the exact version number used and any parameter settings; if default parameter settings were used, please state this. At the moment, these details are given for some of the software, but not all. For example, no version number is given for BWA, PoPoolation and qvalue; autotrim and samtools are missing details of the parameter settings, etc.

Details added.

Reviewer #3 (Recommendations for the authors):

The genetic analysis suggests important hypotheses regarding the drought resistance or resilience of European beech. However, these hypotheses are neither very precisely articulated, nor tested in this study. Some of the information necessary to better articulate the hypotheses seems to be missing, and to articulate and test them would require substantial additional work. For this reason it seems most suitable as a Tools and Resources article accompanying publication of the dataset and improved reference genome, which could be a very valuable resource as described by the authors and the reviewers.

Resulting hypotheses inserted.

The leaf handling protocol is unusual and it is not very precisely reported: collected fresh from field, 50 deg C in lab 30-90 min, then to salt until storage at -80. It seems like the DNA would be better preserved if leaves were harvested directly to silica gel and not heated (in which case storage at -80 deg C can be done but is generally not required).

See above. It worked.

https://doi.org/10.7554/eLife.65532.sa2

Article and author information

Author details

  1. Markus Pfenninger

    1. Molecular Ecology, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    2. Institute for Organismic and Molecular Evolution, Johannes Gutenberg University, Mainz, Germany
    3. LOEWE Centre for Translational Biodiversity Genomics, Frankfurt am Main, Germany
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Supervision, Funding acquisition, Investigation, Visualization, Writing - original draft, Project administration
    For correspondence
    Markus.Pfenninger@senckenberg.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1547-7245
  2. Friederike Reuss

    Molecular Ecology, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    Contribution
    Supervision, Investigation, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Angelika Kiebler

    Molecular Ecology, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    Contribution
    Investigation, Writing - review and editing
    Competing interests
    No competing interests declared
  4. Philipp Schönnenbeck

    1. Molecular Ecology, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    2. Institute of Human Genetics, University Medical Center, Johannes Gutenberg University, Mainz, Germany
    Contribution
    Software, Formal analysis, Writing - review and editing
    Competing interests
    No competing interests declared
  5. Cosima Caliendo

    1. Molecular Ecology, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    2. Institute of Human Genetics, University Medical Center, Johannes Gutenberg University, Mainz, Germany
    Contribution
    Formal analysis, Writing - review and editing
    Competing interests
    No competing interests declared
  6. Susanne Gerber

    Institute of Human Genetics, University Medical Center, Johannes Gutenberg University, Mainz, Germany
    Contribution
    Formal analysis, Writing - review and editing
    Competing interests
    No competing interests declared
  7. Berardino Cocchiararo

    1. LOEWE Centre for Translational Biodiversity Genomics, Frankfurt am Main, Germany
    2. Conservation Genetics Section, Senckenberg Research Institute and Natural History Museum Frankfurt, Gelnhausen, Germany
    Contribution
    Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  8. Sabrina Reuter

    Ecological Networks lab, Department of Biology, Technische Universität Darmstadt, Darmstadt, Germany
    Contribution
    Formal analysis, Investigation, Writing - review and editing
    Competing interests
    No competing interests declared
  9. Nico Blüthgen

    Ecological Networks lab, Department of Biology, Technische Universität Darmstadt, Darmstadt, Germany
    Contribution
    Conceptualization, Supervision, Writing - review and editing
    Competing interests
    No competing interests declared
  10. Karsten Mody

    1. Ecological Networks lab, Department of Biology, Technische Universität Darmstadt, Darmstadt, Germany
    2. Department of Applied Ecology, Hochschule Geisenheim University, Geisenheim, Germany
    Contribution
    Conceptualization, Investigation, Writing - review and editing
    Competing interests
    No competing interests declared
  11. Bagdevi Mishra

    Biological Archives, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    Contribution
    Resources, Writing - review and editing
    Competing interests
    No competing interests declared
  12. Miklós Bálint

    1. LOEWE Centre for Translational Biodiversity Genomics, Frankfurt am Main, Germany
    2. Functional Environmental Genomics, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    3. Agricultural Sciences, Nutritional Sciences, and Environmental Management, Universität Giessen, Giessen, Germany
    Contribution
    Conceptualization, Formal analysis, Writing - review and editing
    Competing interests
    No competing interests declared
  13. Marco Thines

    1. LOEWE Centre for Translational Biodiversity Genomics, Frankfurt am Main, Germany
    2. Biological Archives, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    3. Institute for Ecology, Evolution and Diversity, Johann Wolfgang Goethe-University, Frankfurt am Main, Germany
    Contribution
    Resources, Writing - review and editing
    Competing interests
    No competing interests declared
  14. Barbara Feldmeyer

    Molecular Ecology, Senckenberg Biodiversity and Climate Research Centre, Frankfurt am Main, Germany
    Contribution
    Data curation, Formal analysis, Writing - review and editing
    Competing interests
    No competing interests declared

Funding

No external funding was received for this work.

Acknowledgements

There was no external funding for this study. We thank the SBiK-F Data- and Modelling Centre (Aidin Niamir) for support in obtaining climate and remote sensing data. CC acknowledges support from the Emergent AI Center funded by the Carl-Zeiss-Stiftung.

Senior and Reviewing Editor

  1. Meredith C Schuman, University of Zurich, Switzerland

Reviewers

  1. Lucienne de Witte, Institute for Applied Plant Biology, Switzerland
  2. Meredith C Schuman, University of Zurich, Switzerland

Publication history

  1. Received: December 7, 2020
  2. Accepted: June 7, 2021
  3. Accepted Manuscript published: June 16, 2021 (version 1)
  4. Version of Record published: July 8, 2021 (version 2)

Copyright

© 2021, Pfenninger et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,265
    Page views
  • 173
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Genetics and Genomics
    2. Microbiology and Infectious Disease
    Jessamyn I Perlmutter et al.
    Research Article

    Wolbachia are the most widespread bacterial endosymbionts in animals. Within arthropods, these maternally-transmitted bacteria can selfishly hijack host reproductive processes to increase the relative fitness of their transmitting females. One such form of reproductive parasitism called male killing, or the selective killing of infected males, is recapitulated to degrees by transgenic expression of the WO-mediated killing (wmk) gene. Here, we characterize the genotype-phenotype landscape of wmk-induced male killing in D. melanogaster using transgenic expression. While phylogenetically distant wmk homologs induce no sex-ratio bias, closely-related homologs exhibit complex phenotypes spanning no death, male death, or death of all hosts. We demonstrate that alternative start codons, synonymous codons, and notably a single synonymous nucleotide in wmk can ablate killing. These findings reveal previously unrecognized features of transgenic wmk-induced killing and establish new hypotheses for the impacts of post-transcriptional processes in male killing variation. We conclude that synonymous sequence changes are not necessarily silent in nested endosymbiotic interactions with life-or-death consequences.

    1. Genetics and Genomics
    Kevin R Costello et al.
    Research Article

    Transposable elements (TEs) are mobile genetic elements that make up a large fraction of mammalian genomes. While select TEs have been co-opted in host genomes to have function, the majority of these elements are epigenetically silenced by DNA methylation in somatic cells. However, some TEs in mice, including the Intracisternal A-particle (IAP) subfamily of retrotransposons, have been shown to display interindividual variation in DNA methylation. Recent work has revealed that IAP sequence differences and strain-specific KRAB zinc finger proteins (KZFPs) may influence the methylation state of these IAPs. However, the mechanisms underlying the establishment and maintenance of interindividual variability in DNA methylation still remain unclear. Here we report that sequence content and genomic context influence the likelihood that IAPs become variably methylated. IAPs that differ from consensus IAP sequences have altered KZFP recruitment that can lead to decreased KAP1 recruitment when in proximity of constitutively expressed genes. These variably methylated loci have a high CpG density, similar to CpG islands, and can be bound by ZF-CxxC proteins, providing a potential mechanism to maintain this permissive chromatin environment and protect from DNA methylation. These observations indicate that variably methylated IAPs escape silencing through both attenuation of KZFP binding and recognition by ZF-CxxC proteins to maintain a hypomethylated state.