低温胁迫下枇杷幼果转录组的De novo组装和功能注释

杨 伟，龚荣高，，石佳佳，赖 静，廖明安，梁国鲁

(1 四川农业大学 园艺学院，四川 雅安625014；2 成都龙泉驿区农村发展局，四川 成都610100；3 西南大学 园艺园林学院，重庆400700)

枇杷是调节市场淡季的重要亚热带水果，其以花和幼果过冬，开花坐果期正值一年中温度最低的时期，因而枇杷的抗寒性已成为影响其产量和品质的重要指标之一。目前，有关枇杷抗寒性方面的研究对象相对单一，大多数集中在叶片冻害后的结构、内含物及光响应等方面，而针对果实的研究相对较少[1]，尚未见有关枇杷幼果转录组测序等方面的研究报道。由于转录组表现的是生物或组织细胞在特定状态下的所有RNA的总和，它包含了时间和空间的限定，是连接基因组遗传信息与生物功能蛋白质组的必然纽带，因而，转录组的测序研究成为探究基因结构及功能的基础。运用新一代的Illumina测序技术，可以得到基因组转录区域信息、鉴定转录发生位点、进行可变剪切等，其精确的计数方法更可对基因进行精确的定量分析，已广泛应用于生物学研究领域[2-3]。本研究通过对低温胁迫下枇杷幼果转录组的分析，探索枇杷幼果低温冻害的机理，并挖掘枇杷抗寒性相关基因，以期为枇杷引种、生态区划和优质高产栽培提供相关的理论和实践依据。

1 材料与方法

1.1 材 料

在四川省成都市龙泉驿区园艺场，选取10～15年生树势健壮、常规栽培管理、生长相对一致的5株“大五星”枇杷(EriobotryajaponicaLindley cv.Dawuxing)作为代表性植株,砧木为“解放钟”(E.japonicaLindley cv.Jiefangzhong)，盛花后开始挂牌，并在易受冻害的花后75 d(75 days after full bloom，75 DAFB；正常发育果的直径约1 cm)后，取树冠中部外围带无病虫害幼果的果枝20支，其中10支置于3 ℃的人工气候箱内低温胁迫处理12 h，另外10支置于12 ℃的人工气候箱内作为对照。

1.2 方 法

1.2.1 RNA 的提取及质量控制 枇杷幼果总RNA的提取参照Garg等[4]的方法。以幼果中心柱为对称轴纵剖，在可食部分称取0.3 g样品置于液氮中，研磨成粉末状。RNA样本的质量和浓度采用Garg等[4]的方法评判，用紫外分光光度计测定OD230、OD260和OD280，要求OD260/OD280和OD260/OD230的比值均≥1.8。运用Agilent 2100 生物分析仪检测RNA的完整性，RNA完整性值在8.6～10.0 的用于进一步试验。即将样品用RNase Free Water 在冰上溶解后，离心并充分混匀，取1 μL 样品于70 ℃变性2 min后，利用Agilent 2100生物分析仪进行检测。

1.2.2 文库制备与检测 经过DNase 消化、mRNA分离、mRNA 打断，进行反转录一链的合成，随后进行反转录二链的合成和末端修复。在末端修复产物cDNA 的3′末端加上A 碱基，让接头与A 碱基相连接。连接产物经胶回收后进行PCR 反应及产物回收。使用Agilent 2100 生物分析仪和ABI实时荧光定量PCR系统对文库质量进行检测。

1.2.3 测序分析及Denovo组装 由深圳华大基因研究中心完成转录组测序及Denovo组装。应用Illumina Hiseq 2000测序平台，转录组文库测序为PE 90，reads为180 bp。Denovo组装方法如下：使用短reads组装软件SOAPdenovo进行转录组从头组装，得到Unigene。将Unigene序列与蛋白数据库nr、Swiss-Prot、KEGG和COG做blastx比对(E-value<0.000 01)，确定Unigene的序列方向。如果不同库之间的比对结果有矛盾，则按nr、Swiss-Prot、KEGG和COG的优先级确定Unigene的序列方向；与以上4个库皆比不上的Unigene，用软件ESTScan预测其编码区并确定序列的方向，给出其从5′到3′方向的序列；对于无法确定序列方向的Unigene，给出组装软件得到的序列。

