冬瓜高通量转录组测序及分析

叶新如，张前荣，陈敏氡，王 彬，刘建汀，朱海生，温庆放

(福建省农业科学院 作物研究所，福建 福州 350013)

冬瓜(BenincasahispidaCogn.)原产于中国南部、东南亚以及印度等地，分布于亚洲的热带、亚热带及温带地区[1]，在中国有1 500多年的栽培历史。其作为我国重要的蔬菜品种，对调节蔬菜淡季品种稀少、保证蔬菜周年供应有着非常重要的作用。冬瓜产量高，花、叶、果皮、果肉、种子等都具有很好的药食兼用价值[2]，能够用于治疗咳嗽[3]、消化性溃疡[4]、高血压、炎症[3]等。前人对冬瓜的研究局限在栽培技术[4-5]、药用开发[6-9]、营养价值[10-11]、遗传育种[12-14]等方面，而关于冬瓜基因组学的研究极少，有关冬瓜重要经济性状分子机制的研究进展缓慢。

转录组是特定的细胞或组织在某一发育阶段或功能状态下转录出的所有RNA总和[15]。目前，基于高通量测序的转录组分析成为功能基因研究的主要方法[16-17]。该技术也被广泛应用于植物上，如拟南芥[18]、玉米[19]、香椿[20]、紫色黄秋葵[21]，但在冬瓜上鲜有关于转录组测序的报道。

本研究以冬瓜嫩叶为材料，利用Illumina HiSeqTM2000技术对其进行转录组测序，获得冬瓜转录组序列、遗传变异等信息，从中挖掘冬瓜的基因数据和SSR分子标记，以期为后续研究提供数据支撑。

1 材料与方法

1.1 试验材料

试验材料为“黑皮冬瓜”，统一播种于福建省农业科学院作物研究所福清市东张镇蔬菜科研基地。采用穴盘育苗，待植株长至三叶一心时，取成型的嫩叶3片，等量混样进行转录组测序，测序由广州基迪奥生物有限公司完成。

1.2 冬瓜总RNA提取

冬瓜总RNA提取采用Trizol Reagent法，并对RNA纯度、品质进行检测。若提取的RNA序列满足如下条件： OD260/OD280在1.8～2.2，OD260/OD230≥2，28S/18S≥1，RNA完整值(RIN)≥7.5，则认为RNA符合建库测序标准。

1.3 冬瓜转录组测序

提取冬瓜总RNA后，用带有Oligo(dT)的磁珠富集冬瓜mRNA。以mRNA为模板，用六碱基随机引物(random hexamers)合成第一条cDNA链；然后加入缓冲液、dNTPs、RNase H和DNA polymerase Ⅰ合成第二条cDNA链；经QiaQuick PCR试剂盒纯化及EB缓冲液洗脱后做末端修复，加poly(A)并连接测序接头；最后进行PCR扩增，建好的测序文库用Illumina HiSeqTM2000进行测序。

1.4 冬瓜转录组序列的拼接与组装

利用Illumina HiSeqTM2000测序平台对建立好的冬瓜转录组文库进行测序，得到原始序列(raw reads)，其中含有低品质、重复测序以及接头污染的序列，因此必须进行严格过滤，去除低品质序列，获得干净序列(clean reads)；统计干净序列的长度、数量、GC含量等。使用组装软件Trinity进行拼接，将有重叠信息(overlap)的序列连接成长片段(conting)，通过不断地延伸，拼接组装成转录本(Transcript)，然后使用Tgicl软件对转录本进行去冗余聚类，最终获得单基因簇(Unigene)。利用RPKM法对单基因簇表达量进行计算分析，其计算公式为：

RPKM=C×106/(NL/103)。

式中:C为比对到基因的reads数，N为比对到所有基因的总reads数，L为基因的碱基数。

1.5 冬瓜单基因簇的基本功能注释

