基于基因分型技术的燕麦SNP标记研究

董艳辉，刘龙龙，温 鑫，于宇凤，杨 方，刘根科，崔 林，曹秋芬，秦永军,

(1.山西省农业科学院 生物技术研究中心，山西 太原 030031；2.山西省农业科学院 农作物品种资源研究所，农业部黄土高原种质资源实验室，山西 太原 030031；3.山西省农业科学院 农业科技信息研究所，山西 太原 030031；4.山西省农业科学院 右玉农业试验站，山西 右玉 037200)

燕麦(AvenasativaL.)是世界各地广泛栽培种植的一种重要的粮食兼饲草、饲料作物，其富含β-葡聚糖，被证明具有降血糖和血脂的功效[1]，是美国FDA和英国JHCI认定的功能性谷物[2-3]，随着燕麦及其加工产品的世界消费量的增长，传统燕麦育种技术已经不能有效满足燕麦生产的需要，并且成为制约燕麦育种进程的重要因素，而分子辅助育种是加速燕麦育种的重要手段。

目前，燕麦的育种多依赖于常规育种和表型鉴定研究，已成功育成了多种适应不同需求的燕麦品种并用于大田生产[4]。早期的表型特征量化研究证明燕麦具有丰富的遗传多样性[5]，之后燕麦种质资源遗传连锁图[6]、皮裸性基因相关联SSR[7]和AFLP标记[8]也被相继报道。单核苷酸多态性(SNP)是目前最具优势的分子标记，被广泛应用于作物遗传图谱构建、QTL定位、全基因组关联分析(GWAS)、群体进化等[9-12]。近年来，基于下一代测序技术(Next generation sequencing,NGS)的基因分型技术(Genotyping by sequencing,GBS)是通过酶切处理去除基因组重复序列，显著降低了测序量和测序成本，且能够覆盖整个基因组，不受参考基因组限制，因而被广泛应用到开发高密度、高精度的农作物SNP分子标记研究中。目前，燕麦全基因组序列还没有发表，基因分型技术对农作物遗传多样性研究、种质资源鉴定和分子辅助育种标记开发具有重要的意义[13-16]。基因分型技术从单酶切技术RAD[17]、2b-RAD[18]、GBS[19]到现在的双酶切技术dd-GBS[20]、SLAF[21]均已被应用到玉米、小麦等[22-30]主要农作物分子标记开发研究。Winkler 等[31]采用GBS技术，证明了燕麦颖片颜色和穗型等表型特征和SNP位点相关。在另一项研究中，通过GBS技术构建的高密度SNP则能精准地区分起源于世界不同地区燕麦种质资源[32]，证实GBS-SNP技术是种质资源鉴定的有力工具。

本研究基于NGS的dd-GBS技术，以经过抗旱性鉴定的燕麦种质为试验材料，选用适合燕麦基因组酶切的酶，运用适合燕麦简化基因组参考序列的分析软件，在构建燕麦简化基因组参考序列的基础上，研究燕麦相关SNP标记，旨在丰富当前燕麦的基因组数据库，提供分子标记(SNP)基础数据，为加速燕麦育种进程提供理论支撑。

1 材料和方法

1.1 试验材料

供试材料由中国农业科学院品种资源研究所麦类室于1981-1983年进行抗旱鉴定的燕麦种质42份，包括11份高抗(High tolerance)、13份中抗(Moderate tolerance)、18份低抗(Low tolerance)和山西省农业科学院农作物品种资源研究所燕麦种质资源圑提供[33]的抗旱性未鉴定(Uncharacterized)燕麦育成品种(品系)12个(表1)，2016年种植于山西省农业科学院智能温室。

1.2 试验方法

GBS文库构建与数据质控：燕麦种子发芽14 d，取鲜嫩叶片提取基因组DNA(Mag-MK Plant Genomic DNA Extraction Kit试剂盒，上海生工)，1.0%琼脂糖凝胶电泳和紫外分光光度计分别检测DNA质量。用限制性内切酶PstⅠ(CTGCAG)和MspⅠ(CCGG)对提取的DNA 进行酶切后，回收大小在220～450 bp的酶切片段，之后按照dd-GBS方法进行建库[20]。采用llumina Hiseq测序平台进行双末端(Paired-end，PE)测序，获得的下机数据以双端FASTQ格式保存。

分别对每个样品的下机数据运用 FastQC软件进行质量控制，主要包括碱基质量分布(Per base sequence quality)、质量值(Per sequence quality scores)、GC分布(Per sequence GC content)，对反映下机数据质量的多项指标进行统计。进一步过滤数据包含一些带接头，双末端reads1 5′端6 bp非酶切位点序列CTGCAG或reads2的5′端4 bp非酶切位点序列CCGG。采用Adapter Removal去除3′端的接头污染[34]。采用滑动窗口进行质量过滤，窗口大小设置为5 bp，步长设置为1 bp，每一次往前移动1个碱基，取5个碱基计算窗口的平均Q值，若最后一个碱基的Q值≤2，则仅保留该位置之前的碱基；若窗口的平均Q值≤20，则仅保留该窗口倒数第2个碱基及之前的碱基。滤除小于50 bp 的reads。

