教程：群体演化方法分析玉米的驯化与改良

一般文章在筛选 正选择区间 时，大多 不考虑 群体的 演化历史，即不考虑 群体大小 的变化，只进行亚群之间各种群体遗传参数的对比，这可能会产生大量的假阳性。另一方面，研究一般也 不考虑 遗传信息的 迁移。所以作者希望将群体演化历史及遗传信息的迁移纳入选择区间的考量因素中，以提高结果的准确性。

值得注意的是，此方法误差的主要来源为 预测 的群体演化历史与 真实 历史间的差异。其中影响预测的因素众多，且较多因素 无法确定，如碱基突变率，有的研究预测结果为 $3.3\times 10^{-8}$，有的研究结果为 $4.8\times 10^{-8}$，而群体迁移率更是难以估计。所以，演化历史是在各种猜测的基础上推演出来，准确性难以保证。演化历史又是后续所有分析的基石，所以本方法建议用于 定性讨论，而非定量分析。

步骤

1. SNP 过滤，本文先使用 LD 进行过滤，然后使用 等间距保留，形成了 密度均一（SNP 间隔约为 2 kb）的 SNP 集合。
2. 样本集的 群体结构 分析（PCA、STRUCTURE）。
3. 样本过滤，提高样本之间的 区分度。本文去除样本集中主成份（main membership probability）占比 < 0.6 的样本。
4. 亚群遗传差异分析，使用 Fst、SNP 数量、SFS、群体特异性 SNP、MAF、LD。结果可用于进一步筛选样本及 SNP。
5. 亚群演化历史分析，使用 MSMC、∂a∂i 和 中性 SNP 推算亚群的演化历史，然后再根据历史模型生成 模拟数据，与真实数据比对，验证 演化历史的准确性。
6. 亚群遗传片段血统分析，通过 IBD 在群体中的频率间接表示驯化、改良过程中片段被选择的强度。IBD 频率越高，选择强度越大，作者使用 Top 5% 的 IBD 区间视为正选择区间。

PS：文章前后内容不一致，正文与附表也不一致，可能是文章投稿期间经历了反复的改动且改动篇幅较大。如 Table S1 中描述 MSMC 的输入样本中不包含 SS 亚群，但在文中图 2 中则包含了 SS 亚群。另一方面，文章各段落使用数据不一致，阅读的连续性较差。如推测群体演化是不考虑 SS，但讨论 IBD 时会考虑 SS。所以，更多的是学习本文的研究思路。

Li, X., Jian, Y., Xie, C. et al. Fast diffusion of domesticated maize to temperate zones. Sci Rep 7, 2077 (2017). https://doi.org/10.1038/s41598-017-02125-0

材料

本文数据使用 Fu 等人测定的 368 个玉米自交系（NSS、SS、TS、Mixed）的 RNA 序列数据，总共包含 102.6 w 个 SNP，样本集内存在清晰的群体结构。

368 群体详情参见：玉米关联群体：155、368、527 自交系群体简介（Yan Jianbing，2010、2013）

结果

1. 亚群遗传差异分析

从 Fst、SNP 数量、SFS、群体特异性 SNP、MAF、LD ，5 个角度观察 SS、NSS、TS 亚群之间差异。为保证数据具有可比性，作者每次从 3 个子群中各抽取 50 个个体，重复 10 次（图 1）。

1. Fst 分布情况显示，温带亚群内 NSS-SS 差异小，温热亚群间 NSS-TS、SS-TS 差异大。
2. SNP 数量 TS > NSS > SS，SS 中 SNP 数量最少，遗传多样性 最低。
3. rare SNP 占比 SS > NSS > TS，可能是 SS 受到的 选择压力 较强，有害 SNP 频率升高较难。
4. LD 衰减 到稳定 $r^2=0.1$ 的速度 TS > NSS > SS，平均距离分别是 170/186/224 bp 。衰减速度可能与 驯化时间、有效群体大小 有关，如热带品系驯化改良完成较早或有效群体较大，历史上重组、突变积累较多。

2. 亚群演化历史分析

使用 MSMC 和 ∂a∂i 两种方法推断玉米从热带分化到温带。

1. 使用 MSMC 与 遗传背景差异清晰 的个体，推断 3 个亚群的群体大小变化。样本选自 Hapmap 3 的 916 个玉米自交系中遗传血统（genetic ancestry）大于 99% 的 TS、NSS、SS 品系。结果显示（图 2a），群体大小先下降后快速上升，其中驯化完成时间顺序为 TS > NSS ≈ SS，群体大小为 TS > NSS > SS；共祖率（coalescence rate，图 2b）在 3000-5000 年前大幅下降，遗传信息差异大幅增加，作者推测此时约为温热分化的时间节点。

   PS：MSMC 中突变率设定为 $\mu =3\times 10^{-8}$；因为是初步判断，所以作者没有对 SNP 数据进行筛选。

