張艷俊，楊毅清，袁心悅，王海燕，劉大群

（1.河北農(nóng)業(yè)大學(xué) 植物保護(hù)學(xué)院/河北省農(nóng)作物病蟲害生物防治工程技術(shù)研究中心/國家北方山區(qū)農(nóng)業(yè)工程技術(shù)研究中心，河北 保定 071001；2.中國農(nóng)業(yè)科學(xué)院 研究生院，北京 100081）

小麥葉銹病是由小麥葉銹菌（Puccinia tirticina）引起的真菌病害，是世界上分布范圍最廣的病害［1］，其破壞性強(qiáng)，嚴(yán)重威脅到了小麥的生產(chǎn)安全。我國是小麥的生產(chǎn)、消費(fèi)和貿(mào)易大國，其產(chǎn)量和消費(fèi)量位居世界第一［2］，同時(shí)也是世界上小麥葉銹病的流行區(qū)之一。2008-2009年，小麥葉銹病發(fā)生較為嚴(yán)重，2010年有所下降，2011年輕度發(fā)生，2012年河南、陜西、甘肅、四川、安徽等地小麥葉銹病大面積發(fā)生，2013年山東、河南和新疆局部地區(qū)發(fā)生嚴(yán)重，其他地區(qū)輕度發(fā)生［3-4］。2015和2016年在黃淮海地區(qū)爆發(fā)流行，且發(fā)病比往年約提前一個(gè)月［5-7］。目前，生產(chǎn)中主要依靠化學(xué)農(nóng)藥和抗病育種防治小麥葉銹病。化學(xué)防治雖然見效快，但是污染環(huán)境，不符合人、自然、社會(huì)和諧發(fā)展這一客觀規(guī)律的要求，影響了生態(tài)文明建設(shè)，同時(shí)多次使用同種藥物還會(huì)增加葉銹菌的抗藥性。因此，種植抗病品種不但是防治小麥葉銹病最經(jīng)濟(jì)有效的措施，也是符合生態(tài)文明建設(shè)的要求，適合社會(huì)發(fā)展規(guī)律。

隨著測序技術(shù)的發(fā)展，RNA測序(RNA-seq)通過提供基因表達(dá)的量化，為抗病候選基因的鑒定和分子標(biāo)記開發(fā)提供了基礎(chǔ)。RNA-seq是轉(zhuǎn)錄組學(xué)最強(qiáng)大的測序方法，它基于深度測序技術(shù)，可以比其他方法更精確地測量轉(zhuǎn)錄水平及其亞型［8］，不需要基因組序列，即可分析寄主與病原物的轉(zhuǎn)錄本信息，從而使我們可以評(píng)估病原物是如何調(diào)節(jié)致病基因的表達(dá)以及如何影響寄主植物的防御反應(yīng)［9］。Manickavelu等［10］通過對(duì)近等基因系小麥抗葉銹病與感葉銹病的cDNA文庫測序，分析了Lr10相關(guān)的抗病基因表達(dá)情況。對(duì)抗病文庫與感病文庫中的ESTs進(jìn)行對(duì)比，發(fā)現(xiàn)抗感反應(yīng)間基因的表達(dá)存在差異。Kumar等［11］利用深度測序在小麥中鑒定了497個(gè)保守miRNAs和559個(gè)新的miRNAs，同時(shí)利用降解物組分析篩選到701個(gè)靶基因，并利用熒光定量PCR分析了miRNAs在與葉銹菌互作中的差異表達(dá)。張芳芳等［12］以RNA-seq測序數(shù)據(jù)為基礎(chǔ)，篩選和鑒定了小麥LRRs類基因，又利用qRT-PCR技術(shù)對(duì)其進(jìn)行了轉(zhuǎn)錄水平的表達(dá)分析，為深入研究小麥抗葉銹分子機(jī)制提供了新思路。

