基于SLAF-seq技术的甘蔗遗传多样性和选择信号分析

林萍萍, 张 慧, 徐良年, 邓祖湖, 王勤南, 赵新旺,3

(1.福建农林大学国家甘蔗工程技术研究中心，福建 福州 350002；2.广东省科学院南繁种业研究所,广东 广州 510310；3.广西大学广西甘蔗生物学重点实验室，广西 南宁 530004)

甘蔗是热带、亚热带重要的经济作物[1]。由于甘蔗是非整倍的多倍体作物，遗传方式复杂且遗传多样性丰富，可为“高贵化”育种提供优良的亲本材料[2]。甘蔗倍性复杂且不易开花，主要以无性繁殖为主，使育种进展缓慢[3]。目前应用于遗传育种的分子标记有很多，单核苷酸多态性(single nucleotide polymorphisms, SNP)标记由于具有数量多、方便基因分型、覆盖范围广、效率高及可以大规模筛查的优势[4-6]，得以普遍应用于探究遗传多样性。简化基因组测序技术(specific-locus amplified fragment sequencing, SLAF-seq)是一种高通量的测序技术，可简便、快速获得覆盖全基因组的SNP标记以及基因型。SLAF-seq技术具有构建文库简单快速、成本较低、准确性高、有效标记多等特点，已应用于甘薯[7]、大豆[8]、水稻[9]等作物遗传多样性分析和性状的QTL定位等研究。

研究甘蔗种质资源的遗传多样性可以挖掘控制优良性状的等位基因，通过人工选择使有利等位基因的基因频率不断提高，甚至固定下来。由于连锁，有利基因周围的基因组区域也随之稳定遗传。甘蔗细茎野生种(SaccharumspontaneumL.)俗称割手密，具有抗逆性强、适应性广等优良性状。推广种植的甘蔗栽培品种大多含有割手密血缘。在全基因组范围检测选择性清除信号，能扫描到与遗传进化有关的受选择基因或基因组区域[10-11]。目前利用SNP标记对甘蔗遗传多样性及选择性清除进行分析的相关研究相对较少。本研究以甘蔗近缘属割手密为参考基因组，对107份甘蔗材料进行简化基因组测序，分析不同群体的遗传多样性，并利用遗传分化指标(Fst)和核苷酸多态性(π)对不同群体进行选择性清除分析，进一步分析甘蔗基因组中受人工选择的区域及优良基因，以期为甘蔗分子育种提供参考。

1 材料与方法

1.1 试验材料

选用107份甘蔗材料进行简化基因组测序，包括83份栽培种及24份细茎野生种。其中，粤甘49在两次引种过程中表型有所差异，但性状都优良，分别编号为粤甘49-1、粤甘49-2。甘蔗叶片主要采自福建、海南、广西等地。取新鲜幼嫩无病害的叶片放于锡箔纸中，置-80 ℃冰箱备用。

1.2 试验方法

1.2.1 DNA的提取 采用2×CTAB法[12]提取甘蔗叶片DNA，并用1.0%琼脂糖凝胶电泳。使用紫外分光光度计测定DNA的纯度和浓度，选择单一、无拖带的DNA样品置于-20 ℃备用。

