大蒜 (Allium Sativum) 具有独特的风味和药用价值，是一种重要的蔬菜作物，广泛种植于世界各地，有数千年的种植历史，然而其进化史在很大程度上仍然是未知的。大蒜虽然为二倍体物种，但核基因组复杂，具有大规模 (16.9Gb) 和大比例的重复序列 (91.3%)，这同样限制了对大蒜种群进化史的解析。但随着大蒜基因组的测序组装的完成，通过开展大蒜种质资源重测序及群体分析，构建大蒜全基因组变异图谱，系统探究大蒜种群进化历史具有重要意义。近日中国农业科学研究院麻类研究所、山东农业大学、长沙大学、山东省大蒜工程研究中心和北京诺禾致源科技股份有限公司联合美国夏威夷大学等国内外多家研究机构，在国际权威期刊 Genome Biology 杂志上发表了题为“Genomic insights into the evolutionary history and diversification of bulb traits in garlic”文章。系统解析了大蒜的种群结构，绘制了其变异图谱，并探索了其进化历史，为基于大蒜基因组学的育种研究提供了宝贵的资源。

研究意义

大蒜是一种不育的作物，可以作为蔬菜、调味品和药物，具有重要的经济价值，然而其进化史在很大程度上仍然是未知的。该研究通过开展大蒜种质资源重测序及群体分析，绘制了全面的大蒜基因组变异图谱，系统的揭示了大蒜种群进化历史，发现大蒜种群至少在十万年前开始发生分化，中国主栽的四六瓣蒜和多瓣蒜是独立驯化的两个大蒜组群。通过与转录组关联分析，发现15.0%和17.5%的基因在两个栽培群体中发生了差异表达，导致其转录组结构重塑。这项研究为基于大蒜基因组学的育种提供了宝贵的资源，并对这种营养体繁殖作物的进化史进行了全面的解析。

研究结果与结论

1.大蒜的种群结构

该研究通过GBS (Genotyping-by-sequencing) 技术对从17个国家收集的233个不同的大蒜种质进行测序分析，其中包括3份近缘种Allium ampeloprosum (图1a)。研究共鉴定到2,036,304个SNPs，平均84.3% Allium ampeloprosum 的reads与大蒜基因组一致，覆盖4.9%的大蒜基因组，覆盖深度约为15.6倍，这表明大蒜和Allium ampeloprosum 亲缘关系较近。通过neighbor joining (NJ) tree，贝叶斯聚类和主成分分析 (PCA) 发现，大蒜种质可分为四个类群: 起源组 (OG)、中国主栽的四六瓣蒜组（CG1）、多瓣蒜组（CG2）以及其余组 (G3) (图1b，c)。其中OG组与Allium ampeloprosum 亲缘关系更接近，推测OG基因组更接近大蒜野生种。

图1 大蒜种质的地理分布、种群结构和遗传多样性。a, 230个大蒜种质的地理分布。 b, PC分析。c, NJ法构建系统发育树。d, 核苷酸多样性 (π) 和种群分化 (F_ST)。e, OG和两个栽培组CG1和CG2的鳞茎性状比较

2.基因组变异图谱

通过对84份大蒜核心种质进行重测序分析探索其基因组变异模式。共获得13.51Tb的测序数据，比对大蒜参考基因组，平均比对率为99.05%，深度为9.36倍。在大蒜基因组中共发现了129,409,768个变异位点，包括120,857,927个SNPs和8,551,841个indels (≤5 bp)，基因组中平均每kb有8.0个变异位点。有3,102,210个SNPs (2.57%) 和224,019个 (2.62%) indels位于编码区，其中211,213个非同义SNPs和7,083个移码indels可能改变蛋白质氨基酸序列。该研究绘制了大蒜种质第一个大规模的全基因组变异图谱，为该重要经济作物的生物学研究和育种研究提供了宝贵的资源。

3.种群多样性和分化

为研究大蒜种群的遗传多样性，对230份大蒜种质核苷酸多样性进行评估，发现在4个类群中，遗传多样性相对较低，范围在1.79×10-5到2.88×10-5之间 (图1d)。通过遗传分化指数 FST对四个类群间的基因组多样性进行评估，发现它们之间具有较高的FST值 (>0.368；图1d)， 特别是OG和CG2之间的FST值高达0.474，与两个水稻亚种籼稻和粳稻的FST值 (0.55) 相似，说明组群间存在高度分化。

4.栽培种 CG1 和 CG2 是独立驯化的

