研究：利用古代玉米基因组推断美国玉米群体演化历史（Science 2017）

Genomic estimation of complex traits reveals ancient maize adaptation to temperate North America，Kelly Swarts ，Rafal M. Gutaker，Detlef Weigel，Edward S. Buckler，Hernán A. Burbano，Science，4 Aug 2017 ，DOI: 10.1126/science.aam9425

本文的特殊材料是 15 个古代玉米样本的基因组数据，从宏观和微观两个方面补充了美国玉米群体可能的演化历史。玉米的遗传多样性极为丰富，存在大量由转座子导致的 indel 与重复区间，所以多地域品种的低覆盖度基因组测序数据比对到参考基因组 B73 上时，必然会伴随大量的缺失。宏观层面，作者从 基因组预测表型、Fst、群体结构分析 三个角度分析了不同品种基因组间的整体相似性。微观层面，作者绘制了几个与周期相关 基因 的 物种树，即分析了不同品种相同基因间的相似性，通过物种聚类的结果进一步证明作者的结论。

由于 样本数量较少（15），可能存在较为严重的 采样偏差。另一方面，因为除了 DNA 数据外无法获取古代样本的其他任何数据，所以许多分析需要 基于现有品种 来完成，无法避免群体结构引起的误差，类似于用 A 群体中发现的规律推广到 B 群体中。尽管存在上述误差，这篇文章的研究结论也是具有重要意义的。

结论

1. 玉米从墨西哥引入美国西南部后，在当地经历了驯化，驯化完成后才在美国其他地区普遍种植。
2. 存在大量对新环境适应的突变来自于已有的常态突变。

背景及材料

图 1A 中展示了前人研究的玉米迁移演化的历史与有效基温（Growing degree days，GDD）的变化，图中黑体数字表示了当地玉米产业建立的大致时间（距今，before the present，B.P.），灰色线条表示国界。图 2 是美国及墨西哥的地理情况，可以发现，图 1A 中有效基温的高低受到 纬度 及 海拔 的影响，纬度和海拔的升高都会降低有效基温。图 1B 中 GDD DTS 表示现代地方品种的玉米生长到吐丝期所需的基温，统计显示随着作物向温带转移，所需基温逐渐下降。

玉米从墨西哥中部的低纬度低海拔高基温处，向南北两个方向传播。在距今约 4000 年时，被引入美国西南部。本研究收集的 21 个考古玉米来自 Turkey Pen Shelter（TPS），位于犹他州东南部科罗拉多高原（温带），图 1A 中红色五角星处，海拔约 1500 米，距今约 2000 年前 有玉米产业建立的证据。图 1C 是样本在 TPS 被掩埋的地质截面图，左侧纵坐标是掩埋深度，右侧数字为样本编号，红色表示内源性玉米 DNA 少于 20%（剔除），蓝色为保留的样本，总计 15 个。放射性碳素断代法确定 15 个样本距今约 1850-1890 年前。

图片来自 Google 地图

步骤

1. 测序与 Call SNP。作者测定了 15 个材料的 DNA 并使用 AID panel 进行了 SNP calling。SNP calling 是通过已有的 Hapmap3.21 完成的，Hapmap3.21 中有 8100w SNP。平均 SNP 覆盖率为 0.55 ± 0.16（表 S1）。
2. 预测表型，分析古代品种与现代品种的表型差异。利用已有的线性无偏模型预测 TPS 15 个古代样本及 80 个现代西南地方品种的开花期（Days to flowering）表型。3 个环境（Ames、Blanding） 2 个表型（DTS，days to silking；DTA，days to anthesis）。由于无法验证预测 TPS 表型的准确性，作者通过比较 80 个现代西南地方品种表型的预测精度来间接证明。
3. 驯化区间比较，分析美国玉米的驯化历程。作者计算了墨西哥中部低海拔地区（高基温）玉米品种与加拿大魁北克地区（低基温）Northern Flint 玉米品种间的 $F_{st}$（fixation index）。选择 Top10k 的 SNP，并计算这些 SNP 在 TPS 与温带、热带地方品种之间的 $F_{st}$。
4. 群体结构分析，进一步证明美国玉米的驯化历程。

PS：AID（Ames Inbred Diversity）panel 是指储存在爱荷华州埃姆斯（Iowa Ames）的由 2648 个收集自全世界的自交系所构成的基因组数据库。

结果

表型预测

从预测结果来看（Fig 2A），TPS 的表型预测值与相同地区的现代种（US.Puebloan）的表型预测值之间，没有显著性差异。因为表型基于基因值，所以表型上的差异不大，就意味着基因型上总体差异也不大。除了开花期外，作者还预测了株高，也与当地的现代种表型相似。作者认为，这种相似意味着此时品种对环境的适应基本完成。