通过Blastx[22]将冬瓜的单基因簇序列与公共数据库(E值<1×10-5，E值为单基因簇与数据库中匹配序列为同源序列的假阳性概率)进行比对，得到与给定单基因簇具有最高序列相似性的蛋白，从而获得该单基因簇的蛋白功能注释信息。公共数据库包括非冗余蛋白数据库(nonredundant protein database，Nr)、蛋白质序列数据库(Swiss Prot protein database，Swiss Prot)、基因本体论数据库(gene ontology，GO)、蛋白质真核同源数据库(eukaryotic orthologous groups，KOG)、东京基因与基金组百科全书(Kyoto encyclopedia of genes and genomes，KEGG)、蛋白质家族域数据库(protein families database，Pfam)。

通过Blastx软件将冬瓜单基因簇序列比对至Nr数据库，可获得冬瓜单基因簇对应的蛋白编码序列。将冬瓜单基因簇与Swiss Prot数据库比对，可确定冬瓜基因簇的功能。利用Blast2 GO和WE GO软件对冬瓜单基因簇序列进行GO注释及功能分类统计，同时与KOG数据库比对注释，可统计冬瓜基因功能和分类信息。在KEGG数据库中对冬瓜单基因簇进行相似性分析，可获得冬瓜单基因簇代谢通路的注释和预测。利用Hmmer 3.0将冬瓜单基因簇与Pfam蛋白结构数据库进行比对，可获得冬瓜单基因簇的注释信息。

1.6 冬瓜转录因子分析及SSR位点分析

将冬瓜蛋白序列与TF 数据库(plant TFdb/animal TFdb)进行 hmmscan 比对，可获得相应转录因子家族，以进一步进行转录因子分析。使用软件 MISA对转录组所有单基因簇中的SSR进行搜索，并对SSR重复基序进行特征分析，为获得合适的引物序列提供数据支撑。

2 结果与分析

2.1 冬瓜转录组测序结果

2.1.1 转录组测序及组装后不同长度单基因簇的分布 用高通量测序平台对冬瓜嫩叶测序，共获得63 525 954条原始序列，碱基数据量9 528 893 100 bp(9.5 Gb)。去除低品质的原始序列后，获得62 021 032条高品质干净序列，占原始序列的97.63%。碱基数据量9 154 299 990 bp(9.1 Gb)，G+C含量为46.88%，碱基Q20(即测序错误率小于1%的比例)和Q30(即测序错误率小于0.1%的比例)分别为98.31%及94.96%。测序结果表明数据合格，可以用于后续分析。

利用Trinity软件对干净序列进行组装，获得40 611条单基因簇，删除7条比对到核糖体的基因后，对剩余的40 604条单基因簇进行统计及注释，G+C含量为41.64%，覆盖50%核苷酸的最大单基因簇长度(N50)为1 591 bp，最短和最长单基因簇分别为201及13 773 bp，平均长955 bp，不同长度单基因簇的分布情况见图1。由表1可知，长200～999 bp的单基因簇占总数(40 604条)的66.96%，长 1 000～1 999 bp的单基因簇占比为20.93%，长2 000～2 999 bp的单基因簇占比为7.95%，长度≥3 000 bp的单基因簇占比为4.16%。总体而言，随着长度增加，单基因簇数量减少。

图1 冬瓜单基因簇的长度分布Fig.1 Length distribution of Benincasa hispida unigenes

2.1.2 单基因簇的总体表达量 测序得62 021 032条干净序列，能比对到单基因簇序列的 reads有56 708 038条，占干净序列总数的91.43%，其中比对到唯一指定单基因簇序列的 reads有 49 628 128条，占80.01%。利用RPKM 法计算单基因簇的表达量，结果表明，冬瓜40 611条单基因簇序列RPKM平均值为16.27，最大值51 049.34(Unigene-0022103)，最小值0.011 4(不包含0，RPKM值为0的单基因簇共300条)。

2.2 冬瓜单基因簇的总体注释结果

