汪平　高建明　杨峰　郑金龙　刘巧莲　陈河龙　易克贤

摘 要 为克隆剑麻受斑马纹病病原菌烟草疫霉侵染后的应答基因，以病原菌侵染前后的剑麻H.11648叶片为材料，构建剑麻转录组Unigene数据库；组装得到非冗余Unigene 131 422个，其中，500～1 000 bp的Unigene占总数65.39%，Unigene的数量随其长度递减。GO分类分析发现：在得到的所有Unigene中，注释到生物过程的基因最多（152 835）；细胞组分其次（129 197）；参与分子功能的最少（47 420）。KEGG代谢通路分析将剑麻转录组unigene分成128类，其中，参与生化代谢通路的Unigene最多，有9 354个，占总数27.09%；参与植物-病原互作代谢通路的基因有1 795个，占总数5.2%。基因表达丰度RPKM值显示上调基因9 415个，下调基因28 980个，其中参与病原互作的基因占680个。

关键词 烟草疫霉；侵染；剑麻；转录组

中图分类号 S563.8 文献标识码 A

剑麻（学名：Agave sisalana；英文名：Sisal）俗称龙舌兰麻[1]，为龙舌兰科（Agavaceae）龙舌兰属（Agave linnaeus）多年生肉质旱生草本植物，是热带、亚热带地区栽培的重要硬质纤维作物及重要的工业原料，应用广泛。除纤维外，其液汁可提取贵重药物生产原料皂素；麻渣含有丰富的营养物质，是良好的饲料和肥料[2]。中国是世界剑麻主要生产国之一，产地主要分布在广东、广西、海南、云南、福建等热带和亚热带地区，主栽品种为剑麻H.11648（Agave hybrid No.11648）[3]。该品种叶片宽厚，产量高，但容易感染由烟草疫霉（Phytophthora nicotianaeBreda.）引起的斑马纹病（Zebra disease，ZD）[4]。该病1961年首次在坦桑尼亚发生，后传播到中国等世界各地，对剑麻产业造成极大的损失，为剑麻周期产量影响最大的病害之一[5]。由于形态和真菌很相似，过去曾将烟草疫霉归为真菌，实为色菌界（Chromista）二倍体真核微生物[6]，其区别于真菌的显著之处是细胞壁不含几丁质成分，因而传统的针对几丁质的杀菌剂，对杀灭烟草疫霉等卵菌几乎无效[7-8]。防治烟草疫霉的常用药剂多为苯基酰胺类及乙磷铝等同一种或作用机制相同的几种内吸性杀菌剂，很容易产生耐药性，从而引起药效下降[9]。

转录组研究是一种发掘功能基因的重要途径，是基因功能及结构研究的基础和出发点。了解转录组是解读基因组功能元件和揭示细胞及组织中分子组成所必需的，并且对理解机体发育和疾病具有重要作用[10]。转录组学相对于基因组学而言，只研究被转录的基因，研究范围缩小，针对性更强[11]。转录组测序（RNA-seq）是利用大规模测序技术直接对cDNA序列进行测序，产生数以千万计的reads片段，从而使得一段特殊的基因组区域的转录水平可以直接通过比对该基因组区域的reads数来衡量[12]。通过转录组分析，可高通量地获得基因表达的RNA水平，从而揭示基因与生命现象之间的内在联系。据此可以了解细胞生理活动规律，确定细胞代谢特性，并进而对细胞进行修饰改造[13-14]。目前，转录组学用于植物和病原菌研究的很多。但是，用转录组学研究剑麻与烟草疫霉互作机理的文章还未见报道。本研究在RNA-seq技术的基础上，通过对一系列数据的分析、比对，得到剑麻对烟草疫霉侵染的应答基因，以期为将来转基因抗病育种奠定基础。

1 材料与方法

1.1 材料

1.1.1 植物材料 实验材料取自中国热带农业科学院（海口）院内本实验室基地大棚，选取株高50～70 cm无病害、长势良好的剑麻H.11648植株，移栽入营养盆后置于室内控温28 ℃左右培养备用；