小麥抗葉銹病基因Lr19來源于長穗偃麥草，定位在7DL染色體上，為全生育期抗性基因，至今在我國很少發(fā)現(xiàn)對(duì)Lr19有毒性的菌株，具有很好的開發(fā)應(yīng)用價(jià)值，是一個(gè)應(yīng)用潛力很大的抗葉銹病基因。本實(shí)驗(yàn)室的張楠、王珅和高琳在小麥‘TcLr19’中克隆了TaPR1、TaPR2和TaTLP1等3個(gè)PR基因，并明確其表達(dá)受葉銹菌的誘導(dǎo)，利用基因槍法和VIGS技術(shù)對(duì)3個(gè)PR基因進(jìn)行了功能驗(yàn)證［13-15］。本研究在對(duì)接菌后小麥材料TcLr19轉(zhuǎn)錄組測序的基礎(chǔ)上，依據(jù)基因差異表達(dá)分析和基因功能注釋結(jié)果篩選與抗葉銹菌相關(guān)的候選基因，采用qRT-PCR技術(shù)對(duì)RNA-seq數(shù)據(jù)進(jìn)行了驗(yàn)證，為進(jìn)一步鑒定和篩選抗病相關(guān)基因提供了可靠信息。

1 材料與方法

1.1 試驗(yàn)材料及處理

植物材料為小麥抗葉銹病近等基因系材料‘TcLr19’,其攜帶有小麥抗葉銹病基因Lr19，通過6代回交獲得。

采用撒粉接種法接種葉銹菌07-10-421-3(FHJT)的夏孢子（孢子萌發(fā)率在80%以上）于一葉一心期的小麥材料‘TcLr19’上，分別剪取接菌0、1和6 d后的葉片，用液氮將其迅速冷凍，儲(chǔ)藏于-80℃冰箱中，備用。

1.2 RNA-seq測序及原始數(shù)據(jù)處理、分析

將準(zhǔn)備好的樣品送上海派森諾生物科技股份有限公司進(jìn)行轉(zhuǎn)錄組測序?；贗llumina NextSeq500測序平臺(tái)，采用第2代測序技術(shù)（Next-Generation Sequencing，NGS）進(jìn)行測序。用Sanger質(zhì)量值Q來評(píng)估原始數(shù)據(jù)的測序質(zhì)量，保證后續(xù)定量分析的準(zhǔn)確性；利用Tophat2（http://tophat.cbcb.umd.edu/）將原始數(shù)據(jù)過濾后的Reads與已發(fā)表的中國春小麥基因組序列草圖（the bread wheat cultivar Chinese Spring，IWGSC）［16］進(jìn)行比對(duì)和定位分析，獲得序列的來源基因以及表達(dá)產(chǎn)物的結(jié)構(gòu)。

1.3 差異基因分析

利用HTSeq 0.6.1p2（http://www-uber.embl.de/users/anders/HTSeq）軟件和RPKM（Reads Per Kilo bases per Million reads）方法，得出各樣品中每一個(gè)基因的表達(dá)量，使用DESeq軟件對(duì)基因表達(dá)進(jìn)行差異分析。設(shè)置篩選兩兩比較組樣本間顯著差異表達(dá)基因條件為：|fold change| ＞ 2，P-value ＜ 0.05。將得到的差異基因進(jìn)行GO功能富集分析和KEGG Pathway富集分析，以P-value≤0.05為閾值。

1.4 熒光定量PCR驗(yàn)證

選取差異表達(dá)基因，利用primer5.0軟件設(shè)計(jì)引物（表1），以測序?qū)嶒?yàn)同批次采集樣品的cDNA為模板，GAPDH（GenBank登錄號(hào)：AF251217）為內(nèi)參基因進(jìn)行qRT-PCR。設(shè)3次生物學(xué)重復(fù)，應(yīng)用LightCycler? 96熒光定量PCR儀進(jìn)行擴(kuò)增，采用2-△△Ct計(jì)算基因的相對(duì)表達(dá)量。

