黑果枸杞果实发育过程中转录组测序分析

孟小伟，牛 赟，2，马彦军

（1.甘肃农业大学 林学院，甘肃 兰州 730070；2.淮阴师范学院，江苏 淮安 223300）

高通量测序技术（Illumina测序技术）具有数据量大、准确性好、快速便捷、运行成本低等优点[1]。对于缺乏全基因组信息的非模式物种，利用Illumina测序技术能对研究对象在某一特定状态下基因转录表达信息进行全面地分析[2]，有利于挖掘重要功能基因，而且还可以更加便利地揭示特定生物学过程的分子机制[3]，这已成为研究植物优良性状的重要技术手段[4]。近年来，利用这一技术已成功地对多种药用植物在生理活动、生长发育以及次生代谢产物生物合成的相关基因进行了发掘和鉴定[5]。

黑果枸杞Lycium ruthenicum为茄科Solanceae枸杞属Lycium多年生灌木，是典型的盐生植物，也是一种药食两用的优良野生植物资源，在我国藏药和维药等传统民族医药中广泛应用[6]，主要分布于中国西北地区盐渍化土壤或荒漠环境区域，在荒漠区对维持生态平衡方面发挥着重要的作用，是盐碱地治理的先锋树种[7]。近年来黑果枸杞受到人们极大的关注，主要是由于黑果枸杞成熟果实中含有大量的花色素，其稳定性好，着色性强，是理想的食用天然花色苷[8]，同时黑果枸杞果实多糖，具有抗疲劳、降血糖等作用[9]，因此黑果枸杞是一种既有生态价值又有经济价值的树种。目前对黑果枸杞果实研究主要在黑果枸杞花色苷组成[10]、花青素含量及功能[11]、黑果枸杞多糖的功效[12]等方面，对黑果枸杞果实转录组及有关基因方面的研究较少[13]。本研究通过高通量测序技术，对黑果枸杞果实发育过程中青果期、变色期及成熟期进行测序，通过对黑果枸杞果实不同发育阶段测序数据进行分析，以期从整体水平了解黑果枸杞果实发育过程中基因表达模式和功能，为分析与黑果枸杞多糖和花青素等代谢相关的结构基因和调控基因及SSR分子标记的引物设计提供基础研究资料。

1 材料与方法

1.1 试验材料

供试材料为黑果枸杞不同发育阶段的果实，采自甘肃农业大学校内试验地（38°28′N，106°16′E），树龄为5 a。试验始于2017年6月26日，选取并标记长势优良、无病虫害的黑果枸杞植株（扦插繁育），采摘3个不同发育时期的果实，分别为青果期（S1）、变色期（S2）、成熟期（S3）。每3棵枸杞植株的果实为一个重复，共重复3次，果实采集后用铝箔纸包装，置于液氮中速冻，并立即放置到-80 ℃冰箱中保存，备用。

1.2 试验方法

1.2.1 总RNA的提取、文库构建和转录组测序

用mirVana™ miRNA ISOlation Kit (Ambion-1561)试剂盒提取和纯化黑果枸杞果实总RNA。分别用NanoDrop 2000（美国，赛默飞）与Agilent 2100 Bioanalyzer（美国，安捷伦）进行总RNA质量和纯度检测。将符合质量要求的黑果枸杞果实总RNA送样至上海欧易生物医学科技有限公司进行转录组测序。

1.2.2 De novo拼接

De novo拼接是指在无参考基因组的情况下，将有overlap的reads连接成一个更长的序列，通过不断的延伸，最后拼接出transcript。本研究采用Trinity（version:trinityrnaseq_r20131110）[14]软件paired-end的拼接方法得到Transcript序列，依据序列的相似性及其长度，筛选出最长的一条Transcript作为Unigene，以此作为后续分析的参考序列[15]。

1.2.3 Unigene功能注释

