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
3 figures, 1 table and 2 additional files


Figure 1 with 2 supplements
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

Figure 1—figure supplement 1
Climate change dynamics.

Mean change of sampling site values along principal component analysis (PCA) axis 1 (temperature, red) and PCA axis 2 (precipitation, blue) among decades.

Figure 1—figure supplement 2
Plot of ‘relative leaf area index’ change relative to 2014 values against ‘cumulated evatransporation potential’ change during the growth season relative to 2014 values for all 1 × 1 km plots encompassing the 27 sampling sites.

The overall correlation is r = 0.695.

Figure 2 with 2 supplements
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.

Figure 2—figure supplement 1
Distribution of pairwise distances between the paired trees.
Figure 2—figure supplement 2
Exemplary pictures of damaged and healthy beech tree pairs from several sampling sites.
Figure 3 with 7 supplements
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).

Figure 3—figure supplement 1
Genome-wide linkage disequilibrium (LD) and principal component analysis on genome-wide single-nucleotide polymorphism (SNP) data.

(A) Decay of genome-wide linkage disequilibrium, measured as r² on allele frequencies gained from individual resequencing, with distance from focal SNP in base pairs. (B) Plot of the first two principal component axes of LD-pruned SNP data from individually sequenced beech individuals from the North population. Healthy trees are indicated by a green dot, damaged ones by ochre. Individuals sampled from the same site are grouped by convex hulls, limited with dotted lines.

Figure 3—figure supplement 2
Genome-wide FST distributions in 1 kb windows.

Comparisons between all pools. Please note that the standard deviation included zero in all cases.

Figure 3—figure supplement 3
Genomic similarity among individuals within and among phenotypic classes.

Permutation ANOSIM shows that the differences are not significant (p=0.76).

Figure 3—figure supplement 4
Uncorrected Manhatten-plot.

(A) Manhattan plot of uncorrected p values from Cochran–Mantel–Haenszel test and (B) corresponding QQ-plot. Single-nucleotide polymorphisms on different chromosomes are alternatingly coloured blue and black.

Figure 3—figure supplement 5
Increase of linear discriminant analysis prediction success with the number of loci involved.

Loci were added according to their decreasing contribution in the final analysis.

Figure 3—figure supplement 6
Comparison of observed heterozygosity between lower and upper half of predictive values ‘healthy’ in LDA.

The medians are not significantly different (Mann–Whitney U = 293.5, p same median = 0.72).

Figure 3—figure supplement 7
Histogram of linear discriminant analysis results.

Individuals where the predicted phenotype did not match the observed phenotype are shown in grey, and individuals with matching observed/predicted phenotype in green (healthy) or in ochre (damaged).


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
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
7331100007 .g3816.t1NoneGAL--
745048137 .g552.t1NoneGAG--
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

Additional files

Supplementary file 1

Tables with detailed information.

(A) Sampled Fagus sylvatica individuals. (B) Results from Mann–Whitney U-test on difference between damaged and healthy trees for various parameters. (C) List of genes with significant single-nucleotide polymorphisms (SNPs). (D) List of genes closest to significant SNPs. (E) List of 20 most informative SNPs as selected by the entropy-based Scalable Probabilistic Analysis (eSPA) allowing for 85% correct classification. (F) Software pipeline and commands used for PoolSeq analysis. (G) Workflow individual reseq GWAS.
Transparent reporting form

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)

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

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

  1. Markus Pfenninger
  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
Genomic basis for drought resistance in European beech forests threatened by climate change
eLife 10:e65532.