杜鹃花叶片转录组测序数据组装及功能注释

王 华，汪王微，王冬良，张石虎，胡新芳，卢诗雨，龚雪梅

(1.安徽农业大学 园艺学院 观赏园艺系，安徽 合肥 230036； 2.阜阳职业技术学院 生化工程系，安徽 阜阳 236031)

杜鹃属 (Rhododendron) 植物大多花美色艳，是中国乃至世界著名的观赏植物。但是，时至今日，杜鹃花在中国园林中的应用依然相对较少，其主要原因在于野生杜鹃引种栽培后生长或开花不良，景观效果不理想。为此，国内外学者开展了杜鹃对环境需求的研究，如光照[1-2]、温度[3-5]、水分[6]、土壤[7]等，这些研究为制定杜鹃花合理引种栽培策略提供了理论依据，但都仅从形态学、生理生化或叶片超微结构方面进行探讨，而对逆境条件下杜鹃花代谢相关基因及其功能研究鲜有报道。利用Illumina sequencing进行第二代高通量测序可在较短的时间内获取海量数据，并对数据进行处理分析。近年来，利用这一技术已成功地对多种植物的逆境胁迫进行了研究。杨伟等[8]以枇杷幼果为试材，研究低温胁迫下幼果转录组Denovo组装和功能注释，应用Illumina高通量测序技术得到64 814条Unigenes，注释到了4 017个差异表达基因，为枇杷的抗旱基因挖掘提供了资料。韩超等[9]以高温和适温下的梭枝条为实验材料，利用Solexa测序技术，获得了162 504个差异表达基因，为研究高温胁迫下梭梭的分子机理奠定了基础。黄玉兰等[10]以薏苡苗期叶片为实验材料，通过转录组测序，获得了92 865条Unigenes，通过与基因本体(gene ontology，GO)数据库比对，发现可能有1 053条Unigenes参与干旱胁迫。彭振等[11]以棉花幼苗真叶为材料，利用转录组测序(RNA-Seq)技术分析其在盐胁迫下转录因子家族及转录因子的表达情况，结果发现了大量的盐胁迫响应转录因子，且随着胁迫时间的延长，棉花对盐响应的转录因子急剧增加。

本研究利用高通量测序技术对杜鹃花白凤4号中部功能叶进行转录组测序，并对数据进行注释和分析，可以有效地发掘和鉴定杜鹃花抗逆代谢相关调控基因，为克隆开发杜鹃花重要抗逆功能基因提供有力的研究基础，为杜鹃花选育、扩大园林应用提供理论基础和技术支撑。

1 材料与方法

1.1 实验材料

于2016年7月，以安徽农业大学实验基地所栽种的长势良好的2年生杜鹃花白凤4号为材料，取从上往下第2-4片成熟功能叶用于总RNA的提取。

1.2 样品处理及转录组数据库分析流程

采用 TIANGEN 公司生产的多糖多酚植物总RNA提取试剂盒(RNAprep Pure Plant Kit Polysaccharides﹠Polyphenolics-rich) 进行杜鹃花叶片的总RNA提取。取叶片100 mg加入液氮研磨至粉末状，总RNA的提取流程参照RNAprep Pure Plant Kit说明书(离心柱型)(天根生化科技(北京)有限公司)方案进行。杜鹃花转录组序列的测定和组装分析委托广州基地奥公司进行。利用Oligo(dT)磁珠从总RNA中纯化mRNA，将mRNA打断为片段，构建杜鹃花叶片转录组的cDNA文库，将通过质量检测后所构建的文库用第二代测序仪IlluminaHi Seq TM 2000进行测序，并对所获得的序列运用Trinity软件进行Denovo拼接组装。随后，将组装得到的杜鹃花转录本数据与公共数据库(NR、Swiss-Prot、KOG、KEGG和GO)进行比对，对基因注释、基因功能分类及转录因子等进行分析。

2 结果与分析

2.1 杜鹃花白凤4号转录组组装质量分析

