基于GBS测序的全基因组SNP揭示大凉山玉米地方品种资源的亲缘关系与遗传分化

冷益丰,罗 樊,陈从顺,丁 鑫,蔡光泽

(1.西昌学院 农业科学学院,四川 西昌 615013; 2.攀西特色作物研究与利用四川省重点实验室,四川 西昌 615013)

玉米(ZeamaysL.)是禾本科玉蜀黍属一年生草本植物,又名苞谷、苞米棒子、玉蜀黍、珍珠米等,原产于中美洲和南美洲,是世界上重要的粮食作物,广泛分布于美国、中国、巴西和其他国家。玉米是我国种植面积最大、产量最多的第一大作物,具有重要的粮用、饲用与工业用价值。西南山地玉米区作为我国五大玉米产区之一,玉米单产仅占全国平均水平的80%。2022年6月,农业农村部印发的《“十四五”全国种植业发展规划》提出,加强粮食作物种质资源普查收集,开展种源关键核心技术攻关,加快突破性新品种培育推广。近年来,随着育种技术的发展,种质资源研究成为育种工作的重要任务。

基因分型测序(genotyping-by-sequencing, GBS)技术是在第二代测序技术基础上发展起来的一种简化基因组测序技术,该技术具有成本低廉、操作简单、方便快捷和高通量等特点[1]。GBS技术通过酶切加标签的方法,测序获得酶切位点附近的基因序列信息,进而检测大量高准确性的单核苷酸多态性(single nucleotide polymorphism, SNP)变异位点信息,是一种识别大规模单核苷酸多态性、降低基因组复杂性和基因分型最经济有效的方法。SNP作为第三代分子标记技术,在植物群体进化、种群基因交流、构建遗传图谱和全基因组关联分析(genome wide association study, GWAS)等方面有着广泛的应用[2-6]。SuperGBS是一种改进型的GBS测序技术,主要利用甲基化敏感的限制性核酸内切酶将基因组DNA进行酶切,选择性回收富含基因编码区的酶切片段,然后对其进行高通量测序,通过分析获得SNP信息[7],有或无参考基因组,都可以开发SNP标记[8]。

大凉山地处西南高原,光、热资源丰富,其境内地形分布错综复杂,有河谷平坝、低山、中山、高山等。由于独特的地理条件和生态环境,多数玉米杂交种难以适应高海拔玉米生态区,导致当地玉米生产仍沿用白马牙、白鹤、大黄等农家常规种,尤其是部分高山区域,老百姓的自留种成为玉米生产的主导品种。近年来,随着杂交玉米的推广,人们盲目追求高产,加上不断变化的环境给玉米地方品种保存带来了前所未有的威胁,导致玉米地方品种资源逐年减少,这严重影响优异种质资源的挖掘与利用。到目前为止,尚未有人对大凉山玉米地方品种进行系统收集并对其群体结构、混杂程度、遗传多样性等开展相关研究。因此,本研究广泛收集和整理了来自大凉山26个县(市)的360份玉米地方品种,利用GBS简化基因组测序技术分析品种间的亲缘关系和群体结构,揭示不同品种间的进化关系,完善本区域的玉米分类系统,以了解大凉山玉米地方品种资源的遗传背景和系统演化,为下一步进行育种开发与利用提供科学理论依据。

1 材料与方法

1.1 试验材料

2019—2020年广泛收集并整理来自凉山州全境17个县(市),甘孜州稻城、乡城、得荣、巴塘,西藏察隅、芒康,云南巧家、永德,攀枝花米易的玉米地方品种360份(表1)。2021年种植于凉山州西昌市阿七镇大田村,单行区,行长5 m,株距0.25 m。每667 m2种植4 000株。

1.2 试验方法

1.2.1 基因组DNA提取

苗期剪取各样本植株幼嫩叶片200 g,放入事先标号且预冷的离心管,液氮速冻。待取完所有植株叶片后带回实验室,在-80 ℃超低温冰箱存放。用天根生物科技有限公司的DNAsecure新型植物基因组DNA提取试剂盒(DNAsecure Plant Kit)提取DNA。通过1%琼脂糖凝胶电泳和Nanodrop 2000分光光度计检测提取样品DNA的质量和浓度。

1.2.2 GBS文库构建

质检合格的样品DNA由青岛欧易生物科技有限公司利用SuperGBS技术[7]进行测序文库构建。采用PstI-HF/MspI对DNA进行酶切,用T4连接酶对酶切后的片段两端加接头和样品标签(barcode),用磁珠回收系统通过调整磁珠溶液与连接产物的体积比回收300～700 bp的片段,使用高保真酶对回收片段进行PCR扩增,PCR产物浓度经Qubit测定(质量浓度需大于5 ng·μL-1),将混好的文库上机(Illumina Nova,PE150)测序。

