大麥在網斑病菌侵染后的轉錄組學分析

李 桃,汪軍成,姚立蓉,李葆春,馬小樂, 楊 柯,司二靜,尚勛武,王化俊,孟亞雄

(1.甘肅農業大學農學院，甘肅蘭州 730070; 2.甘肅省干旱生境作物學重點實驗室/甘肅省作物遺傳改良與種質 創新重點實驗室，甘肅蘭州 730070; 3.甘肅農業大學生命科學技術學院，甘肅蘭州 730070)

栽培大麥(HordeumvulgareL.)是隸屬禾本科的一年生作物，主要用途是飼用、啤酒釀造及食品加工等。由子囊菌亞門真菌網斑病菌(Pyrenophorateres，Ascomycotina)引起的大麥網斑病在國內大麥主要產區開始流行，并有逐年加重的趨勢[1-2],其在世界其他國家如澳大利亞[3]、新西蘭和丹麥[4]、南非[5]、美國[6]等地也均有發生。網斑病可抑制大麥生長，嚴重時會造成大麥減產、品質變劣，甚至會造成絕收，是限制大麥生產的主要原因之一。為解決網斑病對大麥的危害，大麥抗網斑病機制研究成為近年來被關注的重點，而通過基因組學方法鑒定和分離新基因是研究大麥抗病機制的有效方法之一。

有關大麥抗病基因分離和鑒定的報道并不多，且大多與大麥抗逆境脅迫有關，如HvNAC1[8]、Mlo[9]、ICS和PAL[10]等家族的基因。黃德華等[7]證實了大麥GER4基因簇的啟動子表達與白粉病真菌病原體攻擊有關，但有關大麥抗真菌類病害基因的研究并不多見。大麥基因組測序的完成[11]以及高通量測序技術的日益成熟，對大麥抗病基因的深入研究提供了方便。

植物的抗逆境脅迫不僅受多個信號傳導途徑的調控，還受多個基因協同表達和轉錄因子的調控[12]。有研究表明，與差異表達基因趨勢相同的蛋白大多與植物抗逆相關；乙烯、赤霉素、苯丙烷類代謝參與甘蔗黑穗病侵染的響應過程[13]。因此，掌握關鍵基因的調控機制是實現促進多基因表達以提高植物綜合抗病性的有效途徑。研究發現，次生代謝產物[14]、植物激素信號轉導[15]、MAPK信號傳導[16]、鞘脂[17]等代謝途徑在植物的抗病防御反應中發揮著重要作用，如NAC[8]、NF-Y[18]、HD-Zip[19-20]和WRKY[21]等家族的基因。為了解已知抗逆境基因是否與抗網斑病有關，本研究擬比較抗性和敏感型大麥材料對網斑病侵染后響應基因的表達特性，并通過比較其轉錄組(RNA-sep)探索大麥中有關抗病代謝途徑，以期獲得與抗病性有關的基因，了解大麥對網斑病菌侵染的響應機制，為大麥抗網斑病的進一步深入研究提供參考。

1 材料與方法

1.1 材料及設計

由甘肅農業大學麥類課題組提供抗病材料BYT-CYA3(B)和敏感型材料美41/I(M)。種子用72%乙醇處理 15 s，10%次氯酸鈉處理 10 min，用蒸餾水沖洗6次后，置于鋪有雙層無菌濾紙的無菌培養皿里發芽，噴霧保持濾紙濕潤；一周后，將發芽種子移栽入花盆。每盆種植5 株，每個品種共30盆。待幼苗長至兩葉一芯期，對25盆大麥噴施由甘肅農業大學麥類課題組提供的活化網斑病病原菌菌株(LQA)懸浮液(濃度約106·mL-1)，其余5盆為空白對照。接菌完成后置于21 ℃黑暗培養，于0、3、6、12、24和72 h后分別剪取接菌組與對照組葉片3～5 g，液氮速凍，于 -80 ℃保存備用。3次生物學重復。

1.2 RNA提取和cDNA 文庫構建

