油茶根系与内生细菌枯草芽孢杆菌互作早期的转录组分析*

李梓杨 陈晓琳 李丽丽 许诗萍 何苑皞

(1.南方人工林病虫害防治国家林业和草原局重点实验室 森林有害生物防控湖南省重点实验室 经济类培育与保护教育部重点实验室 中南林业科技大学 长沙 410004)

内生细菌广泛存在于健康植物体内，具有促进植物生长、提高作物产量、提高植物存活率、抑制病原菌生长、修复环境污染、溶解磷酸盐、促进氮吸收等作用(Bulgarellietal., 2012),通过长期的进化与寄主植物形成了互利互惠和的关系。内生细菌在互作早期，利用自身的趋化性、运动性、黏附作用以及分泌一些小分子物质与植物建立互作关系实现定殖并发挥作用；植物通过激素应答、免疫反应启动、次生代谢物合成、碳代谢等一系列的变化来响应内生细菌的定殖。植物根系或根系分泌物接触后内生菌转录组变化明显，其中氨基酸转运、抗压、IV型菌毛、粘附、植物来源化合物的利用等相关基因的表达变化显著(Shidoreetal., 2014; Yietal., 2014; Pankieviczetal., 2014; Sheibani-Tezerjietal., 2015)。植物根系在接触内生菌后，受体激酶、植物激素信号通路、木质素积累、细胞壁合成相关基因的表达变化明显(Drogueetal., 2014; Irizarry, 2018; Galambosetal., 2020)。目前研究主要以互作后期为主，而互作早期的研究较少，而互作早期植物各项生理变化较为显著，是互作的关键时间，以互作早期的植物根系为研究对象。目前，内生细菌与植物互作研究主要集中在拟南芥(Arabidopsisthaliana)、水稻(Oryzasativa)、棉花(Anemonevitifolia)等模式生物上，对于木本植物则关注较少，而内生细菌与油茶(Camelliaolefolia)互作的研究还未见报道。

笔者课题组前期从油茶内分离出1株内生细菌——枯草芽孢杆菌(Bacillussubtilis1-L-29gfpr)，该菌株可成功在油茶根系内定殖，并具有促进生长、抗病作用(Xuetal., 2020)。在本研究将该菌株接种油茶根系构建油茶-内生细菌活体互作体系，利用转录组测序技术研究不同时间(0、6、12、24 h)油茶根系转录组的变化，旨在解析油茶根系响应内生细菌侵染的变化规律，为后期深入研究内生细菌与油茶的互作机制提供参考。

1 材料与方法

1.1 试验材料

将油茶种子表面消毒后，放入无菌组培瓶中，17 ℃保湿培养，直至幼茎5 cm左右。内生细菌B.subtilis1-L-29gfpr为中南林业科技大学微生物菌种保藏室提供。

1.2 试验方法

1.2.1 油茶接种内生菌 将培养18 h的内生细菌1-L-29gfpr离心后重悬于PBS缓冲液中，浓度调至1×108cfu·mL-1，用注射器吸取20 mL缓冲液注射到组培瓶中，将油茶根系浸泡于缓冲液中， 20 min后将根系取出以无菌水冲洗3遍后置于无菌组培瓶中继续保湿培养。分别于0 h(不接种对照)、接种后6 h(6 hpi)、12 h(12 hpi)、24 h(24 hpi)截取油茶根系根尖部位。取样后将样品液氮速冻，于-80 ℃保存。每个处理3个生物学重复。

1.2.2 RNA提取、文库构建及测序 采用TRIZOL试剂盒(Invitrogen，美国) 分别提取0 、6 、12 、24 h共 12 个样品的总RNA，按照试剂盒说明书进行RNA提取。利用 Agilent 2100 Bioanalyzer (Agilent，美国)检测RNA质量，用带有 Oligo(dT) 的磁珠富集真核生物Mrna； 使用缓冲液将富集的mRNA片段化为短片段，并使用随机引物将其逆转录为cDNA； 通过DNA聚合酶I，RNase H，dNTP和缓冲液合成第二链cDNA； 用QiaQuick PCR提取试剂盒(Qiagen，荷兰)纯化cDNA片段，黏性末端修复，添加poly(A)，并连接测序接头。通过琼脂糖凝胶电泳选择连接产物的大小，PCR扩增，委托广州基迪奥生物科技有限公司使用Illumina Novaseq 6000进行测序。