1.2.3 数据质控与比对

基于标签(barcode)与酶切位点信息,使用Stacks软件[9]对从测序仪下机的高通量测序数据进行拆分,得到每个样本的raw reads。使用软件fastp[10]对raw reads进行质控获得clean reads。使用BWA软件[11]将clean reads比对到B73参考基因组(B73_RefGen_v4)。

1.2.4 SNP检测与注释

基于样品与参考基因组的比对结果,使用GATK软件[12]的Haplotype Caller程序生成每个样品中的g VCF文件,再通过GATK软件的Genotype GVCFs程序进行群体SNP检测,使用GATK软件的Select Variants程序对得到的群体SNP预测结果进行筛选,得到初步SNP和InDel结果。使用VCFtools软件[13]对获得的SNP和InDel分型结果进行过滤,过滤条件为：(1)reads支持数不低于4;(2)剔除最小等位基因频率(MAF)小于0.01的位点;(3)剔除SNP或者InDel分型缺失率高于20%的位点。使用SnpEff软件[14]对得到的SNP进行注释,以确定SNP在基因元件的位置、对氨基酸的变化影响等。

1.2.5 系统关系分析

采用邻接法(neighbor-joining, NJ)[15]对样本的SNP位点进行建树,SNP位点中有缺失则用“-”代替。在treebest软件[16]中用p-distance方法计算距离矩阵,进化树的可靠性通过bootstrap法[17]进行重复1 000次检验。使用plink2[18]软件对获得的SNP标记进行主成分(principle component analysis, PCA)分析。利用ADMIXTURE软件[19]按照K=2到K=10进行群体结构分析,选取10个不同的seed进行10次重复分析。根据交叉验证误差(cross validation error, CVE)确定最优K值,对每个K值10次重复CVE值绘制盒形图。

1.2.6 群体分化指数(fixation index,Fst)和核苷酸多样性(π)的计算与选择性清除分析

基于高质量的SNP,按照100 kb的窗口、10 kb的步长对染色体进行选择性清除区域检测。使用R语言的PopGenome软件包,计算Fst、π和多样性变化倍数(θπratio)等,并利用Fst值绘制曼哈顿图。在筛选选择性清除区域时,分别以Fst和π前5%分位数对应的数据作为阈值,取该区域的交集为选择性清除区域。

1.2.7 候选基因的检测与注释功能富集分析