采用Plant RNA Kit試劑盒(OMEGA，中國上海)按照說明提取大麥葉片總的RNA。在RNA樣品中加DNaseI于 37 ℃水浴30 min，防止DNA污染。用超微量分光光度計檢測總RNA的濃度。按照試劑盒說明富集mRNA和去除rRNA；用DNA探針雜交rRNA，再用RNaseH和DNaseI純化RNA；用打斷buffer使RNA片段化；反轉錄合成cDNA，組建cDNA文庫。

1.3 轉錄組測序和數據過濾

在BGISEQ-500測序平臺對構建好的cDNA文庫進行雙末端測序，所有測序工作由華大基因生物公司(深圳)完成。

測序所得的數據為原始數據(raw reads)。使用在線的SOAPnuke軟件(https://github.com/BGI-flexlab/SOAPnuke)過濾原始數據，去除接頭污染、未知堿基N過高的reads和質量低的reads，獲得clean reads。利用在線軟件HISAT(Hierarchical Indexing for Spliced Alignment of Transcripts)將clean reads與大麥參考基因組序列IBSC_PGSB_v2.39進行比對。使用 Bowtie2 (http://bowtie-bio.sourceforge.net/Bowtie2/index.shtml)將clean reads比對到參考基因序列得到比對結果[22]，對比對結果使用RSEM(http://deweylab.biostat.wisc.edu/rsem/rsem-calculate-expression.html)[23]計算差異表達基因。分析轉錄組測序數據，包括轉錄本的覆蓋度及分布。

1.4 基因的功能注釋

將轉錄組測序鑒定到的所有基因合并為一個基因集數據庫，使用getorf檢測所有基因的 ORF；使用hmmsearch將ORF與轉錄因子蛋白結構域(數據來源于TF)進行比對，根據植物轉錄因子數據庫(PlantTFDB)描述的轉錄因子家族特征對所有基因進行功能鑒定。使用HMMER程序的 hmmpfam與Pfam 23.0對轉錄因子與PlantTFDB進行比對，搜索蛋白結構域以鑒定編碼轉錄因子的基因，使用DIAMOND(https://github.com/bbuchfink/diamond)軟件將轉錄組測序鑒定的所有基因比對至植物抗病基因數據庫(plant resistance gene database，PRGdb)進行注釋，依據比對覆蓋度(query coverage)和一致性(identity)等條件對注釋結果進行篩選。對差異表達基因進行GO和KEGG功能注釋富集分析(Wiki https://en.wikipedia.org/wiki/Hypergeometric_distribution)。本研究的轉錄組原始數據可以從NCBI數據庫獲取，登陸號SUB6851949。

1.5 網斑菌侵染大麥葉片響應的差異表達基因的鑒定和qRT-PCR驗證

試驗以接菌后0 h樣品為對照(B-CK，M-CK)，比較抗病和敏感型材料接種不同時間的樣品(B-3h，M-3h，B-6h，M-6h，B-12h，M-12h，B-24h，M-24h，B-72h，M-72h)與相應對照材料基因的表達水平,參照Audic等[24]的方法，進行基因表達量分析，滿足錯配率[25](false discovery rate，FDR) ≤0.001且表達差異在2倍以上的基因為差異表達基因。

利用Primer-NCBI在線設計軟件(https://www.ncbi.nlm.nih.gov/tools/primer-blast/#)設計引物(表1)，對差異表達基因進行qRT-PCR驗證分析。將樣品3次生物學重復的 RNA 采用第1鏈cDNA合成試劑盒(TIANGEN，中國北京)按照說明合成反轉錄cDNA，以反轉錄cDNA為模板，以大麥Actin為內參，采用SYBR Green RT-PCR 試劑盒(TIANGEN)，按照說明使用ABI 7500 Real Time PCR擴增儀進行熒光定量檢測，結果按照2-ΔΔCT計算相對表達量。

2 結果與分析

2.1 測序結果質量評估

為了研究大麥在不同時間段響應網斑病菌侵染的轉錄組學變化特征，采集大麥接菌處理0 h、3 h、6 h、12 h、24 h和72 h的幼苗葉片，提取總RNA，構建了12個cDNA文庫，用BGISEQ-500平臺對其進行轉錄組學測序，平均每個文庫測序產生10.63 Gb的數據，共獲得1 401 716 667個raw reads，去除含接頭的reads、未知堿基N含量大于5%的reads和低質量的reads，共得到 127 601 667條clean reads，過濾后12個cDNA文庫的clean reads占raw reads的比例均在80%以上；clean reads的Q20比率均在97%以上，其 Q30比率均在86%以上，每個文庫clean reads數均在101.61×106條以上(表2)，獲得的測序數據與質量符合下一步分析要求。

樣品比對到基因組的平均比對率為 87.61%，比對到基因集的平均比對率為 74.93%；共檢測到基因數為35 545個，其中已知基因30 030個，新基因為5 515個。

表1 實時熒光定量PCR所用的引物Table 1 Primers for real-time quantitative PCR

表2 供試材料Reads 的質量統計Table 2 Quality statistics of transcriptome sequencing after filtering

2.2 基因注釋分析

對檢測到的30 030個已知基因進行功能注釋分析，發現只有1 977個基因在TF數據庫有注釋，可歸類到58個基因家族，注釋到基因數最多的是MYB家族(228)，其次是AP2-EREBP(170)、bHLH(157)、FAR1(146)、NAC(145)、ABI3VP1(125)和WRKY家族(107)，注釋到基因數最少的有CAMTA(1)、HRT(1)、S1Fa-like(1)(圖1 A)。在PRGdb數據庫中，有3 427個基因被注釋到13個結構域中，其中注釋到RLP結構域中的基因數目最多，其次是NL、CNL、N和TNL，注釋到的基因數分別有789、508、449和352個，注釋到的基因數目最少的結構域是TNL-OT(圖1 B)。

A：在TF數據庫中的比對結果；B：在 PRGdb數據庫中的比對結果。

2.3 網斑病菌侵染后差異表達基因(DEGs)分析

與網斑病菌侵染大麥0 h處理對照相比，在網斑病菌侵染大麥3 h、6 h、12 h、24 h和72 h后，抗病材料中分別有9 842、12 073、9 065、4 674和5 939個差異表達基因；敏感型材料中分別有10 149、11 322、11 142、2 305和5 023個差異表達基因。隨著侵染時間的持續，抗病材料與敏感型材料上調和下調的差異表達基因(DEGs)數整體呈先增后減的趨勢。抗病材料中的上調差異表達基因數目均大于敏感型材料(圖2)。

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

對差異表達基因進行GO富集分析，發現這些差異表達基因在2種材料不同階段的GO富集結果相似。在生物過程(biological process)中，差異表達基因主要富集的條目有生物調節(biological regulation)、細胞成分組織或生物發生(cellular component organization or biogenesis)、細胞過程(cellular process)、定位(localization)、代謝過程(metabolic process)、對刺激的反應(response to stimulus)。在細胞組分(cellular component)中，差異表達基因主要富集的條目為細胞(cell)、大分子復合物(macromolecular complex)、細胞膜(membrane)、細胞膜片段(membrane part)、細胞器(organelle)、細胞器的部分(organelle part)。在分子功能(molecular function)中，差異表達基因主要富集的條目有整合(binding)、催化活力(catalytic activity)、轉運活動(transporter activity)。

在抗病材料BYT-CYA3中，網斑病菌侵染大麥接觸階段到侵入階段(接種3～6 h)，富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目均明顯增加；在網斑病菌侵染的擴展階段(接種12～24 h)，富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目較侵入階段均驟減；而在大麥發病階段(接種72 h)，與擴展階段相比，富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目均明顯增加。在敏感型材料美41/I中，從網斑病菌侵染大麥到擴展前期階段(3～12 h)，富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目沒有明顯變化；而從擴展后期到發病階段，富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目均明顯減少。

圖2 不同時間的差異表達基因數Fig.2 Number of differentially expressed genes at different time

在網斑病菌侵染的接觸階段(3 h)，抗病材料中富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目明顯少于敏感型材料，其中，最為明顯的條目有細胞成分組織或生物發生(cellular component organization or biogenesis)和細胞(cell)。在網斑病菌侵入階段，抗病材料中富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目明顯多于敏感型材料，最為明顯的條目有細胞成分組織或生物發生(cellular component organization or biogenesis)、細胞過程(cellular process)、代謝過程(metabolic process)、細胞(cell)、細胞器(organelle)、細胞器的部分(organelle part)、整合(binding)、催化活力(catalytic activity)。在網斑病菌侵染大麥的擴展階段(12～24 h)，抗病材料中富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目明顯少于敏感型材料，最為明顯的條目有生物調節(biological regulation)、代謝過程(metabolic process)、細胞(cell)和整合(binding)；在大麥發病階段(72 h)，抗病材料中富集到生物過程(biological_process)、細胞組分(cellular component)和分子功能(molecular function)三大功能條目中的差異表達基因數目均明顯多于敏感型材料。

2.5 差異表達基因的KEGG富集分析

通過KEGG數據庫對差異表達基因進行富集分析，發現這些差異表達基因在2種材料的不同階段的KEGG富集結果比較相似。差異表達基因主要富集在12條通路上，即環境適應(environmental adaptation)、核苷酸代謝(nucleotide metabolism)、脂質代謝(lipid metabolism)、全球概覽地圖(global and overview maps)、碳水化合物代謝(carbohydrate metabolism)、其他次生代謝產物的生物合成(biosynthesis of other secondary metabolites)、氨基酸代謝(amino acid metabolism)、遺傳翻譯(translation)、轉錄(transcription)、折疊、分類和降解(folding,sorting and degradation)、信號轉導(signal transduction)、運輸和分解代謝(transport and catabolism)(圖3)。

在抗病材料BYT-CYA3中，從網斑病菌侵染大麥接觸階段到侵入階段(3～6 h)，富集到各代謝通路上的差異表達基因數目均明顯增加，最明顯的通路有全球概覽地圖(global and overview maps)、碳水化合物代謝(carbohydrate metabolism)、其他次生代謝產物的生物合成(biosynthesis of other secondary metabolites)、氨基酸代謝(amino acid metabolism)、遺傳翻譯(translation)、轉錄(transcription)、折疊、分類和降解(folding,sorting and degradation)、信號轉導(signal transduction)、運輸和分解代謝(transport and catabolism)；在網斑病菌侵染的擴展階段(12～24 h)，富集到各代謝通路上的差異表達基因的數目較侵入階段均劇減；而在大麥的發病階段(72 h)，與擴展階段相比，富集到各代謝通路上的差異表達基因數目均明顯增加(圖3 A)。在敏感型材料美41/I中，從網斑病菌侵染大麥接觸階段到侵入階段(3 ～6 h)，發現富集到各代謝通路上的差異表達基因數目均明顯增加；而網斑病菌侵染的擴展階段到大麥發病(12～72 h)，富集到各代謝通路上的差異表達基因數目均較前階段明顯減少(圖3 B)。

在網斑病菌侵染大麥的接觸階段(3 h)，抗病材料中富集到各代謝通路上的差異表達基因多于敏感型材料，差異較大的代謝通路有轉錄(transcription)和全球概覽地圖(global and overview maps)；從網斑病菌侵染大麥侵入階段到擴展階段(6～24 h)，抗病材料和敏感型材料富集到各代謝通路上的差異表達基因數目均不斷減少；在大麥發病階段(72 h)，抗病材料中富集到各代謝通路上的差異表達基因多于敏感型材料，差異較大的代謝通路有轉錄(transcription)和全球概覽地圖(global and overview maps)(圖3)。

2.6 差異表達基因功能分析

在抗病材料中，共獲得差異表達基因435個(圖4A)。對435個差異表達基因進行KEGG途徑富集分析，發現差異表達基因較多的通路有代謝途徑(metabolic pathways；ko01100)、MAPK信號傳導途徑-植物(MAPK signaling pathway-plant；ko04016)、植物-病原體相互作用(plant-pathogen ；ko04626)、植物激素信號轉導(plant hormone signal transduction；ko04075)和次生代謝產物的生物合成(biosynthesis of secondary metabolites；ko01110)(圖4B)。其中有19.77%的差異表達差因基因通路產物是白介素-1受體相關激酶(interleukin-1 receptor-associated kinase)(圖4B)。對435個差異表達基因進行GO途徑富集分析，在生物過程(biological_process)中，差異表達基因主要富集的條目有細胞過程(cellular process)、代謝過程(metabolic process)。細胞組分(cellular component)中，主要富集的條目為細胞(cell)、細胞器(membrane)、細胞器的部分(membrane part)。在分子功能(molecular function)中，主要富集的條目為整合(binding)和催化活力(catalytic activity)(圖4C)。推測這幾條代謝通路與大麥抗病性有關。

為了分析差異表達基因的表達模式，對182個差異表達基因進行KEGG和GO途徑分析，發現差異表達基因主要富集在過氧化物酶體(peroxisome；ko04146)、代謝途徑(metabolic pathways；ko01100)、MAPK信號傳導途徑-植物(MAPK signaling pathway - plant；ko04016)、植物-病原體相互作用(plant-pathogen；ko04626)、植物激素信號轉導(plant hormone signal transduction；ko04075)、代謝途徑(biosynthesis of secondary metabolites；ko01110)。根據182個差異表達基因的功能在代謝通路上的定位，發現這些基因參與了維持細胞死亡的防御反應、維持病原體入侵和維持植物防御中活性氧的積累(圖5B)。這182個差異基因可以作為大麥抗病基因進一步研究的候選基因。

A：網斑病菌侵染抗病材料的差異表達基因韋恩圖分析；B：網斑病菌侵染抗病材料的差異表達基因的KEGG分類注釋；C：表示網斑病菌侵染抗病材料的差異表達基因的GO分類注釋。

灰框是通過轉錄組序列檢測到的KEGG酶編碼基因。

2.7 實時定量熒光PCR驗證

為驗證轉錄組測序的結果，從篩選的182個抗病相關基因中選取信號轉導(signal transduction)、MAPK信號傳導途徑-植物(MAPK signaling pathway-plant)、其他次生代謝產物的生物合成(biosynthesis of other secondary metabolites)和植物-病原體相互作用(plant-pathogen)代謝通路的10個基因(表1)，利用實時熒光定量PCR對其網斑病菌侵染前后的表達量進行分析，發現這10個基因的熒光定量與轉錄測序結果基本一致，且均上調表達。

3 討 論

有研究表明，植物在與病原微生物的長期互作、協同進化過程中逐漸形成一系列的防衛機制，這些機制在誘導后才能快速、充分地表達，從組織結構機制、生理機制和分子機制看，與脅迫相關的基因在病原微生物侵染后的不同階段具有不同的表達模式[26]。這現象可能是大麥調節體內生理代謝，以形成相應的防御機制。

本試驗發現，大麥網斑病病菌侵染后無論是抗病材料還是敏感型材料，有大量不同家族的基因出現表達量的變化，兩種材料在不同侵染階段的差異表達基因的數目不同，這可能是其抗病性不同的根本原因。比較2個材料在網斑病菌侵染3 h、6 h、12 h、24 h和72 h后基因表達的異同，推測不同基因在網斑病菌侵染后，其轉錄水平調控方式在不同階段存在差異。證明了王海波等[27]和翁 宇等[28]的研究結果：植物通過改變基因表達量來響應脅迫。

抗病材料和敏感型材料在病菌侵染后的差異表達基因有所相同，這表明不同抗病材料對逆境的響應機制不同。錢 前等[29]提出，植物有天然免疫抗性。本研究發現，即便在相同時間2個材料中存在的共同差異表達基因，其上調、下調表達程度也存在差異，推測這可能與大麥的基礎抗性有關，同類型材料對相同脅迫有某種相同的響應機制，但響應程度不同。

本研究對差異表達基因做了PRGdb和TF注釋，發現差異表達基因大多富集在CNL、N、NL、RLP和TNL等結構域中。有研究表明，富集在CNL、N、NL、RLP和TNL等結構域的基因與植物防御系統的抗病性有關[30-31]。目前，在植物中發現的與抗病性有關的基因家族有b ZIP、AP2/EREBP、WRKY、MYB、NAC等[32]。本研究中，差異表達基因在AP2/EREBP、NAC、WRKY和bHLH家族中分布數目較多。植物中特有轉錄因子家族是NAC，其參與植物對逆境脅迫的響應和，調控植物的抗病性[33-35]。王立國等[36]、董 帥等[8]和張立全等[37]完成了在陸地棉、大麥、黃花苜蓿中對NAC轉錄因子相關基因的克隆。證明了NAC參與調控植物的抗逆性。除了NAC，植物中還有一種特有的家族是鋅指蛋家族WRKY。研究發現，WRKY類轉錄因子參與植物的形態發育與對病原體的防御響應[38]。目前，已有對WRKY基因的分析和結構域特點及其生物學功能的研究[39-40]和其在植物防御反應中的作用[41-42]方面的報道。本研究發現，MYB家族被注釋到最多的基因數目。在植物中，MYB通過有效調控類黃酮物質的生物合成以調節機體免疫力、抗氧化、抗衰老及抵抗病毒能力[43]。已有研究表明，MYB與小麥對紋枯病、赤霉病抗性有關[44]。對于bHLH家族的研究表明，植物bHLH基因家族參與調控植物的生長和發育過程[45]。而對于AP2/EREBP家族的研究表明，在擬南芥中得到AP2/EREBP家族的獨特成員是At4g13040基因，其參與SA介導的疾病防御1-APD1的Apetala 2家族蛋白合成，使病原體接種后無法誘導表達[46]。

本研究篩選出來182個可能與大麥抗病性有關的基因，其中，35個差異表達基因參與植物激素信號轉導(plant hormone signal transduction；ko04075)代謝途徑的基因(圖3B)，28個差異表達基因富集到信號轉導(signal transduction)，其中有7個差異表達基因參與合成LRR受體樣絲氨酸/蘇氨酸蛋白激酶FLS2(LRR receptor-like serine/threonine-protein kinase FLS2)(附表1)。苗麗麗等[47]研究發現，蔗糖非發酵相關蛋白激酶2(SnRK2)是植物特有的一類絲氨酸/蘇氨酸蛋白激酶，SnRK2作為典型的多功能調節因子參與響應各種逆境脅迫；郭艷玲等[14]提出次生代謝產物參與植物的抗病防御反應，作為生化壁1壘防御病原物侵染，和信號物質共同參與植物的抗病反應。本研究中發現,45個差異表達基因參與其他次生代謝產物的生物合成(biosynthesis of other secondary metabolites)代謝途徑的基因，其中9個差異表達基因參與合成過氧過氧化物酶體(peroxidase；K00430)，而過氧化物酶是病原真菌抗氧化防御系統的主要組分[48]； 36個差異表達基因富集到環境適應(environmental adaptation)，其中5個差異表達基因是抗病蛋白RPM1(disease resistance protein RPM1)，8個差異表達基因參與合成LRR受體樣絲氨酸/蘇氨酸蛋白激酶FLS2(LRR receptor-like serine/threonine-protein kinase FLS2)。這些基因的結構與功能需要更多試驗與探究。

本試驗從分子水平上初步探索了大麥對網斑病侵染的響應機制，提供了抗病基因深入研究的候選基因，為大麥抗病基因的研究奠定基礎。