通过Blastx将测序获得的40 604条冬瓜单基因簇与相关数据库进行比对，结果(表1)显示，Nr数据库注释到最多的单基因簇，共27 474条，占总单基因簇的67.66%；GO数据库注释到的单基因簇最少，共10 659条，占26.25%；Swiss Prot、KOG、KEGG、Pfam分别注释到19 572，16 452，10 799和17 990条单基因簇。

表1 冬瓜单基因簇在6大数据库中的注释结果Table 1 Statistics of database annotations of Benincasa hispida unigenes

2.3 冬瓜单基因簇的Nr功能注释

冬瓜40 604条单基因簇在Nr数据库中的比对结果(图2-A)显示，其中27 474条单基因簇在Nr数据库中比对到相似序列， 0≤E≤10-150的数量最多，有10 367条，占27 474条单基因簇的37.73%，其中E值为0的单基因簇有8 752条(占31.85%)，10-150

图2 冬瓜单基因簇在Nr数据库中的E值分布(A)和物种分布(B)Fig.2 E-value (A) and species distribution (B) of Benincasa hispida unigenes in Nr database

2.4 冬瓜单基因簇的Swiss Prot统计分析

在Swiss Prot数据库中，共有19 572条单基因簇成功注释到，占40 604条单基因簇总数的48.20%，匹配序列相似度(identity)≥80%的单基因簇有14 456条，匹配序列相似度在30%～80%的单基因簇有222 482条，匹配序列相似度≤30%的单基因簇有37 282条。

2.5 冬瓜单基因簇的GO功能分类

GO功能注释是根据Nr注释信息得到的，是对基因功能进行分类的体系。冬瓜单基因簇的GO功能分类结果见表2。

表2 冬瓜单基因簇GO功能分类Table 2 GO functional categories of Benincasa hispida unigenes

从表2可以看出，冬瓜40 604条单基因簇中有10 659条单基因簇被注释到GO数据库内，分别匹配到生物功能、分子功能和细胞组分3个本体的47个功能组中。根据GO功能分类，注释到生物功能的有25 093条单基因簇(占3个本体总数(52 211条单基因簇)的48.06%)，注释到分子功能的有11 319条单基因簇(占21.68%)，注释到细胞组分的有15 799条单基因簇(占30.26%)。其中生物功能包含19个功能组，以代谢进程(6 156条)、细胞进程(5 711条)和单一有机体进程(4 344条)占比较大，分别占25 093条生物功能单基因簇的24.53%，22.76%和17.31%，而节律进程(12条)、脱毒(11条)、细胞定位(3条)注释到的单基因簇较少；分子功能包含12个功能组，催化活性(6 029条)和结合活性(4 314条)占比较大，分别占11 319条单基因簇的53.26%和38.11%，电子载体活性(12条)、转录因子活性(11条)、翻译调控活性(1条)比例较小；细胞组分包含16个功能组，细胞、细胞部分和细胞器占比较大，细胞、细胞部分各注释到3 541条单基因簇，细胞器注释到2 655条单基因簇，胞外基质、胞外区域部分和超分子纤维注释的单基因簇较少，均仅为2条。

2.6 冬瓜单基因簇的KOG功能分类

将冬瓜单基因簇比对到KOG数据库中发现，共有16 452条单基因簇注释到KOG数据库，共获得27 944条单基因簇，被注释到25种KOG分类中，功能区中单基因簇数量分布有所差异(表3)，其中一般功能基因注释到的单基因簇数量最多，包含5 975条单基因簇，占27 944条单基因簇的21.38%；其次是翻译后修饰、蛋白质折叠和分子伴侣共包含3 112条单基因簇，占11.14%；信号传导机制(2 773条)、RNA加工和修饰(1 605条)、转录(1 469条)、膜泡运输和胞内分泌(1 460条)注释到的单基因簇也较多，最少的是细胞活性，仅有11条单基因簇(占0.04%)。

表3 冬瓜单基因簇KOG功能分类Table 3 KOG functional categories of Benincasa hispida Cogn unigenes

2.7 冬瓜单基因簇的KEGG代谢通路分析

