广丰千金薯和铁棍山药脱毒微型块茎的转录组分析

尹明华，白 丽，陈舒敏，程佳慧，冯丽文

(1.上饶师范学院 生命科学学院，江西 上饶 334001； 2.上饶农业技术创新研究院，江西 上饶 334001； 3.上饶市药食同源植物资源保护与利用重点实验室，江西 上饶 334001； 4.上饶市薯芋类作物种质保存与利用重点实验室，江西 上饶 334001)

广丰千金薯(Turczaninow. cv. Guangfeng Qianjin)是江西省上饶市广丰区传统名优地方山药品种，为薯蓣科薯蓣属多年生草本植物。广丰千金薯药食部位为其地下块茎，药用食用效果俱佳，肉质绵密粉嫩、软糯爽口，具有健脾、助消化、补虚劳、益气、祛痰、降脂、抗肿瘤、增加免疫作用等功效，被誉为“广丰人参”。广丰千金薯更是一种天然的纤体美食，纤维较多，食后易产生饱胀感，有控食减肥作用，闻名省内外，深受消费者青睐。有研究表明，广丰千金薯尤以广丰少阳地产的千金薯质量上乘、口感非凡、与众不同。本课题组采用随机扩增多态性DNA-聚合酶链式反应(RAPD-PCR)技术对来自江西的11份山药种质资源进行遗传多样性分析，证实广丰千金薯和铁棍山药(Turczaninow. cv. Tiegun，产地为河南温县)可以聚为一类，但广丰千金薯和铁棍山药之间的分子差异尚无相关报道。随着广丰千金薯需求量剧增，其规模化生产得到极大重视，而传统广丰千金薯多采用块茎的栽子进行营养繁殖，导致其病毒感染严重，产量逐年下降，品质不断退化。因此，脱除广丰千金薯病毒，改善广丰千金薯品质，提高广丰千金薯产量，已成为广丰千金薯生产中亟待解决的问题之一。利用植物茎尖培养技术进行脱毒快繁是培育脱毒苗的有效途径，但通过此种方法获得的脱毒苗植株细弱，移栽成活率偏低，且不利于包装和长途运输，这给脱毒苗的规模化应用带来了极大的局限。微型块茎是山药组织培养过程腋芽形成的变态小块茎，可在黑暗低温下长期保存，且体积小便于运输，可用于山药田间的规模化繁殖；利用脱毒苗获得的脱毒微型块茎更是山药优质种苗的重要来源。微型块茎的化学成分和山药地下块茎大致相同，药理与临床效果也相似，可抗氧化、降血糖、调节免疫、抗肝损伤、抗病毒等，具有重要的药用价值。

转录组测序分析已成为开展基因功能研究的重要方法。关于微型块茎的研究主要集中在酶基因克隆和序列分析、遗传转化体系建立、诱导形成生理生化、萌发因素、诱变处理等方面，而关于薯蓣类植物品种间微型块茎的转录组分析尚无报道。本研究以广丰千金薯和铁棍山药的脱毒微型块茎作为试验材料，对其进行转录组测序分析，探究广丰千金薯脱毒微型块茎生化成分代谢的分子机制，旨在为进一步利用广丰千金薯脱毒微型块茎进行遗传育种和品种鉴定提供参考。

1 材料与方法

1.1 材料

脱毒微型块茎是广丰千金薯(以下简称QJS，上饶市广丰区少阳乡)和铁棍山药(以下简称TGSY，河南温县)脱毒苗(由上饶师范学院上饶市薯芋类作物种质保存与利用重点实验室提供)在统一时间(60 d)诱导形成的。

1.2 方法

1.2.1 RNA提取和检测

采用植物总RNA 提取试剂盒(Invitrogen Ambion)提取QJS和TGSY脱毒微型块茎的RNA，用Agilent 2100 bioanalyzer精确检测RNA完整性。

1.2.2 转录组测序与分析

QJS和TGSY脱毒微型块茎的RNA交由北京诺禾致源生物科技有限公司利用Illumina HiSeqTM2500测序平台进行转录组测序和分析。转录组测序进行3次生物学重复。首先通过Oligo(dT)磁珠富集带有polyA尾的mRNA，随后在NEB Fragmentation Buffer中用二价阳离子将得到的mRNA随机打断，以片段化的mRNA为模版，以随机寡核苷酸为引物，在M-MuLV逆转录酶体系中合成cDNA第一条链，随后用RNaseH降解RNA链，并在DNA polymerase Ⅰ体系下，以dNTPs为原料合成cDNA第二条链。纯化后的双链cDNA经过末端修复、加A尾并连接测序接头，用AMPure XP beads筛选250～300 bp的cDNA，进行PCR扩增并再次使用AMPure XP beads纯化PCR产物，最终获得文库。文库构建后，用Illumina平台进行测序。测序的基本原理是边合成边测序(sequencing by synthesis)。经过原始数据过滤、测序错误率检查、GC含量分布检查，获得后续分析使用的clean reads再进行转录组重构，得到全部的转录本。