在对杜鹃花的叶片样品进行转录组测序后，通过严格质量控制和数据过滤，共检测得到213 723 424个高质量Clean reads，占总数的95.05%，过滤后质量阈值≥20。采用Trinity软件，对这些高质量数据进行Denovo组装，得到N50为1 600 nt，平均长度为930 nt，GC核苷酸的含量(GC%)为48.05%的杜鹃花Unigene 53 568 条(不含N的组装片段)，其最大长度为13 121 nt，最小长度为201 nt(表1)。对组装出来的Unigene做长度分布统计(图1)，在组装得到的Unigene中，片段长度大于3 000 nt的共有2 331个。

2.2 杜鹃花白凤4号基因注释信息

将所得到的Unigene序列利用BLAST工具分别与NR、Swiss-Prot、KOG和KEGG进行比对，共有28 877条Unigenes (阈值e-value<1E-5)得到注释，占基因总数的53.90%，其中比对到Nr库的有28 626个Unigenes，占比对上基因的99.13%，且有14 465个Unigenes的长度≥1 000 nt。在Swiss-Prot、KOG和KEGG数据库中共有21 548，18 283和11 409条Unigenes获得注释信息(表2)，而且在四大库中都比对上的基因有9 162个，而在四大库中没有任何比对信息的Unigenes有24 691条(占基因总数的46.09%)。

2.2.1 杜鹃花白凤4号转录组Unigenes在Nr数据库的比对结果分析

进一步对在Nr中比对上的基因进行分析，结果如图2。从图2-A可以看出，有59.65%的转录物序列具有极强的同源比对信息(e-value<1E-50)，而且有40.34%的转录物序列相似性E值在1.0E-5～1.0E-50，说明同源性比对效果较好。

图1 杜鹃花白凤4号转录组Unigenes长度分布Fig.1 Assembled Unigenes length distribution of Rhododendron pulchurum cv. Baifeng 4 transcriptome

对杜鹃花转录组所有Unigene的CDS进行预测，共通过BLAST比对获得CDS序列28 455个，利用ESTscan数据库中进行分析，可获得CDS序列1 897个。

收集有Nr注释的28 626条Unigene并与其他物种进行比对，统计其物种信息如图2-B。杜鹃花在葡萄(Vitisvinifera)、可可(Theobromacacao)、芝麻(Sesamumindicum)、莲(Nelumbonucifera)、麻风树(Jatrophacurcas)、梅(Prunusmume)、绒毛状烟草(Nicotianatomentosiformis)、林烟草(Nicotianasylvestris)、亚洲棉(Gossypiumarboreum)和甘蓝型油菜(Brassicanapus)中都有同源序列分布，与葡萄有较高的同源性，有16.94%的同源序列。除此以外，与其他9个植物种类的同源性均相对较低，同源序列都不足10%，还有45.18%的Unigenes属于其他序列。

A，Nr库注释结果E值分布图；B，Nr库注释结果(E值最低)物种分布统计图。A， E value distribution map of Unigenes annotation in Nr database; B， Statistics of species distribution of Unigenes annotation in Nr database( The E values were minimum).图2 杜鹃花白凤4号Nr库比对Unigenes分析Fig.2 Analysis of Unigenes annotation of Rhododendron pulchurum cv. Baifeng 4 in Nr database

2.2.2 杜鹃花白凤4号转录组Unigenes KOG功能分类

通过将转录组的Unigenes与KOG数据库进行比对，可注释得到30 508个有功能信息的Unigenes。依照其注释功能可划分为25类(图3-A，W，Y，Z)，通过对每一类KOG注释的基因功能进行分类统计，可获得杜鹃花转录组中各类功能基因的数量分布信息(图3)。其中，一般功能预测所拥有的Unigenes最多，有6 692条，所占比例最大，达到21.94%。其他所占比例达到5%以上的基因分类还有蛋白质翻译后修饰与转运、分子伴侣(3 433，11.25%)，信号传导机制(3 257，10.68%)，翻译、核糖体结构与生物合成(1 803，5.91%)，RNA加工与修饰(1 774，5.81%)和转录(1 601，5.24%)，而核苷酸运输与代谢(251，0.82%)、辅酶运输与代谢(227，0.74%)、防御机制(171，0.56%)、胞外结构(90，0.30%)、核酸结构(83，0.27%)和细胞运动(10，0.03%)所占比例都不足1%。其他包括脂类转运与代谢(963，3.16%)、无机离子转运与代谢(578，1.89%)和胞壁/细胞膜/细胞被膜再生(325，1.07%)在内的13个基因分类所占比例均在1%～5%(图3)。