将测序得到的40 604条冬瓜单基因簇比对到KEGG数据库，得到单基因簇的Pathway注释(表略)。结果表明，注释匹配到的10 799条冬瓜单基因簇，可分为5个大类、19个亚类、125条代谢途径。5个大类中代谢注释到的单基因簇最多，共5 868条，占带有通路注释的9 623条单基因簇的60.98%；遗传信息处理次之，注释匹配到2 572条单基因簇；生物系统最少，仅211条。进一步分析冬瓜125条代谢途径发现，核糖体的单基因簇数量最多(363条)，碳代谢和植物激素信号转导居次，分别有337条和311条单基因簇。在KEGG数据库功能注释中预测到的125条代谢途径中，含有氨基酸合成、不饱和脂肪酸的生物合成、维生素B6代谢、类黄酮生物合成等相关的代谢途径分别有295，47，20和18条单基因簇，其中部分单基因簇涉及不饱和脂肪酸、类黄酮、维生素等生理活性物质合成，这些物质可能与冬瓜的营养保健功效有关。

2.8 冬瓜单基因簇的Pfam统计分析

利用Pfam数据库对冬瓜蛋白结构域进行预测，对蛋白家族进行分析。将冬瓜40 604条单基因簇与Pfam数据库比对，其中有17 990条单基因簇成功注释，所有注释到的单基因簇编码蛋白序列在数据库中匹配结构的数目均为1个，分成369个类群，三角状五肽重复结构域(1 096条)、三角状五肽蛋白家族(907条)、富亮氨酸重复(769条)、WD重复域的G-beta蛋白(641条)、蛋白激酶结构域(568条)在整体中占比较大。

2.9 冬瓜单基因簇的转录因子分析

对冬瓜的单基因簇进行转录因子分析(表4)，在40 604条单基因簇中，具有转录因子序列特征的有1 099条，分为55个转录因子家族。家族中所占比例最大的是bHLH(81条，占1 099条单基因簇的7.37%)，然后依次是ERF(74条，占6.73%)、C2H2(73条，占6.64%)、MYB_related(70条，占6.37%)。转录因子是读取并解释DNA中遗传“蓝图”的蛋白质组之一，其与DNA结合，可以帮助启动一个增加或减少基因转录的程序，因此转录因子对于许多细胞过程是至关重要的。前人研究表明，bHLH、ERF在植物抗逆境胁迫和生长发育、生物合成及信号转导方面均有作用[23-26]。

表4 冬瓜转录因子家族的单基因簇数量及种类Table 4 Number and species of unigenes in the transcription factor family of Benincasa hispida Cogn

2.10 冬瓜转录组的SSR位点分析

对冬瓜转录组进行测序分析，共获得40 604条单基因簇，其中5 086条单基因簇获得5 474个SSR位点(表5)。SSR的发生频率为12.53%，出现频率13.48%，平均7.08 kb出现一个SSR。冬瓜转录组SSR重复单元次数在4～30次，类型丰富，其中4～10次重复的SSR有5 127个，占5 474个位点的93.66%；11～30次重复的SSR有347个，占6.34%。

由表6可以看出，本研究中SSR有二、三、四、五、六核苷酸5种重复类型，共含170种重复基序，二、三、四、五、六核苷酸重复基序数量分别为4，10，24，42和90个。二核苷酸重复基序中AG/CT发生频率最高，占5 474个位点的26.21%，占二核苷酸重复总数(2 047个)的70.10%；三核苷酸重复基序中AAG/CTT占比较大，在重复基序中占24.08%；而五核苷酸和六核苷酸重复基序数量较少，总体占比也较小。

表5 冬瓜叶转录组的SSR类型与数量Table 5 Number of SSR types in the transcription group of Benincasa hispida

表6 冬瓜转录组的SSR分类重复类型频率Table 6 Frequency of SSR categories repetition in the transcription group of Benincasa hispida

3 讨 论

