轉(zhuǎn)GhB301基因棉花響應(yīng)枯萎病菌侵染的轉(zhuǎn)錄組分析

劉戈輝 韓澤剛 孫士超 張 薇,*

(1 石河子大學(xué)農(nóng)學(xué)院，新疆 石河子 832000；2 南京農(nóng)業(yè)大學(xué)作物遺傳與種質(zhì)創(chuàng)新國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210095)

棉花枯萎病(Fusariumoxysporumf.sp.Vasinfectum)作為棉花的土傳病害之一，其病菌可以從根部侵染棉花，毒害棉花維管束，造成棉花葉片青枯、落葉，甚至死亡。其病菌孢子在土壤中的休眠期長(zhǎng)達(dá)10年之久，農(nóng)藥防治對(duì)病情控制的作用有限，并且會(huì)對(duì)土壤水源產(chǎn)生污染，因此培育棉花抗病品種是解決棉花枯萎病最為經(jīng)濟(jì)有效的措施。近年來(lái)，較多的研究集中于棉花黃萎病，且主要集中在抗病種質(zhì)資源的篩選、抗性基因的定位以及基因的遺傳分析等方面[1-2]，而對(duì)棉花枯萎病的研究相對(duì)較少。

隨著高通量測(cè)序技術(shù)的迅猛發(fā)展，其具有時(shí)間短、通量大、靈敏度高等優(yōu)點(diǎn)，在研究差異表達(dá)基因(Differentially expressed genes,DEGs)方面顯示出強(qiáng)大的優(yōu)勢(shì)[3-4]。第二代高通量轉(zhuǎn)錄組測(cè)序(RNA-seq)已被廣泛應(yīng)用于植物基因表達(dá)分析研究。韓澤剛等[5]通過(guò)Solexa高通量測(cè)序技術(shù)構(gòu)建棉花抗、感枯萎病品種的基因表達(dá)譜，發(fā)現(xiàn)1 080個(gè)DEGs在棉花抗病和感病品種中差異表達(dá)，并發(fā)現(xiàn)在植物激素信號(hào)轉(zhuǎn)導(dǎo)通路中，26個(gè)已知功能的基因可能參與植物抗病反應(yīng)。王琛等[6]通過(guò)分析比較接種和未接種枯萎病菌的棉花葉片轉(zhuǎn)錄組測(cè)序結(jié)果，發(fā)現(xiàn)枯萎病菌侵染后，5 062個(gè)基因的表達(dá)量發(fā)生顯著變化，其中2 091個(gè)基因的轉(zhuǎn)錄水平升高，2 971個(gè)基因的轉(zhuǎn)錄水平下降。杜榮光[7]通過(guò)對(duì)抗感枯萎病海島棉材料接菌前后進(jìn)行轉(zhuǎn)錄組測(cè)序分析，發(fā)現(xiàn)眾多已知抗病相關(guān)基因均顯著高表達(dá)。而這些研究均未對(duì)比分析轉(zhuǎn)基因材料和對(duì)照材料。

本試驗(yàn)在前期研究中獲得了1個(gè)ERF轉(zhuǎn)錄因子家族基因GhB301，在受到枯萎病菌誘導(dǎo)后，該基因的相對(duì)表達(dá)量明顯增加，并且在抗病品種中的表達(dá)量顯著高于感病品種；在煙草中過(guò)表達(dá)該基因，可以顯著增加防御保護(hù)酶活性以及抗病相關(guān)基因的表達(dá)，暗示該基因可能在棉花對(duì)抗枯萎病菌脅迫響應(yīng)中發(fā)揮重要作用[8]。在此基礎(chǔ)上，本研究將GhB301轉(zhuǎn)入陸地棉中并獲得純合株系，研究轉(zhuǎn)基因棉花株系與野生型對(duì)照的抗病性變化，包括分析枯萎病菌誘導(dǎo)前后轉(zhuǎn)基因系和對(duì)照系中的差異表達(dá)基因，篩選與棉花抗枯萎病菌相關(guān)的候選基因，在轉(zhuǎn)錄水平上分析GhB301基因在棉花抗枯萎病中的功能。

1 材料與方法

1.1 棉花枯萎病菌和植物材料的制備

棉花枯萎病菌7號(hào)生理小種F430菌株由石河子大學(xué)棉花分子育種實(shí)驗(yàn)室保存。棉花枯萎病菌在固體馬鈴薯葡萄糖瓊脂(potato dextrose ager,PDA)培養(yǎng)基上培養(yǎng)5 d后(28℃)，轉(zhuǎn)移至查氏液體培養(yǎng)基中，于200 r·min-1的孵育振蕩器(哈爾濱東安電子)中，在28℃黑暗培養(yǎng)7～10 d。最后用無(wú)菌水將孢子懸浮液的濃度調(diào)整至1×107個(gè)·mL-1，備用。

