基于RNA-seq的番茄青枯病响应基因鉴定与表达分析

史建磊，熊自立，苏世闻，王克磊，宰文珊

(1.温州科技职业学院，浙江 温州 325006；2.福建农林大学 农学院，福建 福州 350002)

青枯病是一种毁灭性的土传维管束病害，又称细菌性枯萎病，病原为青枯雷尔氏菌(Ralstoniasolanacearum，Rs)，具有明显的生理分化和寄主多样性。作为世界上最重要的植物病原细菌之一，青枯菌可从植物伤口或次生根的根冠侵入，一定条件下大量繁殖并快速传播，破坏并阻塞维管束，影响水分运输，从而导致植株萎蔫死亡。该病主要发生在热带亚热带气候湿润、土壤呈酸性的地区，在我国主要发生在长江流域及以南的广大地区，但受全球气候变暖等的影响，有向北蔓延的趋势，是造成多种作物减产、甚至绝收的主要原因，特别是对茄科作物番茄的危害尤为严重。番茄(Solanumlycopersicum)是全球重要蔬菜，但遗传基础狭窄，经常受到各种病害的影响。由于连作重茬、土壤富集病菌和高温高湿气候，我国番茄产区青枯病日趋严重，已成为番茄生产的主要障碍。

随着二代测序技术(Next generation sequencing，NGS)的快速发展和广泛应用，转录组测序(RNA-seq)已成为研究基因表达变化和表达调控网络的有效手段，为抗病候选基因的鉴定和分子标记开发提供了基础。目前，RNA-seq技术已在番茄多种生物或非生物胁迫研究中得到广泛应用[1-8]。关于青枯病响应基因的转录组鉴定，Ishihara等[9]利用基因芯片技术对茎刺法接种Rs的番茄茎进行转录组分析，发现致病、激素信号和木质素生物合成相关的140个基因在抗病番茄LS-89中上调表达，β-1，3-葡聚糖酶基因尤为明显；French等[10]对伤根法接种Rs的番茄根进行了转录组分析，发现抗病番茄Hawaii7996根部防御基因更早更强被激活，同时生长素信号途径被抑制；Zuluaga等[11]对伤根法接种Rs的野生马铃薯根进行转录组分析，在抗感材料中分别发现221，644个差异表达基因，前者诱导产生大量新的转录本，后者激素相关基因更突出；Chen等[12]对伤根法接种Rs的茄子茎进行转录组分析，发现抗病材料分别有1 137，9 048个基因上调或下调表达，而感病材料对应数目是6 087，5 832个；谢冰等[13]叶转录组分析发现，Y87作为砧木能够提高接穗HD PAL活性和木质素、多酚及病程相关蛋白等多种抗性基因表达量，从而提高烟草青枯病抗性。

面对各种病原微生物伤害，植物在漫长的进化过程中形成了完善的自我保护机制，抗病基因(R基因)在识别和抵御病原菌入侵方面具有重要作用。尽管前人已在番茄青枯病响应方面进行了相关研究，但不同材料、不同接种方法、不同取样部位等获得的研究结果不尽一致，有必要拓展和丰富这一课题。

本研究利用高通量测序技术对伤根法接种Rs后的抗感番茄自交系地上部进行转录组分析，结合基因功能注释，挖掘青枯病响应相关基因并进行RT-qPCR验证，以期为候选基因抗源的有效利用提供理论指导。

1 材料和方法

1.1 试验材料

番茄AH13112111为温州农业科学院番茄课题组2013年春从荷兰瑞克斯旺公司引进品种爱好与自育AH1杂交选育的高代自交系，黄色小果型，苗期人工接种鉴定为高抗青枯病；G149351121为以色列海泽拉公司同步引进的品种3951与自育G14杂交选育的高代自交系，红色普通型，对青枯病表现为高感。青枯菌株(生理小种1生化变种Ⅲ)从田间发病番茄植株上自行分离。

1.2 试验方法