1.2.4 基因功能注释 由深圳华大基因研究中心完成基因测序及其功能注释。功能注释信息给出Unigene的蛋白功能注释、COG功能注释、Pathway注释，以及GO(GO；http://www.geneontology.org)分类。

1.2.5 低温胁迫下基因的表达 由深圳华大基因研究中心完成差异基因的表达分析。方法如下：Unigene表达量的计算采用FPKM法，其表达量用来比较低温胁迫处理后的基因表达差异。随后进行差异Unigene的GO和Pathway分析。

2 结果与分析

2.1 枇杷幼果cDNA序列分析与组装

运用Illumina测序平台对枇杷幼果cDNA样品进行测序，以弄清幼果的转录组及基因状况。由表1可知，枇杷幼果中总体产生30 911 414条原始序列数据(reads)，经过严格的质量评估和数据过滤，得到26 024 270条clean reads。质量要求Q(质量值)>20的碱基百分比超过95%且G+C含量为35%～65%；数据过滤方法为去除测序接头序列、重复冗余序列(含N比例大于10%)和质量差的序列(Q≤10 的碱基数占整条read 的50%以上)。过滤后转录物Q>20的碱基百分比为97.47%，对其作进一步分析。过滤后碱基G和C的数量占总碱基数的比例为20.3%～70.71%，平均为48.94%。

表 1 枇杷幼果转录组测序数据质量统计结果

从表2可知，对reads进行拼接，共组装 116 723 条Contigs，Contig的平均长度为312 nt(核苷酸)，N50的长度则为487 bp(碱基数)。组装的Contigs长度为200～7 653 nt(图1A)。

表 2 枇杷幼果基因序列组装统计结果

图 1 枇杷幼果基因序列组装的Contigs和Unigenes的长度分布

由表2还可见，对Contigs进一步拼接后，共组装获得转录物Unigenes序列64 814条，Unigenes的平均长度为541 nt，其N50的长度则是797 bp。Unigenes的长度为200～8 895 nt，总长度为35 036 073 nt(图1B)。此外，长度大于1 000 nt的Contigs和Unigenes序列分别占到23.32% 和 13.79%。

2.2 枇杷幼果基因序列功能注释信息

为了整体上了解转录物(Unigenes)的序列信息，将所有Unigene与nr、Swiss-Prot、KEGG、COG等蛋白或基因数据库进行blastx比对,以得到基因的方向、阅读框及基因的功能注释，并给出对应的与数据库中基因同源程度的衡量标准(E-value值)。利用nr、Swiss-Prot等蛋白数据库对64 814条转录物数据进行blastx比对的结果表明，有42 830条转录物序列被NCBI nr数据库功能注释。从图2A的E-value分布图可知，有42.5%的转录物序列具有极强的同源比对信息，而且有57.5%的转录物序列相似性E-value值在1.0E-5～1.0E-60。由图2B的Identity分布图可见，与具体蛋白匹配相似性超过80%的转录物占29.7%，而与具体蛋白匹配相似性范围在19%～80%的转录物占70.3%，相似性低于40%的转录物占5.9%。由图2C物种分布图可见，有32.5%的转录物序列与葡萄序列高度匹配，与蓖麻匹配的占18.4%，与毛果杨匹配的占16.4%，与大豆匹配的占11.3%，与苹果匹配的占4.1%，与蒺藜苜蓿匹配的占3.2%，与鼠耳芥匹配的占1.1%。

图 2 枇杷幼果组装转录物与nr数据库比对后的分类分布图

根据nr注释、GO分析对枇杷幼果基因功能进行可能的分类。由图3可知，在GO功能分类体系中，有34 200条转录物具有具体功能定义，并获得了273 645条转录物的功能注释，其中有143 799条转录物注释归类为生物学过程，有90 020条归为细胞组分，有39 826条归为分子功能。在转录物的生物学过程功能类型中，细胞过程和代谢过程占绝大部分，其次是刺激反应和生物调节过程。而在转录物的细胞组分功能类型中，细胞和细胞分离所占比例最高，其次是细胞器和生物膜。蛋白结合和催化活性是转录物分子功能类型中的主体部分。

为了进一步预测和分类转录物，运用COG库比对注释了枇杷幼果的14 538条转录物，证明其有蛋白功能。COG功能注释了至少25个功能类别的具体蛋白功能，它包括了细胞结构、生物化学代谢、分子过程、信号传递等生理生化过程(图4)。其中，在COG转录物功能注释分类中，一般功能基因所占比例最大(4 622，15%)，所占比例5%以上的有转录(3 181，10.32%)，复制、重组与修复(2 422，7.86%)，蛋白质翻译后修饰与转运、分子伴侣(2 360，7.66%)，信号传导机制(2 060，6.69%)，翻译、核糖体结构与生物合成(1 969，6.42%)，碳水化合物运输与代谢(1 889，6.13%)和功能未知(1 855，6.02%)，而防御机制(321，0.96%)、染色体结构与变化(314，0.94%)、核苷酸运输与代谢(274，0.82%)、RNA加工与修饰(189，0.57%)、胞外结构(16，0.05%)和核酸结构(3，0.01%)所占比例则均在1%以下。

图 3 根据GO分析功能注释枇杷幼果转录物

2.3 低温胁迫下枇杷幼果的差异表达基因

为探索低温胁迫对枇杷幼果转录物转录水平的影响，对低温胁迫处理的2个样本的基因差异表达，按照Audic等[5]根据转录物测序的基因表达差异方法进行筛选，结果总共检测到57 792条转录物与对照相比在表达上发生明显变化，这些转录物定位到4 017个基因中，有1 619个基因表达上调，有2 398个基因表达下调。

2.4 低温胁迫下枇杷幼果差异表达基因的功能注释

通过对枇杷幼果低温胁迫处理后的转录物的分析和定位，找到差异表达的基因。随后通过对这些差异基因进行GO分类和KEGG通径分析，弄清低温胁迫差异表达基因的功能，特别是在细胞的生物代谢途径方面的基因功能。总体上看，枇杷幼果低温胁迫诱导产生了2 007个差异表达的转录物。这些差异表达转录物可以被归于113条生化代谢或信号转导途径中，但最主要的途径是参与生化代谢和信号转导。这113条途径绝大部分是代谢途径、生物合成次生代谢物、细胞的内吞作用、甘油代谢、植物激素信号传导、醚酯代谢、淀粉和糖代谢、植物病原菌相互作用等途径。本试验结果中，值得注意的是枇杷幼果低温胁迫产生油菜素内酯生物合成和磷脂酰肌醇信号系统的几个基因差异表达上调，其中，共有6个参与调控油菜素内酯生物合成的基因差异表达上调(图5、表3)。而PLC(磷脂酶 C)基因则在磷脂酰肌醇信号系统中表达上调(图6、表3)。

图 4 根据COG分析功能注释枇杷幼果转录物

3 讨 论

转录组是生物或细胞特定状态下转录的所有RNA，含mRNA及非编码RNA。转录组的测序逐渐成为研究基因结构及功能的基础。根据转录测序信息可以弄清许多生物过程，比如基因表达情况[6]、试验处理或病害后的基因表达谱、基因调控、组织中生物标识物、基因挖掘、基因含量、分子系统学中的保守同源基因分离等[7-9]。因此，通过转录组的测序技术可以为研究对象提供许多候选功能基因。高通量技术是目前转录组测序技术中一种快捷而又可靠的方法，许多研究通过Denovo组装进行无参考基因组的非模式生物转录组分析[4]。然而，尚无有关枇杷果实转录组信息方面的研究报道。本研究中，运用Illumina高能量技术，对无参考基因组的枇杷果实通过Denovo测序仪进行了组装转录分析，鉴定和分析了大量表达的基因，得到了2千多万条序列数据，用多种互补方法组装了42 830个转录物。功能注释为枇杷研究提供了有价值的特定细胞过程、分子功能和生物代谢途径。本试验结果表明，高能量测序技术是研究无参考基因组生物的有效平台。这些枇杷果实转录组数据为今后开发枇杷分子和细胞标记提供了理想的序列资源。

图 5 通过KEGG 推定低温胁迫下枇杷幼果油菜素内酯的生物合成途径

图 6 通过KEGG 推定低温胁迫下枇杷幼果磷脂酰肌醇的信号系统途径

表 3 低温胁迫下枇杷幼果油菜素内脂生物合成和磷脂酰肌醇信号系统途径中的差异基因

低温胁迫影响植物的生长发育，但植物在低温胁迫下会产生一个调节系统，有千百个基因参与到胁迫下的信号传递和下游的胁迫反应过程中[10-12]。本研究通过枇杷低温胁迫转录组分析，试图鉴定可能参与到低温胁迫的一组基因，经鉴定发现有4 017个基因有可能参与到枇杷低温胁迫中，其中有1 619个基因表达上调，而2 398个基因表达下调。这些基因促进了人们对枇杷幼果耐低温胁迫分子机制方面的认识。

通过分子功能分析可全面探讨枇杷幼果对低温胁迫的反应。本研究结果表明，有6个差异基因表达上调，其参与调控油菜素内酯的生物合成。油菜素内酯(BRs)是植物甾醇类激素，在植物的正常生长发育中是不可缺少的。目前许多植物生理学家已将其列为植物的第6大类激素。而且，油菜素内酯还能增加植物的环境抗逆性。有研究发现，油菜素内酯能提高拟南芥幼苗忍耐低温的能力[13]。本试验结果表明，油菜素内酯参与了枇杷幼果低温胁迫调节系统。

本试验也发现，PLC(磷脂酶 C)的基因差异表达上调参与调控了磷脂酰肌醇信号系统，它是磷脂酰肌醇信号转换系统非常重要的组成部分。磷脂酰肌醇信号转换产生磷脂酶 C和Ca2+，而它们使胞外信号转换为胞内信号，并产生各种非生物胁迫响应。肌醇磷脂被磷脂酶 C水解，在受体刺激的同时，打开钙库释放Ca2+到细胞质中与钙调蛋白结合，随后参与一系列的反应。研究发现，AtPLC1(类磷脂酶 C突变体)很可能参与环境胁迫信号转换途径，提高AtPLC1水平可增强信号，这样有助于提高植物对环境胁迫的适应能力[14]。当PI-PLC的协同作用受抑制后，RD29A和COR47(冷应答信号通路的基因)2个基因表达下调，渗透胁迫诱导基因也出现下调[15]。因而，由本试验结果可以推测磷脂酶 C很可能参与调节枇杷幼果低温抗性反应。

本研究首次报道了运用高能量测序技术并通过Denovo组装，进行枇杷幼果转录组分析和低温胁迫下基因差异表达分析。本研究是在没有参考基因的基础上完成转录本的组装，得到了高质量的枇杷转录组序列，并进一步通过对枇杷转录组的注释，深入了解了枇杷幼果基因及其功能。这些发现将更有利于探索枇杷基因组资源，也能促进今后枇杷基因组和育种方面的研究。此外，本研究运用高能量测序技术分析低温胁迫下枇杷幼果基因的表达情况，表明低温胁迫下功能基因的研究方法是可行的，这也必然会促进对低温胁迫下枇杷幼果分子机理的深入了解。

[参考文献]

[1] 郑国华,张贺英,钟秀容.低温胁迫下枇杷叶片细胞超微结构及膜透性和保护酶活性的变化 [J].中国生态农业学报，2009,17(4):739-745.

Zheng G H,Zhang H Y,Zhong X R.Changes in cell ultra-structure,membrane permeability and protective enzyme activity inEriobotryajaponicaLindl.leaves under cold stress [J].Chinese Journal of Eco-Agricultre,2009,17(4):739-745.(in Chinese)

[2] Krost C,Petersen R,Schmidt E R.The transcriptomes of columnar and standard type apple trees (Malus×domestica):A comparative study [J].Gene,2012,498:223-230.

[3] Varshney R K,Hoisington D A,Tyagi A K.Advances in cereal genomics and applications in crop breeding [J].Trends Biotechnol,2006,24:490-499.

[4] Garg R,Patel R K,Tyagi A K,et al.Denevoassembly of chic-kpea transcriptome using short reads for gene discovery and marker identification [J].DNA Res,2011,18:53-63.

[5] Audic S,Claverie J M.The significance of digital gene expression profiles [J].Genome Res,1997,7:986-995.

[6] Torres T T,Metta M,Ottenwalder B,et al.Gene expression pr-ofiling by massively parallel sequencing [J].Genome Res,2008,18:172-177.

[7] Hegedus Z,Zakrzewska A,Agoston V C,et al.Deep sequencing of the zebrafish transcriptome response to mycobacterium infection [J].Mol Immunol,2009,46:2918-2930.

[8] Hughes J,Longhorn S J,Papadopoulou A,et al.Dense taxono-mic EST sampling and its applications for molecular systematics of theColeoptera(beetles) [J].Mol Biol Evol,2006,23:268-278.

[9] He R.Next-generation sequencing-based transcriptomic and pr-oteomic analysis of the common reed,Phrafmites australis (Poaceae),reveals genes involved in invasiveness and rhizome specificity [J].Amer J Bot,2012,99:232-247.

[10] Kreps J A,Wu Y,Chang H S,et al.Transcriptome changes forArabidopsisin response to salt,osmotic,and cold stress [J].Plant Physiol,2002,130:2129-2141.

[11] Shinozaki K,Yamaguchi-Shinozaki K,Seki M.Regulatory network of gene expression in the drought and cold stress responses [J].Cur Opin Plant Biol,2003,6:410-417.

[12] Xiong L,Schumaker K S,Zhu J K.Cell signaling during cold,drought,and salt stress [J].Plant Cell,2002,14:165-183.

[13] Divi U K,Krishna P.Overexpression of the brassinosteroid biosynthetic gene AtDWF4 inArabidopsisseeds overcomes abscisic acid-induced inhibition of germination and increases cold tolerance in transgenic seedlings [J].Plant Growth Reg,2010,29:385-393.

[14] Hirayama T,Ohto C,Mizoguchi T,et al.A gene encoding a phosphatidylinositol-specific phospholipase C is induced by dehydration and salt stress inArabidopsisthaliana[J].PNAS,1995,92:3903-3907.

[15] Takahashi S.Hyperosmotic stress induces a rapid and transient increase in inositol 1,4,5-trisphosphate independent of abscisic acid inArabidopsiscell culture [J].Plant Cell Physio,2001,42:214-222.