2.2.3 杜鹃花白凤4号转录组Unigenes的GO功能分类分析

通过GO数据库进行杜鹃花转录组功能分析，共有17 218个Unigenes注释到生物过程、细胞组分类和分子功能类三个方面的功能，共计55个不同的功能分支(表3)。其中有38 434个Unigenes注释到生物过程，其中代谢过程、细胞过程占绝大多数，其次是单一的生物过程和刺激反应；在细胞组分功能中注释到33 059个Unigenes，其中占比例最高的是细胞和细胞部分，其次是细胞器和细胞膜；与分子功能相关的Unigenes有19 752个，其中具有催化活性和绑定的Unigenes最多，其次是输送分子和结构分子。

进一步对GO注释中与抗逆相关的注释进行梳理，共发现11 928条注释(表4)，其中注释最多的为参与各类环境应激响应的注释(6 896条)，其他依次为蛋白激酶(1 454条)、激素代谢(59条)相关、渗透调节(460条)、氧化还原反应(2 004条)、受体(469条)有关，以及离子转运(586条)相关。

2.2.4 杜鹃花白凤4号Unigenes的KEGG代谢途径分类分析

在杜鹃花的转录组数据中共有6 475个基因可以通过Blast比对注释到132类KEGG的代谢通路中，属于基因信息处理的有21条代谢通路，且Unigenes超过200个的有核糖体(775)，内质网蛋白加工(354)，剪接体(348)和RNA转运(287)属于新陈代谢的有99条代谢途径，含有200个以上Unigenes的有碳代谢(389)，氨基酸生物合成(311)，氧化磷酸化(235)，淀粉与蔗糖代谢(221)，糖酵解途径(213)和嘌呤代谢(202)；属于细胞过程的有4个代谢通路，其中(细胞)内噬作用的Unigenes为265，其他的Unigenes都不足200个；属于环境信息处理的有3个代谢通路，分别是植物激素信号转导(220)，磷脂酰肌醇信号系统(77)和ABC转运蛋白(51)；属于其他的代谢途径有5个，其中植物-病原菌的交互含有297个Unigenes，其余的Uningenes都不足100个(表5)。除此以外，在杜鹃花的转录组数据中也注释到了植物内源激素代谢的相关途径，分别是与乙烯(ethylene，ETH)相关的半胱氨酸和蛋氨酸代谢(cysteine and methionine metabolism，121)，与生长素(auxin，IAA)相关的色氨酸代谢途径(tryptophan metabolism，62)，与水杨酸(salicylic acid，SA)相关的苯丙氨酸代谢途径(phenylalanine metabolism，55)，与茉莉酸类物质(jasmonic acid，JAs)相关的α-亚麻酸的代谢途径(alpha-linolenic acid metabolism，47)，与脱落酸(abscisic acid，ABA)相关的类胡萝卜素生物合成途径(carotenoid biosynthesis，40)，与赤霉素(Gibberellin，GAs)相关的二萜类化合物的生物合成途径(diterpenoid biosynthesis，21)，与细胞分裂素(cytokinine，CTKs)相关的玉米素合成途径(zeatin biosynthesis，19)，以及与新型植物激素油菜素内酯(brassinolide，BR)相关的油菜素内酯的生物合成途径(brassinosteroid biosynthesis，19)(表5)。

