基于GBS技术的233份大蒜种质资源群体进化分析

管梦娇，谢亚楠，顾启玉，孙秀东，高荣广

（山东农业大学园艺科学与工程学院，山东 泰安 271018）

大蒜（Allium sativum L.）为百合科（Liliacaea）葱属（Allium）两年生草本植物，其幼苗（蒜苗）、花茎（蒜薹）和鳞茎（蒜头）甚至青蒜和蒜黄都能作为蔬菜来食用，又因为大蒜拥有独特的香气，也是人们日常生活所必备的重要香辛料之一［1］。大蒜富含大蒜素、蛋白质、糖类、抗坏血酸和多种还原性硫化物等物质，具有非常高的烹饪价值、功能性价值和医学价值［2，3］。大蒜原产自西亚地区的高原地带，在中国已有2 000多年的种植历史，我国已成为世界上大蒜种植面积和产量最多的国家。

大蒜是我国传统的出口创汇蔬菜产品，属于我国的名特蔬菜种类。由于我国大蒜种质资源保护工作起步晚、收集保存资源不足、遗传多样性及亲缘关系缺乏系统研究，再加上环境和气候变化使大蒜不断发生变异及各地间频繁引种，导致大蒜种质名称混乱，种系不清，同名异物或异名同物的情况非常严重。因此，亟需对我国大蒜种质的遗传进化特性进行鉴定，构建种质资源库，为探索种质创新方法、研究配套关键技术以促进大蒜的丰产丰收奠定基础。

简化基因组测序（GBS）技术是指通过测序进行基因分型，是通过选取合适的限制性内切酶结合高通量群体测序来构建SNP分子标记的方法，可应用于超高密度遗传图谱构建［4，5］、群体GWAS分析［6，7］、群体遗传分析［8］等领域。GBS避免了传统方法的一些不足，虽不能达到覆盖全基因组深度的程度，但由于其省时、价廉、有效等优点而得到大家的青睐［9］。GBS的发展有助于发现植物性状及其相关变异来源，对今后品种选育、分子育种都有一定的指导意义。

群体进化分析是对同一物种的不同亚群或不同地理分布的品种进行全基因组测序或基于酶切的简化基因组测序，通过与参考序列比对，获得大量高准确性的SNP、InDel等变异信息，研究生物群体中的基因频率和基因型频率，探讨群体的遗传结构、遗传平衡及影响群体遗传平衡的因素，从分子层面揭示该物种的进化机制、环境适应性等系列问题［10］。现已有黄瓜基因组的驯化和多样性研究［11］、黑脉金斑蝶迁徙和色变［12］、番茄的选育史［13］、大非洲稻的驯化［14］、豆的驯化［15］等模式生物的群体遗传进化研究报道。

目前，已有多人使用不同的分子标记技术对大蒜进行了多样性分析，如陈昕等利用RAPD和ISSR两种分子标记技术对中国10个不同地区的大蒜进行种质资源多样性研究［16］；韩曙等利用RAPD标记对18个大蒜品种进行遗传多样性和亲缘关系分析［17］；陈书霞等利用SSR技术对40个大蒜品种进行聚类分析、主成分分析和遗传多样性评价［18］；王海平等用AFLP、SSR和InDel三种分子标记对212份大蒜资源进行分析［19］；孙亚丽等利用SSR分子标记技术对55份大蒜品种进行多样性分析［20］。

本研究在前人研究的基础上，应用GBS技术获得了高质量的SNP位点信息，进而构建系统进化树，并进行主成分分析，以期为大蒜种质资源的收集和保护提供依据，发掘优良种质、优良基因，为日后解析鳞茎等性状的遗传基础、基因克隆奠定基础，为大蒜性状分子改良提供资源。

1 材料与方法

1.1 试验材料

试验于2018年9月至2019年6月在山东农业大学园艺作物生物学实验室进行。所用233份大蒜种质材料及来源见表1，来自于不同时期的资源收集活动。其中，170份种质材料收集于国内，59份来源于国外，剩余4份来源不详。

表1 试验用大蒜种质材料及来源

表1（续）

表1（续）

1.2 群体进化分析

1.2.1 DNA提取和GBS文库构建 DNA提取：2018年12月采集所有大蒜资源的嫩叶，利用植物基因组DNA提取试剂盒（Magen Hi Pure Plant DNA Mini Kit，Guangzhou，China）提取DNA。将提取的样品DNA送诺禾致源进行DNA质检、建库和测序，使用Qubit®2.0荧光测定计（Invitrogen，Carlsbad，USA）检测核酸浓度，使用Nanodrop检测DNA的纯度（OD260／280），使用Qubit对DNA浓度进行精确定量。选用高质量的DNA用于GBS文库构建和测序。