1.2.3 差异表达基因的筛选和富集分析

差异表达基因的筛选标准为差异倍数(FC)≥2，且错误发现率(FDR)<0.01。筛选出的差异表达基因的富集分析：GO(gene ontology)富集分析是基于描述基因功能的综合性数据库(http：//www.geneontology.org/)的分析；KEGG(kyoto encyclopedia of genes and genomes)富集分析是基于有关pathway主要公共数据库(http：//www.genome.jp/kegg/kegg2.html)的分析。

1.2.4 SNP、InDel和SSR分析

通过samtools和picard-tools等工具对比对结果进行染色体坐标排序、去掉重复的read等处理，最后通过变异检测软件GATK3分别进行SNP Calling和InDel Calling，并对原始结果进行过滤(过滤掉质量值小于40，距离小于2的SNP)；采用MISA软件(1.0版)对差异表达基因进行SSR检测和分析。

1.3 差异表达基因的实时荧光定量PCR(qRT-PCR)验证

根据1.2.3节的结果，找出与淀粉合成、块茎膨大、黄酮类代谢物相关的差异表达基因，以为内参基因，设计基因引物(表1)，进行qRT-PCR验证。PCR反应程序：95 ℃ 30 s；95 ℃ 10 s，60 ℃ 30 s，40个循环。使用2-△△C法计算基因相对表达水平。实验重复3次。

表1 差异表达基因的qRT-PCR引物

所有数据表示为平均值，并使用SPSS19.0软件进行统计分析，应用单因素方差分析(One-way ANOVA)检验与淀粉合成、块茎膨大、黄酮类代谢物相关的差异表达基因在QJS和TGSY脱毒微型块茎中表达的差异显著性。

2 结果与分析

2.1 测序数据产出统计

经过原始数据过滤、测序错误率检查、GC含量分布检查，获得后续分析使用的clean read(原始数据过滤后的read数)(表2)。由表2可知，QJS组raw read(原始数据中的read数)为22 686 888，clean read(原始数据过滤后的read数)为21 821 043，clean bases(原始数据过滤后的碱基数)为6.55 G，error rate(数据整体测序错误率)为0.03%，Q20和Q30(Phred数值大于20或大于30的碱基占总碱基的百分比)分别为98%和94.1%，GC含量为46%；TGSY组raw read为22 152 964，clean read为21 470 765，clean bases为6.44 G，error rate为0.02%，Q20和Q30分别为98.08%和94.3%，GC含量为45.59%。

表2 QJS组和TGSY组的测序数据产出统计

表3 QJS组和TGSY组的测序数据产出统计

2.2 拼接转录本长度分布

将Trinity拼接得到的转录本序列作为后续分析的参考序列。以Corset层次聚类后得到的最长cluster序列进行后续分析。转录本与聚类序列长度的统计结果见表3。由表3可知，300～500 bp的transcript数量和unigene数量分别为21 413和11 661，500～1 000 bp的transcript数量和unigene数量分别为25 033和9 924，1 000～2 000 bp的transcript数量和unigene数量分别为23 152和8 117，>2 000 bp的transcript数量和unigene数量分别为11 940和4 365，transcript和unigene的总数量分别为81 538和34 067。

2.3 参考序列比对

QJS组比对到参考序列上的read总数为43 642 086，能定位到参考序列上的测序序列数量和百分比为33 797 934、77.44%；TGSY组比对到参考序列上的read总数为42 941 530，能定位到参考序列上的测序序列数量和百分比为32 366 580、75.37%。

2.4 样本间相关性

样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接近1，表明样品之间表达模式的相似度越高。经过分析可知，QJS组和TGSY组样本间相关系数为0.42，表明QJS组和TGSY组2个样本间的相关性较低。

2.5 表达量统计

由图1-A可知，QJS组和TGSY组表达量FPKM的对数值在0～1.5，其中，QJS组表达量FPKM的对数值为0.2～1.2，TGSY组表达量FPKM的对数值为0.3～1.3，表明TGSY组基因的表达量比QJS组高。由图1-B可知，QJS组和TGSY组表达量密度在0～1.0，其中，TGSY组表达量密度为0～0.6，QJS组表达量密度为0～0.9，表明QJS组基因的表达量密度比TGSY组高。

图1 QJS组和TGSY组FPKM盒形图和FPKM密度分布图Fig.1 FPKM distribution Boxplo and FPKM density distribution map in QJS and TGSY group

2.6 差异表达基因数量统计

与TGSY组相比，QJS组有4 765个基因下调，有5 112个基因上调。图2为QJS组和TGSY组差异表达基因Venn图。Venn图中重叠区域表示2个品种间共有的差异表达基因数量。从图2可知，QJS组和TGSY组表达的共有基因数为25 207，QJS组单独表达的基因数为5 261，TGSY组单独表达的基因数量为3 571。

图2 QJS组和TGSY组差异表达基因Venn图Fig.2 Venn map of differentially expressed genes in QJS group and TGSY group

2.7 差异表达基因GO富集分析

GO(gene ontology)是描述基因功能的综合性数据库(http：//www.geneontology.org/)，可分为生物过程(biological process)、细胞组成(cellular component)和分子功能(molecular function)3个部分。共有4 528个差异表达基因富集到3 892个GO term，上调基因富集到3 194个GO term，下调基因富集到3 181个GO term(图3)，富集显著的前20个GO term见表4，包括核糖体的结构成分(structural constituent of ribosome)、结构分子活性(structural molecule activity)、细胞成分组织或生物发生(cellular component organization or biogenesis)、核糖体(ribosome)、细胞成分生物发生(cellular component biogenesis)、核糖体生物发生(ribosome biogenesis)等。

1，细胞成分组织；2，细胞成分与生物发生；3，核糖体；4，核糖体结构成分；5，结构分子活性。1， Cellular component organization； 2， Cellular component biogenesis； 3， Ribosome； 4， Structural constituent of ribosome； 5， Structural molecule activity.图3 QJS组和TGSY组差异表达基因GO富集分析柱形图Fig.3 Column chart of GO enrichment analysis of differentially expressed genes in QJS group and TGSY group

表4 部分差异表达基因GO富集结果

2.8 SNP和InDel分析

从表5可知，QJS组和TGSY组共发现611 498个SNP位点。SNP的碱基替换方式分为转换和颠换2种类型，其中碱基转换位点数量(237 083个，6.56%)显著少于颠换位点(374 415个，10.37%)。在转换类型中，C/T发生频率(30.783%)大于A/G(3.28%)；4种颠换类型中，A/T的发生频率最高，为1.15%；其次是A/C和T/G，分别为0.95%、0.97%；C/G发生频率最低，为0.73%。

表5 SNP部分分析结果

2.9 SSR分析

采用MISA软件对差异表达基因进行SSR检测，检测结果见表6。同时还对不同SSR类型在基因转录本的密度分布进行统计，结果如图4所示。从表6和图4可知，大多数基因具有单碱基重复(p1)、2个碱基重复(p2)、3个碱基重复(p3)和4个碱基重复(p4)，共发现12 056个SRR，其中包含1 179个混合型SSR和8 663个完整型SSR位点，表明广丰千金薯SSR标记的数量极其丰富。

表6 部分SSR检测结果

图4 SSR在转录本上的密度分布Fig.4 Density distribution of SSR on transcripts

2.10 差异表达基因的qRT-PCR验证

为了验证转录组测序差异表达基因表达的可靠性，选取与淀粉合成、块茎膨大、黄酮类代谢物相关的差异表达基因进行qRT-PCR验证。由表7可知，qRT-PCR测定的10个DEGs在QJS组和TGSY组中的表达水平趋势与转录组测序基因的表达趋势一致，表明本研究中得到的DEGs是可信的。从表7还可以看出，与TGSY组相比，QJS组的淀粉合成酶3、α-淀粉酶3、β-淀粉酶、类黄酮3′-单加氧酶、类黄酮3′-羟化酶、花青素合酶基因上调，而颗粒结合淀粉合酶2、蔗糖合酶1、扩展蛋白2和苯丙氨酸解氨酶等基因下调。

表7 差异表达基因的qRT-PCR验证

3 讨论

山药脱毒微型块茎可代替脱毒苗应用于大田生产，也可用于食疗和功能性食品，对一些疾病的治疗起到重要的辅助作用。研究表明，SNP和Indel位点广泛分布在基因组中。小果甜柿(Thunb.)的转录组检测到1 070 614个SNP位点。本试验共发现611 498个SNP位点。从理论上讲，转换的SNP位点与颠换的SNP位点比例一般为0.5，在本试验中，转换的SNP位点与颠换的SNP位点比例为0.63，与理论值较为接近，而小果甜柿转换的SNP位点与颠换的SNP位点比例为1.69，属于转换偏差，可能与进化选择相关，并非随机出现。SSR序列一般可以体现物种的进化水平，并影响基因的功能。小果甜柿果实的转录组中发现SSR种类十分丰富，以单个碱基、2个碱基和3个碱基为主。本试验结果与此一致。李府贡枣(Mill.)转录组的42 570条基因组数据库序列中发现18 016个SSR，其中包含1 442个混合型SSR和13 033个完整型SSR位点。在本试验中，广丰千金薯转录组共发现12 056个SSR，其中包含1 179个混合型SSR和8 663个完整型SSR位点，表明广丰千金薯 SSR标记的数量较为丰富。

目前发现的淀粉合成酶基因主要有3种，包括淀粉合成酶1、淀粉合成酶2和淀粉合成酶3，三者高度同源，区别在于N-端氨基酸长度存在差异。淀粉合成酶起到延伸支链淀粉和直链淀粉葡萄糖链的作用，颗粒结合淀粉合酶主要涉及直链淀粉的延长。在本试验中，广丰千金薯脱毒微型块茎的淀粉合成酶3基因上调、颗粒结合淀粉合酶2下调，表明与铁棍山药脱毒微型块茎相比，广丰千金薯脱毒微型块茎的淀粉合成酶3活性较强，有利于支链淀粉的合成，颗粒结合淀粉合酶2活性较低，直链淀粉的合成稍微偏弱，这可能是广丰千金薯脱毒微型块茎肉质绵密粉嫩、软糯爽口的内在原因。蔗糖合酶可参与蔗糖的运输和分配、淀粉和纤维素的合成、细胞壁组成、生物与非生物逆境等生理过程。在本试验中，广丰千金薯脱毒微型块茎的蔗糖合酶1基因下调，表明与铁棍山药脱毒微型块茎相比，广丰千金薯脱毒微型块茎的蔗糖合酶1活性较弱，不利于其蔗糖的运输和分配。α-淀粉酶是一种内切酶，可以将淀粉分解为麦芽糖、麦芽三糖、葡萄糖和少量大分子极限糊精，积极参与植物体的多种生命活动，响应外界的生物与非生物胁迫。β-淀粉酶负责储藏型淀粉的水解和过渡型淀粉的降解，产生β-极限糊精和β-麦芽糖，参与植物的生长发育和对非生物胁迫的响应。在本研究中，广丰千金薯脱毒微型块茎的α-淀粉酶3和β-淀粉酶基因上调，表明与脱毒微型块茎铁棍山药相比，广丰千金薯脱毒微型块茎的α-淀粉酶3和β-淀粉酶活性较强，有利于其淀粉转化为麦芽糖，在调控广丰千金薯脱毒微型块茎生长发育和非生物胁迫耐受过程中起到重要作用。油类黄酮3′-单加氧酶和类黄酮-3′-羟化酶均是类黄酮化合物生物合成途径中的关键酶，花青素合酶可将无色花青素转化为有色花青素，是花青素合成途径的关键酶之一，可参与低温、紫外线胁迫和植物病害防治。在本试验中，广丰千金薯脱毒微型块茎的类黄酮3′-单加氧酶、类黄酮-3′-羟化酶和花青素合酶基因上调，表明与铁棍山药脱毒微型块茎相比，广丰千金薯脱毒微型块茎的类黄酮3′-单加氧酶和类黄酮-3′-羟化酶活性较强，有利于其类黄酮化合物和花青素合成，在调控广丰千金薯脱毒微型块茎生长发育和防御过程中起到重要作用。植物扩展蛋白主要与植株抗逆性相关，在高温、高盐、干旱等非生物胁迫，以及病害、虫害等生物胁迫应答中具有重要作用。在研究中，广丰千金薯脱毒微型块茎的扩展蛋白2基因下调，表明与铁棍山药脱毒微型块茎相比，广丰千金薯脱毒微型块茎的扩展蛋白2活性较弱，打断半纤维素和纤维素氢键、疏松细胞壁结构能力不足，导致纤维素含量较丰富，这也是广丰千金薯脱毒微型块茎是一种天然的纤体美食、食后易产生饱胀感、可控食减肥的内在原因。苯丙氨酸解氨酶可以催化-苯丙氨酸或酪氨酸脱氨，产生反式肉桂酸，最终进入类黄酮代谢、木质素代谢、香豆素代谢等，可以抵御植物抗病原菌的入侵，免除病原生物的侵染。在本研究中，广丰千金薯脱毒微型块茎的苯丙氨酸解氨酶基因下调，表明与铁棍山药脱毒微型块茎相比，广丰千金薯脱毒微型块茎的苯丙氨酸解氨酶活性较弱，木质素和酚类等次生物质的合成下降，导致次生壁中的木质素含量下降，广丰千金薯脱毒微型块茎耐病程度降低。