1.1.2 试验试剂及药品 康为世纪RNApure Plant Kit（DNaseⅠ）（Cat.No.：CW0559）试剂盒；Thermo Scientific RevertAid First Strand cDNA Synthesis Kit （Cat.No.：#K1622）；其它试剂均为国产分析纯。

1.1.3 试验仪器 Eppendorf 5810R台式高速冷冻离心机；琼脂糖凝胶电泳仪；Tanon-4500凝胶成像系统；美国Thermo Nanodrop 2000微量紫外分光光度计等。

1.2 方法

1.2.1 菌株准备 用于接种的烟草疫霉菌株为本实验室前期分离获得，经PDA平板活化后保存备用。

1.2.2 接种 选取剑麻叶片基部幼嫩部位进行接种，先用75%酒精对其进行消毒，再用灭菌水清洗，待干燥后接种。实验设处理（CL）与对照（CK）两组，每组选5株作为接种对象，处理组在叶片基部幼嫩部位刺伤并接种烟草疫霉，对照组只做刺伤，不接种烟草疫霉，做好上述处理后，在植株叶面上方喷无菌水，并用透明薄膜将植株从底部套住以达到保湿的效果，以后每天喷水2次。

1.2.3 总RNA提取 RNA提取方法参照康为世纪RNApure Plant Kit（DNaseⅠ），接种后1、2、3、4、5 d的感病叶片等量混合后提取总RNA（CL组），同时取对照（CK组）。提取后RNA进行琼脂糖凝胶电泳检测及Thermo Scientific NanoDrop 2 000分光光度计分析。

1.2.4 转录组测序 检测合格的RNA样品委托深圳华大基因科技服务有限公司测序（项目编号：F13FTSSCKF0046）。

2 结果与分析

2.1 产量统计

对CK和CL 2组样品进行转录组测序后共获得8G原始数据。其中，CK组获得原始片段（reads）68、802、552条，去除带接头，重复，测序质量低片段后clean reads数59、305、846；CL组reads数74、975、368，clean reads 66、778、374个。CL获得的测序片段数略多于CK，这一结果与预期相符，因为CL中除剑麻表达基因外还有烟草疫霉表达的基因，理论上比CK表达的基因多。两者的GC含量介于45%～50%之间，较为一致。测序的碱基质量参数（Q20）均在97%左右，不确定碱基比例0.00%（表1）。

2.2 组装质量统计

组装质量评估通过组装后所得片段的长度分布及其N50值来反应，N50指从组装最长的Unigene 依次向下求长度的总加和，当这个加和达到组装长度的一半时，对应的Unigene的长度就是N50的长度。N50值越大，反映组装得到的长片段越多，组装效果越好。CK组Unigene 总数127、361，Contig总数253、520；CL组Unigene 总数153、951，Contig总数320、529。All-Unigene的长度分布大致在300～3 000之间， 500～1 000 bp的Unigene最多，占总数65.39%，大于3 000 bp的大片段较少。Unigene数量随长度成递减关系见图1。

CK组装转录组N50=275，Unigene N50=638；CL组装转录组N50=299，Unigene N50=667。将上述2个处理的Uingene数据进行比对聚类合并，作为剑麻Uingene数据库。结果显示，该Unigene库共有131、422个Uingene，N50=861，组装结果较好，能满足后续分析需要。

2.3 Uingene的覆盖度、深度和表达量分析

对Unigene进行覆盖度分析发现，131、422个Unigene的覆盖范围介于5.94%～96.91%之间；测序深度从0.059 4%～9.166 7%不等；样品中能唯一比对到指定unigene序列的reads数（Unique-mapped-Reads）从1～299个不等；能比对到多个unigene序列的reads数（Multiple-mapped-Reads）从1～385个不等。

2.4 Unigene的功能注释

将All-Unigene分别注释到NR、NT、Swiss-Prot、KEGG、COG、GO 6个数据库中，对注释到每个库的Unigene数目统计结果见表2。

将剑麻数据库中Unigene比对到氨基酸序列数据库（Non-redundant sequences，Nr），结果显示，注释到葡萄（vitis vinifera）的Uingene最多，粳稻（oryza sativa japonica group）其次，具体物种分布见图2。