1.2.3 差异表达基因 qRT-PCR 验证 测序获得的原始序列去除低质量序列得到Clean reads。使用对比软件bowtie2将测序数据与枯草芽孢杆菌ASM211716V1菌株的基因组进行mapping(Langmeadetal., 2012)，mapping后的枯草芽孢杆菌reads剔除掉，将所有油茶样本剩余的reads进行混合组装，得到油茶的参考序列信息。然后，使用序列聚类软件 TGICL 做进一步序列拼接和去冗余处理(Perteaetal., 2003)，得到Unigene。FPKM (fragments per kilobases per millionreads)值用于评估每个基因的表达水平。差异表达基因的筛选条件为P-adjust<0.05 且∣log2FC∣>=1。采用 NR、Swiss-Prot、Pfam、EggNOG、GO 和 KEGG 等数据库进行基因功能注释。通过 GO 和 KEGG 富集分析，确定差异表达基因主要富集的通路及其可能的生物学功能。

1.2.4 差异表达基因 qRT-PCR 验证 为了验证RNA-Seq 数据的可靠性，随机选取5个差异表达的基因进行 qRT-PCR表达量验证。利用软件Primer5.0设计引物，内参基因选用EF1α，引物列表见表1，随后采用 ABI 7000 系统进行 PCR扩增和分析。每个基因设 3 个重复。采用2-ΔΔCt法计算不同侵染时间点差异基因相对表达量。

表1 差异表达基因 qRT-PCR 验证使用的引物Tab.1 Primers used for qRT-PCR gene expression analysis

2 结果与分析

2.1 测序数据分析

油茶和内生细菌互作的文库0 (不接种对照)、6 、12 、24 h共产生序列52 958 922条，约46.1 Gb数据，差异表达基因12 561条，GC含量变化为43.14%～44.47%。过滤后Q30比例为91.94%以上。组装的Unigene序列N50为1 175，Unigene平均长度为866 bp，数据的比对率为71.13%。通过blastx将Unigene序列比对到蛋白数据库nr、SwissProt、KEGG 和 COG/KOG(evalue<0.000 01)发现有67 700获得了注释(44.5%)(表2)。对样品的生物学重复和样品组间的差异表达情况进行评估，相关系数均大于0.8，这表明重复样品间的相关性较强。

表2 内生细菌接种油茶不同时间点油茶转录组测序数据Tab.2 Summary of De novo data generated for C. oleifera transcriptomes inoculated with endophytic bacteria at different time points

2.2 差异表达基因分析

通过对转录组测序后的基因进行注释和表达量(RPKM 值)计算可知，3个时间点共有相同差异表达基因84个，内生菌接种油茶根部6 h时的差异表达基因数量最多为5 715个，其次为24 h，差异表达基因数量为1 871个(图1)。

图1 3个不同侵染时间点的油茶差异表达基因分析Fig. 1 DEGs of three different infection time points of C. oleiferaA： 不同时间点基因表达差异； B： 特异差异基因和共有差异基因的维恩图分析。 A： Differences in gene expression at different time points； B： Venn diagrams representing the numbers of DEGs and the overlaps of sets obtained across four comparisons.

2.3 差异表达基因的 GO 富集分析

以 correctedP-value ≤0.05 为筛选阈值，对不同时间获得的差异表达基因进行 GO 注释分析，这些差异表达基因被富集注释到3个 GO 分类，分别是分子功能(molecular function)、细胞组成(cellular component)、生物过程 (biological process)。差异表达基因生物学过程中有2 188 个差异表达基因；细胞组分有3 489 个； 分子功能有2 490 个。

由图2可知，油茶根系接种内生细菌6 h时生物学过程中磷酸代谢途径、含磷化合物代谢途径、响应激素途径、茉莉酸代谢途径、响应细菌途径等得到富集，这些途径均与植物响应细菌刺激有关； 分子功能中的萜烯合酶活性，双加氧酶活性等得到富集； 细胞组分中膜、膜固有成分等得到富集。接种12 h时,富集的生物学过程主要有单个有机体代谢过程、小分子代谢过程、羧酸代谢过程、酮酸代谢过程、有机酸代谢过程； 富集的分子功能主要有氧化还原酶活性、四吡咯绑定等； 富集的细胞组成有膜结合细胞器、胞内膜细胞器等。接种24 h时，代谢途径、单有机体代谢途径、小分子代谢途径等生物学过程得到富集； 富集的分子功能主要有催化活性； 细胞组分中膜结合细胞器、细胞等得到富集(图2)。

