小麥條銹菌鑒別寄主中4抗條銹性的轉錄組學分析

湯亞琪，蘇文文，李 瑩，王保通，李 強

(西北農林科技大學植物保護學院/旱區作物逆境生物學國家重點實驗室，陜西楊凌 712100)

小麥條銹病是世界范圍內小麥主產區最重要的病害之一[1]，其病原菌條形柄銹菌小麥?；?PucciniastriiformisWest. f. sp.triticiEriks & Henn，Pst)為專性寄生菌[2]，生理小種致病性變異迅速且頻繁，容易造成主栽小麥品種抗條銹的頻繁“喪失”[3]，導致小麥條銹病的爆發流行，嚴重威脅著全球小麥的安全生產[4]。

中4是以普通小麥(TriticumaestivumL.)(2n=6x=42，AABBDD)與中間偃麥草(Thinopyrumintermedium)(2n=6x=42，E1E1E2E2XX)經有性雜交、選育獲得的異源八倍體小偃麥(AABBDDXX)[5]。自1983年列為中國小麥條銹菌鑒別寄主至今，對我國現有小麥條銹菌流行小種均表現為抗病[6]，其抗性穩定且持久，具有較高的研究和利用價值。但由于中4是異源八倍體小偃麥，與六倍體普通小麥經常規雜交后獲得穩定遺傳后代的難度較大，導致其抗性研究嚴重滯后[7]。藺瑞明等[7]將中4分別與中國春phb1突變體和條銹病感病品種銘賢169雜交，對其雜交后代進行遺傳分析，發現其對條銹菌小種CYR31和CYR32的抗病性由1對顯性基因控制。曹世勤等[8]利用基因推導法發現中4含有未知抗條銹病基因。甘肅省天水市農業科學研究所以中4為親本，先后選育出中梁22號[9]、中93444[10]、中93447[11]等多個具有較好豐產性和抗病性的品種(系)材料。但是其抗條銹病機制目前尚未見報道。

本研究利用轉錄組測序，篩選中4接種條銹菌后的差異表達基因，以期為進一步解析中4抗條銹性分子機制提供科學依據。

1 材料與方法

1.1 供試材料

供試小麥品種為我國小麥條銹菌鑒別寄主中4和條銹病感病品種銘賢169，小麥條銹菌為我國小麥條銹菌流行小種CYR32，均由西北農林科技大學植物保護學院小麥病原真菌監測及抗病遺傳實驗室提供。

1.2 樣品采集

挑選籽粒發育良好、顏色形狀一致的中4種子，催芽后播種在7 cm × 7 cm的小花盆內。待小麥幼苗長至第一葉完全展開時，接種條銹菌CYR32，接種0、12、24和48 h后，分別采集接種葉片樣品，以未接種葉片(0 h)為對照，每個時間點采集3個生物學重復，用于轉錄組數據質量測定。接種48 h后的材料在人工氣候箱(溫度15～17 ℃、濕度90%、16 h光照/8 h黑暗、光照強度10 000 lx)繼續培養16 d，用于觀察表型，以銘賢169作為感病對照。

1.3 RNA提取及文庫構建

樣品總RNA的提取按照Plant RNA Midi Kit(OMEGA，美國)試劑盒說明書進行，獲取總RNA后，檢測RNA樣品的純度、濃度和完整性，用帶有Oligo(dT)的磁珠富集真核生物mRNA，加入Fragmentation buffer將mRNA隨機打斷，以First Strand cDNA Synthesis Kit(OMEGA，美國)合成第一條cDNA鏈，然后加入緩沖液、dNTPs、RNase H和DNA polymerase I合成第二條cDNA鏈，利用AMPure XP beads純化雙鏈cDNA，純化后的雙鏈cDNA須先進行末端修復、加A尾并連接測序接頭，再用AMPure XP beads進行片段大小選擇，最后通過PCR富集得到 cDNA文庫。

1.4 轉錄組測序及功能注釋

