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

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

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

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

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

1 材料與方法

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

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

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

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

1.2 表型鑒定與病情指數分析

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

1.3 RNA提取和轉錄組測序

參照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，美國)試劑盒構建cDNA文庫。使用Qubit2.0 Fluorometer(賽默飛，美國)對文庫進行初步定量，稀釋文庫至1.5 ng·μL-1。使用Agilent 2100 bioanalyzer對文庫的insert size進行檢測。庫檢合格后，委托北京諾禾致源生物信息科技有限公司利用Illumina HiSeq 4000進行測序。

1.4 數據指控和參考基因組比對

對原始數據進行過濾，去除帶接頭、無法確定堿基信息以及低質量reads，得到clean reads。對clean reads分別進行Q20、Q30和GC含量計算。從CottonGene數據庫(http://www.cottongen.org)下載陸地棉TM-1參考基因組[11]和基因模型注釋文件，利用軟件Tophat2[12]將高質量的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)驗證，反應體系及程序參考北京全式金生物技術有限公司Trans start Tip Green qPCR Super Mix說明書進行，試驗設3個重復。采用2-ΔΔCt法[17]計算基因相對表達量，所用引物見表1(用Primer5軟件設計引物，由北京華大基因公司合成)。

2 結果與分析

2.1 表型鑒定及病情指數調查

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

2.2 RNA-seq分析

轉錄組分析結果表明，共獲得279 114 320個raw reads，過濾后保留273 111 170個clean reads，數據量為81.95 Gb，每個文庫平均8.20 Gb，如表2所示，WT棉花測序共得到131 456 222個raw reads，對低質量reads進行過濾后仍有128 669 540個clean reads，其中clean reads數最多的是N6，最少的是N48。而N材料測序后共得到147 658 098個raw reads，過濾后剩余144 441 630個clean reads，其中clean reads數最多的是WT24，最少的是WT48。從測序結果可以看出，各樣本文庫的Q30值均超過87.64%，GC含量均大于43.06%，其平均值分別為89.86%和43.28%。將clean reads比對到參考基因組中，發(fā)現比對率的均值是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ù)分析奠定了堅實的數據基礎。

注：A：3次重復的病級統(tǒng)計；B：病情指數分析；**表示與WT相比，N的病情指數在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篩選與注釋

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

組間比較結果顯示，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組內比較發(fā)現共存在694個共同的DEGs(圖2-B)，轉N組內比較則有829個共同的DEGs(圖2-C)。組間比較結果顯示，存在41個共同的DEGs(圖2-D)，推測這些共同的DEGs可能與棉花枯萎病脅迫響應具有密切關系。其中，值得特別關注的是，Gh_A03G0673(MYB1R1)和Gh_A06G1626(抗病蛋白RPS5)在N植株中均不表達，而在WT植株中表達，推測轉入GhB301基因可能促進ERF轉錄因子的表達而抑制MYB轉錄因子和PR基因的表達。Gh_A02G1734(戊二酸/鐵依賴性雙加氧酶)基因在N植株中的表達量均高于WT，推測轉入GhB301基因可以提高氧化還原酶相關基因的表達。Gh_A13G2026(生長素應答蛋白SAUR32)只在N株系中表達，推測轉入GhB301基因合影響棉花植株的生長發(fā)育。

注：A：總差異表達基因GO富集分析；其中1氧化應激反應；2：氨基多糖代謝過程；3：氨基多糖分解代謝過程；4：幾丁質代謝過程；5：幾丁質分解代謝過程；6：氨基糖代謝過程；7：細胞壁大分子分解代謝過程；8：細胞壁大分子代謝過程；9：氨基糖分解代謝過程；10：含谷氨酸的復合代謝過程；11：細胞外區(qū)；12：細胞外圍；13：質外體；14：細胞壁；15：外部封裝結構；16：外囊腫；17：細胞皮層；18：細胞皮層部分；19：細胞質區(qū)；20：光系統(tǒng)II；21：氧化還原酶活性，作用于過氧化物；22：過氧化物酶活性；23：抗氧化活性；24：銅離子結合；25：氧化還原酶活性作用于單個供體；26：氧化還原酶活性作用于供體；27：雙加氧酶活性；28：幾丁質酶活性；29：十葡聚糖：木糖基轉移酶活性；30：轉移酶活性；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富集分析及不同樣本間的基因數目比較Fig.2 The GO enrichment analysis of the total differentially expressed genes and the number of genes between different samples