表1 不同抗旱性的燕麦试验材料Tab.1 Sequenced oats germplasm of different drought tolerance

1.3 数据分析

1.3.1 燕麦简化基因组参考序列组装与群体SNP检测 应用Stacks软件包[35]中的ustacks对每个样品的reads进行聚类，同一个stack代表1个酶切位点(loci)，聚类参数-m设置为4，对每一个样品的loci及loci的测序深度进行统计。下一步用cstacks将所有样品的loci合并，不同样品loci之间最多允许2个错配，获得每个loci的catalog consensus序列。

采用sstacks将每个样品的loci序列与catalog consensus序列比对后populations过滤获得群体SNP。主要参数包括：1个位点最少要在1个群体中出现；1个群体中检测到同一位点的个体最低百分数50%(当群体中该位点的缺失率超过 50%，则去除该位点)；1个位点的最小等位基因频率0.05。

1.3.2 燕麦SNP分析 采用 GCTA软件利用群体SNP数据进行主成分聚类分析；对高、低抗材料的SNP数据进行Fisher Test统计分析获取相关SNP标记；采用结合校正标准错误发现率(False Discovery Rate，FDR)小于0.05和0.001确定差异显著性标记和极显著性标记。

1.3.3 燕麦SNP功能预测 进一步对差异极显著(FDR<0.001)的SNP所在源序列和小麦(TriticumaestivumL.)基因组序列进行Blast比对，筛查SNP相关基因。

2 结果与分析

2.1 燕麦GBS测序结果

本研究共测序54个样本，注释reads数在4 111 218～21 782 382；质控后保留的高质量数据(reads数4 111 218～19 019 296)用于后续的数据分析(图1)。

图1 燕麦GBS-SNP注释序列Fig.1 Oat GBS-SNP annotation sequence

2.2 燕麦参考序列与群体SNP位点

采用Stacks数据包对燕麦reads进行聚类，结果显示：每个loci平均测序深度为20.61×，合并所有样本同一个loci，共获得753 325个参考基因组序列(catalog consensus)。每个样品的loci与其对应的catalog比对，去除群体位点缺失(./.)率超过50%和最小等位基因频率小于0.05的SNP位点，共获得74 657个群体SNP位点。由图2可知，所有测序样本中纯合状态的SNP位点(0/0或1/1)占少数，大部分呈现不同程度的杂合(0/1)或缺失状态(./.)。

0/0 和1/1.纯合体；0/1.杂合体；./..缺失。0/0 和1/1. Homozygote; 0/1.Heterozygote;./.. Deletion.

2.3 燕麦SNP研究

根据群体SNP数据，对54份燕麦种质的PCA结果显示，SNP基因型明显聚类为2簇，全部11份高抗、4份抵抗、8份中抗种质和6个未鉴定品种聚为一簇；而另外一簇则包含14份低抗和5份中抗种质；其余6个未鉴定的燕麦品种则分布在2簇之外(图3)。

基于PCA基因型聚类结果，进一步对高抗(n=11)和低抗(n=18)种质2组进行SNP基因型差异显著分析，共获得2 937个燕麦抗旱性相关SNP(Fisher Test，FDR<0.05)。其中，差异极显著(FDR<0.001)的55个SNP位点坐落于41个燕麦简化基因组参考序列(同一参考序列包含2或3个SNP位点)。这55个相关SNP可分为3类，23个SNP 位点在高抗种质中缺失； 28个SNP位点在低抗种质中缺失；其余4个SNP位点在2种种质中都存在(表2)。

表2 高抗与低抗燕麦的差异显著SNP位点Tab.2 Oat significant SNPs of high tolerance and low tolerance

表2(续)

注：▲.SNP位点缺失。

Note: ▲. Represents site deletion of SNP.

2.4 燕麦SNP的Blast功能预测

进一步对差异极显著的55个相关SNP的源序列与小麦基因组序列进行了Blast比对，结果显示，10个源序列(包括14个SNP位点)和多个基因转录序列联配(表3)。其中，参与线粒体活性氧信号转导、植物激素信号转导的调控蛋白、自噬基因表达与植物抗逆、抗旱水平有关。

表3 燕麦SNP源序列联配基因(比对小麦基因组)Tab.3 Oat SNPs alignment genes (against wheat genome)

3 结论与讨论

3.1 GBS技术开发高密度SNP分子标记