GBS文库构建：为了得到适合的marker密度，首先用限制性内切酶对0.1～1.0μg基因组DNA进行酶切；将P1和P2 Adapter接头加在酶切得到的片段两端后，再对样品进行扩增，分别对含有P1和P2接头的tag序列进行PCR两端扩增，DNA片段pooling，然后混合样品，电泳回收需要区间的DNA，选择所需的片段进行GBS文库构建。利用IlluminaHiSeq测序平台，进行双末端（Paired-End）150测序。

1.2.2 SNP检测与注释 用碱基识别（base calling）测序得到的原始图像数据文件，再经过分析转化为原始测序序列（sequenced reads）。分析前需将sequenced reads中可能包含的接头信息、低质量碱基以及未测出的碱基（以N表示）等干扰信息去掉，最终得到的数据即为有效数据（clean data或clean reads）。原始数据过滤方法如下：首先过滤掉含有接头序列的序列；然后检查是否有单端测序序列中未测出碱基的含量超过该条序列长度10%的现象，如果有，也需要去除此paired reads；最后，若单端测序序列中含有的低质量（≤5）碱基数超过该条序列长度50%时，也需要去除此paired reads。经过严格过滤的测序数据就是高质量的有效数据。将高质量、有效的测序数据通过BWA软件（参数：mem-t4-k32-M）比对到参考基因组（http：／／gigadb.org／dataset／view／id／100209／File_page／1）。

进行群体SNP检测，使用的是SAMTOOLS等软件。利用贝叶斯模型来检测待测群体中的SNP位点。再利用ANNOVAR软件对得到的高质量SNP进行群体SNP注释，并对由多个基因组检测出的基因变异进行功能注释。

1.2.3 群体进化分析 SNP检测之后，得到的个体SNP可以用于计算种群之间的距离（p-距离）。两个体i和j之间的p-距离通过如下公式计算：

公式中dij表示两个个体为不同基因型的赋值表，L为高质量SNPs区域长度。若在位置1的等位基因为A／C，那么：

运 用 TreeBest（http：／／treesoft.sourceforge.net／treebest.shtml）软件计算距离矩阵，以此为基础，通过邻接法（neighbor-joining method）构建系统进化树。引导值（bootstrap values）经过1 000次计算获得。

1.2.4 主成分分析 主成分分析（PCA）是一种纯数学的运算方法，可以将多个相关变量经过线性转换选出较少个数的重要变量。PCA方法被广泛运用于多学科，在遗传学当中，主要用于聚类分析，它是基于个体基因组SNP差异程度，按照不同性状特征将个体按主成分聚类成不同的亚群，同时用于和其它方法做相互验证。本试验通过GCTA（http：／／cnsgenomics.com／software／gcta／pca.html）软件计算特征向量以及特征值，并利用R软件绘制PCA分布图。

2 结果与分析

2.1 测序质量

2.1.1 酶切效率 选用MseⅠ+Eco RⅠ限制性内切酶组合进行酶切，统计clean reads两端为MseⅠ捕获的reads数、酶捕获率（捕获的reads数占clean reads数的比率），用来评估酶切效率。酶切效率平均在99.0%左右，说明MseⅠ+Eco RⅠ限制性内切酶的组合酶切效率很高。

2.1.2 测序数据统计与质量评估 233个大蒜样本总的测序数据量为1 486.6871 Gb，去除低质量序列后剩下的高质量有效数量为1 486.5896 Gb，平均每个样本6.3802 Gb。测序质量高（Q20≥93.27%、Q30≥83.50%），GC分布正常，说明233个大蒜样本都没有被污染，此次建库测序成功。

2.1.3 测序深度及覆盖度 平均测序深度为7.68X～17.11X，1X覆盖度（至少有一个碱基的覆盖）在4.21%以上。各个样本的比对结果显示，它们与参考基因组的相似度达到重测序分析的要求，同时又有非常不错的覆盖深度和覆盖度。

2.1.4 SNP的检测 经SAMTOOLS软件检测，共获得了67 817 466个SNPs位点。用个体深度DP＝4、缺失率Miss＝0.2、最小等位基因频率MAF＝0.01等为条件进行过滤，最后获得了高质量的SNPs位点2 036 116个（表2），可用于后续分析。

表2 SNPs统计信息

2.2 群体进化树分析