2.4 氧化還原過程中相關DEGs

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

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

呼吸爆發(fā)氧化酶同系物(respiratory burst oxidase homologs，RBOHs)，又被稱作植物NADPH氧化酶，在細胞生長、植物發(fā)育、植物對非生物環(huán)境約束以及與生物相互作用的反應中作為ROS生產者發(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種材料中均呈現先上升后下降的趨勢，在N株系中，其表達量在6 h即達到峰值，WT中則在12 h達到峰值。而RBOHC、RBOHD、RBOHF在N植株中呈現先下降后上升的趨勢，在WT中則呈現持續(xù)下降的趨勢。

對ROS清除劑抗壞血酸過氧化物酶(ascorbic acid peroxidase，APX)，過氧化氫酶(catalase，CAT)，過氧化物酶(peroxidase，POD)，谷胱甘肽過氧化物酶(glutathione peroxidase，GPX)等進行表達模式分析，共檢測到5種APX蛋白質基因(Gh_D13G2402、Gh_A06G0383、Gh_A08G1745、Gh_A02G1648、Gh_A01G1388)，這些基因在N株系中表達量均明顯高于WT。共鑒定到2種CAT基因(Gh_D03G0021和Gh_A05G1539)，在2種材料中均表現出先上調后下調的趨勢；2種POD基因(Gh_A08G0714和Gh_D08G0832)，在N株系中上調表達，在WT中則先上調后下調；鑒定出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產生與清除過程中，眾多DEGs參與其中，其表達量復雜多樣的動態(tài)變化也反映出棉花在抵抗枯萎病感染的動態(tài)平衡中，氧化還原過程相關基因起著重要作用。

2.5 植物防御反應相關DEGs

在植物與病原物長期演化過程中，植物產生了一系列防御反應機制來阻止病原物的侵害[20]。對所有DEGs進行GO富集，發(fā)現共有67個與防御反應(GO:0006952)相關的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種材料中均呈現先上調后下調的趨勢，而Gh_D04G1394在N株系和WT中的表達模式不同，其在N株系先上調后下調，而在WT中則先下調后上調。

圖4 與植物防御反應相關基因的熱圖Fig.4 The heat map of genes related to plant defense response

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

另外，還發(fā)現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個材料中均呈現下調趨勢，且在N株系中的表達量均低于在WT中的表達量。PR-4B基因在2個材料中均呈現逐漸升高的趨勢，特別是在48 h，在N株系中的表達量遠高于WT。4個NIMIN-1基因則在2個材料中表現出相反的表達趨勢，在N株系中表現為先降低后升高的表達趨勢，在WT中則相反。6個HSPRO2基因均大致表現出先降低后升高的趨勢，在N中大部分的基因均在24 h降到最低點，而在WT中則在6 h降到最低點。在48 h，N中所有基因的表達量均高于WT中對應基因的表達量。

2.6 苯丙烷代謝途徑相關DEGs

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

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

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

圖5 與苯丙烷代謝途徑相關基因的熱圖分析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中均表現出先升高后降低的趨勢，在WT中則表現出逐漸下降趨勢；且該基因在N中的表達量高于在WT中的表達量。

在苯丙烷代謝通路中，還發(fā)現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株系中的響應時間早于WT，并且表達量也更高，推測接種枯萎病菌后，N株系苯丙烷代謝途徑激活響應的速度更快，這可能與其抗病性顯著增強有密切關系。

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

2.7 qRT-PCR檢驗

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

3 討論

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

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

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

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

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

4 結論

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