表1 用于qRT-PCR表達(dá)驗(yàn)證的基因及其引物Table 1 The qRT-PCR primers for the selected genes

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測序質(zhì)量檢測和RAW數(shù)據(jù)整理分析

通過對(duì)各樣品的轉(zhuǎn)錄組測序數(shù)據(jù)的分析，分別獲得Raw reads 48 258 292、47 346 190和56 341 268條，數(shù)據(jù)經(jīng)過去除接頭和質(zhì)量過濾后，與原始數(shù)據(jù)相比，有效數(shù)據(jù)量均高于99%。各樣品中GC含量均為58%以上，Q20值均高于98%（表2）。以上分析結(jié)果表明，各樣品的RNA 樣本測序所得原始數(shù)據(jù)的堿基組成基本平衡，且大部分reads 的堿基質(zhì)量值大于20，符合數(shù)據(jù)質(zhì)量錯(cuò)誤率低于1%的要求，說明測序質(zhì)量較好，數(shù)據(jù)具有可靠性。

表2 數(shù)據(jù)統(tǒng)計(jì)Table 2 The statistics of data

2.2 基因比對(duì)分析

通過Bowtie2建立參考基因組索引，然后使用 Tophat2（http://tophat.cbcb.umd.edu/）將過濾后的Reads比對(duì)到參考基因組上。0、1和6 d樣本比對(duì)上參考基因組的序列總數(shù)分別為34 275 708、33 853 145和37 746 006，比對(duì)上參考基因組的序列所占百分比分別為71.66%、72.14%和83.01%（表3）。以上分析結(jié)果表明，試驗(yàn)所產(chǎn)生的測序序列的Mapping比例都高于70%，Mapping成功，轉(zhuǎn)錄組測序參考基因組選擇合適，且相關(guān)實(shí)驗(yàn)不存在污染。

表3 RNA-seq Map統(tǒng)計(jì)Table 3 The statistics of RNA-seq Map

2.3 差異表達(dá)基因的篩選分析

利用DESeq軟件，根據(jù)read count值分別篩選接種小麥葉銹菌后TcLr19在不同時(shí)間的差異表達(dá)基因。為了進(jìn)一步獲得與抗葉銹相關(guān)的抗病基因，對(duì)不同時(shí)間點(diǎn)數(shù)據(jù)的兩兩對(duì)比分析，在TcLr19中共篩選到差異表達(dá)基因8 970個(gè)，其中上調(diào)差異表達(dá)基因4 953個(gè)，下調(diào)差異表達(dá)基因4 017個(gè)。在1 d vs 0 d組合中差異基因最少，在6 d vs 1 d組合中差異基因最多（圖1）。

圖1 不同侵染時(shí)間TcLr19的差異表達(dá)基因Fig.1 DEGs of TcLr19 with different infection time

2.4 差異表達(dá)基因的Gene Ontology富集分析

對(duì)差異基因進(jìn)行GO富集分析，結(jié)果顯示：差異基因被GO富集到分子功能、生物過程和細(xì)胞組分等3大類和1 790個(gè)小類中，其中分子功能部分包括489個(gè)、生物進(jìn)程部分包括1 140個(gè)、細(xì)胞組分部分包括161個(gè)。選取GO富集基因中富集程度高和差異基因占比大的term作圖（圖2）。差異基因主要富集到催化活性和基礎(chǔ)代謝term且富集程度高，6 d vs 1 d的富集term中差異基因數(shù)最多。

圖2 差異表達(dá)基因的GO富集Fig.2 GO enrichment of DEGs

2.5 差異表達(dá)基因的KEGG富集分析