由GO富集分析中可知，油茶根系接种枯草芽孢杆菌1-L-29gfpr6 h后，生物学过程分子功能变化最明显。

图2 油茶-内生细菌互作早期差异表达基因GO富集分析Fig. 2 GO enrichment analysis of DEGs at early stage of C. olefolia-endophytic bacteria

2.4 差异表达基因的 KEGG 代谢通路富集分析

对不同时间点的差异表达基因进行KEGG pathways富集分析表明： 差异表达基因被富集到碳水化合物代谢、氨基酸代谢、信号转导、环境适应性等通路。

油茶根系接种内生细菌1-L-29gfpr6 h，植物-病原互作通路、植物激素信号转导、MAPK信号途径、芪类化合物、二芳基庚烷类化合物及姜辣素生物合成通路都有一定程度的提高，这些途径均参与植物对病原菌的响应。接种12 h，玉米素生物合成、核黄素代谢、双萜生物合成有一定程度的提高。接种24 h，酪氨酸代谢、硫代谢等有一定程度的提高。这些代谢途径可能参与植物与内生菌的互作(图3)。由此可见，互作过程中6 h 富集的代谢通路与12 h 和24 h的差别较大。

图3 油茶-内生细菌互作早期差异表达基因的 KEGG 富集分析Fig. 3 KEGG enrichment analysis of DEGs at early stage of C. olefolia-endophytic bacteria

2.5 植物激素信号转导相关的差异基因

调节植物激素水平是内生细菌重要的特点之一。在内生细菌与宿主植物互作过程中生长素、赤霉素、脱落酸、茉莉酸、水杨酸等植物激素信号转导相关基因均发生变化(Pinskietal.， 2019)。本研究中植物激素转导途径差异表达基因占总差异表达基因的比例随时间递减，6 、12 、24 h分别为6.41%、 3.96%、1.98%。由图4A和表3可知，生长素响应基因(SAUR15/32/50)、生长素早期响应基因(IAA)表达量在6 h或12 h都有一定程度的提高，表明内生细菌1-L-29gfpr在一定程度可以刺激油茶根系的生长。木葡聚糖内糖基转移酶/水解蛋白酶基因(XTH22/23/25)、TIFY家族基因(TIFY9/10B)、转录因子AP2/ERF(ERF1B)表达量在6 h显著上调随后下降。转录因子PIF(PIF3)、转录因子bZIP、DELLA蛋白编码基因(GAI)表达量呈现低-高-低的状态，在12 h最大。生长素响应蛋白编码基因(GH3.1)表达量则在24 h达到最大值。

2.6 植物病原互作相关差异基因

油茶根部接种内生细菌1-L-29gfpr后，植物病原互作通路检测到差异基因100个，其差异基因占总差异基因比例6、12、24 h分别为9.61%、3.96%、3.06%。油茶根部多个参与植物病原互作的基因均有上调，但到达峰值时间不同。由图4B和表4可知，钙调蛋白编码基因(CAM2)、呼吸爆发氧化酶同源基因(RBOH)以及RPM1互作蛋白编码基因(RIN4)表达量在6 h达到峰值，随后一直下降。转录因子MYB、MAPK激酶(NTF4)表达量呈现高-低-高的状态，在12 h达到最低值。抗病蛋白相关基因(RPM1、RPS5)表达量则在12 h达到最大值。由此可见，内生菌接种12 h后植物启动抗病蛋白相关基因，LysM类受体激酶基因(CERK1)则持续性的下调。

2.7 苯丙素生物合成相关差异基因

苯丙素参与植物对生物和非生物刺激的的反应(Vogt， 2010)。接种内生细菌后6 、12 、24 h，苯丙素生物合成途径差异基因占总差异基因比例分别为5.74%、3.96%、2.52%。在接种内生细菌后的不同时间点，检测到差异基因51个，6 h上调基因35个，下调基因16个； 12 h上调基因33个，下调基因18个； 24 h上调基因21个，下调基因30个。由图4C和表5可知，过氧化物酶基因(pod，PER5/44/58，PNC2)表达量在6 h达到峰值。糖苷水解酶(BGLU12、AA7GT、BGLU11)表达量在12 h达到峰值。肉桂醇脱氢酶(CAD1)表达量在24 h达到峰值。其中，大部分的基因都参与羟基苯基木质素、愈创木基木质素、紫丁香木基木质素以及香豆素的合成。