以轉(zhuǎn)GhB301基因棉花株系(N)與野生型受體YZ-1(WT)為材料，經(jīng)硫酸脫絨的種子用10%的雙氧水浸泡2～3 h，無(wú)菌水沖洗5～6遍后，種植于無(wú)菌的營(yíng)養(yǎng)缽(花土∶蛭石=2∶1)中。待棉苗長(zhǎng)至第1片真葉完全展開時(shí)，進(jìn)行傷根接菌處理。

從WT和N中分別采集枯萎病菌誘導(dǎo)后0、6、12、24和48 h的棉花根部組織，液氮速凍后，保存于-80℃冰箱，用于轉(zhuǎn)錄組測(cè)序分析。

1.2 表型鑒定與病情指數(shù)分析

接菌處理20 d后，對(duì)N與WT進(jìn)行抗病性鑒定(3次重復(fù)的植株分別為22、50和38株)。根據(jù)子葉和真葉的癥狀將棉花幼苗劃分為5個(gè)等級(jí)(0、1、2、3、4級(jí))，計(jì)算棉花植株的病情指數(shù)(disease index,DI)：DI=[(∑病級(jí)數(shù)×各級(jí)病株數(shù))/(總株數(shù)×4)]×100[9-10]。利用SPSS和Excel 2010軟件對(duì)數(shù)據(jù)進(jìn)行處理分析。

1.3 RNA提取和轉(zhuǎn)錄組測(cè)序

參照RNA prep Pure Plant Kit(北京天根生化科技有限公司)說(shuō)明書提取各處理樣品總RNA。通過(guò)1%瓊脂糖凝膠電泳，分析RNA的完整性及是否存在DNA污染，使用Qubit2.0 Fluorometer(賽默飛，美國(guó))對(duì)RNA濃度精確定量，使用Agilent 2100 bioanalyzer (Nanodrop,美國(guó))精確檢測(cè)RNA的完整性。

使用Illumina的NEBNext?UltraTM RNA Library Prep Kit(NEB，美國(guó))試劑盒構(gòu)建cDNA文庫(kù)。使用Qubit2.0 Fluorometer(賽默飛，美國(guó))對(duì)文庫(kù)進(jìn)行初步定量，稀釋文庫(kù)至1.5 ng·μL-1。使用Agilent 2100 bioanalyzer對(duì)文庫(kù)的insert size進(jìn)行檢測(cè)。庫(kù)檢合格后，委托北京諾禾致源生物信息科技有限公司利用Illumina HiSeq 4000進(jìn)行測(cè)序。

1.4 數(shù)據(jù)指控和參考基因組比對(duì)

