甘蔗葉片響應褐銹病菌(Puccinia melanocephala)侵染的轉錄組分析

高小寧，吳自林，黃詠虹，劉 睿，齊永文

（廣東省科學院生物工程研究所/廣東省甘蔗遺傳改良工程技術中心，廣州510316）

0 引言

甘蔗(Saccharum spp.)是中國最重要的糖料作物，蔗糖占到中國食糖總產量的90%左右，在保障食糖安全方面具有重要地位。由黑頂柄銹菌(P.melanocephala)引起的甘蔗褐銹病(brown rust)于1890年在爪哇首次發現，1981年在中國福建蔗區發現其危害[1]，現已遍布世界各甘蔗種植區[2]。褐銹病主要危害甘蔗葉片，初為黃色小斑點，后擴展為褐色短條斑，病斑處的葉背面產生紅褐色夏孢子堆。發病嚴重時蔗葉提早枯死，分蘗減少，蔗莖變細，可造成甘蔗減產15%～40%，蔗糖含量降低10%～36%[3]。已有的研究結果表明不同甘蔗品種及親本對褐銹病的抗性差異顯著，而感病品種的大面積種植成為病害流行的重要原因[4]。由于褐銹病的常發態勢，致使中國主栽的一批豐產高糖甘蔗品種因高感褐銹病而面臨淘汰，從而極大影響了中國甘蔗產業的持續健康發展。因此，開展甘蔗與褐銹病互作研究，揭示甘蔗對褐銹病的抗性規律進而選育抗病品種就成為主要的研究方向和目標。研究人員通過構建遺傳圖譜，發現并確認了甘蔗抗褐銹病基因Bru1和Bru2；并建立了PCR、SSR、AFLP和RFLP的分子標記技術用于甘蔗對褐銹病的抗性分析[5-9]。李文鳳等[10-13]結合抗褐銹病基因Bru1的分子檢測和抗病性鑒定，系統評價了中國甘蔗栽培原種、主要育種親本及品種對褐銹病的抗性特點。近年來轉錄組測序技術(RNA-seq)已經越來越廣泛地用于植物與病原菌互作的研究中，挖掘出了大量與抗病相關的新基因及代謝途徑，為深入研究植物與病原菌互作的分子機制和植物抗病機理提供了物質基礎和科學依據。在甘蔗病害研究中，研究人員利用RNA-seq技術，分析了甘蔗與黑穗病菌(Sporisorium scitamineum)[14]、甘蔗花葉病毒(Sugarcane mosaic virus)[15]、赤 條 病 菌 (Acidovorax avenae subsp.avenae)[16]、白 條 病 菌 (Xanthomonas albilineans)[17]互作過程中的基因表達特征。然而，至今未見從轉錄水平分析甘蔗響應褐銹病菌侵染的相關報道。為此，本研究以Illumina HiSeq 2500高通量測序技術為基礎，對受褐銹病菌侵染以及健康甘蔗葉片進行轉錄組測序，找到二者在轉錄水平基因表達的差異，分析甘蔗響應褐銹菌侵染的相關基因和途徑，對進一步研究甘蔗抗褐銹病機理和抗病育種提供一些理論依據。

1 材料與方法

1.1 試驗材料

以廣東省科學院生物工程研究所科研試驗基地內種植的高感褐銹病甘蔗品種‘粵糖60號’為材料，測序樣品來自田間管理方式完全相同的同一地塊，根據葉片的感病狀況，試驗材料分為感病組YTCL（感病初期甘蔗葉片）和對照組YTCK（健康的甘蔗葉片）。感病組為受褐銹病感染，初發病的黃色小斑點且無其他病斑及害蟲危害的甘蔗葉片，取病斑處；對照組取與感病組同等嫩度且無任何損傷的面積相等的健康葉片；所取樣品液氮速凍，每處理組3個重復。試驗在廣東省科學院生物工程研究所，于2019年10月—2020年5月份進行。

1.2 總RNA提取、cDNA文庫構建與測序

使用Trizol法提取各樣品的總RNA，1.0%瓊脂糖電泳、NanoDrop 2000、Agilent Bioanalyzer 2100 system(Agilent Technologies,CA,USA)檢測RNA純度、濃度和完整性。檢測合格的RNA樣品，使用NEBNext UltraTM RNA Library Prep Kit試劑盒標準流程構建文庫，委托北京百邁客生物科技有限公司使用高通量測序平臺Illumina HiSeq 2500進行測序，得到原始測序序列。

1.3 數據的過濾與組裝

將測序得到的raw reads，去除帶接頭的reads、未知序列比例大于10%的reads和低質量reads后得到clean reads，作為后續試驗的基礎。