通过分析OG，CG1和CG2的有效群体大小，发现三个大蒜类群在50万年前就已出现分化。此后，这三个种群大小呈现不同的变化趋势，其中， CG1在1.5万年出现种群瓶颈， CG2在4万年前发生种群扩张，紧接着种群持续减少，在约7000年前达到一个瓶颈期 (图2a)。由于人类仅有约一万年的作物驯化历史，考古发现最早在公元前2600年才有大蒜相关记载。据此，研究人员推测CG1和CG2在驯化之前就已发生分化，即CG1和CG2为独立驯化。D-statistic分析发现三个种群间并无基因流同样验证这一推测 (图2b)。研究者还观察到CG1和CG2的基因组表现出明显的分化，其中30.1% (总长度为4.89Gb) 的FST值超过0.5 (图2c)。驯化选择信号分析发现，在57,561个大蒜基因中只有105个在两个栽培组CG1和CG2中受到共同选择 (图2d)，该比例甚至低于两个独立驯化的甜瓜种melo 和agrestis。这些结果进一步证实了CG1和CG2是独立驯化的。

图2 CG1和CG2独立驯化证据。a, OG、CG1和CG2的有效种群数量。b, Patterson's D检验大蒜种间基因流。c, CG1和CG2的F_ST分布。d, 基因组中putative selective sweeps 的分布。

5.CG1和CG2的转录组存在显著的差异

该研究同时对OG，CG1和CG2三个组群的81份种质进行了RNA-sequencing，并分析群体中大蒜基因的转录本丰度。香农-维纳指数（Shannon-Wiener index）分析发现大蒜基因在三个类群间的表达多样性呈现明显差异 (图3a)。与OG相比，15.0%的大蒜基因  (8,621个基因) 在CG1中存在差异表达，而在CG2中这一比例达到17.5% (100,96个基因；图3b)。此外，研究者发现16.3%的大蒜基因 (9400个基因) 在CG1和CG2之间表现出差异性表达，包括3795个在CG1或CG2中特异性表达的基因。将转录组数据与鳞茎性状进行进一步关联分析发现，这些差异表达的基因在与鳞茎性状相关的功能组中显著富集 (P < 0.001；图3c)。这一结果表明，转录组的差异可能与CG1和CG2中鳞茎性状的差异有关。综上所述，经过分化和各自的进化，两个栽培大蒜群体的转录组结构表现出明显的差异，并且这些差异与CG1和CG2中鳞茎性状的差异有关。

图3 大蒜基因的表达演变。a, 三个大蒜组的基因表达特征。b,  OG、CG1、CG2组间差异表达基因数量。c, 差异表达基因的富集分析。d, Asa2G02184.1和Asa4G00255.1在三个大蒜组中的表达存在明显差异。e, Asa2G02184.1和Asa4G00255.1的转录本丰度与蒜瓣数量性状的关联分析。

6. 有害突变在CG1和CG2中经历了单独的积累和选择清除

对84个大蒜品种的突变负荷进行评估，发现在487,215个非同义突变中，28.64% (139，518) 为有害变异。探究任意两个种质间有害突变的重复率，发现个体间有害突变重复率变异幅度很大（2%-35%），但三组之间的个体有害突变重复率均很小（<9%）(图4)。这些结果表明，三个大蒜组群有害突变是单独积累形成的。大蒜是不育的，一般通过鳞芽进行无性繁殖，这给通过重组消除其有害突变带来了挑战。研究人员比较了驯化选择区间和非选择区间的有害突变，发现CG1和CG2的受选择区域内有害突变率仅为非选择区域的1.57%和1.34%，并且在CG1的选择性区域有21.02%的有害突变被移除，这高于CG1全基因组的平均比例 (10. 32%；图4)，在CG2中也发现了类似的结果。这一结果说明驯化选择有效清除了大蒜的有害突变负荷。

图4 大蒜基因组中的有害突变负荷。a,  84个大蒜品种中突变负荷的高度多样性。b,  OG、CG1和CG2之间有害突变数量的比较。c,  全基因组和选择性区域去除有害突变的百分比。d, 选择性区域和非选择性区域之间的突变负荷比较 。

7.CG1和CG2的差异化选择产生鳞茎性状的多样性