高通量测序技术的面世是现代分子生物学技术的革命性飞跃，NGS技术促使生命科学研究进入基因组学时代，NGS在农作物上的重要应用之一就是GBS技术。GBS技术已经成为加速农作物种质资源鉴定筛选以及分子辅助育种的重要技术之一[16,19,46]，燕麦基因组相对复杂(2n=6x=42)，具有丰富的遗传多样性。本研究采用新一代ddGBS技术，成功组装了基于54个样本的燕麦简化基因组参考序列。质控显示，获得的高质量reads通量和SNP杂合状态与抗旱性没有明显关联。GBS通过酶切过滤重复序列，降低基因组复杂程度，虽不代表全基因组序列结构，但非常适合于农作物分子辅助育种SNP标记开发。本研究之所以选择经过多年抗旱性鉴定、遗传稳定、抗旱梯度差异明显(高抗、中抗、低抗)的燕麦种质资源进行测序，目的是通过抗旱性表型与SNP基因型关联分析，开发可能与抗旱性有关联的SNP。因本研究样本数量相对偏小(n=54)，在构建简化基因组参考序列中没有采用报道的针对大样本TASSEL分析流程[33]，而应用更适合于小群体的Stacks分析流程[35]。虽技术路线不同，但每项研究贡献的简化基因组参考序列无疑丰富了目前的燕麦基因组数据库。

3.2 抗旱性和SNP基因型关联分析

本研究的主要目的是挖掘可能与燕麦抗旱性有关联的SNP，为燕麦种质资源鉴定和分子辅助育种服务。这些试验材料根据多年的抗旱性鉴定划分为高抗、中抗和低抗3组，但根据SNP基因型聚类分析则大致分为2簇，高抗材料SNP基因型和大多数低抗材料基因型各聚类为一簇，明确显示抗旱性和SNP基因型是关联的，也就是说SNP基因型很大程度上能反映抗旱性。中抗材料SNP基因型没有聚类成簇，而是呈分散状态，也从侧面说明表型是基因型与环境共同作用的结果，同时作物抗旱性又是一个多指标的综合性状表现，但基因型的选择是改良作物性状的有效途径[47-48]。燕麦抗旱表型是抗旱基因决定的性状在环境作用下的具体表现，其抗旱性鉴定结果可能在不同的环境条件下会有所偏差。

因高、低抗材料的抗旱性和SNP基因型相关联，本研究针对这2组材料，开发可能与燕麦抗旱性相关的SNP，大部分SNP属于杂合或缺失状态，只有极少数呈现纯合状态，可能是因为燕麦基因组高度杂合(六倍体)以及测序技术的原因。下一步仍需要应用精度更高的数字化 PCR 技术在更大群体验证本研究发现的SNP。

利用GBS技术开发可能与燕麦抗旱性相关SNP标记在其他作物育种研究中已经报道[49-53]。对于燕麦来讲，本研究采用GBS技术，第1个报道可能与燕麦抗旱性有关联的SNP标记，以期为燕麦种质资源精准筛查和分子辅助育种提供数据基础。

3.3 相关SNP生物学功能

基因组技术已被广泛应用到抗旱基因挖掘[54-56]。GRLs[57]、OsMYB55[58]和GIF[59]被相继报道与不同作物抗旱机制有关。本研究经Blast比对获得的相关SNP源序列和不同生物通路基因转录子联配。其中cat 186303的第83个碱基在高抗材料中位点缺失，而在低抗材料中则呈现GT杂合状态。cat 186303源序列Blast 比对小麦基因组数据库显示，其序列与小麦硫氧还蛋白Trx序列相似。Trx蛋白参与细胞生长凋亡和基因转录[60]，多项研究证明，Trx基因和作物抗旱、抗逆功能有关[16, 61-65]。后续功能验证包括研究 cat 186303对燕麦Trx基因表达调控，是解释功能SNP的重要步骤。此外，其他抗旱性相关SNP源序列联配的基因包括植物激素信号转导的DUF221、双特异蛋白磷酸酶和核孔蛋白，与植物抗逆有关。

3.4 本研究的局限性

燕麦是六倍体作物，本研究所开发的燕麦抗旱相关联SNP不能定位，但证明不同抗性材料的燕麦在SNP存在表达差异，是作为燕麦材料分子筛选的潜在标记。因为抗旱性状是一个综合指数，单一参数不能全面代表抗旱能力。本试验得出的特异性SNP 需要今后在大量的群体材料中进行验证。此外，燕麦抗旱指标有待量化也是今后挖掘抗旱SNP 的重要基础。本研究结果支持表型-基因型互作关系的同时也指出，抗旱性是基因型和环境因子综合体现，初期开发的可能与燕麦抗旱性有关联的SNP必须在更大种质群体中验证以加强其SNP的广谱性和实用性。

常规育种主要依赖于表型数据，而分子辅助育种则侧重于基因型数据，构建燕麦抗旱性表型-基因型关联则是精准育种的必经之路。应用GBS技术，本研究成功组装了燕麦简化基因组参考序列并初步研究了可能与燕麦抗旱性有关联的SNP标记，研究结果不仅丰富了燕麦基因组数据库，同时为燕麦种质资源早期筛选和燕麦基因数据库提供数据基础。

致谢：感谢荷兰莱顿大学统计系Goeman J J教授的统计学方法指导；感谢美国加州戴维斯分校植物系罗明诚教授对本试验给予的建议。