利用KEGG 數(shù)據(jù)庫對(duì)上面已篩選出的差異基因進(jìn)行了功能分類和Pathway 注釋（表4），結(jié)果顯示：3個(gè)時(shí)間對(duì)比組合中的差異基因定位到了84個(gè)代謝通路中，顯著富集到了23個(gè)通路中，如氨基酸的生物合成（Biosynthesis of amino acids）、硫胺素代謝（Glycerolipid metabolism）、半胱氨酸蛋氨酸代謝（Cysteine and methionine metabolism）、甘氨酸、蘇氨酸和絲氨酸代謝（Glycine, serine and threonine metabolism）等。

表4 差異基因KEGG富集分析Table 4 KEGG pathway enrichment of DEGs

2.6 差異基因Blast比對(duì)和在小麥與葉銹菌互作過程中的表達(dá)模式分析

通過以上差異基因的GO和KEGG分析，從轉(zhuǎn)錄組數(shù)據(jù)庫中分別獲得了Traes_5DL_D1D0CA79E、Traes_2AL_3D7807781、Traes_2BS_9931EE821和Traes_2BL_71ED6B5BD的序列。BLAST同源序列比對(duì)得知這些序列分別為NAC類轉(zhuǎn)錄因子、bZIP型轉(zhuǎn)錄因子、絲氨酸/蘇氨酸激酶和酪氨酸蛋白激酶等與抗病相關(guān)的基因。又根據(jù)其序列設(shè)計(jì)qRTPCR引物，驗(yàn)證上面篩選出來的候選基因是否受葉銹菌的誘導(dǎo)及轉(zhuǎn)錄組分析結(jié)果的可靠性，結(jié)果顯示：4個(gè)候選基因在0、1和6 d的表達(dá)量存在差異，且qRT-PCR與RNA-seq差異倍數(shù)分析結(jié)果雖在數(shù)值上不一致，但其差異基因的表達(dá)趨勢(shì)與RNA-seq分析結(jié)果是一致的（圖3），以上結(jié)果分析說明篩選出的差異基因受葉銹菌的誘導(dǎo)及利用RNA-seq技術(shù)分析基因表達(dá)的結(jié)果是可靠的。

圖3 基因IDFig.3 Gene identify

3 討論與結(jié)論

植物受到生物或者非生物的傷害時(shí)，植物體內(nèi)的某些防御反應(yīng)基因的表達(dá)就會(huì)發(fā)生變化，主要體現(xiàn)在時(shí)間、空間及產(chǎn)物含量的變化，從而起到抵抗外界物質(zhì)的侵染及擴(kuò)展的作用。本研究根據(jù)小麥近等基因系TcLr19與葉銹菌生理小種FHJT互作過程中可能會(huì)引起的變化，對(duì)接菌后0、1和6 d的樣品，進(jìn)行了 RNA-Seq測序。這與 Dobon等［17］采用RNA-Seq技術(shù)對(duì)條銹菌侵染小麥過程樣品測序時(shí)間相似，能夠徹底地記錄這些防御調(diào)控基因在表達(dá)上的變化情況。Khalil-Ur-Rehman等［18］在RNA-seq分析的基礎(chǔ)上，比較了夏芽(C1)、半休眠(C2)和內(nèi)休眠(C3)時(shí)的赤霉素(GA)和脫落酸(ABA)相關(guān)轉(zhuǎn)錄本的水平，從而推測ABA和GA通路是調(diào)節(jié)葡萄芽休眠的調(diào)節(jié)開關(guān)。Yang等［19］對(duì)衰老后期的玉米葉片RNA-seq分析，發(fā)現(xiàn)了4 013個(gè)DEGs，且植物體內(nèi)對(duì)氮運(yùn)輸和儲(chǔ)藏起作用的天門冬酰胺合成酶(ZmAS3和ZmAS4)的2個(gè)基因也在干旱處理中被顯著誘導(dǎo)，表明了在干旱誘導(dǎo)葉片衰老過程中，這些參與碳和氮代謝的基因表達(dá)譜發(fā)揮著調(diào)節(jié)作用。本研究通過對(duì)RNA-Seq數(shù)據(jù)的分析，獲得差異基因8 970個(gè)，主要富集到了催化活性和代謝term中，參與到了硫胺素代謝、半胱氨酸蛋氨酸代謝和甘氨酸、蘇氨酸和絲氨酸代謝等代謝通路，為抗葉銹基因的挖掘提供了大量的信息。