将鉴定的选择性清除区域通过Interproscan软件[20]进行GO注释分析。在网站(http：//geneontology.org)上进行GO富集分析,并与Swiss-Prot[21]、GO[22]和KEGG[23]等数据库进行对比得到注释信息。参考NCBI数据库(https：//www.ncbi.nlm.nih.gov/)和玉米数据库(https：//www.maizegdb.org),对选择区域的基因进行功能注释,对候选基因进行GO功能富集分析。

2 结果与分析

2.1 基因组DNA提取效果

部分样品的DNA电泳结果如图1所示,条带清晰明亮。同时,经Nanodrop 2000分光光度计检测,样品的DNA质量浓度均大于50 ng·μL-1,满足GBS文库构建要求。

M,DL2000 DNA分子量标准;1～24,不同样本DNA。M, DL2000 DNA marker; 1-24, DNA samples from different samples.图1 部分样品的基因组DNA电泳条带Fig.1 Genomic DNA electrophoresis of partial samples

2.2 GBS测序质量

360个样本因低质量(low_quality)被过滤掉的reads占比0.945%～5.180%,因含N数量≥5(too_many_N)被过滤掉的reads占比0.003%～0.056%,因质控后R1长度短于142 bp或者R2长度短于150 bp(too_short_reads)被过滤掉的reads占比0.314%～3.139%,质控后得到的clean reads占比92%以上(92.90%～98.36%)。各样本的单碱基平均错误率均在0.000 3左右,最高在0.000 5左右。所有样本的A、T碱基比例均在30%左右,G、C碱基比例均大于20%,平均GC含量为48.20%(47.18%～49.57%)。

对360份玉米地方品种基因组DNA进行GBS测序,共产生250.99 GB有效数据,获得1 659 033 712个clean reads,平均reads数量为4 608 427;平均每份样品的raw base为0.70 GB,clean base为0.67 GB。各样品的Q20平均值大于等于94.57%,Q30平均值大于等于87.14%。表明测序的质量高,可以用于后续的信息分析。

2.3 参考基因组比对

使用BWA软件将clean reads比对到B73参考基因组上,通过比对率检测样品与参考基因组的相似程度,一般比对率越高,相似程度越高。根据比对结果统计样品深度信息并计算测序数据对参考基因组的覆盖度。全部样本比对率为93.55%～98.74%,平均测序深度为14.79×,覆盖度范围为1.31%～3.02%。使用Qualimap 2[24]统计各样本的插入片段长度分布,预估的插入片段长度符合正态分布(图2),峰值位置接近建库时DNA片段的平均长度。

图2 丹红_1的插入片段长度分布Fig.2 Distribution of insert size of Danhong_1

2.4 突变(SNP/InDel)检测

经过突变位点检测并过滤,共获得124 342个SNP位点,32 063个InDel位点(其中插入15 738个位点、缺失16 325个位点),总共检测到156 405个变异位点。发生碱基转换位点12 955 634个、碱基颠换位点7 955 116个,碱基转换/碱基颠换比例为1.628 591。在染色体上通过滑窗(窗口和步长均为500 kb)统计,得到各样本SNP/InDel在染色体上的分布,如图3所示。

2.5 SNP注释

SNP注释分析表明,124 342个SNP变异位点中对蛋白质具有高效(破坏性)影响,可能引起蛋白质截断、功能丧失等的变异(high)有1 674个,占0.174 9%;低影响,可能不影响蛋白质变异的(low)有91 776个,占9.587 2%;中等影响,非破坏性变异,可能影响蛋白质功效的变异(moderate)有63 458个,占6.629 0%;修饰,通常为非编码区变异,只影响非编码基因(难以预测其影响程度)的变异(modifier)有800 364个,占83.608 8%(表2)。错义突变有63 654个,占44.469 8%;无义突变有836个,占0.584 0%;同义突变,有78 650个,占54.946 2%(表2)。

表2 变异影响程度与功能级别影响统计

变异分布最多的是内含子区域的内含子变异,占32.793 6%,其次是基因下游区域的基因下游变异,占22.508 0%,第三是外显子区域与基因上游变异,分别占14.878 9%和14.673 4%(表3、表4)。

表4 变异类型分布

2.6 系统关系分析

从系统进化树(图4)可以看出,整个群体被分成2大类群,类群之间分层明显。各类群可进一步分离出若干亚群,且存在离群个体。

图中的品种编号与表1一致。下同。The variety codes in the image are consistent with those in Table 1. The same as below.图4 360份大凉山玉米地方品种的系统进化树Fig.4 Phylogenetic tree of 360 maize landraces in Daliangshan Mountain area

2.7 主成分分析

对获得的124 342个SNP标记进行PCA分析,得到影响最大的3个特征向量,3个主成分变异贡献率分别为16.13%、2.57%和2.10%,对3个主成分作图,结果见图5。由图5可知,360份大凉山玉米地方品种明显分成两部分,各形成一个紧密类群,说明2个类群遗传背景相差较大;其中一个类群比较分散,而另一个类群相对集中,没有明显的遗传差异,呈混杂在一起的趋势。主成分分析结果与系统进化树分析结果一致。

图5 主成分分析Fig.5 PCA analysis of three-dimensional graph

2.8 群体遗传结构分析

最佳K值分析结果表明(图6),随着K值的增大,CVE值顺势减小,当K=9时,CVE值最小,随着K值再增大,CVE值也随之增加,表明K=9时是最合适的分组,K=9时的分组情况见图7。群体遗传结构指遗传变异在物种或群体中的一种非随机分布。按照地理分布或亲缘关系可将一个群体分为若干亚群,处于同一亚群内的不同个体亲缘关系较近,而亚群与亚群之间则亲缘关系稍远。从图7可以看出,360份大凉山玉米地方品种被分成9个亚群。

图6 不同K值下群体结构的CVE值曲线Fig.6 The CVE curve of the population structure with different K values

图7 360份大凉山玉米地方品种的群体遗传结构Fig.7 Genetic structure of 360 maize landraces in Daliangshan Mountain area

2.9 群体遗传分化分析

估算群体间的Fst可有效区分群体间相对遗传变异大小,是解释群体遗传变异程度的主要来源。从图8可见,类群A与类群B的Fst值大于0.462 2,说明类群间有很大的遗传分化。类群A的核苷酸多样性πA高于类群B的核苷酸多样性πB。可能是因为类群A遗传背景较为复杂,群体内玉米材料血缘丰富,导致类群A遗传多样性高于类群B,深层次原因需进一步研究。

2.10 差异化的候选区域与候选基因

以玉米B73基因组为参考,分析360份玉米差异基因组区域中的基因。从图9可见,在A群体和B群体中共检测到96个受选择区域,其中418个基因具有强烈选择信号。进一步对选择性扫描的区域进行基因功能注释,发现富集的基因组区域中包含寒冷响应基因Zm00001d027468(OJ2056_H01.33)、Zm00001d006508(YUP8H12.19)、Zm00001d020961(T24D18.5)、Zm00001d048168(F18A8.1),缺水响应基因Zm00001d033107(OsJ_16068)、Zm00001d013399(F2H15.17)、Zm00001d022002(OsJ_11395)、Zm00001d048165(F13H10.2),病菌响应基因Zm00001d028273(T9L24.32)、Zm00001d005727(T8O11.18)、Zm00001d006507(F8K4.6)、Zm00001d047690(MZN1.6)、Zm00001d025318(F1N19.11)等与逆境应答相关的基因。

lg(π ratio)和Fst分别对应上面的频率分布图和右侧的频率分布图,中部的点图则代表不同窗口内的相应的Fst和π比值。其中最上方绿色和红色区域为π选择出来的top 5%区域,右侧橘色区域为Fst所选择top 5%区域,中间红色和紫色区域为Fst和π的交集,即为候选的位点。lg(π ratio) and Fst values, corresponding to the frequency distribution chart above and on the right, the point plot in the middle represents the corresponding ratio of Fst to π in different windows. The green and red areas at the top are the top 5% areas selected by π, and the orange area on the right is the top 5% area selected by Fst, and the red and purple areas in the middle are the intersection of Fst and π, which are the candidate loci.图9 类群A和B的差异基因组区域Fig.9 Different genomic regions between group A and B

3 结论与讨论

本研究采用GBS简化基因组测序技术对收集自大凉山26个县(市)的360份玉米地方品种进行了基因分型,获得了大量的SNP分子标记,经过SNP calling并过滤,共获得124 342个SNP位点、32 063个InDel位点(其中插入15 738个位点、缺失16 325个位点),总共检测到156 405个变异位点。基于获得的124 342个SNP位点,系统进化树聚类分析和PCA主成分分析将群体分成A和B两大类群,类群间遗传分化明显(Fst大于0.462 2),类群A的核苷酸多样性高于类群B。玉米地方品种在经历了长时间的人工选择和对大凉山地区的气候环境产生高度适应性后,类群间已经产生了明显的遗传差异,形成了独立的遗传背景,且在进化过程中受到的选择强度较大,推测可能是由于长期的自然选择和人工选育造成的。这些数据将为大凉山玉米地方品种的表型差异和复杂性状的遗传分析提供重要的遗传基础依据。研究结果将为本区域玉米重要目标性状基因挖掘、农艺性状QTL定位以及玉米新品系、新品种培育提供重要支持。

玉米种质类群的划分问题一直存有争议。我国玉米种质资源类群划分研究始于20世纪50年代中后期,主要为瑞德群、兰卡斯特群、塘四平头群和旅大红骨群4个类群。殷洪达等[25]根据1年2点产量和相关性状的配合力与杂种优势表现,将40份中晚熟、耐密、抗病自交系分为Reid、Lancaster、塘四平头和旅大红骨4个杂种优势群。吕学高等[26]利用40个SSR分子标记将浙江省90份玉米地方种质资源划分为塘四平头、Reid、旅大红骨和Lancaster等类群。吴承来等[27]利用112对SSR标记将我国97个玉米自交系划分为PB、Reid、塘四平头和旅大红骨4个类群。类群多而复杂并不利于提高育种效率。张世煌等[28]认为杂种优势群的划分应尽量向两群靠拢,并最终归纳为两个杂种优势群,故他提出了两群划分策略。本研究通过GBS测序技术进行系统发育与群体结构分析,把360份大凉山玉米地方品种分为两大支系9个亚群,这与Mumm等[29]利用46个RFLP标记将148个美国常用自交系划分为Lancaster Sure Crop(LSC)和Iowa Stiff Stalk Synthetic(BSSS)两大类群11个亚群,孙友位等[30]利用70个SSR分子标记将85份早熟玉米自交系划分A(Reid、热带种质、旅大红骨)和B(Lancaster、塘四平头、农家种)两大类群6个亚群,张伟玮等[31]利用43个SSR标记将56份玉米自交系划分为SS(旅大红骨、PA、Reid)和NSS(PB、四平头、Lancaster)两大种质5个亚群,王文斌等[32]利用2 846个高质量SNP标记将陕西31份玉米自交系划分为A(丹340、掖478等)和B(Mo17、黄早四等)两个类群,吴迅[33]利用56 110个SNP标记将367玉米自交系分为外引种质和本地种质两大类群5个亚群(瑞德、兰卡斯特、塘四平头、温热I和P群)是基本一致的。本研究的系统进化树结果支持玉米杂优类群的两群学说,但类群间是否具有育种上的杂交组配优势,还需要对具体材料的配合力进行分析。本研究的群体结构分析进一步将两大类群分成9个亚群,与上述前人研究结果不同,这可能与研究材料本身的血缘有很大关系,导致分析所用的数据集不同而呈现差异。

大凉山地处西南高原,是玉米传播到中国后较早的种植区之一。境内大山密布,交通极不方便。在海拔1 800～2 800 m的山地和山谷旱地之间,因杂交种难以适应其独特的地方性气候和生态条件,至今依然保留了较多的地方品种资源。这些品种由于长期所处气候、生态条件的差异,形成了各自鲜明的特点,如抗旱、抗寒、耐贫瘠、抗病等,籽粒颜色丰富多彩(红色、白色、紫色、黄色等),其中蕴含了复杂多样的基因资源,使其成为我国玉米地方种质最为丰富的地区。此外,由于相对落后的文化、交通和经济条件,使得该区域内种质资源交流较少,同时外来种质也难以适应局部小气候,进一步导致种质交流和遗传渗透远远落后于其他玉米生态区,从而使得该区域玉米地方种质资源能够得到很好的保存和利用,遗传多样性也远远高于其他玉米生态区。大凉山的玉米发展离不开地方特有玉米种质资源的挖掘,尤其是耐寒、耐旱和耐瘠群体、白粒适口玉米等自交系材料或群体的改良和扩增工作。本研究中,通过选择信号分析在大凉山地方玉米品种群体中发现96个基因组区域受到选择,检测到受选择区域的418个候选基因,部分基因已被证明与玉米的抗寒、抗旱和抗病等有关。杨宇昕等[34]利用群体分化指数和群体间扩展单倍型纯合度(cross population extended haplotype homozygosity, XP-EHH)法分析温带和热带玉米群体间的选择信号分布情况,选择Fst和XP-EHH值的top 1%为阈值,分别鉴定到557和1 913个候选基因,多个候选基因与玉米的开花调控密切相关,包括ZmCCT9、COLl、GRMZM2G387528。周玲等[35]基于187份种质资源材料的全基因组重测序数据开发了120 583个高质量SNP标记,通过整合选择信号检测和GWAS分析结果,共识别出与1 153个候选位点紧密连锁的324个候选基因,涉及氮化合物代谢、叶酸代谢、糖酵解、发育过程的负调控与转录调控等重要生物学途径。

西南山地玉米区是我国三大玉米主产区之一,全区玉米种植面积500万hm2左右,约占全国玉米生产面积的1/5[36]。该区地形地貌十分复杂,山地丘陵和高原土地总面积占90%以上,从海拔250 m河谷平坝直至2 500 m的高山均有玉米种植。生态条件区域差异很大,垂直分布明显,农业立体性强。长期以来,西南地区玉米育种使用的自交系绝大多数为二环系,主要选育来自美国玉米带杂交种和我国北方玉米优良杂交种,种质基础极为狭窄,导致所选育的杂交种产量没有突破、稳定性差,抗病性、适应性不强。这主要是由于急功近利的商业育种热衷于抄近道,追求多(育种材料多)、短(组配时间短)、快(上市品种快),盲目跟风育种,却未能静下心来认真研究本土种质。西南地区类似大凉山这样的山地玉米生态区还有很多,区内地形地貌复杂,生态条件多种多样,种植玉米历史悠久,加之长期以来部分区域交通不便造成相对闭塞,形成了许多适应性强、独具特色的玉米地方种质资源。这部分地方玉米种质资源具有很大的开发利用潜力,目前存在的主要问题是鉴定、筛选和研究利用严重滞后。因此,对于西南玉米杂交育种,面对多种多样的杂交种类型需求,必须把种质资源的拓展与创新作为今后重要的攻关方向;建议加强对本土地方玉米优良种质的评价与遗传多样性研究,从特异抗逆基因的发掘与利用着手,简化杂优模式,分类合成与改良育种群体。