鳞芽数和鳞芽重量决定了大蒜鳞茎产量，它们在三个组群间存在显著差异。其中鳞芽数，与OG相比，在CG1中明显下降，但在CG2中略有增加 (图5a)，表明在两个栽培组中性状演化的潜在差异。为了确定与鳞芽数量相关的基因，对230个大蒜材料进行全基因组关联分析，并确定了22个显著关联信号。鉴定到一个TB1 同源的Asa2G00245.1 和一个AXR 同源的Asa2G00503.1 为鳞芽数性状相关的候选基因，并且其在CG1中受到强烈选择 (图5c)。变异分析发现Asa2G00245.1 基因中的一个突变产生了一个终止密码子，携带该突变位点的种质会产生较少的鳞芽 (图5d，e)。此外，73.5%的CG1在Asa2G00245.1的启动子区域带有等位突变，导致该基因表达量降低并产生较少的鳞芽 (图5d-f)。Asa2G00503.1 的表达分析显示，在CG1中的表达量较高，其表达与大蒜群体的鳞芽数量呈负相关，表明Asa2G00503.1 是鳞芽数量的候选抑制因子。Asa2G00503.1 的三个单倍型在84份大蒜种质中均被鉴定出来，但单倍型1和单倍型2只出现在OG和CG1中，CG2只携带单倍型3 (图5h)；相应地，携带单倍型2的种质生产的鳞芽数量较少 (图5i)。

图5 CG组群中鳞芽数量性状的选择印记。a, πOG/πCG1值最高的5%的基因组区域。b, 大蒜TB1同源基因的系统发育树。c, 2号染色体核苷酸多样性 (πOG/πCG1) 和F_ST值的分布 。 d, 三个大蒜组中Asa2G00245.1的重要变异及其等位基因频率。e, 具有Asa2G000245.1不同等位基因的大蒜材料中的鳞芽数量。f, 三个大蒜组中Asa2G000245.1的表达水平。g, 2号染色体160-170Mb区域曼哈顿图。h, Asa2G00253.1 基因三种单倍型的频率。 i, Asa2G00253.1不同单倍型鳞芽数量。

鳞芽重量是鳞茎产量的另一个关键组成部分。为了研究与鳞芽重量驯化有关的选择特征，研究者对鳞芽膨大前后的鳞茎部位进行转录组分析，发现4658个基因在鳞茎膨大期发生了差异表达。多组学关联分析选着差异表达的Asa6G06199.1（Flowering locus T 同源基因），DELLA基因Asa2G00237.1和GA20氧化酶基因 Asa4G02474.1 为鳞芽重量性状候选基因。同样，这些基因在不同的组群中经历了不同的选择 (图6)。值得注意的是，研究人员发现在CG1和CG2中分别有125个和248个基因经过选择，在不同发育阶段的鳞茎中呈现出不同的表达，表明它们可能参与了鳞茎的发育。在这些差异表达的基因中，只有8个基因在CG1和CG2中受到共同的选择，多数鳞茎性状候选基因要么仅在CG1中受到选择，要么仅在CG2中具有选择信号，这些结果表明，鳞茎性状在两个组群间的进化选择是独立进行，这进一步支持了CG1和CG2是独立驯化的。

图6 CG2组群中鳞芽重量的选择印记。 a,  πOG/πCG2值前5%的基因组区域。b和d, 分别为鳞芽重量候选基因Asa06G6199.1 和Asa04G02474.1 的单倍型。 c和e, Asa06G6199.1和Asa04G02474.1 在三个大蒜组中的表达水平

该研究由国内外多家单位研究人员共同协作完成。山东农业大学李宁阳教授，孙秀东副教授，已毕业博士生高松（扬州大学），中国农业科学院麻类研究所博士生张雪钰，朱四元研究员，程毅博士和北京诺禾致源科技股份有限公司刘孟为本文第一作者。中国农业科学研究院麻类研究所、扬州大学园艺与植物保护学院特聘教授刘头明，中国农业科学研究院麻类研究所粟建光研究员，西北农林科技大学程智慧教授和北京诺禾致源科技股份有限公司研发与合作中心首席科学家田仕林研究员为论文通讯作者。 早在2020年刘头明研究团队就携手北京诺禾致源科技股份有限公司破译了大蒜栽培种二水早的基因组，组装了一个16.3 Gb的染色体水平大蒜基因组，使大蒜成为第一个完成基因组测序的葱属植物，为基于大蒜的全基因组的研究奠定了基础。在该研究在此基础上，本研究再次开展大蒜种质资源重测序及群体分析，系统构建大蒜基因组变异图谱、揭示大蒜种群进化历史，为葱属植物的生物学和育种研究提供宝贵的遗传学资源。