煙草疫霉侵染前后劍麻葉片轉錄組學研究

汪平　高建明　楊峰　鄭金龍　劉巧蓮　陳河龍　易克賢

摘 要 為克隆劍麻受斑馬紋病病原菌煙草疫霉侵染后的應答基因，以病原菌侵染前后的劍麻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.