通過Illumina高通量測序獲得12個cDNA文庫，用FastQC對cDNA文庫的堿基質量值(Q30)、堿基含量分布(GC含量)和重復相關性(r2)等指標進行檢測，以判斷文庫數據質量并評估其是否滿足后續分析。用HISAT2對reads進行高效比對，用StringTie進行從頭組裝，基于中國春全基因組(https://wheat-urgi.versailles.inra.fr/Seq-Repository)進行有參組裝。

將完成組裝的轉錄本數據用非冗余(non-redundant，NR)數據庫和基因本體(gene ontology，GO)數據庫進行比對注釋，用真核生物同源蛋白簇(eukaryotic orthologous groups，KOG)、蛋白質直系同源簇(cluster of orthologous groups of proteins，COG)以及基因的進化譜系：非監督直系群(evolutionary genealogy of genes:non-supervised orthologous groups，eggNOG)等數據庫進行基因富集聚類，用京都基因與基因組百科全書(kyoto encyclopedia of genes and genomes，KEGG)數據庫進行代謝通路的富集分析。

1.5 實時熒光定量PCR

選取轉錄組測序得到的差異表達基因，根據其全長cDNA序列，使用Primer 5設計特異性引物(表1)，以小麥TaEF為內參基因。提取中4各個樣品的總RNA并反轉錄為cDNA，按照ChamQTMSYBR○RqPCR Master Mix(Vazyme，中國)說明書的反應體系和條件，使用QuantStudio 5(ABI，美國)進行qRT-PCR，按照2-△△Ct計算基因相對表達量，重復3次。

表1 熒光定量PCR所用引物Table 1 Primers used for fluorescence quantitative PCR

2 結果與分析

2.1 中4抗條銹菌鑒定結果

從圖1可以看出，中4葉片只有少許壞死斑且不連片，表現為高抗，而銘賢169產生大量夏孢子堆。

2.2 轉錄組數據質量

取接種條銹菌0、12、24和48 h后的中4葉片樣品，每個處理3個重復，共計12個樣品進行轉錄組測序分析，共獲得130.51 Gb Clean Data，各樣品Clean Data均達到8.42 Gb；堿基質量值Q30的堿基百分比為92.83%～94.14%，表明堿基識別可靠，準確度高；G和C、A和T的含量為58.33%～60.83%。分別將各樣品的Clean Reads與中國春全基因組IWGSC RefSeq v1.0進行序列比對，比對效率為60.29%～69.03%，由于中4為異源八倍體小偃麥(AABBDDXX)[6]，其X染色體組與普通小麥親緣關系較遠，所以轉錄組的比對效率沒有出現高度匹配。基于比對結果進行基因表達量分析、功能注釋、富集分析以及新基因的發掘[8]，發掘新基因20 788個，其中 15 069個得到功能注釋。關于樣品相關性，r2值越接近1，說明兩個重復樣品相關性越強且樣品重復性良好，測序結果顯示，中4接種條銹菌CYR32 0 h后，3個樣品的r2值為0.960～ 0.983；接種條銹菌CYR32 12 h后，3個樣品的r2值為0.958～0.975；接種條銹菌CYR32 24 h后，3個樣品的r2值為0.972～ 0.982；而接種條銹菌CYR32 48 h后，3個樣品的r2值為0.824～0.932，表明樣品間整體相關性較強，重復性良好。

2.3 差異表達基因數目及染色體定位

以Fold Change≥2且FDR<0.01為篩選標準[12]，對不同時間點差異表達基因的分布進行調查，發現0 h VS 12 h(表示兩個侵染時間點差異表達基因的比較，下同)、0 h VS 24 h、0 h VS 48 h分別篩選出7 816、10 412、9 563個上調表達差異基因和7 209、10 159、9 612個下調表達差異基因；12 h VS 24 h、12 h VS 48 h、24 h VS 48 h分別篩選出937、1 150、57個上調表達差異基因和 2 866、 1 643、229個下調表達差異基因。結果表明，以0 h為對照，24 h獲得的差異表達基因最多，48 h次之，12 h最少；以12 h為對照，24 h獲得的差異表達基因多于48 h；48 h較24 h新增加了286個差異表達基因。

為獲得潛在的抗病相關基因，本研究重點關注上調表達差異基因，從圖2可以看出，共獲得13 674個上調表達差異基因，12 h與24 h、48 h共有的上調表達差異基因分別有6 241、5 611個，24 h與48 h共有的上調表達差異基因7 416個，0～48 h篩選到持續上調表達的差異基因5 151個，這些基因可能與中4的的抗條銹性有關。對5 151個上調差異基因進行染色體定位，其中 4 482個獲得定位信息，其余669個無法確定其染色體位置，為新基因。進一步分析發現，差異表達基因在B染色體組上分布最多，為1 581個；D染色體組次之，為1 523個；A染色體組最少，為 1 378個。上調表達差異基因在2B(7.12%)、5B (6.38%)、3B(6.34%)、2D(5.85%)、7D(5.42%)染色體上所占比例較高。

2.4 共同上調表達差異基因的分類注釋及富集

2.4.1 共同上調表達差異基因的GO功能注釋

中4轉錄組測序數據中有58 059個轉錄本在GO數據庫中被注釋，5 151個共同上調表達差異基因中3 411個基因獲得GO注釋。獲得注釋的轉錄本可分為分子功能、細胞組分和生物學過程3大類，46個亞類(圖3)。分子功能包含13個亞類，其中蛋白結合、催化活性通路富集到的上調表達差異基因數目較多，此外，在轉運活性、受體活性、酶調節活性、抗氧化活性、蛋白結合轉錄因子活性等通路也富集到上調表達差異基因；細胞組分包含13個亞類，其中細胞部件、細胞、細胞器等通路富集到的上調表達差異基因數目較多，在大分子復合物、細胞連接、超分子復合物等通路也富集到上調表達差異基因；生物學過程包含20個亞類，其中細胞過程、代謝過程、單有機體過程等通路富集到的上調表達差異基因數目較多，此外，在生殖過程、多有機體過程、信號轉導等通路也富集到上調表達差異基因。

2.4.2 共同上調表達差異基因的功能富集分析

利用COG、KOG、eggNOG數據庫對5 151個上調表達差異基因分別進行直系同源分類，以期獲得更準確的分類及功能注釋。COG數據庫中有2 638個上調表達差異基因獲得注釋分類，最大類別為R類(僅一般功能預測)，包含555個基因，占21.04%；其次是K類(轉錄)，包含306個基因，占11.60%(圖4A)；KOG數據庫中有 2 486個上調表達差異基因獲得注釋分類，最大類別為R類，包含436個基因，占17.54%，其次是O類(翻譯后修飾，蛋白質轉換，分子伴侶)，包含271個基因，占10.90%(圖4B)；eggNOG數據庫中有3 906個上調表達差異基因獲得注釋分類，最大類別為S類(未知功能)，包含859個基因，占21.99%，其次是R類，包含666個基因，占 17.05%(圖4C)。此外，COG、KOG、eggNOG三個數據庫分別注釋到與防御相關的V類(防御機制)上調表達差異基因，分別為47、28、25個。雖然eggNOG數據庫獲得的注釋基因最多，但是注釋到與防御相關的基因最少，在未知功能和一般功能預測富集到的基因最多，為1 525個，占 39.04%；KOG數據庫在所占比例為 5.83%～ 5.39%處出現分類富集，且數據庫25個分類中均注釋到上調表達差異基因；COG和KOG數據庫雖然在最大類別注釋相同，但是其余類別差異很大。三個數據庫注釋分類結果有較大差異，因此，應綜合分析三個數據庫彌補單一數據庫的不足。

2.4.3 共同上調表達差異基因的KEGG通路富集分析

5 151個差異表達基因在KEGG數據庫中注釋到431個上調表達差異基因，主要分為5大類20個小類，共涉及109個代謝通路。其中，細胞過程、環境信息處理、遺傳信息加工、代謝和有機系統5大類注釋到的上調表達差異基因數目分別為39、24、146、203和19。有機系統類別注釋到14個參與植物與病原菌互作通路的差異表達基因。將顯著性設為qvalue<0.05，對上調差異基因進行富集分析，發現差異基因在剪接體、苯丙氨酸代謝、過氧化物酶體、糖和氨基酸的代謝等通路顯著富集 (圖5)。

2.5 R基因和新基因的發掘

基于所選參考基因組序列的比對結果，共注釋到包括7種類型的33個R(resistance)基因，其中，12個與富含半胱氨酸類受體蛋白激酶有關，10個與抗病性蛋白有關(表2)。此外，還篩選到β-半乳糖苷酶、纖維素合酶、過氧化物酶等病程相關蛋白。

表2 中4抗條銹性上調差異基因中R基因的功能分析Table 2 Functional analysis of R gene in up-regulated differential genes of Zhong 4 resistance to stripe rust

將篩選到的669個上調表達的差異新基因在7種數據庫注釋分析，結果顯示，有335個獲得GO注釋，123個獲得COG注釋，79個獲得KOG注釋，511個獲得eggNOG注釋；其中，eggNOG數據庫最大的兩類注釋分別為功能未知和一般功能預測，占48.57%，KOG數據庫中細胞周期控制、細胞分裂和染色體分裂通路中獲得注釋信息的基因比例較高。151個上調表達差異新基因獲得KEGG注釋，且其中85個差異新基因共涉及89個信號通路；562個獲得NR注釋。其中，苯丙素的生物合成和苯丙氨酸代謝通路分別富集到7個上調表達差異基因，谷胱甘肽代謝通路分別富集到6個，內質網蛋白加工和碳代謝通路分別富集到5個，植物與病原菌互作、剪接體、糖酵解和糖異生以及淀粉和蔗糖代謝通路分別富集到4個，說明顯著富集的差異新基因主要與糖的合成與代謝、氨基酸的代謝以及植物與病原菌互作有關。上調表達差異新基因還注釋到R基因、WRKY 和NAC轉錄因子等應答非生物脅迫的基因。

2.6 qRT-PCR對RNA-seq數據的驗證

為驗證轉錄組測序RNA-seq數據的準確度，在差異基因中隨機選取與NB-ARC、鉀離子轉運蛋白和熱激蛋白相關的3個基因進行qRT-PCR，其中，轉錄本數據經NCBI Blast確定其特異序列位置，用Premier 6設計與內參引物TaEF-1α 退火溫度值相近的特異引物。以0 h為對照，計算基因相對表達水平與轉錄組數據表達量并進行對比。結果(表3)表明，轉錄組測序RNA-seq數據與qRT-PCR數據整體趨勢一致，轉錄組測序數據準確可信。

表3 熒光定量PCR結果與轉錄組測序結果對比Table 3 Comparison of real-time PCR results with RNA-seq results

3 討 論

條銹病的流行嚴重威脅小麥的產量和品質，培育抗病品種是防治該病最有效、經濟和環保的措施[1-4]。中4自育成至今，對我國現有條銹菌小種均表現為抗病，然而其抗條銹病機制尚不明確。

本研究利用轉錄組測序，分析了條銹菌CYR32侵染中412、24和48 h后的差異表達基因，獲得 5 151個持續上調表達差異基因。經多個數據庫注釋及富集分析，篩選得到R基因、PR基因、蛋白酶體、WRKY和NAC轉錄因子等與調控小麥條銹菌有關的基因。其中33個上調表達差異基因參與植物與病原菌互作過程；大量上調表達差異新基因獲得注釋，在植物與病原菌互作過程、糖的合成與代謝以及氨基酸的代謝通路顯著富集。

Chen等[13]通過轉錄組學分析小種專化性抗性、全生育期抗性和高溫成株期抗性，發現高溫成株期抗性中基因多樣性最高，但是基因誘導表達緩慢，解釋了其抗性持久性的原因。Coram等[14]分析Yr39介導的與條銹菌親和與不親和的兩個F7重組自交系的轉錄本，發現條銹菌侵染48 h后時，轉錄本數量均達到峰值；對轉錄組數據進行注釋后發現，50%以上的轉錄本參與防御和信號轉導過程，例如R基因同源物、PR蛋白、苯丙氨酸的生物合成和蛋白激酶等。Coram等[15]分析Yr5介導的與條銹菌親和與不親和的相關轉錄本，發現在侵染24 h時轉錄本積累量出現峰值，注釋后發現被誘導的轉錄本均參與基礎防御相關過程，Yr5快速響應的信號通路和防御反應有R基因參與，包括蛋白激酶信號轉導和活性氧導致的過敏性壞死反應。本研究結果發現，條銹菌侵染24 h時，差異表達基因數目出現峰值，表明小麥在應對條銹菌的誘導下，大量與寄主抗性及相關抗性基因轉錄表達來調控免疫反應。Hao等[16]研究發現，在條銹菌侵染24 h時，KEGG顯著富集的通路有苯丙氨酸的生物合成和代謝，該途徑可以產生水楊酸、類黃酮和木質素等化合物，并且苯丙氨酸代謝途徑在條銹菌侵染初期參與植物防御反應；在條銹菌侵染120 h時，KEGG富集到與淀粉、蔗糖和碳水化合物代謝相關的基因，推測小麥光合作用的強度影響其對條銹菌的抗性。本研究KEGG通路分析發現，上調表達差異基因顯著富集在蛋白質加工、苯丙氨酸代謝、過氧化物酶體、糖和氨基酸代謝等通路，與前人研究結果一致。Garnica等[17]研究發現，ABC轉運蛋白基因在孢子萌發中差異表達，并在附著期高度富集；MFS轉運蛋白優先在孢子萌發期表達，主要影響其生長過程；己糖轉運蛋白和氨基酸轉運體在吸器中高度表達，條銹菌的生活方式對轉運體類型起決定作用。葡萄糖轉運體在發芽夏孢子中高表達，說明條銹菌侵染小麥后消耗寄主碳水化合物吸收營養。

PR蛋白在系統獲得抗性和植物防御中尤為重要。研究表明，PR1蛋白能夠抑制小麥條銹菌的生長[18]。另外，TaWRKY49負調控小偃6號對小麥條銹菌的抗性，并誘導水楊酸、茉莉酸和活性氧響應相關基因的表達，同時抑制乙烯相關基因的表達；TaWRKY45和aWRKY62正調控小偃6號對小麥條銹菌的抗性，且水楊酸、乙烯、茉莉酸甲酯、過氧化氫和脫落酸均能夠誘導TaWRKY62基因表達[19-20]。NAC基因家族編碼的蛋白質是目前發現的最大的植物特異性轉錄因子之一[21-23]，在植物應答非生物和生物脅迫以及調控植物發育等方面發揮著重要的作用，TaNAC4[21]、TaNAC21/22[22]和TaNAC30[23]負調控小麥對條銹菌的抗性，且茉莉酸甲酯、乙烯和脫落酸能夠誘導TaNAC4基因的表達，而水楊酸對TaNAC4基因的表達無影響[21]；TaNAC21/22作為tae-miR164的靶標在小麥對條銹菌的抗性方面發揮著重要作用[22]；沉默TaNAC30可顯著增加H2O2的積累量，說明TaNAC30可增加小麥對條銹菌的抗性[23]。

本研究分析了條銹菌CYR32侵染中4的轉錄組數據，通過基因組序列比對獲得7種類型R基因并得到大量的新基因，新基因中也預測到R基因和抗條銹性相關基因，表明小麥品種中4中蘊含豐富的抗性基因。