1.4 轉錄組數據分析

使用Hisat2將得到的clean reads比對到割手密參考基因序列(S.spontaneum)[18]，堿基錯配值設為2。使用 Nr(NCBI non-redundant protein sequences)、Nt(NCBI nucleotide sequences)、Pfam(Protein family)、KOG/COG(Clusters of orthologous groups of proteins)、Swiss-prot(A manually annotated and reviewed protein sequence database)，KEGG (Kyoto encyclopedia of genes and genomes)、GO(Gene ontology)等數據庫進行基因功能注釋。

基因表達水平通過RSEM軟件計算FPKM(fragments perkilobase oftranscriptpermillion fragments mapped reads)。使用DESeq2進行差異基因的篩選，篩選條件為差異倍數(Fold Change)≥4且錯誤發現率(FDR,False Discovery Rate)＜0.01。通過GOseq對差異基因進行GO功能富集分析，利用KOBAS進行KEGG通路富集分析。

1.5 差異表達基因的qRT-PCR驗證

以甘蔗的GAPDH為內參基因，選取5個差異表達基因進行qRT-PCR驗證。試驗材料和取樣方法同1.1，采用Trizol法提取總RNA。按照SMARTer PCR cDNA Synthesis Kit試劑盒的方法合成cDNA。使用Primer 5.0軟件參考轉錄組測序序列設計實時熒光定量PCR引物（表1）。按照BBI公司SG Fast qPCR Master Mix說明書進行實時熒光定量PCR。反應體系(20 μL)：2×SYBR Green PCR Master Mix 溶液 10 μL（TOYOBO，日本），10 μmol/L上下游引物各0.5 μL，模板 1 μL，超純水補足 20 μL。PCR 反應程序：95℃5 min，95℃ 10 s，60℃ 45 s，40個循環。3次試驗重復，基因相對表達量的計算采用2-ΔΔCt法。

表1 實時熒光定量PCR的基因及引物

2 結果與分析

2.1 測序數據質量分析及比對統計

經過測序質量控制，6個樣品共得到42.71 Gb Clean Data，各樣品Clean Data都達到6.43 Gb，Q30堿基百分比均高于93%，對照葉片和感病葉片GC含量均大于54%。分別將各樣品的Clean Reads與割手密(S.spontaneum)的基因組進行序列比對，比對效率在69.67%～80.75%；基于比對結果，進行可變剪接預測分析、基因結構優化分析以及新基因的發掘，發掘新基因33,788個，其中27,999個得到功能注釋。

2.2 基因表達分析

在本研究的RNA-seq分析中，通過最大流量算法，采用FPKM作為衡量基因表達水平的指標。對不同樣品的FPKM做箱線圖（圖1），FPKM密度分布顯示，感褐銹病葉片(YTCL)和健康葉片(YTCK)的基因總體表達量在離散度和總體分布度上存在一定差異，表明甘蔗葉片受到褐銹病菌侵染后，部分基因的表達發生了改變。

圖1 FPKM分布箱線圖

2.3 差異表達基因篩選

使用DESeq2進行對照組與處理組間的差異表達分析，將Fold Change≥4且FDR＜0.01作為篩選標準，統計分析所得的差異基因個數，共鑒定出差異基因4716個，其中上調表達2979個，下調表達1737個（圖2a）。進一步對差異表達基因進行聚類分析（圖2b），發現處理組與對照組差異明顯，受到褐銹病菌侵染的甘蔗葉片與健康葉片相比，上調表達基因數量高于下調表達基因。結果表明，大量基因在甘蔗葉片受到褐銹病菌侵染時呈現差異表達模式。

圖2 基因差異表達分析火山圖和熱圖

2.4 差異表達基因GO富集分析

GO功能分析將差異表達基因分為生物過程、細胞組分和分子功能3個功能注釋本體。圖3展示了每類功能之下Top10的功能亞分組。分析發現受褐銹菌侵染的甘蔗葉片涉及生物過程中代謝過程、細胞進程、單體代謝過程出現了大量差異表達基因；分子功能的催化活性和結合2個條目中差異基因數量最多；細胞、細胞組分和細胞膜是細胞組分中最大的3個條目。

圖3 差異表達基因的GO功能分析

為了進一步了解甘蔗葉片對褐銹菌侵染的響應，通過topGO有向無環圖對上調差異表達基因分析富集程度，結果顯示（表2）：分子功能中有較多的GO term顯著富集，包括ATP結合、谷胱甘肽轉移酶活性、解旋酶活性、黃嘌呤脫氫酶活性、L-乳酸脫氫酶活性、纖維素合成酶活性、轉座酶活性、查耳酮合成酶活性；生物過程中的谷胱甘肽代謝過程、毒素分解過程、纖維素生物合成過程、黃嘌呤分解過程、RNA聚合酶II啟動子轉錄調控、氧化還原過程顯著富集；細胞質在細胞組分中顯著富集。