E-Value值是序列比对到该物种的可靠性评价，它表明在随机情况下，其它序列与目标序列相似度要大于这条显示的序列的可能性，因而它的值越低越好。当E值小于10-5时，表明两序列有较高的同源性，本次转录组测序All-Unigene注释到Nr的E-value值分布图见图3。从图3中可以看出，比对到的物种序列同源性均小于10-5，其中小于1e-30的序列占比对序列总数超过50%，说明比对结果可靠性很高。

对All-Unigene的注释结果进行物种相似性分析发现：多数Unigene序列比对相似性在40%～80%之间，序列相似性高于95%的很少，占All-Unigene总数不到4%，但是从整个序列相似性分布来看，对无参考基因组的剑麻All-Unigene的功能注释结果较好。剑麻All-Unigene物种相似性分布图见图4。

对剑麻Unigene数据库做蛋白相邻类的聚簇分析COG（Cluster of Orthologous Groups of proteins）显示：注释到通用功能（R）的Unigene最多；转录组相关功能（K）其次；胞外结构（W）和核酸结构相关功能（Y）最少；值得一提的是，有507个Unigene被注释到防卫反应（V），为后续的病原互作基因研究打下了很好的基础；此外，还有3 679个未知功能的Unigene。COG分类注释结果见表3。

2.5 Unigene的GO分类

Gene Ontology（简称GO）是一个国际标准化的基因功能分类体系，对剑麻Unigene数据库做GO分类分析发现：共有42 987个Unigene注释到GO数据库，且较多的单个基因可以与多种基因相对应，最多的一个基因“CL4909.Contig1_All”可以同时与GO数据库中69个基因相似。建立了42 987个对应关系，从而得到更多的分类和注释信息。剑麻转录组中Unigene根据GO功能分类可分为3个ontology，其中生物过程（biological process）类中基因最多（152 835）；细胞组分（cellular component）其次（129 197）；分子功能（molecular function）最少（47 420）。各分支中基因表达丰度不尽相同。

2.6 Unigene的KEGG代谢通路分析

根据KEGG数据库代谢通路可以将剑麻转录组数据分成128类，包括生化代谢途径、次生代谢产物合成、内吞作用、环甘油磷脂代谢、植物-病原互作、植物激素信号转导、苯丙氨酸代谢、萜类化合物生物合成、过氧化物酶体、类黄酮生物合成、基底转录因子、泛素介导的蛋白质水解等。其中，参与生化代谢通路的Unigene最多，有9 354个（占总数的27.09%）；参与植物-病原互作代谢通路的基因有1 795（5.2%）个；参与苯丙素生物合成的基因有431（1.25%）个；参与吞噬体的基因有412（1.19%）个；参与过氧化物酶体的基因有284（0.82%）个；参与黄酮类生物合成的基因有227 （0.66%）个；参与苯丙氨酸代谢的基因有224（0.65%）个；参与萜类生物合成的基因有122（0.35%）个；参与黄酮和黄酮醇类生物合成的基因有107（0.31%）个；参与天然杀伤细胞介导的细胞毒作用的基因有106（0.31%）个；参与异黄酮合成的基因有47（0.14%）个；参与倍半萜和三萜生物合成的基因有37（0.11%）个。

2.7 差异基因分析

对剑麻数据库进行差异分析发现：转录本中有38 395个差异表达Unigene，其中上调基因9 415个，下调基因28 980个；对Unigene的表达量计算采用RPKM算法，计算公式为：RPKM=（1 000 000×C）/（N×L×1 000），设RPKM为Unigene A的表达量，则C为比对到Unigene A的reads数，N为比对到所有Unigene的总reads数，L为Unigene A的碱基数；用log2（CL_RPKM/CK_RPKM）来表示2个处理所表达基因的差异倍数；用FDR（False Discovery Rate）值来进行多重检验控制错误率，本研究取FDR≤0.001。对处理后的数据进一步分析，筛选植物-病原互作相关的基因发现：在所有680个植物-病原互作相关基因中，有309个为上调表达基因，其log2（CL_RPKM/CK_RPKM）值自12.151 1至1.003 8不等，下调基因371个，log2（CL_RPKM/CK_RPKM）值自-1.0005至-14.119 7不等。