對(duì)原始數(shù)據(jù)進(jìn)行過(guò)濾，去除帶接頭、無(wú)法確定堿基信息以及低質(zhì)量reads，得到clean reads。對(duì)clean reads分別進(jìn)行Q20、Q30和GC含量計(jì)算。從CottonGene數(shù)據(jù)庫(kù)(http://www.cottongen.org)下載陸地棉TM-1參考基因組[11]和基因模型注釋文件，利用軟件Tophat2[12]將高質(zhì)量的clean reads比對(duì)至參考基因組。

1.5 基因表達(dá)量與注釋

使用FPKM(expected number of fragments per Kb of transcript sequence per millions base pairs sequenced)來(lái)表示RNA-seq的基因表達(dá)值[13]。使用edgeR軟件包[14](3.18.1)進(jìn)行DEGs分析。當(dāng)FPKM>1時(shí)認(rèn)為該基因表達(dá)。以校正后的P值≤0.05以及|log2foldchange|≥2作為顯著差異表達(dá)的閾值。利用cluster Profiler R 軟件(v3.10.1)進(jìn)行DEGs的GO富集分析[15]和 KEGG 通路富集分析[16]。

1.6 通過(guò)qRT-PCR驗(yàn)證RNA-seq

隨機(jī)選取15個(gè)基因進(jìn)行實(shí)時(shí)熒光定量PCR(quantitative real-time PCR，qRT-PCR)驗(yàn)證，反應(yīng)體系及程序參考北京全式金生物技術(shù)有限公司Trans start Tip Green qPCR Super Mix說(shuō)明書進(jìn)行，試驗(yàn)設(shè)3個(gè)重復(fù)。采用2-ΔΔCt法[17]計(jì)算基因相對(duì)表達(dá)量，所用引物見表1(用Primer5軟件設(shè)計(jì)引物，由北京華大基因公司合成)。

2 結(jié)果與分析

2.1 表型鑒定及病情指數(shù)調(diào)查

N和WT接種枯萎病菌15 d后，WT株系開始發(fā)病，N株系無(wú)明顯病癥；接菌第20天后WT和N植株均發(fā)病，對(duì)發(fā)病情況進(jìn)行調(diào)查(圖1-A)。結(jié)果表明，N株系中病級(jí)多分布在1級(jí)或2級(jí)，而WT株系中病級(jí)則多分布在3級(jí)及4級(jí)，且明顯多于N株系，說(shuō)明WT株系比N株系感染枯萎病菌的情況更為嚴(yán)重。對(duì)3次病級(jí)統(tǒng)計(jì)進(jìn)行病情指數(shù)計(jì)算(圖1-B)，發(fā)現(xiàn)WT棉花接種枯萎病菌后其病情指數(shù)平均值為37.50。而N株系接種枯萎病菌后其病情指數(shù)平均值僅為14.77。這說(shuō)明GhB301基因的轉(zhuǎn)入可以顯著提高棉花對(duì)枯萎病菌的抗性。

2.2 RNA-seq分析

轉(zhuǎn)錄組分析結(jié)果表明，共獲得279 114 320個(gè)raw reads，過(guò)濾后保留273 111 170個(gè)clean reads，數(shù)據(jù)量為81.95 Gb，每個(gè)文庫(kù)平均8.20 Gb，如表2所示，WT棉花測(cè)序共得到131 456 222個(gè)raw reads，對(duì)低質(zhì)量reads進(jìn)行過(guò)濾后仍有128 669 540個(gè)clean reads，其中clean reads數(shù)最多的是N6，最少的是N48。而N材料測(cè)序后共得到147 658 098個(gè)raw reads，過(guò)濾后剩余144 441 630個(gè)clean reads，其中clean reads數(shù)最多的是WT24，最少的是WT48。從測(cè)序結(jié)果可以看出，各樣本文庫(kù)的Q30值均超過(guò)87.64%，GC含量均大于43.06%，其平均值分別為89.86%和43.28%。將clean reads比對(duì)到參考基因組中，發(fā)現(xiàn)比對(duì)率的均值是94.37%，其中最小的是N12，為92.94%，最大的是WT24，為94.95%，85.26%～87.56%的reads能夠比對(duì)到參考基因組的唯一位置，但是仍有7.27%～8.10%的reads比對(duì)到參考基因組的多個(gè)位置。綜合分析10個(gè)樣本材料可以看出，平均有27 311 117個(gè)clean reads 以94.37%的平均比對(duì)率比對(duì)至參考基因組上，并且86.78%的clean reads可以比對(duì)至參考基因組的唯一位置，為后續(xù)分析奠定了堅(jiān)實(shí)的數(shù)據(jù)基礎(chǔ)。

注：A：3次重復(fù)的病級(jí)統(tǒng)計(jì)；B：病情指數(shù)分析；**表示與WT相比，N的病情指數(shù)在0.01水平上差異顯著。Note:A:Three replicated disease statistics.B:Disease index statistics.** indicate that N′s disease index is significantly different compared with WT at 0.01 level.圖1 棉花抗病性檢測(cè)Fig.1 Detection of disease resistance in Gossypium hirsutum

2.3 DEGs篩選與注釋

為了研究轉(zhuǎn)基因棉花和野生型棉花接種枯萎病菌后不同時(shí)期的轉(zhuǎn)錄組情況，利用EdgeR進(jìn)行差異表達(dá)基因的篩選，以P≤0.05，|log2foldchange|≥2為篩選條件，共發(fā)現(xiàn)14 021個(gè)DEGs(其中1 290個(gè)新基因)能夠響應(yīng)枯萎病菌誘導(dǎo)。對(duì)這些DEGs進(jìn)行GO富集分析(圖2-A)，共富集注釋到30個(gè)亞類，涉及生物過(guò)程，細(xì)胞組成和分子功能三大類，主要富集在抗氧化活性、氧化還原酶活性、過(guò)氧化物酶活性、氧化應(yīng)激的反應(yīng)、光系統(tǒng)Ⅱ、碳水化合物結(jié)合以及細(xì)胞壁等功能類別，其中富集DEGs最多的是氧化還原酶活性(GO:0016684)，共注釋到135個(gè)基因，最少的是光系統(tǒng)Ⅱ，僅注釋到24個(gè)基因。成對(duì)比較結(jié)果顯示(圖2-B、C、D)：野生型樣本W(wǎng)T6、WT12、WT24、WT48 與對(duì)照WT0相比，DEGs的數(shù)目分別為3 494、4 807、3 137、3 814個(gè)，WT12 vs WT0的DEGs最多，WT24 vs WT0的DEGs最少；而轉(zhuǎn)基因樣本N6、N12、N24、N48與對(duì)照N0相比，DEGs的數(shù)目分別為6 019、5 220、4 859、2 620個(gè)，N6 vs N0的DEGs最多，N48 vs N0的DEGs最少。比較發(fā)現(xiàn)，轉(zhuǎn)基因株系在病原菌侵染后6 h內(nèi)大量的DEGs即開始做出響應(yīng)，而野生型株系則需要12 h。由此推測(cè)，轉(zhuǎn)基因棉花在應(yīng)對(duì)枯萎病菌脅迫時(shí)響應(yīng)時(shí)更迅速，差異表達(dá)基因更多。

組間比較結(jié)果顯示，N0 vs WT0的DEGs有790個(gè)，N6 vs WT6的DEGs最多，共有935個(gè)，N12 vs WT12存在DEGs 610個(gè)，N24 vs WT24的DEGs最少，僅有607個(gè)，N48 vs WT48存在918個(gè) DEGs。綜上所述，WT組內(nèi)比較發(fā)現(xiàn)共存在694個(gè)共同的DEGs(圖2-B)，轉(zhuǎn)N組內(nèi)比較則有829個(gè)共同的DEGs(圖2-C)。組間比較結(jié)果顯示，存在41個(gè)共同的DEGs(圖2-D)，推測(cè)這些共同的DEGs可能與棉花枯萎病脅迫響應(yīng)具有密切關(guān)系。其中，值得特別關(guān)注的是，Gh_A03G0673(MYB1R1)和Gh_A06G1626(抗病蛋白R(shí)PS5)在N植株中均不表達(dá)，而在WT植株中表達(dá)，推測(cè)轉(zhuǎn)入GhB301基因可能促進(jìn)ERF轉(zhuǎn)錄因子的表達(dá)而抑制MYB轉(zhuǎn)錄因子和PR基因的表達(dá)。Gh_A02G1734(戊二酸/鐵依賴性雙加氧酶)基因在N植株中的表達(dá)量均高于WT，推測(cè)轉(zhuǎn)入GhB301基因可以提高氧化還原酶相關(guān)基因的表達(dá)。Gh_A13G2026(生長(zhǎng)素應(yīng)答蛋白SAUR32)只在N株系中表達(dá)，推測(cè)轉(zhuǎn)入GhB301基因合影響棉花植株的生長(zhǎng)發(fā)育。

注：A：總差異表達(dá)基因GO富集分析；其中1氧化應(yīng)激反應(yīng)；2：氨基多糖代謝過(guò)程；3：氨基多糖分解代謝過(guò)程；4：幾丁質(zhì)代謝過(guò)程；5：幾丁質(zhì)分解代謝過(guò)程；6：氨基糖代謝過(guò)程；7：細(xì)胞壁大分子分解代謝過(guò)程；8：細(xì)胞壁大分子代謝過(guò)程；9：氨基糖分解代謝過(guò)程；10：含谷氨酸的復(fù)合代謝過(guò)程；11：細(xì)胞外區(qū)；12：細(xì)胞外圍；13：質(zhì)外體；14：細(xì)胞壁；15：外部封裝結(jié)構(gòu)；16：外囊腫；17：細(xì)胞皮層；18：細(xì)胞皮層部分；19：細(xì)胞質(zhì)區(qū)；20：光系統(tǒng)II；21：氧化還原酶活性，作用于過(guò)氧化物；22：過(guò)氧化物酶活性；23：抗氧化活性；24：銅離子結(jié)合；25：氧化還原酶活性作用于單個(gè)供體；26：氧化還原酶活性作用于供體；27：雙加氧酶活性；28：幾丁質(zhì)酶活性；29：十葡聚糖：木糖基轉(zhuǎn)移酶活性；30：轉(zhuǎn)移酶活性；B：不同時(shí)期的WT樣本比較；C：不同時(shí)期的N樣本比較；D：相同時(shí)期的兩樣本間的比較。Note:A:The GO enrichment analysis of the total differentially expressed genes;1 RESPONSE to oxidative stress.2:Aminoglycan metabolic process.3:Aminoglycan catabolic process.4:Chitin metabolic process.5:Chitin catabolic process.6:Amino sugar metabolic process.7:Cell wall macromolecule catabolic process.8:Cell wall macromolecule metabolic process.9:Amino sugar catabolic process.10 Glucosamine-containing compound metabolic process.11:Extracelular region.12:Cell periphery.13:Apoplast.14:Cell wall.15:External encapsulating structure.16:Excyst.17:Cell cortex.18:cell cortex part.19:Cytoplasmic region.20:Photosystem II.21:Oxidoreductase activity.acting on peroxide.22:Peroxidase activity.23:Antioxidant activity.24:Copper ion binding.25:Oxidoreductase activity acting on single donors.26:Oxidoreductase activity acting on donors.27:Dioxygenase activity.28:Chitinase activity.29:Xyloglucan:xyloglucosyl transferase activity.30:Transferase activity.B:Represents the comparisons of WT samples at the different stages.C:Represents the comparisons of N samples at the different stages.D:Represents the comparisons of two lines at the different stages.圖2 差異表達(dá)基因GO富集分析及不同樣本間的基因數(shù)目比較Fig.2 The GO enrichment analysis of the total differentially expressed genes and the number of genes between different samples

2.4 氧化還原過(guò)程中相關(guān)DEGs

活性氧(reactive oxygen species，ROS)作為植物活性氧產(chǎn)生與清除的氧化還原產(chǎn)物，能夠參與植物的生長(zhǎng)發(fā)育及應(yīng)激反應(yīng)[18]。對(duì)DEGs進(jìn)行GO富集分析，發(fā)現(xiàn)共有135個(gè)DEGs參與氧化還原過(guò)程(GO:0016684)(圖3)。

圖3 與氧化還原過(guò)程相關(guān)基因的熱圖Fig.3 The heat map of genes related to oxidation-reduction process

呼吸爆發(fā)氧化酶同系物(respiratory burst oxidase homologs，RBOHs)，又被稱作植物NADPH氧化酶，在細(xì)胞生長(zhǎng)、植物發(fā)育、植物對(duì)非生物環(huán)境約束以及與生物相互作用的反應(yīng)中作為ROS生產(chǎn)者發(fā)揮作用，是重要的分子“中樞”[19]。本研究共鑒定出9個(gè)RBOHs參與氧化還原過(guò)程，即1個(gè)RBOHA(Gh_D08G1257)，2個(gè)RBOHB(Gh_D11G2743和Gh_A11G2426)，2個(gè)RBOHC(Gh_D07G0463和Gh_A07G0398)，3個(gè)RBOHD(Gh_A05G2211、Gh_D01G0990和Gh_D05G2471)，1個(gè)RBOHF(Gh_A12G2653)。其中RBOHB在2種材料中均呈現(xiàn)先上升后下降的趨勢(shì)，在N株系中，其表達(dá)量在6 h即達(dá)到峰值，WT中則在12 h達(dá)到峰值。而RBOHC、RBOHD、RBOHF在N植株中呈現(xiàn)先下降后上升的趨勢(shì)，在WT中則呈現(xiàn)持續(xù)下降的趨勢(shì)。

對(duì)ROS清除劑抗壞血酸過(guò)氧化物酶(ascorbic acid peroxidase，APX)，過(guò)氧化氫酶(catalase，CAT)，過(guò)氧化物酶(peroxidase，POD)，谷胱甘肽過(guò)氧化物酶(glutathione peroxidase，GPX)等進(jìn)行表達(dá)模式分析，共檢測(cè)到5種APX蛋白質(zhì)基因(Gh_D13G2402、Gh_A06G0383、Gh_A08G1745、Gh_A02G1648、Gh_A01G1388)，這些基因在N株系中表達(dá)量均明顯高于WT。共鑒定到2種CAT基因(Gh_D03G0021和Gh_A05G1539)，在2種材料中均表現(xiàn)出先上調(diào)后下調(diào)的趨勢(shì)；2種POD基因(Gh_A08G0714和Gh_D08G0832)，在N株系中上調(diào)表達(dá)，在WT中則先上調(diào)后下調(diào)；鑒定出3種GPX基因(Gh_A09G2408、Gh_D12G2260和Gh_A08G0673)，其在2種材料中的表達(dá)無(wú)明顯變化。另外，鑒定得到2個(gè)1-Cys(1-Cys peroxiredoxin)基因(Gh_D10G1825和Gh_A10G1567)和2個(gè)DOX(alpha-dioxygenase)基因(Gh_A05G3137和Gh_D09G1660)，在N株系中的表達(dá)量均低于WT。

綜上，在ROS產(chǎn)生與清除過(guò)程中，眾多DEGs參與其中，其表達(dá)量復(fù)雜多樣的動(dòng)態(tài)變化也反映出棉花在抵抗枯萎病感染的動(dòng)態(tài)平衡中，氧化還原過(guò)程相關(guān)基因起著重要作用。

2.5 植物防御反應(yīng)相關(guān)DEGs

在植物與病原物長(zhǎng)期演化過(guò)程中，植物產(chǎn)生了一系列防御反應(yīng)機(jī)制來(lái)阻止病原物的侵害[20]。對(duì)所有DEGs進(jìn)行GO富集，發(fā)現(xiàn)共有67個(gè)與防御反應(yīng)(GO:0006952)相關(guān)的DEGs(圖4)，共注釋到33個(gè)PR10基因。這些PR10基因中存在20個(gè)MLP(major latex protein)基因、7個(gè)BETV1(major pollen allergen bet v 1)基因、4個(gè)NCS2(s-norcoclaurine synthase 2)基因、2個(gè)STH-2(pathogenesis-related protein STH-2)基因。Gh_D02G1678、Gh_A03G1240和Gh_D04G1394均注釋為MLP423基因，前2個(gè)基因在2種材料中均呈現(xiàn)先上調(diào)后下調(diào)的趨勢(shì)，而Gh_D04G1394在N株系和WT中的表達(dá)模式不同，其在N株系先上調(diào)后下調(diào)，而在WT中則先下調(diào)后上調(diào)。

圖4 與植物防御反應(yīng)相關(guān)基因的熱圖Fig.4 The heat map of genes related to plant defense response

MLO(mildew resistance locus O)蛋白的表達(dá)受到生物與非生物脅迫的影響[21]。在這些與防御反應(yīng)相關(guān)的DEGs中共存在9個(gè)MLO基因，整體來(lái)說(shuō)，在N株系中的表達(dá)量均低于在WT中的表達(dá)量。其中Gh_D07G1738和Gh_A07G1067在N和WT中皆下調(diào)表達(dá)，而Gh_D02G0670、Gh_A10G1852和Gh_D05G0753在2種材料中均顯示先上調(diào)后下調(diào)的趨勢(shì)。

另外，還發(fā)現(xiàn)1個(gè)編碼HEV1(pro-hevein)蛋白的基因(Gh_D11G2595)，2個(gè)編碼PR-4B(pathogenesis-related protein)基因(Gh_D13G1816,Gh_D13G1817)，以及4個(gè)編碼NIMIN-1(protein NIM1-INTERACTING 1)基因(Gh_D11G1125、Gh_D12G1232、Gh_A11G0975、Gh_A12G1108)，6個(gè)HSPRO2(nematode resistance protein-like)基因(Gh_A05G2623、Gh_D05G2913、Gh_D13G0174、Gh_A10G2093、Gh_D10G2363、Gh_A13G0153)。HEV1基因在2個(gè)材料中均呈現(xiàn)下調(diào)趨勢(shì)，且在N株系中的表達(dá)量均低于在WT中的表達(dá)量。PR-4B基因在2個(gè)材料中均呈現(xiàn)逐漸升高的趨勢(shì)，特別是在48 h，在N株系中的表達(dá)量遠(yuǎn)高于WT。4個(gè)NIMIN-1基因則在2個(gè)材料中表現(xiàn)出相反的表達(dá)趨勢(shì)，在N株系中表現(xiàn)為先降低后升高的表達(dá)趨勢(shì)，在WT中則相反。6個(gè)HSPRO2基因均大致表現(xiàn)出先降低后升高的趨勢(shì)，在N中大部分的基因均在24 h降到最低點(diǎn)，而在WT中則在6 h降到最低點(diǎn)。在48 h，N中所有基因的表達(dá)量均高于WT中對(duì)應(yīng)基因的表達(dá)量。

2.6 苯丙烷代謝途徑相關(guān)DEGs

本研究還鑒定到577個(gè)DEGs富集在KEGG途徑，主要涉及植物激素信號(hào)轉(zhuǎn)導(dǎo)、植物病原體相互作用、苯丙烷生物合成、MAPK信號(hào)傳導(dǎo)途徑-植物等代謝途徑(表3)。

表3 差異基因KEGG通路富集統(tǒng)計(jì)分析Table 3 Analysis of KEGG function annotation of differential expression genes

研究發(fā)現(xiàn)苯丙烷代謝途徑以及與這些途徑相關(guān)的基因在棉花與病原菌的互作中發(fā)揮著重要作用[22]，因此對(duì)苯丙烷類生物合成途徑相關(guān)基因進(jìn)行重點(diǎn)分析。結(jié)果表明，31個(gè)DEGs參與該代謝途徑(圖5)，其中包含15個(gè)PER(peroxidase)基因。在PER基因中，Gh_A05G1452的表達(dá)量最高，在N和WT株系中均表現(xiàn)出先降低后升高的趨勢(shì)，但在N中的表達(dá)量顯著高于WT。此外，編碼PER的基因Gh_D04G1101在2個(gè)材料中均表現(xiàn)出先升高后降低的趨勢(shì)，并且在N中的表達(dá)量是在WT的2倍。編碼PER的基因Gh_A10G2288在N中6 h的表達(dá)量是WT的2倍，但在其他時(shí)間點(diǎn)均顯著低于WT。

圖5 與苯丙烷代謝途徑相關(guān)基因的熱圖分析Fig.5 The heat map analysis of genes related with phenylpropane metabolic pathway

鑒定出2個(gè)β-葡萄糖苷酶(beta-glucosidase，BGLU)基因(Gh_A06G2098、Gh_D13G0844)，其中Gh_A06G2098基因的WT中，枯萎病菌侵染初期該基因基本不表達(dá)，直至48 h才逐漸發(fā)揮作用，然而在N株系中，該基因侵染初期即開始表達(dá)。Gh_D13G0844基因在N中均表現(xiàn)出先升高后降低的趨勢(shì)，在WT中則表現(xiàn)出逐漸下降趨勢(shì)；且該基因在N中的表達(dá)量高于在WT中的表達(dá)量。

在苯丙烷代謝通路中，還發(fā)現(xiàn)3個(gè)CAD(cinnamyl alcohol dehydrogenase-like protein)基因(Gh_D12G2787、Gh_D11G1980和Gh_A02G1320)在N中的表達(dá)量顯著高于WT。3個(gè)4CL(4-coumarate-CoA ligase)基因(Gh_A08G1375、Gh_A02G1344、Gh_D09G1372)，其中Gh_A08G1375在N中的表達(dá)更早且表達(dá)量更高。1個(gè)P450(CytochromeP450)基因(Gh_D01G1930)，2個(gè)PNC(Cationic peroxidase)基因(Gh_Sca007747G01、Gh_D02G1998)。6 h時(shí)P450(Gh_D01G1930)基因在N中的表達(dá)量是WT的3倍。0 h時(shí)，PNC基因Gh_D02G1998在N中的表達(dá)量高于在WT中的表達(dá)量，然而在其他時(shí)間點(diǎn)未顯著變化。

綜上所述，大量DEGs參與苯丙烷代謝途徑，其差異表達(dá)基因在N株系中的響應(yīng)時(shí)間早于WT，并且表達(dá)量也更高，推測(cè)接種枯萎病菌后，N株系苯丙烷代謝途徑激活響應(yīng)的速度更快，這可能與其抗病性顯著增強(qiáng)有密切關(guān)系。

注：柱狀圖表示qRT-PCR的結(jié)果，折線圖表示轉(zhuǎn)錄組測(cè)序的結(jié)果。Note:The histogram shows the results of qRT-PCR,and the line graph shows the results of transcriptome sequencing.圖6 通過(guò)qRT-PCR驗(yàn)證RNA-seq數(shù)據(jù)Fig.6 Validation of RNA-seq data by qRT-PCR

2.7 qRT-PCR檢驗(yàn)

為了進(jìn)一步確定轉(zhuǎn)錄組數(shù)據(jù)的可靠性，從差異表達(dá)基因中隨機(jī)抽取15個(gè)基因進(jìn)行qRT-PCR檢驗(yàn)。結(jié)果表明(圖6)，轉(zhuǎn)錄組得到的FPKM值與qRT-PCR得到的基因相對(duì)表達(dá)量在0～48 h的變化趨勢(shì)基本一致，說(shuō)明RNA-seq數(shù)據(jù)分析具有可靠性。

3 討論

研究表明ERF轉(zhuǎn)錄因子與植物抗病響應(yīng)密切相關(guān)[23]。ERF1作為ET/JA響應(yīng)應(yīng)答調(diào)控因子起作用，過(guò)表達(dá)ERF1的擬南芥對(duì)枯萎病菌抗性增強(qiáng)，ERF2也具有類似的功能[24]。在擬南芥中過(guò)表達(dá)B3亞組的AtERF14基因能夠增強(qiáng)防衛(wèi)基因的表達(dá)，并且ERF14突變體對(duì)枯萎病菌表現(xiàn)為更易感病[25]。本研究通過(guò)抗病性鑒定，發(fā)現(xiàn)接種枯萎病菌后過(guò)表達(dá)GhB301轉(zhuǎn)基因棉花與WT相比病情指數(shù)顯著降低，表明轉(zhuǎn)GhB301基因可以增強(qiáng)棉花的抗病性。

Fradin等[26]研究發(fā)現(xiàn)真菌孢子的萌發(fā)以及向表皮細(xì)胞的滲透一般發(fā)生在接菌后的0～12 h。而植物的應(yīng)激反應(yīng)、信號(hào)轉(zhuǎn)導(dǎo)和信息交換一般發(fā)生在真菌感染早期的40 h左右[18,27]。本研究通過(guò)比較差異表達(dá)基因，發(fā)現(xiàn)接菌后6 h時(shí)，轉(zhuǎn)基因棉花中大量的差異表達(dá)基因開始響應(yīng)，而WT棉花則在12 h開始大量表達(dá)差異表達(dá)基因，說(shuō)明轉(zhuǎn)GhB301棉花能夠更快響應(yīng)枯萎病菌的脅迫，這些差異表達(dá)基因的提前響應(yīng)可能會(huì)增強(qiáng)棉花的抗病性。

在這些差異表達(dá)基因中，本研究重點(diǎn)關(guān)注了參與氧化還原過(guò)程、植物防御反應(yīng)以及苯丙烷代謝途徑的基因，其中氧化還原過(guò)程主要是植物產(chǎn)生和清除ROS的過(guò)程[19]。本研究鑒定了9個(gè)呼吸爆發(fā)氧化酶同系物(RBOHs)，并且發(fā)現(xiàn)在N株系中的表達(dá)早于在WT中的表達(dá)。當(dāng)植物受到病原菌等非生物脅迫時(shí)，植物在短時(shí)間內(nèi)產(chǎn)生大量的ROS，如果無(wú)法及時(shí)清除ROS，則會(huì)傷害植株。本研究鑒定出多個(gè)調(diào)控APX、CAT、POD和GPX等防御保護(hù)性酶合成的基因，并發(fā)現(xiàn)在N株系中的表達(dá)量高于WT。研究表明，這些酶可以清除植物體內(nèi)多余的ROS，從而保護(hù)植株免受枯萎病侵害[28]，本研究結(jié)果與其一致。

在植物防御反應(yīng)中共鑒定到33個(gè)PR10蛋白，而PR10蛋白是植物抵御外界生物脅迫與非生物脅迫產(chǎn)生的一類病程相關(guān)蛋白。植物在受到尖孢鐮刀菌、黃萎病、青枯菌或者外源激素侵染后均能夠誘導(dǎo)PR10蛋白的表達(dá)[29-31]。本研究還發(fā)現(xiàn)了3個(gè)MLP423基因，而前人研究發(fā)現(xiàn)煙草NbMLP423通過(guò)內(nèi)源激素信號(hào)途徑參與調(diào)控?zé)煵輰?duì)赤星病的抗性及其他生物過(guò)程[32]。推測(cè)GhB301基因轉(zhuǎn)入激活了植物防御反應(yīng)從而提高了棉花對(duì)枯萎病的抗性。本研究還發(fā)現(xiàn)9個(gè)MLO基因，在2個(gè)材料中的表達(dá)量差異較大，說(shuō)明轉(zhuǎn)入ERF轉(zhuǎn)錄因子能夠影響MLO基因的表達(dá)，而MLO基因在植物激素信號(hào)轉(zhuǎn)導(dǎo)中發(fā)揮著重要作用[33]。此外，Wang等[34]研究表明，NIMIN1.2可以通過(guò)介導(dǎo)在擬南芥中的氧化還原穩(wěn)態(tài)和JA/ET途徑來(lái)積極調(diào)控灰霉病的抗性。推測(cè)本研究發(fā)現(xiàn)的4個(gè)NIMIN-1基因可能通過(guò)茉莉酸/乙烯途徑影響棉花對(duì)枯萎病的抗性。HSPRO2基因在植物防御反應(yīng)中起重要作用，能夠響應(yīng)丁香假單胞菌的防御反應(yīng)[35-36]。本研究檢測(cè)到6個(gè)HSPRO2蛋白在N株系中的表達(dá)量高于WT，說(shuō)明HSPRO2蛋白在棉花抵御枯萎病中發(fā)揮作用。