通过blastx将获得的Unigene序列分别与非冗余数据库（non-redundant，NR）、SwissProt和真核生物蛋白相邻类的聚簇（Clusters of orthologous groups for eukaryotic complete genomes，KOG）库进行比对，取e＜1e-5的注释基于BLAST（Basic Local Alignment Search Tool）寻找蛋白质或核酸的相似序列。有两个值，一个是S值，一个是E值。S值表示两序列的同源性，分值越高表明它们之间相似的程度越大。E值就是S值可靠性的评价。它表明在随机的情况下，其它序列与目标序列相似度要大于S值的可能性，所以它的分值越低越好[16]。得到跟给定Unigene具有最高序列相似性的蛋白，通过利用KAAS（http://www.genome.jp/kaas-bin/kaas_main）得到Unigene的京都基因与基因组百科全书（KEGG）注释信息，基于SwissProt的注释结果，根据Uniprot ID映射GO term，从而得到该Unigene的基因本体（gene ontology，GO）蛋白功能注释信息。

2 结果与分析

2.1 总RNA质量检测和组装结果分析

由表1可看出，黑果枸杞果实在3个发育时期9个样本，通过转录组测序共获得49019994～50537322原始读序（raw reads），各样本的clean data均达到6 Gb以上，Q30（测序错误率小于0.1%）均在89.85%以上，GC含量均在40%以上。利用Trinity软件对测序所得数据进行合并组装（表2），共获得43 573个Unigene，其中长度在1 kb上的Unigene有19 453条，这些Unigene可作为后续实验研究的重点对象；Unigene平均长度为1 262.65 bp。一般N50长度超过800 bp就可认为组装序列完整性较好，本实验所组装得到的序列N50为1 743 bp，这说明本次转录组测序数据量和质量都较高，可以用于后续分析。

表1 有效数据评估统计Table 1 Valid data evaluation statistics

表2 组装结果统计分析Table 2 Statistical analysis of assembly results

2.2 基因功能注释

由表3可知，黑果枸杞果实发育过程中43 573个Unigene中有23 723个Unigene在不同数据库中得到了注释，占总Unigene的54.44%，还有19 850个Unigene在这些数据库中没有得到注释，其中注释到NR数据库的Unigene最多，达到23 559个，占54.07%。共有3 726个Unigene在所有数据库中都得以注释（图1）。

表3 单基因序列注释统计Table 3 Statistics of single gene sequence annotation

2.2.1 Unigene的Nr及SwissProt数据库比对分析

由表3～4可看出，黑果枸杞果实转录组所得Unigene有23 559个Unigene在NR数据库中得到注释。在NR数据库注释中与马铃薯Solanum tuberosum同源序列最多，为4 586个，占注释Unigene的19.47%；与拟南芥Arabidopsis thaliana同源序列最少，为248个，占注释Unigened的1.05 %；与其他物种的为2 154个，占注释Unigene的9.14%。E值等于0的Unigene有2 993个，占注释到NR数据库Unigene的12.70%，E值在0～1e-100的Unigene最多，为9 024，占注释到NR数据库Unigene的38.30%。

图1 Unigene中NR,SwissProt,KOG,GO,KEGG注释交并集Fig.1 Unigene NR, SwissProt, KOG, GO, KEGG annotation intersection

由表3～5可看出，黑果枸杞果实转录组所得Unigene在SwissProt数据库中得以注释的Unigene有17 212个，其中E值为0注释的Unigene有3 393个，占19.71%；E值 介 于0到1e-100的Unigene有2 970个，占注释Unigene的17.26%。黑果枸杞果实转录组所得Unigene在SwissProt数据库注释中与拟南芥同源序列最多，为11 159个，占注释Unigene的64.83%；与爪蟾Xenopus laevis同源序列最少，为195，占注释Unigene的1.13%。

2.2.2 GO注释和分类

通过GO数据库对比分析，共有15 064个Unigene在GO数据库中得到了注释，获得的GO数据库注释的Unigene 可分为分子功能（Molecular function）、细 胞 组 分（Cellular component）和生物过程（Biological process）3大类别，进一步可细分为58个功能组。分子功能涉及20 197个Unigene，分为16个功能组，其中结合（Binding，8 959个）和催化活性（Catalytic activity，7 557个）含Unigene较多；细胞组分涉及51 077个Unigene，分为19个功能组，其中细胞（Cell，11 571个）、细胞部分（Cell part，11 561个）、细胞器（Organelle，8 998）涉及的Unigene较多；生物过程涉及56 167个Unigene，分为23个功能组，其中细胞进程（cellular process，10 157）和代谢进程（metabolic process，8 695）涉及的Unigene较多。

表4 黑果枸杞Unigene的NR比对分析Table 4 Nr mapping analysis of Lycium ruthenicum Unigene