表2 轉錄組測序分析中上調差異表達基因的GO富集分析

2.5 差異表達基因的Pathway分析

利用KEGG數據庫可以分析差異表達基因參與的主要生化代謝途徑和信號轉導途徑。本研究中差異表達基因的KEGG注釋顯示，1919個差異表達基因被注釋到114個通路；進一步按照KEEG的通路類型進行分類，結果如圖4所示：參與苯丙烷生物合成的差異表達基因最多，其次為谷胱甘肽代謝途徑、植物激素信號轉導、淀粉和蔗糖代謝、類黃酮生物合成、碳代謝、糖酵解/糖異生、氨基酸生物合成、甘油磷酸脂代謝和植物與病原菌互作等通路。

圖4 差異表達基因KEGG分類圖

以KEGG數據庫中Pathway為單位，應用超幾何檢驗，與整個基因組背景相比，進行Pathway顯著性富集分析。本研究中上調和下調差異表達基因KEGG通路富集分析結果如表3所示：從上調差異表達基因的富集程度可以看出，當甘蔗葉片受到褐銹菌侵染時，苯丙烷生物合成、類黃酮生物合成、谷胱甘肽代謝等植物抵御生物脅迫有關的代謝途徑顯著富集，大量差異基因上調表達。

表3 差異表達基因KEEG通路富集分析

2.6 甘蔗響應褐銹菌侵染的防御基因

本研究針對植物與病原菌互作途徑的富集基因，共篩選到13個差異表達基因，其中10個上調表達，3個下調表達。植物與病原菌長期互作過程中形成了包括PTI(PAMP-triggered immunity)和ETI(effectortriggered immunity)的免疫系統。在PTI系統中，鈣依賴蛋白激酶(CDPK)和呼吸爆發氧化酶(Rboh)協同作用產生活性氧并通過過敏性壞死反應提高植物的抗病性。真菌細胞壁成分如寡糖和幾丁質真菌PAMP，可以被寄主受體識別并激活防御反應。在甘蔗受到褐銹菌侵染后，編碼幾丁質激發子受體激酶(CERK1)、鈣依賴蛋白激酶(CDPK)、呼吸爆發氧化酶(Rboh)、鈣調蛋白(CaM)、PR1等基因上調表達。此外，在甘蔗響應褐銹病菌侵染時，RPM1、RPS2、RPS4等相關抗病蛋白基因同樣出現上調表達。

2.7 差異表達基因的qRT-PCR分析

以甘蔗的GAPDH基因為內參基因，對選取的PR1、PR10、CDPK、NPR1、CTR1等5個上調表達基因進行qRT-PCR驗證，結果顯示這5個差異表達基因的表達趨勢與轉錄組數據一致，說明轉錄組測序分析結果是可靠的。

3 討論

甘蔗褐銹病自20世紀80年代初期在中國甘蔗種植區發現以來，因其適應性廣、產孢量大、傳播速度快、流行性強，加之單一甘蔗品種的大面積種植，致使病害迅速蔓延，已發展成為甘蔗種植中重要的葉部病害，而篩選和培育抗褐銹病的甘蔗種質和品種是防控該病害最為經濟和有效的途徑[19]。在前期的研究中，研究者基于cDNA-AFLP和SSH方法分析了甘蔗與褐銹病互作過程中的差異表達基因[20-21]，但是由于研究技術的局限，在一定程度上限制了對甘蔗響應褐銹菌侵染分子機制的認知水平。本研究首次利用高通量測序技術(RNA-seq)分析了感病品種粵糖60受到褐銹菌侵染的轉錄組，更為全面地揭示了基因表達特征，鑒定到4716個差異表達基因，并分析了差異表達基因的功能及其參與的代謝通路。取得的研究結果為揭示甘蔗與病菌互作的分子機制提供了豐富的數據資源。

