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

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

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

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

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

1 材料與方法

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

棉花枯萎病菌7號生理小種F430菌株由石河子大學棉花分子育種實驗室保存。棉花枯萎病菌在固體馬鈴薯葡萄糖瓊脂(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。最后用無菌水將孢子懸浮液的濃度調(diào)整至1×107個·mL-1，備用。

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

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

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

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

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

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

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

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

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

1.5 基因表達量與注釋

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

1.6 通過qRT-PCR驗證RNA-seq

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

2 結(jié)果與分析

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

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

2.2 RNA-seq分析

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

注：A：3次重復的病級統(tǒng)計；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 棉花抗病性檢測Fig.1 Detection of disease resistance in Gossypium hirsutum

2.3 DEGs篩選與注釋

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

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

注：A：總差異表達基因GO富集分析；其中1氧化應(yīng)激反應(yīng)；2：氨基多糖代謝過程；3：氨基多糖分解代謝過程；4：幾丁質(zhì)代謝過程；5：幾丁質(zhì)分解代謝過程；6：氨基糖代謝過程；7：細胞壁大分子分解代謝過程；8：細胞壁大分子代謝過程；9：氨基糖分解代謝過程；10：含谷氨酸的復合代謝過程；11：細胞外區(qū)；12：細胞外圍；13：質(zhì)外體；14：細胞壁；15：外部封裝結(jié)構(gòu)；16：外囊腫；17：細胞皮層；18：細胞皮層部分；19：細胞質(zhì)區(qū)；20：光系統(tǒng)II；21：氧化還原酶活性，作用于過氧化物；22：過氧化物酶活性；23：抗氧化活性；24：銅離子結(jié)合；25：氧化還原酶活性作用于單個供體；26：氧化還原酶活性作用于供體；27：雙加氧酶活性；28：幾丁質(zhì)酶活性；29：十葡聚糖：木糖基轉(zhuǎn)移酶活性；30：轉(zhuǎn)移酶活性；B：不同時期的WT樣本比較；C：不同時期的N樣本比較；D：相同時期的兩樣本間的比較。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 差異表達基因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ān)DEGs

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

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

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

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

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

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

在植物與病原物長期演化過程中，植物產(chǎn)生了一系列防御反應(yīng)機制來阻止病原物的侵害[20]。對所有DEGs進行GO富集，發(fā)現(xiàn)共有67個與防御反應(yīng)(GO:0006952)相關(guān)的DEGs(圖4)，共注釋到33個PR10基因。這些PR10基因中存在20個MLP(major latex protein)基因、7個BETV1(major pollen allergen bet v 1)基因、4個NCS2(s-norcoclaurine synthase 2)基因、2個STH-2(pathogenesis-related protein STH-2)基因。Gh_D02G1678、Gh_A03G1240和Gh_D04G1394均注釋為MLP423基因，前2個基因在2種材料中均呈現(xiàn)先上調(diào)后下調(diào)的趨勢，而Gh_D04G1394在N株系和WT中的表達模式不同，其在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)蛋白的表達受到生物與非生物脅迫的影響[21]。在這些與防御反應(yīng)相關(guān)的DEGs中共存在9個MLO基因，整體來說，在N株系中的表達量均低于在WT中的表達量。其中Gh_D07G1738和Gh_A07G1067在N和WT中皆下調(diào)表達，而Gh_D02G0670、Gh_A10G1852和Gh_D05G0753在2種材料中均顯示先上調(diào)后下調(diào)的趨勢。

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

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

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

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

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

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

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

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

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

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

2.7 qRT-PCR檢驗

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

3 討論

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

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

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

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

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

4 結(jié)論

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