X轴上的字母代表不同的功能类别: A，RNA加工与修饰；B，染色质结构与动力；C，能量生产与转化；D，细胞周期调控、细胞分裂、染色体分裂；E，氨基酸转运与代谢；F，核苷酸转运与代谢；G，糖转运与代谢；H，辅酶转运与代谢；I，脂类转运与代谢；J，翻译、核糖体结构和生物合成；K，转录；L，复制、重组和修复；M，细胞壁、膜、包体生物合成；N，细胞运动性；O，蛋白质翻译后修饰与周转；P，无机离子转运与代谢；Q，次生代谢产物的生物合成；R，一般功能预测；S，未知功能；T，信号传导机制；U，胞内运输、分泌和囊泡运输；V，防御机制；W，胞外结构；Y，核酸结构；Z，细胞骨架。The capital letters in X-axis indicate the KOG categories， and the Y-axis means the number of Unigenes. A， RNA processing and modification；B， Chromatin structure and dynamics；C， Energy production and conversion；D， Cell cycle control， cell division， chromosome partitioning；E， Amino acid transport and metabolism；F. Nucleotide transport and metabolism；G， Carbohydrate transport and metabolism；H， Coenzyme transport and metabolism；I， Lipid transport and metabolism；J， Translation， ribosomal structure and biogenesis；K， Transcription；L， Replication， recombination and repair；M， Cell wall/membrane/envelope biogenesis；N， Cell motility；O， Posttranslational modification， protein turnover， chaperones；P， Inorganic ion transport and metabolism； Q， Secondary metabolites biosynthesis， transport and catabolism；R， General function prediction only；S， Function unknown；T， Signal transduction mechanisms；U， Intracellular trafficking， secretion， and vesicular transport；V， Defense mechanisms；W， Extracellular structures；Y， Nuclear structure；Z， Cytoskeleton.图3 杜鹃花白凤4号Unigenes 在KOG 功能分类Fig.3 Clusters of orthologous group (COG) classification of Unigenes of Rhododendron pulchurum cv.Baifeng 4

表3 杜鹃白凤4号GO注释结果

表4 杜鹃白凤4号转录组中与逆境胁迫相关的基因

表5 杜鹃花白凤4号转录组序列的KEGG通路分析

2.3 杜鹃花转录组转录因子注释

杜鹃花转录组中1 062个Unigenes被注释为转录因子，这些转录因子分布在55个转录因子家族中：富集到20个以上Unigenes的家族共有18个，其中富集基因数量在50个以上的有bHLH、MYB-related、FAR1、ERF、WRKY、NAC和C2H2家族，富集基因数量在20～49个的有MYB、bZIP、Nin-like、C3H、GRAS、HD-ZIP、B3、Trihelix、G2-like、Dof及GATA家族(图4)；其他37个转录因子家族共注释到Unigenes 239个，但每个家族富集的基因最多17个，最少只有1个。

2.4 杜鹃花白凤4号转录组SSR分析

利用MISA 软件分析筛选得到的所有53 568个Unigenes中的SSR，发现包含 SSR 的序列数目为8 738个条，识别的SSR总数为11 108条。其中双碱基重复SSR 7 805条，占70.26%；三碱基重复SSR 2 469条，占22.23%；四碱基重复SSR 385条，占3.47%；五碱基和六碱基重复SSR分别有240条(2.16%)和209条(1.88%)(表6)。双碱基SSR中，发生频率最高的是AG，其次是AC和AT。三碱基重复中共有10个类别的SSR，其中发生频率较高的是AAG，AGG和ACC，发生频率最低的是AAT；四碱基重复中AAAT是发生频率最高的，其次是AAAG。

图4 杜鹃花白凤4号转录组中转录因子分类Fig.4 Classification of transcription factors of the assembled of Rhododendron pulchurum cv. Baifeng 4

表6 杜鹃花白凤4号SSR分析结果统计

3 讨论