2. 使用 ∂a∂i、中性、密度均一 的 SNP、MSMC 结果，推断群体的演化历史，筛选 SNP 过程如下。最终，作者因为 SS 有效群体太小而将 SS 排除，仅使用 NSS 代表温带玉米群体。如果 SS 与 NSS 合并，可能会由于两者遗传多样性差异较大，统一分析产生大量的假阳性。

   1. 筛选 fastSTRUCTURE 中主要关系概率（main membership probability）> 0.6 的样本，同时，作者使用 R 中 SNPRelate 包的 snpgdsIBS 函数计算 IBS，进而估算个体之间的亲缘关系，随机剔除了亲缘关系大于 0.5 组合中的 1 个样本，；
   2. 使用 SnpEff 鉴定出的 24.7w 同义 SNP 进行分析，以排除驯化、改良等选择造成的影响；
   3. 使用 plink 的 --thin 参数随机采样 SNP 以确保它们之间 相距 > 2 kb，筛选后有 3.36 w（13.6%）SNP 被保留。
   
   

   PS：群体大小 $Na$ 使用方程 $\theta =4\mu LNa$ 计算得到，其中 $\mu$ 是突变率，$L$ 是有效序列长度。$L$ 的总长度是包含 SNP 的外显子长度的总和，$L=3.8Mb$ ；突变率参考文献，$\mu =4.8\times 10^{-8}$，由于非同义词突变与同义突变的比值为 1.14，所以同义突变率 ${\mu }_s=2.2\times 10^{-8}$ 。

3. 作者使用 3 种群体演化模型（中性平衡、两时期、三时期）与 ∂a∂i 的结果来模拟数据，通过标准似然比检验（standard likelihood-ratio test）评估 模拟群体 与 真实群体 之间的差异，进而判断 ∂a∂i 结果的准确性。结果显示，两时期、三时期模型之间没有显著差异，都优于中性平衡模型，推测 TS 和 NSS 群体分别在 4000 和 3000 年前发生了一次群体收缩。

   PS：两时期是指群体发生了 1 次群体数量的瞬时改变，三时期是指发生了 2 次，一收一放，类似瓶颈模型。

4. 因为考古和分子数据都表明玉米是 单一驯化（single domestication），所以作者假设 NSS 是直接从 TS 驯化得到。使用 ∂a∂i 预测群体演化历史。作者构建了 4 种模型（图 3）。作者以 AIC（Akaike information criterion）值作为模型评价标准，AIC 值越低，模型拟合的越好。因为模型 3、4 拟合相似，作者选择了较为简单的模型 3（迁移率对称，图 3c、图 4a），并展示了模型 3 下 ∂a∂i 推测的各参数的最适值（图 4b）。

PS：Na（祖先种群大小）、nu1（分裂后 TS 种群大小）、nu2（分裂后 NSS 种群大小）、T1（驯化）、T2（亚群分离）、T3（亚群分化）和 m（种群之间的迁移率）。

3. 亚群遗传片段血统分析

作者计算了亚群之间 IBD 的情况，研究亚群之间的基因流动与选择区间。

1. 作者使用 fastIBD 计算 IBD，然后剔除 IBD < 1cM 的片段，因为较小的 IBD 容易受到背景 LD 影响。结果显示， TS 群体内 IBD 片段较少，可能是 TS 驯化完成较早，群体内 IBD 片段被重组和突变破坏而无法识别。
2. 作者讨论了人工选择使温热品系之间的基因 相互渗透 的情况。比如 Mixed 可能就是温带、热带相互杂交的结果。
3. 作者认为具有广泛 IBD 共享的区间与驯化、改良中的正选择有关，IBD 共享比例越高，选择强度越大。据此，作者绘制了 IBD 频率图（图 5），选择频率 Top 5% 的 IBD 区间作为 驯化、改良 的候选区间。作者使用蛋白质结构域比对（pfam HMM），举例研究了 bin 3.08、bin 8.06 中的基因功能。

附录

硬粒玉米、马齿玉米、NSS、SS

Flint，硬粒玉米，籽粒 顶部 及 四周 胚乳都是 角质 淀粉，仅 中心 近胚部分为 粉质 淀粉，外表半透明有光泽、坚硬饱满，主要用作粮食或饲料。

Dent，马齿玉米，籽粒 两侧 胚乳为 角质 淀粉，中间和顶部 为 粉质 淀粉，籽粒成熟时顶部干燥凹陷呈马齿状。马齿玉米单产高，但对肥水要求也高，是大面积生产应用的主要品种。籽粒食味较差，多作 饲料 或作为淀粉、酒精等工业的 原料。

温带玉米中 SS 与 NSS（Stiff Stalk Synthetic，SS，硬秆综合种；non-stiff stalk，NSS，非硬秆综合种）是被广泛用于杂交育种的两个群体。SS 和 NSS 起源相似，都是 Northern Flint（硬粒玉米）、Southern Dent（马齿玉米）、Tropical highland 的杂交后代。

群体演化历史推断方法简介

目前，各类推断群体演化历史的算法包括：

1. 溯祖模拟（Coalescent）模拟
2. 近似贝叶斯计算（approximate Bayesian computation，ABC）
3. 复合似然法（composite likelihood approach），如 Diffusion Approximation for Demographic Inference（∂a∂i），扩散近似群体推断方法。
4. 马尔可夫溯祖分析（Markovian coalescent analysis，MSMC）

相比于 ABC 方法，MSMC 与 ∂a∂i 更适合处理基因组数据，计算速度更快，本文使用这两种方法研究玉米群体演化历史，因为 ∂a∂i 需要提供一个预制的群体演变历史，而 MSMC 不需要，所以作者先用 MSMC 分析，MSNC 推断 TS 和 NSS 的分离时间约 3000-5000 年前。然后将 MSMC 输出的群体演变模型作为预制模型输入 ∂a∂i，推断瓶颈发生的持续时间与强度。最后挖掘可能与人工选择有关的同源片段（identical-by-descent tracts，IBD tracts），作为后续育种研究的目标。