冬瓜营养丰富，含有多种生理活性成分。近年来有研究表明，冬瓜的多种提取物具有抗氧化[27]、治疗荨麻疹[28]、降压利尿[29]、降糖[30]等功效。冬瓜优良的药用价值和食疗保健价值使其市场需求量旺盛。然而，目前关于冬瓜相关代谢途径的机理研究还未见报道，本研究利用Illumina HiSeqTM2000技术进行转录组测序，通过拼接获得了大量冬瓜单基因簇，为冬瓜后续相关功能基因研究奠定了基础。

本研究通过高通量测序平台对冬瓜转录组进行测序及分析发现，覆盖50%核苷酸的最大单基因簇长度为1 591 bp，表明组装质量较好，适合进行下一步分析。组装后获得40 611条单基因簇，平均长955 bp，对其中40 604条单基因簇进行注释，仍有13 004条单基因簇未成功注释，占总数的32.03%，这可能与序列片段过短、注释信息缺乏以及初次从冬瓜中获得新基因等有关[31]。用GO数据库对冬瓜单基因簇进行功能分类，匹配到3个本体的47个功能组中，其中生物功能中代谢过程注释到的单基因簇数量最多，共6 156条。从功能注释结果可以看出，冬瓜在基因表达水平上与代谢过程相关的表达较多，表明植物存在着不同代谢途径，且不同的代谢途径同时在植物体内与相应的生理过程联系密切。KEGG代谢通路分析得到单基因簇的Pathway注释，显示成功注释了10 799条冬瓜的单基因簇，其中代谢途径占5 868条，代谢途径中次生代谢途径有268条。同时，KEGG数据库注释将冬瓜单基因簇分为125条代谢途径。用Pfam数据库对冬瓜蛋白结构域进行预测、对蛋白家族进行分析，结果显示注释到的17 990条单基因簇可分成369个类群，匹配结构数目均为1个。

转录因子能够调控基因的表达水平及调控网络。本研究与TF数据库进行hmmscan比对，注释到1 099条单基因簇分为55个转录因子家族，占比最大的是bHLH家族，其次是ERF、C2H2、MYB_related、NAC、MYB家族。bHLH是真核生物中存在最广泛的一大类转录因子，Cui等[23]发现其在花药发育调控方面具有重要作用，何洁等[24]认为其在西瓜抵御低温、ABA和盐胁迫中有可能发挥作用，Gonzalez等[32]认为TTG1/bHLH/Myb转录因子对拟南芥幼苗花青素生物合成途径具有调控作用。另外，Jung等[33]研究表明，bHLH转录因子可能与水稻雄性不育有关。转录因子在植物的生长发育[34]、进化[35]、抗逆境[36]等方面具有重要的调控作用，这为冬瓜后续在生长和抗逆性的研究上提供了有力的依据。

冬瓜资源较为丰富，利用表型观察虽然能够直观地确定作物的产量、品质以及相关的性状特征表现，但同时也受到很多条件的影响与限制，因而需要进一步结合分子标记进行亲缘关系鉴定。本研究通过转录组测序得到冬瓜5 474个SSR位点，SSR发生频率为12.53%，出现频率13.48%，低于丝瓜[37]和莲雾[38]，但高于南瓜[39]和黄秋葵[40]，这可能与物种、单基因簇总量及长度、筛选条件等因素不同有关[41]。冬瓜SSR有二、三、四、五、六核苷酸5种重复类型，其中二核苷酸重复类型和三核苷酸重复类型占比较大，而五核苷酸和六核苷酸重复占比较小。这与丝瓜[37]、莲雾[38]、南瓜[39]和黄秋葵[40]的结果相似，可能与SSR搜索标准有关。冬瓜中AG/CT和AAG/CTT为主要重复基序，AG/CT为主要基序，这与前人对平欧杂交榛[42]和菜薹[43]的研究结果一致；AAG/CTT在双子叶植物中是常见基序，与Morgante等[44]的研究结论相似。本研究结果可为后续的冬瓜研究提供参考数据。