转GhB301基因棉花响应枯萎病菌侵染的转录组分析

刘戈辉 韩泽刚 孙士超 张 薇,*

(1 石河子大学农学院，新疆 石河子 832000；2 南京农业大学作物遗传与种质创新国家重点实验室，江苏 南京 210095)

棉花枯萎病(Fusariumoxysporumf.sp.Vasinfectum)作为棉花的土传病害之一，其病菌可以从根部侵染棉花，毒害棉花维管束，造成棉花叶片青枯、落叶，甚至死亡。其病菌孢子在土壤中的休眠期长达10年之久，农药防治对病情控制的作用有限，并且会对土壤水源产生污染，因此培育棉花抗病品种是解决棉花枯萎病最为经济有效的措施。近年来，较多的研究集中于棉花黄萎病，且主要集中在抗病种质资源的筛选、抗性基因的定位以及基因的遗传分析等方面[1-2]，而对棉花枯萎病的研究相对较少。

随着高通量测序技术的迅猛发展，其具有时间短、通量大、灵敏度高等优点，在研究差异表达基因(Differentially expressed genes,DEGs)方面显示出强大的优势[3-4]。第二代高通量转录组测序(RNA-seq)已被广泛应用于植物基因表达分析研究。韩泽刚等[5]通过Solexa高通量测序技术构建棉花抗、感枯萎病品种的基因表达谱，发现1 080个DEGs在棉花抗病和感病品种中差异表达，并发现在植物激素信号转导通路中，26个已知功能的基因可能参与植物抗病反应。王琛等[6]通过分析比较接种和未接种枯萎病菌的棉花叶片转录组测序结果，发现枯萎病菌侵染后，5 062个基因的表达量发生显著变化，其中2 091个基因的转录水平升高，2 971个基因的转录水平下降。杜荣光[7]通过对抗感枯萎病海岛棉材料接菌前后进行转录组测序分析，发现众多已知抗病相关基因均显著高表达。而这些研究均未对比分析转基因材料和对照材料。

本试验在前期研究中获得了1个ERF转录因子家族基因GhB301，在受到枯萎病菌诱导后，该基因的相对表达量明显增加，并且在抗病品种中的表达量显著高于感病品种；在烟草中过表达该基因，可以显著增加防御保护酶活性以及抗病相关基因的表达，暗示该基因可能在棉花对抗枯萎病菌胁迫响应中发挥重要作用[8]。在此基础上，本研究将GhB301转入陆地棉中并获得纯合株系，研究转基因棉花株系与野生型对照的抗病性变化，包括分析枯萎病菌诱导前后转基因系和对照系中的差异表达基因，筛选与棉花抗枯萎病菌相关的候选基因，在转录水平上分析GhB301基因在棉花抗枯萎病中的功能。

1 材料与方法

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

棉花枯萎病菌7号生理小种F430菌株由石河子大学棉花分子育种实验室保存。棉花枯萎病菌在固体马铃薯葡萄糖琼脂(potato dextrose ager,PDA)培养基上培养5 d后(28℃)，转移至查氏液体培养基中，于200 r·min-1的孵育振荡器(哈尔滨东安电子)中，在28℃黑暗培养7～10 d。最后用无菌水将孢子悬浮液的浓度调整至1×107个·mL-1，备用。

以转GhB301基因棉花株系(N)与野生型受体YZ-1(WT)为材料，经硫酸脱绒的种子用10%的双氧水浸泡2～3 h，无菌水冲洗5～6遍后，种植于无菌的营养钵(花土∶蛭石=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株系开始发病，N株系无明显病症；接菌第20天后WT和N植株均发病，对发病情况进行调查(图1-A)。结果表明，N株系中病级多分布在1级或2级，而WT株系中病级则多分布在3级及4级，且明显多于N株系，说明WT株系比N株系感染枯萎病菌的情况更为严重。对3次病级统计进行病情指数计算(图1-B)，发现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比对到参考基因组中，发现比对率的均值是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可以比对至参考基因组的唯一位置，为后续分析奠定了坚实的数据基础。

注：A：3次重复的病级统计；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为筛选条件，共发现14 021个DEGs(其中1 290个新基因)能够响应枯萎病菌诱导。对这些DEGs进行GO富集分析(图2-A)，共富集注释到30个亚类，涉及生物过程，细胞组成和分子功能三大类，主要富集在抗氧化活性、氧化还原酶活性、过氧化物酶活性、氧化应激的反应、光系统Ⅱ、碳水化合物结合以及细胞壁等功能类别，其中富集DEGs最多的是氧化还原酶活性(GO:0016684)，共注释到135个基因，最少的是光系统Ⅱ，仅注释到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最少。比较发现，转基因株系在病原菌侵染后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组内比较发现共存在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基因合影响棉花植株的生长发育。

注：A：总差异表达基因GO富集分析；其中1氧化应激反应；2：氨基多糖代谢过程；3：氨基多糖分解代谢过程；4：几丁质代谢过程；5：几丁质分解代谢过程；6：氨基糖代谢过程；7：细胞壁大分子分解代谢过程；8：细胞壁大分子代谢过程；9：氨基糖分解代谢过程；10：含谷氨酸的复合代谢过程；11：细胞外区；12：细胞外围；13：质外体；14：细胞壁；15：外部封装结构；16：外囊肿；17：细胞皮层；18：细胞皮层部分；19：细胞质区；20：光系统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)作为植物活性氧产生与清除的氧化还原产物，能够参与植物的生长发育及应激反应[18]。对DEGs进行GO富集分析，发现共有135个DEGs参与氧化还原过程(GO:0016684)(图3)。

图3 与氧化还原过程相关基因的热图Fig.3 The heat map of genes related to oxidation-reduction process