再根据表达倍数差异优先选取前50位上调/下调基因共100个，以进行后续干扰表达和超表达研究，这100个基因是注释到Nr-annotation数据库后，剔除疫霉基因所获。通过对这些基因的初步分析发现：⑴在差异倍数上，上调基因差异倍数最大为log2（CL_RPKM/CK_RPKM）等于15.851 2，最小为13.842 6；下调表达基因中，下调倍数最多为-13.998 5，最少有-12.304 2。⑵注释到Nt-annotation和Nr-annotation数据库后，有67个Unigene可以比对到两个数据库，还有33个无法比对到序列。比对结果如下：葡萄（Vitis vinifera）17；二蕙短柄草（Brachypodium distachyon）7；毛果杨（Populus trichocarpa）6；水稻（Oryza sativa Japonica Group）6；蓖麻（Ricinus communis）5；高粱（Sorghum bicolor）4；橹豆（Glycine max）3；拟南芥（Arabidopsis thaliana）3；苜蓿（Medicago truncatula）2；大麦（Hordeum vulgare subsp. Vulgare）2；玉米（Zea mays）2；蓝桉（Eucalyptus globulus subsp. Globulus）1；闭鞘姜（Costus speciosus）1；红松（Pinus koraiensis）1；黑麦（Secale cereale）1；欧芹（Petroselinum crispum）1；烟草（Nicotiana tabacum）1；兰花（Cymbidium hybrid cultivar）1；沉水樟（Cinnamomum micranthum f. kanehirae）1；芝麻（Sesamum indicum）1；油棕（Elaeis guineensis）1。

经过以上的基础生物信息学分析，再择优选取Unigene构建载体进行RNA干扰和超表达转基因研究。

2.8 植物-病原菌互作通路分析

生物体内不同基因间通过相互协调作用来完成某一生物学功能，Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。以KEGG Pathway为单位，背景序列的KO号通过同blast比对来获取，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的pathway，一般Qvalue值设为≤0.05。基于Pathway的植物-病原物互作通路的分析有助于更进一步了解剑麻与斑马纹病互作的机理，其代谢通路见图5。

图中每一小方框代表一个代谢控制点，方框由黑、红、绿三色组成，红色部分代表此处有上调基因表达，绿色代表下调基因，全红或全绿代表此处基因全部上调或下调，黑色代表此处无基因差异表达。从图中可以看出，病原互作这一通路中，各代谢支路相互交错，基因间通过各自的上调或下调来共同完成这一生命活动。

3 讨论与结论

转录组测序对测序样品RNA纯度、完整度及浓度均要求较高，鉴于剑麻成熟叶片纤维含量丰富[15]，提取叶片RNA时很难达到测序的浓度和纯度要求，取样时选取株高50～70 cm的剑麻基部幼嫩叶片做为RNA提取材料。实验证明，在该部位接种烟草疫霉病原菌提取总RNA能达到测序的A类标准要求，是理想的取材部位。

本研究剑麻All-Unigene数据库有95.31%（60341）的Unigenen被注释到Nr数据库，跟番薯[16]，芝麻[17]，萝卜[18]46.21%，53.91%，85.51%的注释比例相比有显著的提高，究其原因主要跟测序技术的进一步成熟以及拼接水平的提升有关[16]。此外，注释到NT、Swiss-Prot、KEGG、COG、GO数据库的比例也相对较高，分别为66.81%、59.25%、54.53%、32.96%、67.90%，说明此次测序结果相对较好。

实验需用烟草疫霉对剑麻植株进行侵染，这样势必导致转录组测序结果中含有烟草疫霉表达的基因，而干扰实验结果。可以利用本文提到的方法，即生物信息学分析剔除干扰基因，或者根据实验需要设置烟草疫霉侵染其他寄主植株的转录组测序以消除干扰。