1.2.2 酶切建库及SNP标记的开发 以甘蔗近缘属割手密基因组(S.spontaneumAP85-441 genome, http://sugarcane.zhangjisenlab.cn/sgd/html/index.html)为参考，利用北京百迈客生物科技有限公司的酶切预测软件对割手密基因组进行最佳酶切方案预测，并对107份甘蔗DNA样品进行酶切。对得到的SLAF标签进行3′端加A处理，连接Dual-index测序接头，通过PCR目的片段的扩增、纯化、混样，选取目的片段，构建合格的测序文库。以水稻品种‘日本晴’作为对照进行测序，评估酶切的准确性。利用BWA软件[13]将测序后的数据比对到参考基因组上，并使用GATK和SAMtools[13]两种软件检测SNP。两种方法取交集以确认最终的SNP位点，并且以最小等位基因频率(minor allele frequency, MAF)大于0.05以及位点完整度大于0.5为标准进行数据过滤，获得高质量的SNP。

1.2.3 系统发育树构建与群体结构分析 采用MEGA 6软件[14]，基于邻接法、Kimura 2-parameter模型和自展值重复1 000次，构建系统发育树。基于筛选出的SNP分别对材料进行群体结构和主成分(PCA)分析。(1)通过Admixture软件[15]分析群体结构。预先设定亚群数目(K值)为1～10进行聚类，并对聚类结果进行交叉验证。根据交叉验证错误率的谷值确定最优分群数。(2)通过EIGENSOFT软件进行PCA分析，对所有材料进行聚类。

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

1.2.5 候选基因的检测与注释功能富集分析 将鉴定的选择性清除区域通过Interproscan软件[16]进行GO注释分析。在网站(http://plantregmap.cbi.pku.edu.cn/go.php)上进行GO富集分析，并与Swiss-Prot[17]、GO[18]和KEGG[19]等数据库进行对比得到注释信息。参考NCBI数据库(https://www.ncbi.nlm.nih.gov/)以及割手密基因组数据库，对受选择区域的基因进行功能注释，对候选基因进行GO功能富集分析。

2 结果与分析

2.1 测序数据统计与评估

通过对107份甘蔗基因组DNA进行SLAF酶切及建库，总共获得223.63 Mb序列数据。其中，云南82-29获得的数据量(4.66 Mb)最大，湛蔗80-101等样品的数据量(0.27 Mb)最小(表1)。107份材料测序结果显示，GC含量在42.08%～45.69%之间，平均44.11%；测序碱基质量值Q30在89.41%～93.08%之间，平均90.93%，说明碱基测序错误率较低，获得的数据可靠。本研究中第42个样品即新台糖26号的碱基含量分布和质量分布见图1。通过序列分析，从107份甘蔗材料中共获得7 869 726个高质量的SNP。根据SNP在染色体上的分布，绘制SNP在染色体上的分布图(图2)。由图2可知，开发的SNP标记在染色体上分布比较均匀，Chr7D的SNP位点相对较多，Chr8C的SNP位点相对较少。

表1 SLAF-seq基因组测序数据统计表Table 1 Summary of genomic sequences generated by SLAF-seq

续表1

图1 新台糖26号的碱基含量及测序质量分布Figure 1 Base distribution and quality distribution of Xintaitang 26 reads

每个条带代表一条染色体，横坐标为染色体长度，颜色越深代表SNP标记数越多。图2 SNP标签在甘蔗染色体上的分布Figure 2 Distribution of SNPs on the chromosomes of sugarcane

2.2 甘蔗系统发育分析及群体结构划分

为了从基因水平上分析栽培种群体与细茎野生种群体的遗传结构，利用筛选出的7 869 726个高质量的SNP，通过Admixture软件进行107份甘蔗材料群体结构分析(图3)。根据交叉错误率确定最优分群。当K从1到2时，交叉错误率逐渐减小；从2到10时，交叉错误率逐渐增大(图3C)。说明K=2时，交叉错误率最小。因此，107份甘蔗群体可分为2个亚群，一个为栽培种群体，另一个为细茎野生种群体(图3D)。

利用SNP信息，构建107份甘蔗材料的进化树(图3A)。从进化树可见，细茎野生种群体、栽培种群体分别聚在一起，说明栽培品种间的序列相似性较高，与细茎野生种存在一定的差异。利用R软件绘制主成分分析图(图3B)，发现2个亚群在PC1轴上的分布存在差异。该结果与聚类分析结果相一致。周珊等[20]研究表明，甘蔗不同亚群间存在遗传渗透，尤其是栽培种群体中混有细茎野生种的遗传成分。细茎野生种将优良的抗性基因渗入栽培品种中，通过品种间不断地杂交、人工选择及自然选择，丰富了甘蔗遗传背景，培育出抗逆性强和产量高的优良品种[21]，与本研究遗传多样性的分析结果相类似。

A.进化树分析，Ⅰ为栽培种群体、Ⅱ为细茎野生种群体；B.群体主成分分析；C.群体结构交叉验证错误分析；D.群体结构。图3 107份甘蔗材料遗传结构的分析Figure 3 Genetic structure analysis of 107 sugarcane varieties

2.3 甘蔗群体遗传分化分析

为了探究甘蔗群体的分化，本研究分析了栽培种群体和细茎野生种群体的Fst值，并利用100 kb滑动窗口对基因组区域Fst值进行分析(图4)。从图4可见，栽培种群体与细茎野生种群体之间Fst=0.216，栽培种群体的核苷酸多样性(π=8.93×10-6)高于细茎野生种群体(π=6.11×10-6)。可能是因为栽培种遗传背景较为复杂，除了割手密血缘外，还有热带种血缘或者大茎野生种血缘[22]，导致栽培种群体遗传多样性可能高于细茎野生种群体，深层次原因需进一步研究。

2.4 差异化的候选区域及候选基因

利用群体间的群体分化系数以及核苷酸多样性比率进行选择性清除分析。使用100 kb的滑动窗口以及10 kb步长分别计算群体分化系数和核苷酸多样性比率，并取基因组前5%的交集区域作为选择清除区域(图5)。以割手密基因组为参考，分析107份甘蔗差异基因组区域中的基因。从图5可见，在栽培种群体和细茎野生种群体中共检测到72个受选择区域，其中439个基因具有强烈选择信号。进一步对选择性扫描的区域进行基因功能注释，发现富集的基因组区域中包含UBA3、CBP2、GSTU8等与抗性相关的基因(表2、图5)。

红点和蓝点分别表示Fst和多样性变化倍数大于95%群体间筛选出的基因组区域。图4 甘蔗栽培种与细茎野生种群体间的差异基因组区域Figure 4 Different genomic regions between cultivated sugarcane population and wild S.spontaneum population

红色、蓝色线分别表示基因组前1%、前5%水平的阈值线；对受选择区域的基因进行功能注释分析获得箭头所示基因。图5 甘蔗栽培种与细茎野生种群体分化系数在染色体上的整体分布Figure 5 Global Fst distribution between cultivated sugarcane population and wild S.spontaneum population

表2 甘蔗栽培种与细茎野生种群体间高差异基因组区域前5% SNP的部分基因Table 2 Partial candidate genes from highly different genomic regions between cultivated sugarcane population and wildS.spontaneum population (top 5% of SNPs in each region)

通过富集的基因进行GO注释和富集分析发现，两个群体受选择基因主要富集在泛素蛋白转移酶活性过程、核苷酸结合过程、对热的反应过程、磷酸蛋白酶活性过程、复制后修复过程、丝氨酸/苏氨酸蛋白激酶活性过程、谷胱甘肽转移酶活性过程、钾离子跨膜转运过程和丝氨酸型羧肽酶活性过程等通路。其中，位于2号染色体的候选区域锚定在10.63～10.64 Mb之间，基因Sspon.02G0053560-1C与野生土豆抗病蛋白RGA2同源[23]，说明抗性基因在选择性清除中受到选择。

3 讨论

3.1 野生种质是拓宽甘蔗遗传基础的有效途径

本研究通过简化基因组测序对107份甘蔗材料进行全基因组SNP标记的开发，获得7 869 726个高质量的SNP，并进行群体遗传结构分析。利用Admixture软件对甘蔗自然群体的群体结构分析表明，K=2时交叉错误率最小。因此，107份甘蔗材料划分成2个亚群，分别为栽培种群体和细茎野生种群体。该结果与进化树和主成分分析结果一致，说明甘蔗群体结构分析结果较可靠。值得注意的是，栽培种群体结构存在细茎野生种血缘的渗透现象。郎荣斌等[24]研究显示，甘蔗是多倍体植物且为非整倍体，遗传背景复杂，经过多年人工杂交，有较为明显的血缘渗透，大多数推广的栽培种中混有细茎野生种的遗传成分。本研究选用了83份栽培种和24份细茎野生种，之中的栽培品种是我国近50年育成及推广应用的品种，细茎野生种是我国甘蔗育种中应用最多的野生资源[25]。因此，该107份材料在一定程度上能反映我国甘蔗育成品种的遗传多样性。同时也说明利用甘蔗野生种质资源可以拓宽甘蔗遗传基础，是提高遗传多样性的有效途径[26]。

3.2 我国甘蔗栽培品种遗传多样性良好

本研究表明，栽培种群体的核苷酸多样性(π=8.93×10-6)高于细茎野生种群体(π=6.11×10-6)，可能是因为栽培种遗传背景较为复杂[27]。田春艳等[28]研究统计，我国现有甘蔗栽培品种中超过90%的甘蔗含有POJ2878的血缘，栽培品种中大多数的血缘来自于热带种、细茎野生种和印度种。通过甘蔗的“高贵化”育种[29]，将细茎野生种作为父本，与热带种进行杂交并回交以选育出优良的品种，说明大多数甘蔗栽培材料含有细茎野生种血缘，提高了栽培群体的核苷酸多样性。

3.3 抗性基因在甘蔗育种中受到选择

甘蔗栽培品种在选育和生产过程中，经常受干旱[30]、寒害[31]、黑穗病[32]和花叶病[33]等胁迫，一定程度上制约了甘蔗的生产。在甘蔗抗性性状的遗传改良方面，细茎野生种使甘蔗种质的开发取得突破性的进展。利用分子标记挖掘细茎野生种优良的抗性基因有利于提升甘蔗抗逆水平。通过对107份材料进行选择性清除分析发现，抗性基因在育种过程中持续受到选择。如沈怿丹[34]研究表明，GST家族在作物抗旱、耐盐、耐高温等环境胁迫上起重要作用，并且朱丹等[35]发现在干旱、高温胁迫中薄荷GSTU8基因上调表达； 殷龙飞等[36]研究发现，GRAS家族在调控植物信号传导的过程中具有重要作用，在干旱胁迫下作物茎叶部和根部调控SCL9基因的表达，SCL9基因是GRAS家族的成员之一。本研究中抗性基因在甘蔗育种过程中得到了强烈的选择，这能够为甘蔗抗逆分子育种提供候选基因，并为后续关联分析挖掘优良基因提供参考。