1.2.1 建库与测序 番茄幼苗长到4～6片真叶时用伤根法接种青枯菌，接后48 h植株萎蔫前分别取幼嫩组织，锡箔纸包裹，立即置于液氮中，再放入-80 ℃超低温冰箱中保存备用[14]。设置3次生物学重复。参照植物TRIzol总RNA提取试剂盒(Sangon Biotech)说明书提取RNA。用带有Oligo(dT)的磁珠富集mRNA；加入Fragmentation Buffer将mRNA进行随机打断；以mRNA为模板，用六碱基随机引物合成第1条cDNA链，然后加入缓冲液、dNTPs、RNase H和DNA Polymerase Ⅰ合成第2条cDNA链，利用AMPure XP beads纯化cDNA；纯化的双链cDNA再进行末端修复、加A尾并连接测序接头，然后用AMPure XP beads进行片段大小选择；最后通过PCR富集得到cDNA文库。文库质控合格后进行Illumina转录组测序(Biomarker Technologies)。

1.2.2 数据质量控制 去除原始数据中带有接头序列和低质量的reads，得到Clean reads。计算Clean reads的Q30值和GC含量。从茄科基因组数据库(Solanaceae Genomics Network，https：//solgenomics.net/organism/Solanum_lycopersicum/genome)下载番茄参考基因组(ITAG4.0)的序列文件和注释文件。使用HISAT2[15]进行测序序列的比对。比对分析完成后利用StringTie[16]对符合要求的reads进行组装和定量。基于所选参考基因组序列，使用StringTie对mapped reads进行拼接，并与原有的基因组注释信息进行比较，寻找原来未被注释的转录区，发掘该物种的新转录本和新基因。

1.2.3 差异基因筛选 以差异倍数(Fold Change，FC)≥2且错误发现率(False Discovery Rate，FDR)<0.01作为标准，利用DESeq2[17]进行样品组间的差异表达分析，筛选差异表达基因。通过Blast比对工具，将测序基因与蛋白相邻类聚簇数据库(Cluster of orthologous groups of proteins，COG)[18]、基因本体论数据库(Gene ontology consortium，GO)[19]、东京基因与基因组百科全书(Kyoto encyclopedia of genes and genomes，KEGG)[20]、蛋白质真核同源数据库(Eukaryotic orthologous groups，KOG)[21]、非冗余蛋白数据库(Non-redundant protein datebase，NR)[22]、Pfam[23]和蛋白质序列数据库(Swiss prot protein datebase，Swiss-Prot)[24]进行序列比对，获得基因的注释信息。

1.2.4 基因定量验证 以SlRPL2为内参基因，利用Primer Premier 5.0软件设计定量引物(表1)。RT-qPCR反应体系包括2 μL cDNA，0.4 μL PCR primer，10 μL SYBR，7.2 μL ddH2O。扩增程序为95 ℃ 3 min；95 ℃ 5 s，60 ℃ 30 s，45个循环，整体反应在StepOne Plus型荧光定量PCR仪(ABI，Foster，CA，美国)中进行，采用2-ΔΔCt法计算基因相对表达量。用SPSS 17.0对数据进行单因素方差分析，用最小显著极差法(LSD)进行多重比较和差异显著性分析。

表1 部分RT-qPCR引物Tab.1 Part of the primers used for RT-qPCR

2 结果与分析

2.1 测序数据质量评估

经检测，样本RNA的OD260/280均为2.2，OD260/230≥1.8，28S/18S ≥2.0，RIN≥8，表明所提RNA满足后续分析需求。抗感番茄接种前后12个样本原始测序数据去除接头序列和低质量序列，共获得75.02 Gb Clean reads，各样本净碱基数均不小于5.73 Gb，GC含量在43.27%～44.17%，Q30碱基百分比均超过了94%(表2)。将各样本的Clean reads与指定参考基因组进行序列比对，比对效率为95.00%～96.31%，且比对到唯一位置的reads均超过了90%(表3)。基于基因表达量，重复样本间相关系数在0.947～0.993(表4)，生物学重复良好。

表2 番茄各样本RNA测序数据统计Tab.2 RNA sequencing data statistics of tomato libraries

表3 番茄各样本Clean reads与参考基因组比对Tab.3 Comparing Clean reads onto the reference genome in each tomato library

表4 生物学重复样本的相关性Tab.4 Correlation between biological replicates

2.2 差异表达基因筛选

本研究测序共获得35 371个基因，其中新基因1 296个。经log2FC和FDR过滤，在抗感番茄中分别获得970，695个DEGs，二者重复基因353个，DEG总数1 312个(新基因32个)，占基因总数的3.71%。其中，上调基因分别为457，450个，重复214个，总计693个；下调基因分别为513，245个，重复137个，总计621个(表5)。2种材料接种后均上调或下调的基因分别为214，137个；抗病材料上调而感病材料下调的基因数目是0，反之是2(Solyc04g074440.1和Solyc12g096780.2)；抗病材料上调且感病材料不变的基因是198个，反之是216个；抗病材料下调且感病材料不变的基因是316个，反之是104个，后6种表达模式基因总数为836个(新基因18个)(表6)。

表5 差异表达基因数目Tab.5 Number of differentially expressed genes

2.3 差异表达基因注释

本研究注释基因共31 903个，注释率为90.20%。其中，抗感番茄DEG在GO和KEGG等八大数据库的注释数目分别为951，687个，注释率均超过了98%；总注释DEG是1 290个，注释率是98.32%。这些注释基因包括4个NBS类抗病基因、6个植物-病原互作基因、11个植物激素信号转导基因、22个防御反应基因、32个蛋白激酶、65个转录因子及其他多个重要功能因子。

表6 差异表达基因的模式分布Tab.6 Pattern distribution of differentially expressed genes

1 312个DEGs GO注释基因为973个，注释率74.16%。将这些基因划分成47个功能组，根据功能组的不同，又分为生物过程(Biological process，BP，19个亚类)、细胞组分(Cellular component，CC，15个亚类)和分子功能(Molecular function，MF，13个亚类)三大类(图1)。在生物过程中，主要涉及代谢过程、细胞过程、单生物体过程、生物调控和刺激响应等亚类；对于细胞组分，主要涉及细胞、细胞部分、膜、细胞器、膜部分、细胞器部分和大分子复合物等亚类；结合催化活性在分子功能中非常显著，其次是转运活性和核酸结合转录因子活性，再次是电子载体活动、抗氧化活性、结构分子活性和信号传导活性。

1.代谢过程；2.细胞过程；3.单生物体过程；4.生物调控；5.刺激响应；6.定位；7.细胞组成与发生；8.信号；9.发育过程 ；10.多细胞生物体过程；11.生殖过程；12.生物群体过程；13.生长；14.繁殖；15.免疫系统过程；16.细胞死亡；17.节律过程；18.移动；19.生物粘附；20.细胞；21.细胞部分；22.膜；23.细胞器；24.膜部分；25.细胞器部分；26.大分子复合物；27.膜内腔；28.胞间区；29.细胞连接；30.共质体；31.胞间区部分；32.颗粒；33.颗粒部分；34.类核；35.结合；36.催化活性；37.转运活性；38.核酸结合转录因子活性；39.结构分子活性；40.电子载体活性；41.分子传导活性；42.信号传导活性；43.抗氧化活性；44.转录因子活性-蛋白结合；45.养分存储活性；46.蛋白质标签；47.金属伴侣活性。

本研究209个(15.93%)DEGs被注释到88个代谢通路中，分为细胞过程、环境信息处理、遗传信息处理、代谢和生物系统五大类代谢途径(图2-A)。其中，植物激素信号转导和苯丙烷生物合成相关基因分别为16个，淀粉和蔗糖代谢相关基因13个，植病互作相关基因11个，同时注释还包括核酸、糖、脂肪、蛋白、萜类和黄酮类等的生物合成与代谢相关基因。依据q-value和富集因子大小，取前20个最显著富集的代谢途径(图2-B)，主要包括DNA复制(9个DEGs)、二萜合成(7个DEGs)、组氨酸代谢(6个DEGs)、苯丙烷生物合成(16个DEGs)、类固醇生物合成(5个DEGs)、类黄酮生物合成(6个DEGs)、萜类骨架合成(5个DEGs)、谷胱甘肽代谢(9个DEGs)。另外，植病互作(rank 38)、淀粉和蔗糖代谢(rank 46)、植物激素信号转导(rank 47)相关基因分别为11，13和16个。

图2 差异表达基因KEGG注释(A)与富集(B)Fig.2 KEGG annotation(A)and enrichment(B)of differentially expressed genes

植病互作通路在抗感番茄接种后显著富集，抗病番茄KEGG分析该通路中有7个基因差异表达，感病番茄该通路中有6个基因差异表达，2个材料中共有差异表达基因11个(图3)。该通路中抗病番茄接种后CNGCs(Solyc02g086980.3和Solyc05g050380.4)、EIX1/2(Solyc07g008620.1)、FLS2(Solyc02g031790.3、Solyc02g068820.3和Solyc02g070890.3)和WRKY25/33(Solyc09g014990.4)上调表达；感病番茄接种后CaM/CML(Solyc04g058170.1)、FLS2(Solyc02g068820.3)、PIK1(Solyc06g075550.3)和RPM1(Solyc01g073985.1)上调表达，CNGCs(Solyc05g050380.4和Solyc09g007840.3)既有上调表达又有下调表达。

图3 植病互作通路差异表达基因KEGG富集Fig.3 KEGG enrichment of differentially expressed genes in the plant-pathogen interaction pathway

2.4 差异表达基因启动子分析

基于基因功能注释，选取74个DEGs进行启动子分析。结果发现这些基因启动子主要包括转录起始(TATA-box和CAAT-box)、防御与胁迫(W box、TC-rich repeats和AT-rich sequence)、干旱(MBS)和低温(LTR)、MeJA(CGTCA-motif和TGACG-motif)、ABA(ABRE)、SA(TCA-element)、GA(TATC-box、Pbox和GARE-motif)和AUX(TGA-element、TGA-box、AuxRR-core和AuxRE)、光响应及厌氧诱导等元件(图4)。其中，29个DEGs含42个防御与胁迫响应相关元件。

2.5 差异表达基因的诱导表达

选取上述50个DEGs进行RT-qPCR验证。与转录组测序相比，28个基因青枯菌接种前后的表达变化趋势一致，结合差异显著性分析，11个基因表达调控模式一致，其中6个基因完全一致(图5)。番茄接种青枯菌后，Solyc02g086980.3和Solyc04g011670.3在不同材料中相对表达量均极显著上调(P<0.01)，抗病材料中分别提高2.62，2.36倍，感病材料分别提高1.06，1.34倍，为正调控基因；其他4个基因(Solyc01g073985.1、Solyc09g092580.4、Solyc09g098100.4和Solyc10g081300.1)表达量在感病材料中极显著上调(P<0.01)，分别提高1.16，1.19，1.63，0.87倍，但在抗病材料中变化不显著(P>0.05)，为负调控基因。推测Solyc02g086980.3和Solyc04g011670.3可能参与抗病反应，Solyc01g073985.1、Solyc09g092580.4、Solyc09g098100.4和Solyc10g081300.1可能参与感病反应。

A.已知基因；B.新基因。A.Known genes；B.New genes.

3 结论与讨论

番茄青枯病是由青枯菌引起的一种土传维管束病害，获得抗病相关基因并理解其作用机制非常重要，而RNA-seq技术为这一问题的解决提供了有效途径。本研究测序共获得35 371个基因，其中1 296个新基因，1 312个DEGs，836个材料特异DEGs，并对其进行了功能注释和RT-qPCR验证，这为目的基因筛选、功能验证和相关机制解析提供了基础。

3.1 差异表达基因的筛选与注释

本研究测序共获得75.02 Gb高质量数据，各样本均达到5.73 Gb，Q30碱基占比均大于94%，过滤数据与番茄参考基因组比对率均不小于95%，说明样品不存在污染且参考基因组选择合适，所有样本测序质量良好，数据可靠，可用于后续分析。获得的35 371个基因中注释基因数目是31 903个，略多于Heinz 1706的测序基因数30 855个[25]，可能源自材料的不同和技术手段的改进。基于所选参考基因组序列，共发掘到1 296个新基因，说明RNA-seq技术在新基因挖掘上的可行性。生物学重复可有效消除基因表达在不同个体间的可变性，本研究重复样本皮尔逊相关系数大于0.9，表明生物学重复良好。抗感番茄接种前后DEG数目多于Ishihara等[9]的鉴定结果，但少于French等[10]的鉴定结果，主要归因于接种方法、取样部位和基因表达计算的不同。有研究表明，地上部早期对青枯病的防御反应远没有地下部强烈，但萎蔫病症出现后却有大量基因差异表达[26]。抗病番茄DEG要多于感病番茄，尤其是下调基因，说明抗病性主要与病菌侵染后某些基因下调表达有关。不同材料间DEG数目少于接种前后的1 312个，说明差异表达主要来自接种。感病番茄接种前后DEG少于接种前抗感材料间的DEG，且不同材料接种后DEG变少，说明感病番茄对青枯菌侵染的响应不如抗病番茄剧烈，这和前人研究结果一致[9-10]。抗感材料特异表达的836个DEGs可作为病害响应基因和候选目的基因的筛选基础。

R和S分别表示抗病和感病番茄。不同大小写字母分别表示0.01和0.05水平下的差异显著性。R and S represent resistant and susceptible tomato lines，respectively.Different upper and lowercase letters indicate statistically significant differences at the 0.01 and 0.05 level，respectively.

本研究通过抗感材料转录组分析表明，病菌侵染会启动植物相关基因的表达及多条代谢通路。对这些DEG进行GO富集分析，抗病反应主要发生在膜和细胞内部，同时也发生在细胞连接和细胞外部区域等部位，通过调节催化活性、结合、转运活性、转录因子活性、信号转导活性以及抗氧化活性等分子功能，引起以细胞过程、代谢过程、单生物体过程和生物调控为主，定位、刺激响应、细胞成分组成、信号、免疫系统过程等一系列生物过程的协同变化，抵抗病原菌的入侵。通过KEGG分析表明，DEG在植物激素信号转导、淀粉和蔗糖代谢、DNA复制、嘧啶代谢、嘌呤代谢、植病互作、谷胱甘肽代谢、苯丙烷生物合成、氨基酸合成代谢、玉米素生物合成、泛素化蛋白水解、类黄酮生物合成、萜类物质合成、异喹啉生物碱生物合成、芪类化合物合成、ABC转运蛋白、内吞作用、角质蜡质生物合成、油菜素甾醇生物合成、VB1/6和VC代谢和生理节律等多个代谢通路富集，说明这些代谢基因可能参与了抗病调控。本研究选取番茄幼苗地上部进行分析，其对青枯病的防御反应既与地下部部分重叠，又具有自身独特的代谢途径[10]。

基于基因功能注释，50个DEGs用来RT-qPCR验证(这些基因普遍具有转录起始和胁迫响应相关元件)，结果发现，28个基因在番茄接种青枯菌前后表达变化趋势与转录组数据一致。推测Solyc02g086980.3和Solyc04g011670.3可能参与抗病反应，Solyc01g073985.1、Solyc09g092580.4、Solyc09g098100.4和Solyc10g081300.1可能参与感病反应。2种策略的相关性不强，基因表达结果不完全一致，主要归因于计算原理和表示方法的不同。本研究结合2种方法和统计检验，筛选到6个病害响应基因。基因表达受植物生长发育与环境胁迫响应的动态平衡调控，一方面通过R基因的过量表达抵御外界不良刺激，另一方面通过非编码RNA调控和可变剪接等降低基因表达量[27]。然而，本研究也发现多个感病基因(S基因)，与R基因堆叠相比，破坏单个S基因是实现作物广谱和持久抗性的更有效策略[28]。

3.2 抗病相关分子机制初步解析

上述这些物质有些可能对抗性有直接影响，有些则可能进一步形成抗性物质。其中植物-病原互作代谢通路包括PTI(Pattern-triggered immunity)和ETI(Effector-triggered immunity)，当植物受到胞外病原相关分子模式(Pathogen-associated molecular pattern，PAMP)刺激时，Ca2+浓度增加，激活钙离子蛋白激酶(CDPK)和钙调蛋白(CaM/CML)，将信号传导至NADPH氧化酶活性氧中间体(RBOH)和一氧化氮合成酶(NOS)，激活活性氧(ROS)和一氧化氮(NO)的产生，诱导超敏反应(Hypersensitive response，HR)、细胞壁加固和气孔关闭[29]。当细菌鞭毛蛋白被受体激酶(FLS2)识别，通过胞吞作用将信号传递到细胞内，激活下游的丝裂原活化蛋白激酶(MEKK)MKK1/MKK2→MPK4或MKK4/MKK5→MPK3MPK6。信号传入细胞核后，WRKY22/29激活防御基因或WRKY25/33抑制相关基因[30]。本研究中，环腺苷酸门控通道(CNGCs)基因Solyc02g086980.3在抗病番茄接种后上调表达，Solyc05g050380.4在抗感番茄中均上调，Solyc09g007840.3在感病番茄中下调表达，Ca2+浓度发生变化，激活CaM/CML基因Solyc04g058170.1在感病番茄中的表达。青枯菌侵染后，FLS2基因Solyc02g031790.3和Solyc02g070890.3在抗病番茄中均上调表达，Solyc02g068820.3在抗感番茄中均上调表达，激活下游WRKY25/33基因Solyc09g014990.4在抗病番茄中的表达。植物通过ETI防御机制特异的抵抗病原菌入侵，该通路中RPM1和RPS2是2个能与RIN4基因结合的R基因。当RIN4被磷酸化时，与R基因结合，触发HR反应。蛋白激酶RPS5与PBS1结合，也可触发免疫反应[29]。本研究中，感病番茄接种青枯菌后，RPM1基因Solyc01g073985.1和PIK1基因Solyc06g075550.3上调表达。本研究初步探索了抗感番茄在接种青枯菌后的转录组动态，发现相关抗性涉及一个非常复杂的网络调控系统。

本研究转录因子(Transcription factor，TF)主要涉及AP2/ERF、bHLH、bZIP、Dof、MYB、NAC和WRKY等。这些TFs能激活下游生物胁迫响应相关基因的表达，从而提高茄科作物青枯病抗性[31-40]。32个差异蛋白激酶主要为受体类激酶，其能通过磷酸化在细胞信号转导中发挥作用，从而提高茄科植物抗病性，如AhRLK1、CaCDPK15、CaLRR51和CaRop1等能够与其他基因互作正/负调控植物青枯病抗性[41-43]。4个差异NBS基因分正/负调控。Deslandes等[44]克隆获得了首个抗青枯病基因RRS1。另外，在拟南芥突变体SLH1中发现了另一个编码TIR-NBS-LRR-WRKY类抗病蛋白的基因RPS4。NBS类抗病基因AhRRS5通过多种信号途径交联参与烟草青枯病防御反应[45]。Ben等[46]在苜蓿中发现青枯病抗性位点MtQRRS1，包括15个TIR/NBS类抗病基因。面对病原菌入侵，植物会通过激素信号分子激活抗病基因的表达。本研究植物激素信号转导涉及生长素、细胞分裂素(CTK)、脱落酸(ABA)、水杨酸(SA)和茉莉酸(JA)等通路，SA通路基因通过参与苯丙氨酸代谢过程正调控抗病性，而JA通路基因通过参与α-亚麻酸代谢负调控抗病性。Sánchez-Vallet等[47]和Lim等[48]报道ABA信号途径参与植物抗病免疫，French等[10]也证明生长素信号途径在番茄青枯病抗性上的重要性。

本研究其他次生代谢物和功能分子也能参与植物抗病过程。一些MYB可以通过调控靶基因的表达来调控植物体内次生代谢物质的合成，保护植物抵御生物胁迫，如AtMYB15调控拟南芥防御诱导的木质化和基础免疫[49]；AtMYB11组成表达增强了苯丙醇生物合成途径中基因的表达，导致烟草和番茄植株中类黄酮和绿原酸的积累，AtMYB11和AtMYB12调控番茄和烟草中黄酮类化合物和咖啡酰奎宁酸的生物合成[50]；AtMYB46、AtMYB58和AtMYB63在次生壁形成过程中能够激活木质素生物合成关键基因[51]。PpNAC1在苯丙氨酸生物合成和调控上发挥重要作用[52]。本研究筛选到的许多差异蛋白基因可能是构成抗性的重要组分，有必要对其功能和作用机制进一步解析。