本研究在没有参考基因组基因组测序情况下，以杜鹃花白凤4号叶片为实验材料，通过Denovo组装获得53 568条Unigenes。把这些基因在Nr库中注释到28 626条Unigenes，且通过与其他物种进行同源序列比对发现，这些基因与葡萄的同源性最高。通过Blast对所有Unigenes CDS进行预测，获得了28 455个CDS序列，为后续研究奠定了基础。在GO库注释到的17 218个Unigenes中，我们挖掘到了11 928个与植物抗逆有关的Unigenes，分别隶属于GO注释分类的环境应激响应、蛋白激酶、植物内源激素、渗透调节、氧化还原、受体、离子转运等[10]。

本研究中杜鹃花转录组数据注释到KEGG的有分布在132个代谢途径中的6 475个Unigenes，其中注释到应激激素ABA[12]等植物激素代谢途径中的Unigenes有384个，涉及ABA生物合成和代谢的有40个，占10.42%。其中，注释为ABA合成路径中玉米黄质环氧化酶(zeaxanthin epoxidase，ZEP)[13]和ABA合成限速酶9-顺式-环氧类胡萝卜素双加氧酶(9-cis-epoxycarotenoid di-oxygenase activity，NCED)[14]的各有3个Unigenes，注释为ABA降解关键酶脱落酸8'羟化酶(abscisic acid 8′-hydroxylase，ABA 8′OH)[14]的有9个Unigenes。同时，植物激素信号转导代谢途径在KEGG中注释到220个Unigenes，其中与ABA信号转导有关的有38个，分别是7个ABA受体PYR/PYLs[15]、8个蛋白磷酸酶(type 2C proteinphosphatase，PP2C)[16]、9个SNF1相关激酶2[17](SNF1-related protein kinase 2，SnRK2)以及14个ABFs(ABA response element binding factor)转录因子[18]。植物各生理过程中特定的ABA含量受其合成、代谢及转运等精细调控，以上这些基因的注释表明，当杜鹃花遭受逆境胁迫时会通过调动包括以上基因在内的基因表达来调控ABA在体内的动态变化，从而对逆境做出响应。

转录因子(transcription factor，TF)又称反式作用因子，在植物应对外界环境变化的信号转导或者基因转录水平调控过程中起着重要的作用，并且往往是以基因家族的形式出现[19]。本研究中注释到的转录因子家族有bHLH、MYB-related、ERF和WRKY等，其中WRKY家族基因是近几年研究比较热的与植物胁迫应答相关的转录因子。

WRKY转录因子由WRKY蛋白N端高度保守的WRKYGQK氨基酸序列及其C端特殊的锌指结构而命名，在植物体内当受外界因子诱导后特异性识别W-box([(T) (T) TGAC(C/T)]序列)而调节下游基因的表达，从而在多种生物和非生物胁迫及诸如水杨酸[20]、赤霉素[21]等植物激素信号转导中起重要作用，目前，也有越来越多的证据表明WRKYs是ABA应答信号网络的关键节点[22]。本研究中注释到52个WRKYs，其中4个WRKY4转录因子值得关注。Ren等[23]研究发现烟草WRKY4在叶片生长和生物胁迫方面都具有重要作用，这在WRKYs转录因子中是比较少见的；Zheng等[24]研究发现，刚毛柽柳中WRKY4能被ABA、盐和干旱诱导，同时被ABF和DOF(DNA binding with one finger)转录因子调控，同时WRKY4也能与其他WRKYs交叉调控，过表达WRKY4拟南芥植株能提高其抗盐、抗氧化和抗ABA诱导的能力；徐丽等[25]从核桃中克隆了WRKY4基因，还发现其受低温的诱导。由此可见，WRKY4转录因子在植物体内的作用比较复杂，会参与杜鹃花所的逆境响应，值得深入研究。

简单序列重复SSR(simple sequence repeats)标记是目前最理想的分子标记，广泛用于遗传多样性和进化论研究[26]。但引物开发要建立在已知的DNA序列上，由于之前杜鹃植物DNA序列信息相对缺乏，使得杜鹃SSR标记引物的开发相对较少，本研究在杜鹃花转录组测定的基础上筛选了8 738 SSRs序列，有望会同前人开发的SSR标记[27-30]为杜鹃种质资源的有效分类与鉴定研究奠定基础。