2.8 差异表达基因 qRT-PCR 验证

qRT-PCR 结果表明，随机挑取的 5个基因在不同时间点的差异表达量变化与转录组测序结果的趋势一致(图5)，说明本次转录组测序的结果可靠。

3 讨论

3.1 植物激素信号转导相关的差异基因

内生菌或病原菌的定殖都可引起植物体内激素的变化，植物响应内生菌或病原菌最为重要的是三类激素：水杨酸(SA)、茉莉酸(JA)和乙烯(ET)(Dominiketal., 2012)。本研究中植物激素转导途径差异基因占总差异基因比例从6.41%下降到1.98%，水杨酸、茉莉酸途径产生了变化，但乙烯未发生明显变化。在JA信号转导过程中，ZIM(TIFY)结构域可介导JAZ蛋白与其共抑制因子NINJA相互作用，进而共同抑制JA信号转导(段龙飞等， 2013)。本研究中JAZ编码基因TIFY9、TIFY108在接种6 h表达量上升，随后下降，这表明在内生菌与植物互作开始阶段JA信号转导受到抑制，12 h抑制解除。这与King等(2019)研究相似，接种内生菌后JA的响应会存在延时效应。TGA 基因家族是 bZIP 转录因子家族中十分重要的一组，它可与水杨酸信号途径关键调节因子 NPR1 相互作用，促进 TGA 对下游PR1启动子的结合并上调其表达，提高植株抗病性(Johnsonetal., 2003)。本研究中TGA表达量在接种6 h时下降，但24 h时PR1是有明显的上调表达。由此可见，内生菌接种后，油茶根部的水杨酸(SA)、茉莉酸(JA)开始时受到抑制，但随后抑制解除。这可能是由于内生菌为了定殖成功，在定殖初期(6 h)迫使植物降低了免疫反应，定殖成功后内生菌解除了免疫抑制，从而提高了植物的免疫力，这也证明内生菌可提高植物抗性。关于内生菌是如何迫使植物降低免疫反应还有待进一步研究。

图4 差异基因表达热图Fig. 4 Heat map diagram of the expression levels of DEGsA： 植物激素信号转导相关差异基因表达热图； B： 植物病原互作相关差异基因表达热图； C： 苯丙素生物合成相关差异基因表达热图。A:Heat map diagram of the expression levels of DEGs which are involved in plant hormone signal transduction. B:Heat map diagram of the expression levels of DEGs which are involved in plant-pathogen interaction. C:Heat map diagram of theexpression levels of DEGs which are involved in plant phenylpropanoid biosynthesis.

图5 qRT-PCR表达量分析和转录组数据的验证Fig. 5 Expression analysis of qRT-PCR and verification of transcriptome dataA： TCH4基因； B:RBOHC基因； C:POD基因； D:TIFY基因； E:XTH基因。 A： TCH4 gene； B:RBOHC gene； C:POD gene； D:TIFY gene； E:XTH gene.

表3 植物激素信号转导通路差异基因Tab.3 DEGs of plant hormone signal transduction pathway

内生细菌通过调控植物激素信号传递，起到促进植物生长的作用(Straubetal., 2013)。赤霉素(GA)在植物株高等生长发育过程中起着重要的调控作用(张佳琦等， 2019)。GA 信号通路中的关键组分 DELLA 蛋白，主要起阻遏作用，GA可以解除DELLA 蛋白的阻遏作用。本研究中DELLA 蛋白基因表达量最后降低，表明内生菌定殖后期可以促进油茶生长。同时本研究也发现，DELLA 蛋白基因表达量与JA的相关基因表达量呈现相反的变化趋势，这与Yimer等(2018)的研究结果一致，即GA和JA之间存在对立关系。推测其原因可能是植物将DELLA 蛋白作为调节生长与防御之间平衡的手段之一。

本研究中，接种内生细菌6 h后油茶根系生长素Aux/IAA基因表达量明显上升，但随后开始出现下降，这与Irizarry等(2018)研究结果相反，这可能是由于本研究选取的时间点为侵染前期，而Irizarry的研究选取的是接种10天后的样品。

表4 植物病原互作通路差异基因Tab.4 DEGs of plant-pathogen interaction pathway

表5 苯丙素生物合成通路差异基因Tab.5 DEGs of phenylpropanoid biosynthesis pathway