為了后續(xù)試驗(yàn)順利進(jìn)行下去，大部分研究利用qRT-PCR進(jìn)一步確定RNA-Seq數(shù)據(jù)的可靠性。Liu等［20］在RNA-Seq分析基礎(chǔ)上，利用qPCR驗(yàn)證套袋蘋果梨4個(gè)成熟階段轉(zhuǎn)錄組比較中的32個(gè)被鑒定為花青素生物合成途徑基因的表達(dá)模式，結(jié)果與RNA-Seq結(jié)果高度吻合。其中，部分花青素生物合成通路基因MYB4-like1、MYB4-like2、MYB1R1、WDR在去袋9 d后表達(dá)上調(diào)，還發(fā)現(xiàn)了1 025個(gè)新基因，為揭示袋裝處理蘋果花青素生物合成機(jī)理提供了豐富的信息。Wang等［21］在全基因組高通量RNA測序的基礎(chǔ)上，篩選出9個(gè)可能參與黃瓜性別分化調(diào)控的候選基因，有5個(gè)基因(Cs-MCM6、Cs-MCM2、Cs-CDC45、Cs-Dpri、Cs-CDC20) 參與細(xì)胞周期通路，推測細(xì)胞周期通路可能在黃瓜性別決定中發(fā)揮重要作用。Guan等［22］利用RNA-seq技術(shù)對(duì)冷處理前后的雜草稻(WR 03-35和WR 03-26)和栽培稻(Kongyu 131和9311)進(jìn)行了測定，并對(duì)隨機(jī)選取的12個(gè)DEGs進(jìn)行半定量RT-PCR和qRT-PCR檢測，結(jié)果發(fā)現(xiàn)這些基因半定量RT-PCR和qRT-PCR的表達(dá)模式與RNA-seq表達(dá)一致，進(jìn)一步確定了RNA-seq數(shù)據(jù)的可靠性。本研究從RNA-seq數(shù)據(jù)庫中篩選到了4個(gè)與抗病相關(guān)的候選基因（NAC類轉(zhuǎn)錄因子、bZIP型轉(zhuǎn)錄因子、絲氨酸/蘇氨酸激酶和酪氨酸蛋白激酶），又利用qRTPCR驗(yàn)證了候選基因的表達(dá)模式，結(jié)果其受葉銹菌誘導(dǎo)且表達(dá)趨勢(shì)與RNA-Seq結(jié)果一致，以上研究結(jié)果為獲得小麥抗葉銹基因及其研究其抗病機(jī)制提供了豐富的信息。

通過RNA- seq技術(shù)對(duì)組織或脅迫誘導(dǎo)的差異表達(dá)基因分析，可以研究基因功能、代謝及調(diào)控途徑等［23］。李芳等［24］在 RNA-seq基礎(chǔ)上，通過BLAST同源序列比對(duì)和RACE技術(shù)克隆得到TaGSO1基因，并利用qPCR技術(shù)分析了TaGSO1基因在小麥與葉銹菌互作過程不同組合中的表達(dá)情況，推測TaGSO1基因參與了小麥抵抗葉銹菌侵染過程,為全面闡釋小麥抗葉銹病機(jī)制奠定了基礎(chǔ)。本研究下一步將結(jié)合RNA- seq數(shù)據(jù)的分析結(jié)果，篩選出與抗葉銹菌相關(guān)的候選基因，明確其與小麥葉銹病抗性的相關(guān)性，探索候選基因介導(dǎo)的小麥與葉銹菌互作的分子機(jī)制。