预测误差所致的可能性不大。如线性模型中考虑的 SNP 较多，所以每个 SNP 难以赋予过多的权重，最终表型值的基数取决于大多数与群体结构相关的 SNP 。因为训练群体中包含热带，所以如果与温带表型相似意味着与温带之间总体变异差异较小。虽然存在关键性突变因改良而被排除的情况，即预测是基于当前单倍型预测古代表型，类似与用 A 群体的标记预测 B 群体一样，存在一定误差。但预测结果并没有偏向于其他群体，说明与其他群体间的差异还是较大。驯化基本完成的结论是比较适用的。

Fst

墨西哥中部低海拔地区玉米与加拿大魁北克地区玉米之间 固定系数（fixation index，$F_{st}$） Top10k 的 SNP 被视为代表了玉米从热带到温带的相关驯化基因。 作者比较了 Top10k SNP 在 TPS 与墨西哥、南美、美国西南部地方品种之间的 $F_{st}$。结果显示（Fig 2C），TPS 与墨西哥、南美地方品种之间 $F_{st}$ 相较于背景（随机选取 10k SNP 计算 $F_{st}$）存在显著差异（<0.001），TPS 与大多数美国西南部地方品种之间无显著差异。这说明，玉米从墨西哥引入美国西南部后，在 当地 经历了 驯化，驯化 完成后 才在美国其他地区 普遍种植。

TPS 适应低基温环境的基因型可能是新突变，也可能是常态突变，作者尝试比较 TPS 与高基温野生近缘种（Parviglumis）和低基温野生近缘种（Mexicana）间的序列差异，进而判断突变的来源。大刍草被认为是玉米未驯化的野生近缘种，两种大刍草中 Parviglumis 生长在墨西哥中部平原，Mexicana 生长在墨西哥北部高原，相较于 Parviglumis 具有更高的纬度和海拔，基温较低。作者比较了 TPS 与 Parviglumis、Mexicana 之间 Top10k SNP 和背景 SNP 的 分离 情况（proportion segregating）。结果显示（Fig S8）：① 无论是 Top10k SNP 还是随机挑选的背景 SNP，TPS 与 Mexicana 的分离比例都小于 Parviglumis；② Top10k SNP 间的差异显著大于随机 SNP；③ Mexicana 中 Top10k SNP 与随机 SNP 间分离比例的差异小于 Parviglumis（18.11 < 20.66），说明温热环境相关的 SNP 差异缩小速度快于背景，说明 TPS 与 Mexicana 部分共有突变因自然选择而固定，即存在 大量 对新环境适应的基因型 来自于 已有的常态突变。如果适应的突变全来自新突变，则红线与分布之间的差异应该基本不变。图中红线为 Top10k SNP 的分离比例。

作者绘制了 ZmCCT、ZCN8、Rap2.7 三个已经验证与开花相关的大效应 基因 的 物种树（Fig S9）。除了 ZCN8 外，剩下两个基因的物种树中大刍草群体均 没有 被聚类成为一个 单系群，而是 分散 在整个树中。如果突变是新产生的，新突变固定后群体内每个个体的单倍型都来源于相同的单倍型，会与旧突变的单倍型差异较大。导致携带新突变的现代品种被分在一个大枝下，携带旧突变的大刍草被分在另一个大枝下，而非分散分布结构。所以作者认为物种树的结构也证明了 “存在大量对新环境适应的基因型来自于已有的常态突变” 的结论。

群体结构

作者使用 1316 个样本（大刍草 + 地方品种）的全基因组 SNP 数据，分析了个体之间的 群体结构 关系（STRUCTURE 软件）并绘制了 物种树（TreeMix 软件）。从群体结构的分析结果可以发现，TPS 及美国西南部的地方品种的群体结构与墨西哥、南美等地方品种（无论高低海拔）基本没有重叠，作者认为这进一步证明了美国西南部地方品种在 TPS 发生了独立的驯化，产生了较大规模瓶颈，进而与祖先群体出现了较大的差异。

附录

Corn Belt（玉米带）

Corn 在美国是 Maize 的常用词，即常指玉米而非谷物。玉米带（Corn Belt）是美国 中西部 的一个地区，土地平坦，土壤肥沃，自 1850 年代以来一直主导着美国的玉米生产。

图片源自 https://en.wikipedia.org/wiki/Corn_Belt

GDD

GDD，Growing degree days，有效积温，生长度日，是指在实际环境条件下，完成某一生育阶段所累积的有效积温值，它是基于温度的一个指数，代表着植物生长期积累的热量。公式中 $n$ 为生长天数，$T_{max}$ 为当天的最高温度，$T_{min}$ 为当天的最低温度，$T_b$ 为作物发育基点温度。
$$GDD=\sum _1^n[(T_{max}+T_{min})/2–T_b]$$