表5 黑果枸杞Unigene的SwissProt比对分析Table 5 SwissProt mapping analysis of Lycium ruthenicum Unigene

表6 黑果枸杞的Unigene GO功能分类Table 6 GO functional categories of Lycium ruthenicum Unigene

2.2.3 KOG注释和分类

为了进一步分析黑果枸杞果实发育时期转录组Unigene的功能，进行了KOG功能分类分析，共有13 128个Unigene获得14 487个KOG注释，平均每条1.1个。分类结果如表7所示，共获得25个不同的功能分类。一般功能预测（General function prediction only）的Unigene为4 180个，是最大的功能类群；其次是翻译后修饰、蛋白质翻转和分子伴侣（Posttranslation almodification,protein turnover, chaperones）有1 391个，信号转导类机制（Signal transduction mechanisms）次之，有1 224条，最少的是细胞运动（Cell motility）功能类别，仅7个。黄酮类代谢途径所属的Q类（次生代谢产物生物合成、运输和代谢）共获得了520个Unigene注释。

2.2.4 KEGG代谢通路分析

KEGG是系统分析基因产物在细胞中的代谢途径以及基因产物功能的数据库。根据KEGG数据库的注释信息能进一步得到Unigene的pathway注释。结合KEGG数据库，黑果枸杞果实转录组注释到KEGG的4 951个Unigenes，获得了9 754个KEGG注释。注释的Unigene参与的代谢通路可分为4大类别23个子类。由表8可看出，4大类别代谢通路中，与代谢（Metabolism）相关的通路获得4 351个Unigene注释，遗传信息处理（Genetic information processing）相关的通路获得1 841个Unigene注释，细胞过程（Cellular processes）相关的通路获得1 694个Unigene注释，环境信息处理（Environmental information processing）相关的通路获得1 868个Unigene注释。进一步细分为23个子类代谢通路，其中信号传导（Signal transduction）获得Unigene注释最多，为1 833个，其次为碳水化合物代谢（Carbohydrate metabolism），为1 068个。次生代谢物生物合成（Biosynthesis of other secondary metabolites） 有321个Unigene。

表7 黑果枸杞Unigene的KOG功能分类Table 7 KOG functional categories of Lycium ruthenicum Unigene

以KEGG pathway数据库作为参考，可将注释到KEGG数据库中的9 754个Unigene定位到215个具体的代谢途径分支。表9列出注释到KEEG前10的代谢途径以及类黄酮生物合成、苯丙醇生物合成和花青素生物合成途径。由表9可看出，核糖体（Ribosome）代谢途径注释到的Unigenes数量最多，有263个Unigenes；其次为内质网蛋白加工途径（Protein processing in endoplasmic reticulum），为254个Unigene。注释到类黄酮生物合成途径中的Unigene有45个，苯丙醇生物合成途径的Unigene有133个，花青素生物合成途径的Unigene有1个，黄酮和黄酮醇的生物合成途径（Flavone and flavonol biosynthesis）Unigene有2个。

表8 黑果枸杞Unigene的KEGG功能分类Table 8 KEGG functional categories of Lycium ruthenicum Unigene

表9 黑果枸杞Unigene的KEGG代谢途径分析Table 9 KEGG analysis of metabolic pathway of Lycium ruthenicum Unigene

2.3 SSR信息分析

利用软件MISA对黑果枸杞果实发育过程中转录组测序所获得的Unigene进行SSR预测，结果见表10。由表10可知，共有16 815个SSR位点：单核苷酸SSR最多，为12 262个，占72.92%；6核苷酸SSR最少，为0.04%。重复单元重复出现的次数大于11次以上最多，为5 297，占31.50%；重复单元重复出现9次的最少，为287，占1.71%。

表10 SSR不同重复基序分布Table 10 Distribution of different repeat motifs in SSR

3 讨 论

第二代高通量测序技术因测序时间短、成本低和所获得数据量大等优点，被广泛应用于非模式生物分子生物学研究中。例如对文冠果[17]、杜仲[18]、蓝靛果忍冬[19]等非模式植物果实发育过程中转录组测序分析，全面了解了文冠果、杜仲、蓝靛果忍冬等果实基因表达情况。本研究利用第二代高通量测序技术对黑果枸杞果实发育不同阶段进行转录组测序，建立了黑果枸杞的转录组数据库，获得了大量基础数据。对这些数据进行分析、序列组装，所得序列在不同数据库的功能注释及分类、代谢途径等分析，揭示黑果枸杞果实不同发育阶段整体基因表达特征，为黑果枸杞果实发育过程中的分子生物学研究提供了基础资料。