系统进化树是用来表示群体间进化关系、描述群体间进化顺序的分支树，可以根据群体物理或遗传学特征等方面的共同点或差异推断出其亲缘关系远近，即分析群体中个体间由于共同祖先而产生的相互关系。基于2 036 116个SNPs标记，构建系统进化树（图1），可以看出，233份大蒜自然群体被划分为8个亚群，第一亚群为蓝色部分，共包含10份材料，分别来源于中国的山东、四川、贵州、广东以及埃及、美国、泰国和加拿大等；第二亚群为粉色部分，两份分别来自于中国云南和中国台湾，另一份来自欧洲；第三亚群为黄色部分，7份来源于加拿大，2份来源于乌克兰；第四亚群为橙色部分，12份材料除1份来源于波斯外，其余全部来源于加拿大；第五亚群为棕色部分，18份材料中有6份来源于加拿大，其余12份则来源于中国的新疆、甘肃、青海、山西等地；第六亚群为红色部分，56份材料中，有2份来源于日本，其余54份材料则全部来源于中国，主要来源于山东、河北、新疆、四川等地；第七亚群为绿色部分，12份材料中，3份来源于中国山东和陕西，其余9份则来源于美国、加拿大和澳洲；第八亚群为淡紫色部分，113份材料中，有4份来源不详，其余则大部分来源于中国山东、新疆、河南、湖北等地。

2.3 主成分分析

本试验利用经过滤得到的2 036 116个SNPs进行PCA分析，在忽略错配数据和高于2个等位基因位点的条件下，选取贡献率最大的前两维数据绘制散点图，见图2，其中group1—group 8与进化树依据颜色相对应。可知，PCA将233份大蒜群体分为4大类，将系统进化树划分的8个亚群用不同颜色呈现在PCA分析图中，结果与进化树结果大体一致。类群①对应群体进化树中的蓝色部分，大部分材料来自于中国；类群②对应群体进化树中粉、黄、橙、棕、红色部分，其中约有五分之三的材料来自于中国，另有约五分之一的材料来自于加拿大；类群③对应群体进化树中的绿色部分，约有半数的材料主要来自于美国，其余材料则分别来自于美国周边地区以及中国；类群④对应群体进化树中的紫色部分，几乎全部材料都来自于中国。

图1 233份大蒜种质群体进化树

图2 233份大蒜种质的主成分分析结果

3 讨论

本试验得到高质量的有效数据量为1 486.5896 Gb，平均每个样本6.3802 Gb，与参考基因组（大小为16 557 142 287 bp）相比，群体样本平均比对率为82.94% ～99.23%，说明测序数据与参考基因组的相似程度较高。利用GBS简化基因组测序技术，选用MseⅠ+Eco RⅠ酶切组合，测序深度为7.68X～17.11X，用个体深度DP＝4、缺失率Miss＝0.2、最小等位基因频率MAF＝0.01等为条件进行过滤后，共获得高质量SNPs位点2 036 116个。

233份大蒜资源的系统进化树和PCA分析是基于SNP构建的，可以从中了解到大蒜不同种质之间的遗传距离，进而了解他们之间的亲缘关系远近。在进化树中，来源于同一地区的大蒜大都被聚集在一起，如第六亚群全部来自于国外，且几乎都来自于加拿大，第七亚群大部分来自于北美洲，第五亚群主要来自于加拿大和中国西北部，第六亚群几乎全部来自于中国，这与陈书霞［18］、王海平［19］和孙亚丽［20］等的试验结果相同。不同地域来源的种质也可能被聚在一起，如第一亚群的9个种质来自三个大洲、五个国家、八个地区；第二亚群的3个种质分别来自于中国台湾、中国云南和南欧，虽然三者地理位置相距较远，但均拥有温暖湿润的环境条件，说明相似的地理环境可能使大蒜的遗传多样性趋于相同；第六亚群中仅有1／28的种质来自日本，其余均来自于中国，可能由于两地在地理位置上距离较近、气候相似；第三亚群大蒜种质来自加拿大和乌克兰，这与孙亚丽等［20］的研究结果相同。

4 结论

本研究利用GBS技术对233份大蒜种质资源进行群体进化分析，共得到高质量的有效数据1 486.5896 Gb，平均每个样本6.3802 Gb；共获得高质量的SNPs位点2 036 116个，较均匀地分布在8条染色体上。基于获得的SNPs标记进行PCA与系统进化树分析，PCA将大蒜群体划分为4大类，系统进化树将大蒜群体划分为8个亚群，将8个亚群用不同颜色在PCA图中呈现，两者结果基本一致，表明大蒜种质间的亲缘关系主要受地理起源影响，另外还与各地大蒜种质资源交流频繁、气候环境等有一定关系。本研究结果可从分子水平上为大蒜亲缘关系分析和系统分类等提供辅助作用，并为大蒜栽培和育种提供科学依据。