Phenotypic drought stress prediction of European beech (Fagus sylvatica) by genomic prediction and remote sensing

Current climate change species response models usually do not include evolution. We integrated remote sensing with population genomics to improve phenotypic response prediction to drought stress in the key forest tree species European beech (Fagus sylvatica L.). We used whole-genome sequencing of pooled DNA from natural stands along an ecological gradient from humid-cold to warm-dry climate. We phenotyped stands for leaf area index (LAI) and moisture stress index (MSI) for the period 2016–2022. We predicted this data with matching meteorological data and a newly developed genomic population prediction score in a Generalised Linear Model. Model selection showed that the addition of genomic prediction decisively increased the explanatory power. We then predicted the response of beech to future climate change under evolutionary adaptation scenarios. A moderate climate change scenario would allow persistence of adapted beech forests, but not worst-case scenarios. Our approach can thus guide mitigation measures, such as allowing natural selection or proactive evolutionary management.

Methods
The sites selected for this study were in many regards typical for German beech forests. They comprised a large ecological amplitude with regard to long-term conditions. Within the range from humid-cold to warm-dry climatic conditions, important demographic and genetic processes in beech forests are taking place. Therefore, the aim of selection and stratification was to cover the climatic spectrum of forested areas in central and southern Germany to a large extent. For the characterisation of the climatic niche, we used the climatic marginality towards the rear edge, i.e. the dry and warm border of species distributions.

Construction of population pools, sequencing, population structure and genetic diversity

From each tree, 1-2 buds or 3-4 leaf discs of 0.5 mm diameter (approx. 50 mg of fresh plant material) were dried in Silicagel prior to homogenization. DNA was extracted using an in-house protocol. We constructed DNA pools per population using the same DNA quantity per individual beech. DNA concentration was measured using a Quantus fluorometer (Promega). Library preparation and 150bp paired-end sequencing with 450bp insert was conducted at Novogene.

Reads were trimmed using Trimmomatic v.0.39 and quality controlled with FastQC v.0.11.9. We used BWA mem v.0.7.17 to map the reads onto the newest version of the beech reference genome and Samtools v.1.10 to convert, sort and pile up the bam files. Duplicates were marked and removed with Picard v.2.20.8 (https://github.com/broadinstitute/picard). PoPoolation2 v.2.201 pipeline was used to remove indels, and calculated allele frequencies for every position and pairwise FSTs in non-overlapping 1kb windows. Genetic variation (ϴ) and nucleotide diversity (π) were estimated in the same windows with PoPoolation1 v.1.2.2. We considered only sites within a coverage range of 15-50X.

当前的气候变化物种响应模型通常未纳入进化因素。本研究将遥感技术与群体基因组学相结合，以提升对关键林木树种欧洲山毛榉（Fagus sylvatica L.）干旱胁迫表型响应的预测能力。我们对沿湿润寒冷至暖干气候生态梯度的天然林分的混合DNA进行全基因组测序。针对2016-2022年期间的林分，我们测定了其叶面积指数（LAI）与水分胁迫指数（MSI）。我们通过匹配的气象数据与新开发的基因组群体预测评分，借助广义线性模型对该数据进行预测。模型选择结果显示，加入基因组预测可显著提升模型的解释力。随后，我们在进化适应情景下预测了山毛榉对未来气候变化的响应：中度气候变化情景可使适应后的山毛榉林得以存续，但极端最坏情景下则无法实现。因此，本研究方法可为减缓措施提供指导，例如允许自然选择或开展主动进化管理。

研究方法
本研究选取的样地在诸多方面均具有德国山毛榉林的典型特征，其涵盖了针对长期气候条件的较大生态幅度。在湿润寒冷至暖干的气候条件范围内，山毛榉林正发生着重要的种群动态与遗传过程。因此，样地选择与分层的目标是在很大程度上覆盖德国中部与南部林区的气候谱。为表征气候生态位，我们采用了指向物种分布后缘（即物种分布的暖干边界）的气候边际性指标。

种群混合池构建、测序、种群结构与遗传多样性分析
从每株树木上采集1-2个芽或3-4片直径0.5mm的叶圆片（约50mg新鲜植物材料），置于硅胶中干燥后进行匀浆。采用实验室自研方案提取DNA。我们按照每个山毛榉个体的DNA量一致的原则，构建了种群水平的DNA混合池。使用Quantus荧光计（普洛麦格Promega公司）测定DNA浓度。文库构建与150bp双端测序（插入片段长度450bp）由诺禾致源（Novogene）完成。

使用Trimmomatic v0.39对测序读段进行修剪，并通过FastQC v0.11.9进行质量控制。采用BWA mem v0.7.17将读段比对至最新版本的山毛榉参考基因组，使用Samtools v1.10对BAM文件进行转换、排序与pileup分析。使用Picard v2.20.8（https://github.com/broadinstitute/picard）标记并移除重复序列。采用PoPoolation2 v2.201流程去除插入缺失序列，并计算每个位点的等位基因频率以及非重叠1kb窗口内的两两FST值。采用PoPoolation1 v1.2.2在相同窗口内估算遗传变异（ϴ）与核苷酸多样性（π）。我们仅考虑测序深度介于15-50X范围内的位点。