在苯丙烷代謝途徑中，研究發(fā)現(xiàn)差異表達(dá)基因主要集中在BGLU、CAD和4CL等相關(guān)基因，其中BGLU基因作為香豆素合成的重要前體，其上調(diào)表達(dá)可以增加香豆素的含量。香豆素作為一種植物抗毒素，當(dāng)植物受到病原菌侵染后迅速生成[37]。接種枯萎病菌后，BGLU相關(guān)基因表達(dá)量增加，表明轉(zhuǎn)GhB301基因可能通過(guò)增加植物抗毒素從而使棉花抗病性增強(qiáng)。綜上所述，與WT相比，轉(zhuǎn)GhB301的棉花中存在大量的DEGs，其在表達(dá)量及表達(dá)時(shí)效方面均呈現(xiàn)出較強(qiáng)優(yōu)勢(shì)，能夠參與多種植物抗病反應(yīng)，從而增強(qiáng)棉花對(duì)枯萎病的抗性。

4 結(jié)論

對(duì)轉(zhuǎn)基因棉花株系與野生型對(duì)照進(jìn)行抗病性調(diào)查，發(fā)現(xiàn)過(guò)表達(dá)GhB301轉(zhuǎn)基因的棉花株系對(duì)枯萎病的抗病性顯著增強(qiáng)。對(duì)枯萎病菌侵染初期的棉花材料進(jìn)行RNA-seq分析，發(fā)現(xiàn)了大量表達(dá)有關(guān)ROS產(chǎn)生和清除的氧化還原基因，以及與植物防御反應(yīng)相關(guān)的基因，并篩選出多個(gè)可能與棉花抗枯萎病有關(guān)的候選基因。推測(cè)過(guò)表達(dá)GhB301基因在棉花抵抗枯萎病的動(dòng)態(tài)平衡中發(fā)揮著重要作用。KEGG富集分析發(fā)現(xiàn)差異表達(dá)基因主要參與植物激素信號(hào)轉(zhuǎn)導(dǎo)、植物病原體相互作用和苯丙烷生物合成等途徑。本研究結(jié)果為揭示GhB301基因的功能及其調(diào)控網(wǎng)絡(luò)奠定了基礎(chǔ)。