植物內源激素在植物與病原菌互作中發揮極其重要的作用。當植物受到病原菌侵染時，通過產生小的信號分子并協調激素信號通路來激活信號網絡，從而啟動針對病原菌攻擊的防御機制。水楊酸(SA)、茉莉酸(JA)、乙烯(ET)、油菜素類脂(BR)、脫落酸(ABA)等激素類物質是植物抗病信號轉導途徑中的重要調控因子[22]。SA作為植物防御反應中一個重要的信號分子，其介導的抗病信號通路主要針對活體營養型或半活體營養型病菌的侵染[23]。TGA基因家族是bZIP轉錄因子家族的重要亞家族之一，能夠與靶基因啟動子上的as-1區域相結合，激活或抑制下游靶標基因的轉錄，在SA介導的植物抗病信號轉導中起重要的調控作用[24]。郭雙元等[25]研究發現小麥TaTGA2.2受條銹菌和外源激素SA的誘導表達，TaTGA2.2可能通過與NPR1互作參與SA信號途徑介導的小麥對條銹菌的防御反應。本研究的試驗結果顯示：甘蔗TGA4和PR1在受到褐銹菌侵染后顯著上調表達，表明SA信號途徑在甘蔗與褐銹菌互作中發揮著相似作用。

轉錄因子(Transcription factor，TF)作為一類重要的調控基因，能夠與基因啟動子區域中的順式元件特異性結合，激活或者抑制目的基因的表達。WRKY、MYB、NAC、ERF等轉錄因子被發現是寄主植物響應生物脅迫的重要調控因子。已有的研究結果發現WRKY轉錄因子中的很大一部分都參與了植物對生物脅迫的響應，在植物多種免疫系統中發揮著重要的調控作用，是植物響應生物脅迫尤其是病原微生物的重要轉錄因子家族[26]。研究發現擬南芥AtWRKY33的過表達增強了對死體營養型真菌灰霉病菌(Botrytis cinerea)和黑斑病菌(Alternaria brassicicola)抗性，但是增強了對丁香假單胞桿菌(Pseudomonas syringae)的感病性[27]；進一步的研究發現AtWRKY33直接被MPK3/6磷酸化進而調控植保素合成基因的表達，將AtWRKY33磷酸化位點突變后，轉基因植株的抗病性明顯削弱[28]。本研究的結果顯示甘蔗WRKY33、WRKY22在褐銹菌侵染后顯著上調表達，但是WRKY轉錄因子在甘蔗與褐銹菌互作中的功能及其作用機理還有待進行深入研究。

在植物信號轉導途徑中，鈣離子是普遍存在的第二信使，參與調控了大多數細胞生理代謝的過程。鈣信號在植物抗病信號轉導途徑中起重要的作用，并通過鈣離子結合蛋白調節信號的傳遞[29]。鈣調蛋白(CaM)是植物細胞內最重要的鈣離子結合蛋白，在植物生長發育和抗逆抗病中具有重要的調控作用。鈣依賴的蛋白激酶(CDPK)是細胞內鈣離子信號的受體，同時具有鈣離子受體蛋白和Ser/Thr蛋白激酶的功能，廣泛參與植物逆境脅迫的分子響應。本研究的KEGG代謝通路富集分析結果顯示，在甘蔗受到褐銹菌侵染后CDPK、CaM、CML等與Ca2+相關的基因顯著上調表達，這與孫云南等[30]在茶葉葉片響應茶餅病侵染的轉錄組分析結果基本一致。

植物與病原物的互作是一個復雜且連續的過程，涉及植物抗病相關基因和病原物致病相關基因的遺傳和變異、表達和調控以及它們之間的相互聯系和影響。本研究僅對親和性互作條件下甘蔗葉片響應褐銹菌侵染的基因表達特征進行了分析，在褐銹病菌侵染甘蔗的不同階段以及甘蔗-褐銹病菌互作的非親和關系中，甘蔗在轉錄水平上的表達特征將是后續需要開展的研究。盡管，轉錄組分析作為挖掘功能基因的重要途徑之一，已被廣泛應用于植物與病原菌互作研究中，但是，單個組學數據具有片面性，可能無法完整描繪整個生物學過程。因此，將基因組學、轉錄組學、蛋白組學、代謝組學和代謝調控網絡有機結合，基于各組學的特點及其各自相適應的分析技術，通過多組學的聯合研究將能更全面的闡述甘蔗與褐銹病菌的互作過程及其規律，進而為甘蔗抗褐銹病的種質創新和品種改良奠定基礎。

4 結論

本研究利用Illumina Hiseq 2500高通量測序平臺，在轉錄組水平上分析了甘蔗響應褐銹病菌侵染的基因表達特征，篩選出差異表達基因4716個，其中2979個上調表達，1737個下調表達。GO和KEGG分析結果發現，甘蔗葉片受到褐銹病菌侵染時，差異表達基因主要參與苯丙烷生物合成、谷胱甘肽代謝途徑、植物激素信號轉導、類黃酮生物合成、植物與病原菌互作等代謝途徑。本研究獲得的甘蔗葉片與褐銹病菌互作的轉錄組數據，將為今后開展甘蔗抗/感病基因及其功能，解析甘蔗與褐銹病菌互作的分子機理提供數據支撐。