呼吸爆发氧化酶同系物(respiratory burst oxidase homologs，RBOHs)，又被称作植物NADPH氧化酶，在细胞生长、植物发育、植物对非生物环境约束以及与生物相互作用的反应中作为ROS生产者发挥作用，是重要的分子“中枢”[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中则呈现持续下降的趋势。

对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参与其中，其表达量复杂多样的动态变化也反映出棉花在抵抗枯萎病感染的动态平衡中，氧化还原过程相关基因起着重要作用。

2.5 植物防御反应相关DEGs

在植物与病原物长期演化过程中，植物产生了一系列防御反应机制来阻止病原物的侵害[20]。对所有DEGs进行GO富集，发现共有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种材料中均显示先上调后下调的趋势。

另外，还发现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通路富集统计分析Table 3 Analysis of KEGG function annotation of differential expression genes

研究发现苯丙烷代谢途径以及与这些途径相关的基因在棉花与病原菌的互作中发挥着重要作用[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才逐渐发挥作用，然而在N株系中，该基因侵染初期即开始表达。Gh_D13G0844基因在N中均表现出先升高后降低的趋势，在WT中则表现出逐渐下降趋势；且该基因在N中的表达量高于在WT中的表达量。

在苯丙烷代谢通路中，还发现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基因能够增强防卫基因的表达，并且ERF14突变体对枯萎病菌表现为更易感病[25]。本研究通过抗病性鉴定，发现接种枯萎病菌后过表达GhB301转基因棉花与WT相比病情指数显著降低，表明转GhB301基因可以增强棉花的抗病性。

Fradin等[26]研究发现真菌孢子的萌发以及向表皮细胞的渗透一般发生在接菌后的0～12 h。而植物的应激反应、信号转导和信息交换一般发生在真菌感染早期的40 h左右[18,27]。本研究通过比较差异表达基因，发现接菌后6 h时，转基因棉花中大量的差异表达基因开始响应，而WT棉花则在12 h开始大量表达差异表达基因，说明转GhB301棉花能够更快响应枯萎病菌的胁迫，这些差异表达基因的提前响应可能会增强棉花的抗病性。

在这些差异表达基因中，本研究重点关注了参与氧化还原过程、植物防御反应以及苯丙烷代谢途径的基因，其中氧化还原过程主要是植物产生和清除ROS的过程[19]。本研究鉴定了9个呼吸爆发氧化酶同系物(RBOHs)，并且发现在N株系中的表达早于在WT中的表达。当植物受到病原菌等非生物胁迫时，植物在短时间内产生大量的ROS，如果无法及时清除ROS，则会伤害植株。本研究鉴定出多个调控APX、CAT、POD和GPX等防御保护性酶合成的基因，并发现在N株系中的表达量高于WT。研究表明，这些酶可以清除植物体内多余的ROS，从而保护植株免受枯萎病侵害[28]，本研究结果与其一致。

在植物防御反应中共鉴定到33个PR10蛋白，而PR10蛋白是植物抵御外界生物胁迫与非生物胁迫产生的一类病程相关蛋白。植物在受到尖孢镰刀菌、黄萎病、青枯菌或者外源激素侵染后均能够诱导PR10蛋白的表达[29-31]。本研究还发现了3个MLP423基因，而前人研究发现烟草NbMLP423通过内源激素信号途径参与调控烟草对赤星病的抗性及其他生物过程[32]。推测GhB301基因转入激活了植物防御反应从而提高了棉花对枯萎病的抗性。本研究还发现9个MLO基因，在2个材料中的表达量差异较大，说明转入ERF转录因子能够影响MLO基因的表达，而MLO基因在植物激素信号转导中发挥着重要作用[33]。此外，Wang等[34]研究表明，NIMIN1.2可以通过介导在拟南芥中的氧化还原稳态和JA/ET途径来积极调控灰霉病的抗性。推测本研究发现的4个NIMIN-1基因可能通过茉莉酸/乙烯途径影响棉花对枯萎病的抗性。HSPRO2基因在植物防御反应中起重要作用，能够响应丁香假单胞菌的防御反应[35-36]。本研究检测到6个HSPRO2蛋白在N株系中的表达量高于WT，说明HSPRO2蛋白在棉花抵御枯萎病中发挥作用。

在苯丙烷代谢途径中，研究发现差异表达基因主要集中在BGLU、CAD和4CL等相关基因，其中BGLU基因作为香豆素合成的重要前体，其上调表达可以增加香豆素的含量。香豆素作为一种植物抗毒素，当植物受到病原菌侵染后迅速生成[37]。接种枯萎病菌后，BGLU相关基因表达量增加，表明转GhB301基因可能通过增加植物抗毒素从而使棉花抗病性增强。综上所述，与WT相比，转GhB301的棉花中存在大量的DEGs，其在表达量及表达时效方面均呈现出较强优势，能够参与多种植物抗病反应，从而增强棉花对枯萎病的抗性。

4 结论

对转基因棉花株系与野生型对照进行抗病性调查，发现过表达GhB301转基因的棉花株系对枯萎病的抗病性显著增强。对枯萎病菌侵染初期的棉花材料进行RNA-seq分析，发现了大量表达有关ROS产生和清除的氧化还原基因，以及与植物防御反应相关的基因，并筛选出多个可能与棉花抗枯萎病有关的候选基因。推测过表达GhB301基因在棉花抵抗枯萎病的动态平衡中发挥着重要作用。KEGG富集分析发现差异表达基因主要参与植物激素信号转导、植物病原体相互作用和苯丙烷生物合成等途径。本研究结果为揭示GhB301基因的功能及其调控网络奠定了基础。