对2个处理的样本进行差异表达基因分析旨在找到剑麻抗斑马纹病病原互作相关基因，为转基因抗病育种奠定基础。本研究通过将两个样本比对到KEGG数据库，找到植物-病原互作通路基因680个，并给出了2个样本间差异表达基因的差异倍数关系。此外，本研究还发现其苯丙氨酸解氨酶（Phenylalanine ammonia-lyase，PAL，E.C. 4.3.1.5）Unigene基因均有不同程度上调现象，上调倍数自log21.1654至log26.6458不等，而在植物抗病反应的次生代谢中，苯丙氨酸解氨酶是形成植保素、木质素和酚类化合物等抗病次生物质代谢途径中的关键酶和限速酶[19]，说明剑麻在抗斑马纹病病原菌侵染过程中，苯丙氨酸解氨酶发挥了较大作用。

参考文献

[1] 裴超群.龙舌兰科植物资源调查报告[J]. 广西热作科技， 1997（1）： 15-21.

[2] 谢红辉， 黄兑武， 韦艳明， 等.广西剑麻病虫害发生现状及防治对策[J]. 中国热带农业， 2012（05）： 47-49.

[3] 高建明， 张世清， 陈河龙， 等. 剑麻抗病育种研究回顾与展望[J]. 热带作物学报， 2011（10）： 1 977-1 981.

[4] 陈锦平， 刘清芬. 龙舌兰麻杂种11648斑马纹病的分布及其病原的研究[J]. 热带作物研究， 1980（02）： 122-129.

[5] 赵艳龙， 何衍彪， 詹儒林. 我国剑麻主要病虫害的发生与防治[J]. 中国麻业科学， 2007（06）： 334-338.

[6] 陈孝仁， 王源超. 卵菌基因功能分析方法的研究进展[J]. 农业生物技术学报， 2012（05）： 568-575.

[7] 钟 恒.植物学教学中的卵菌分类地位[J]. 植物学通报， 1995（03）： 59-61.

[8] Agne's Attard. Strategies of attack and defense in plantoomycete interactions， accentuated for Phytophthora parasitica Dastur（syn.P.Nicotianae BredadeHaan）[J]. Journal of Plant Physiology， 2008， 165： 83-94.

[9] 马国胜， 高智谋， 陈 娟. 烟草黑胫病菌研究进展（Ⅲ）[J]. 烟草科技， 2004（02）： 44-48.

[10] 祁云霞， 刘永斌， 荣威恒.转录组研究新技术： RNA-Seq及其应用[J]. 遗传， 2011（11）： 1 191-1 202.

[11] Ansorge W J. Next-generation DNA sequencing techniques[J]. N Biotechnol， 2009， 25（4）： 195-203.

[12] Jewett M C， Oliveira A P， Patil K R， et al. The role of high throughput transcriptome analysis in metabolic engineering[J].Biotechnol Bioproc Eng， 2005， 10： 385-399.

[13] Donson J， Fang Y， Espiritu-Santo G， et al. Comprehensive gene expression analysis by transcript profiling[J]. Plant Mol Biol， 2002， 48： 75-97.

[14] Yan Li， Yiu-Wing Mai， Lin Ye. Sisal fibre and its composites： a review of recent developments[J]. Composites Science and Technology， 2000， 60： 2 037-2 055.

[15] Wang X W， Luan J B， Li J M， et al. Denovo characterization of a whitefly transcriptome and analysis of its gene expression during development[J]. BMC Genomics， 2010， 11： 400.

[16] Wei W， Qi X， Wang L， et al. Characterization of the sesame （Sesamum indicum L.）global transcriptome using Illumina paired-end sequencing and development of EST-SSR markers[J]. BMC Genomics ， 2011， 12： 451.

[17] Shufen Wang， Xiufeng Wang， Qiwei He， et al. Transcriptome analysis of the roots at early and late seedling stages using Illumina paired-end sequencing and development of EST-SSR markers in radish[J]. Plant Cell Reports， 2012， 31： 1 437-1 447.

[18] 欧阳光察， 薛应龙. 植物苯丙烷类代谢的生理意义及其调控[J].植物生理学通迅， 1988（3）： 9-26.