本研究所得到的黑果枸杞果实发育过程中的43 573个Unigene有23 723个Unigene在NR、SwissProt、GO、KOG、KEGG等不同数据库中得到了注释，占总Unigened的54.44%，还有19 850个Unigene在这些数据库中没有得到注释。这一结果在许多非模式生物转录组测序中都存在[20]，这主要是由于非模式生物缺乏基因组方面研究的基础资料，使得部分Unigene在NR、SwissProt、GO、KOG、KEGG等数据库中无法得以注释，而这些未得到注释的Unigene有可能是非模式生物特有的基因，因此需要对这些未得到注释的Unigene进一步从结构、功能等方面来深入研究，从而揭示他们在该生物生长发育过程中所发挥的作用。由于黑果枸杞现有遗传信息量少的原因，测序结果中所得到的Unigene即使能在NR、SwissProt、GO、KOG、KEGG等数据库中得以注释，也有一部分Uuigene是注释到动物和人类基因数据库中。本研究中SwissProt注释时就有2.04%的Unigene与智人同源，1.55%的Unigene与小鼠同源，还有部分Unigene同黑腹果蝇和爪蟾同源。在KEGG代谢通路pathway分析中，有部分Unigenere与人类疾病路径相关。这一结果在许多其他生物测序中都存在[21]，这主要是因为这些非模式生物现有EST数据很少，而为了确定这些非模式生物在测序过程中获得大量Unigene的功能，只能与其他植物甚至动物和人类的基因数据库进行比对，因此产生了这样的结果。

本研究对黑果枸杞果实发育过程中3个阶段的43 573条Unigene进行SSR搜索，得到16 815个SSR位点，出现频率为38.59%，高于野三七（频率为16.86%）[22]、腊梅（频率为12.35%）[23]、红豆杉（频率为2.07%）[24]等。这表明不同植物转录组测序结果中SSR的发生频率差异较大，产生这种差异可能与所研究物种、转录组测序所得数据量以及检索标准等因素有关[25]。在黑果枸杞果实转录组中，SSR以单核苷酸重复基序（72.92%）为主，这与‘芙蓉李’研究中SSR主要类型为单核苷酸重复基序（42.19%）为主一致[26]。目前大量研究结果表明利用转录组数据开发SSR标记是可行的[27]。本研究基于黑果枸杞果实发育过程中转录组的高通量测序数据，从RNA水平上有针对性地进行了黑果枸杞特异性SSR位点的检索和评价，为进一步开发新的黑果枸杞功能基因及调控基因的SSR分子标记、黑果枸杞遗传多样性研究、黑果枸杞遗传图谱的构建、分子标记辅助育种等研究提供了基础资料。

4 结 论

本研究利用Illumina HiSeqTM2500测序仪对黑果枸杞果实发育过程中转录组进行测序，拼接出Unigenes 43 573条。注释到GO、KEGG、KOG、NR和Swiss-Prot数据库上的Unigenes总数为23 723条，占总Unigened的54.44%，共有3 726个Unigene在所有数据库中都得以注释，还有19 850个Unigene在这些数据库中没有得到注释。获得的GO数据库注释的Unigene可分为细胞组分、分子功能和生物过程3大类别，进一步可细分为62个功能组。KOG功能分类中共获得25个不同的功能分类，黄酮类代谢途径所属的Q类（次生代谢产物生物合成、运输和代谢）共获得了520个Unigene注释。黑果枸杞果实转录组注释到KEGG的4 951个Unigenes获得了9 754个KEGG注释，参与的代谢通路可归为4大类别23个子类。将注释到KEGG数据库中的9 754个Unigene定位到215个具体的代谢途径分支，其中注释到类黄酮生物合成途径中的Unigene有45个，苯丙醇生物合成途径的Unigene有133个，花青素生物合成途径的Unigene有1个，黄酮醇的生物合成途径的Unigene有2个。在黑果枸杞果实发育过程的转录组中发现16 815个SSR位点，最多的为单核苷酸SSR，占72.92%，重复单元重复出现的次数大于11次以上最多，占31.50%。