3.2 植物病原互作相关差异基因

植物与微生物共同进化过程中形成了复杂的免疫防卫体系。该免疫防卫体系的正常运行保证了植物的正常生长。由于植物细胞可以通过先天免疫反应对MAMPs进行监测及响应(Machoetal., 2014)，因此内生菌或病原菌的定殖都可引起植物的免疫反应。Drogue等(2014)研究发现,内生菌可以抑制更多与植物防御有关基因的表达。在本研究中，内生菌在6 h引起了植物病原互作通路大量基因的表达，但在12 h表达量明显下降，表明在内生菌侵染的初期油茶起动了防御反应，到后期防御反应减弱，让内生菌更好的实现定殖。

转录因子在调节植物与微生物互作相关基因的过程中起着重要的作用。已有研究证明WRKY转录因子参与大量的植物-病原菌互作反应(Banerjeeetal., 2015; Ishihamaetal., 2012; Pandeyetal., 2009; Ryuetal., 2006)，但是对于内生菌与植物相互作用中的WRKY转录因子知之甚少。本研究中油茶根系的WRKY22/24/27/29/33都在6 h表达量上调，随后出现一定的下调。Jorge对有益真菌与拟南芥互作过程中WRKY基因的表达进行了研究，发现不同的转录因子调节的途径不同，其中WRKY33是受茉莉酸调节，并且在24 h的表达量是下降的(Saenz-Mataetal., 2014)。

钙是植物中PAMP触发的免疫力重要信号(Yuanetal., 2017)。植物环核苷酸门控通道(CNGC)蛋白介导的钙离子内流是识别受体复合物与钙离子依赖性免疫的一个重要关联。Tian等(2019)研究表明，CNGC通道直接被BIK1磷酸化并响应flg22激活PTI信号传导途径。本研究中，内生细菌枯草芽孢杆菌也存在flg22，CNGC编码基因表达量上调可能与此有关。同时本研究中油茶根系的钙调蛋白基因(CAM2)在6 h上调明显，主要是由于钙调蛋白参与了植物免疫反应(Singhetal., 2020)。

3.3 苯丙素生物合成相关差异基因

植物的苯丙素化合物在抵抗病原体侵袭中起重要作用。苯丙素合成途径中的木质素可以使细胞壁增厚，抵抗真菌侵染钉形成的压力，同时木质化的细胞壁可以减弱病原菌细胞壁降解酶的作用(Bechingeretal., 1999; Naoumkinaetal., 2010)。苯丙素合成途径在病原菌与植物互作过程中变化较大，但内生菌与植物互作中的研究报道较少。本研究中3个主要单元： 对羟基苯基(H)、愈创木脂(G)和丁香基(S)的基因表达量在接种6 h和12 h均有上调，而在24 h下调。王恒照等(2019)对促生菌B.pumilusLZP02接种后的水稻根部的研究发现，LZP02可提高苯丙素合成相关基因的表达。这一结果与本研究存在一定的差异，其主要原因是本研究为12 h与24 h的表达量比较，如果将24 h与未接种内生菌的试验组比较，苯丙素合成基因的表达量仍为上升，同时Hong的取样时间为促生菌接种48 h后，也与本研究存在一定不同。

豆香素是苯丙素合成途径中另一类非常重要的植物次生代谢产物。香豆素是木质素合成的前体，并且香豆素本身就是一种植物抗毒素，可以防止病虫的侵害，同时豆香素还对植物根部微生物组的装配起着重要的作用(Stringlisetal., 2018)。本研究中香豆素合成基因表达量也出现了相应的变化，在12 h达到峰值，这表明油茶内生细菌可通过调节宿主信号传导，从而提高宿主对病原菌的抗性。

4 结论

对内生细菌接种后不同时间点的油茶根部转录组进行测序分析，得到12 561个差异表达基因。差异表达基因在6 h最多，主要集中在植物-病原互作通路、植物激素信号转导、苯丙素化合物等通路。其中htpG、钙调蛋白基因(CAM2)、呼吸爆发氧化酶编码基因(RBOH)、转录因子WRKY22、钙调蛋白基因(CAM2)、木糖转移酶基因(TCH4)、茉莉酸受体ZIM结构域包含蛋白编码基因(TIFY10B)、Aux/IAA基因IAA1、过氧化物酶编码基因等上调显著。接种后6 h差异表达基因最多，为互作的关键时间。