短短芽胞杆菌B011基因组中芽胞形成调控基因spo0E的鉴定与表达分析

周向平，王运生，陈 武，刘天波，王凯歌，刘 峰，周 立，袁志辉

（1.湖南农业大学植物保护学院，长沙 410128；2.湖南省烟草公司永州市公司，永州 425100；3.湖南省烟草科学研究所，长沙 410004；4.湖南科技学院化学与生物工程学院，永州 425199；5.湖南省银杏工程技术研究中心，永州 425199）

形成芽胞是产芽胞类细菌在长期的自然进化过程中形成的抵抗逆境的能力体现。当芽胞类细菌感知到逆境胁迫时，如细胞密度过高、营养物质匮乏等，细胞首先一般可通过三种途径来减轻环境压力，一是激活运动装置，如鞭毛蛋白，通过趋化性寻找新的食物；二是分泌抗生素、抗菌肽和其他化学物质，清除或驱除同一生境中的竞争性微生物；三是分泌大量水解酶类，降解胞外蛋白质和糖类[1]。对枯草芽胞杆菌Bacillussubtilis168 的研究显示，如果上述3 种途径没有减轻环境压力，50%～70%的细胞将形成芽胞[2]。芽胞形成属于群体感应行为，即细胞向胞外分泌信号分子，细胞能感应这些信号分子的浓度（与细胞密度成正比），当达到一定阈值后，通过信号转导通路激活或抑制相关基因的表达，使细胞进入不同发育途径[3]。

spo0A是调控芽胞形成的开关基因[4]，调控因子Spo0A 磷酸化为Spo0A～P，Spo0A～P 是Spo0A 的活性形态，是芽胞形成的关键调控因子[5]，一方面，Spo0A～P 可以正调控spoIIE基因，正式启动芽胞形成途径[6]，另一方面Spo0A～P 负调控abrB基因[7]，AbrB 是一个全局性转录负调控因子，同时Spo0A～P 还可以正调控abbA基因，而AbbA 可抑制AbrB 活性[8]，AbrB 可负调控群体感应关键基因comK[9]。细胞内Spo0A～P 的浓度水平决定了细胞的发育途径，Spo0E 是一类磷酸酶，可将活性态Spo0A 即Spo0A～P，转为非活性态Spo0A，从而起到负调控胞内Spo0A～P 的浓度[10,11]，同时spo0E受AbrB 负调控，Spo0A-AbrB-Spo0E 形成一个调控回路[12]。Spo0E在芽胞杆菌基因组中多以基因家族（SMFP,spo0Emultigene family of phosphatases）存在，基因上游一般存在多个启动子（promoter）区域，可受不同转录因子调控[13]。

短短芽胞杆菌BrevibacillusbrevisB011 菌株是本实验室从烟草根系筛选到的内生细菌，其对烟草青枯病菌、烟草黑胫病菌、辣椒白绢病菌等均有较明显的抑制效果，通过形态鉴定、生理生化试验、16S rDNA测序及比对分析，确定该生防菌株为短短芽胞杆菌[14]。本文在对短短芽胞杆菌B011 菌株进行全基因组测序的基础上，对不同生长阶段的B011 菌株进行转录组测序，通过生物信息学工具预测得到了9 个spo0E同源基因，并通过与其他101 个短芽胞杆菌Brevibacillus菌株全基因组进行比较分析，结合转录组数据，推测spo0E同源基因功能，为进一步验证基因功能和B011 菌株的遗传改造提供候选基因。

1 材料与方法

1.1 短芽胞杆菌基因组序列分析

短短芽胞杆菌B011 菌株为本实验室从土壤环境中分离获得的菌株，经鉴定其为短短芽胞杆菌[14]。对其进行了全基因组测序，并提交序列至NCBI，GenBank 登录号为CP041767。从NCBI RefSeq 数据库[15]（https://www.ncbi.nlm.nih.gov/）下载所有短芽胞杆菌已经公布全基因组序列的菌株核苷酸序列和蛋白序列，统计各个菌株基因组的GC 含量、蛋白编码基因个数、基因组大小等信息。

1.2 短短芽胞杆菌基因家族分析

用OrthoMCL[16]构建基因家族，选取已公布完成图的4个短短芽胞杆菌菌株（B011、X23、NBRC 100599、DZQ7）用来做基因家族分析，OrthoMCL 首先通过blastp 对所有蛋白序列进行两两比对，将任何两两比对e-value 小于1E-5 的序列保留下来，然后根据score 值进行MCL 聚类，聚在同一个cluster 的基因被认为是同一个基因家族。

1.3 Spo0E 家族基因预测

从Pfam 数据库[17]（https://pfam.xfam.org）下载SpoOE-like (PF09388)模型文件，然后用hmmsearch[18]搜索上述Brevibacillus蛋白库，取E 值小于0.05 的蛋白为候选基因，然后经BLASTp[19]进一步验证NR 数据库中相似序列是否含有Spo0E功能描述。

1.4 多序列比较和构建进化树

用MAFFT[20]进行多序列比对，参数选用自动模式（-auto），然后用FastTree[21]构建系统进化树，参数用默认参数，进化树上传到在线分析平台iTOL[22]进行进一步分析与显示作图。

1.5 基序（motif）分析

基序采用MEME Suite 软件包[23]进行分析，选取输出前5 个基序。MEME Suite 是集合众多预测和注释motif 工具的在线网站，其中MEME 算法是基于最大期望值（EM）算法来识别motif。保守基序用seqLogo[24]作图。

1.6 B011 菌株中spo0E 类似基因的表达情况分析

B011 菌株在LB 培养液中，37 ℃、200 r/min 条件下分别振荡培养12、18、30、36 h 后，送交测序公司进行转录组测序。芽胞形成调控相关基因相对表达量采用RNA-Seq 的TPM（Transcripts Per Million）值进行表征，TPM 与RPKM/FPKM 类似，是将原始reads 除去基因长度的影响后标准化到每百万条reads，相当于重新标准化的文库，保证每个样本中所有TPM 的总和是相同的[25]。RNA-Seq 分析流程可简单描述为：首先将质控合格的reads 回帖到参数基因组，然后用TPMCalculator[26]计算每个基因的TPM 值。

2 结果与分析

2.1 短芽胞杆菌基因组序列及特征

截止2020 年8 月1 日，在NCBI 的RefSeq 数据库中已收录并公布全基因组序列的短芽胞杆菌菌株有101 个，加上本实验室的B011 菌株，一共包含了102 个菌株的全基因组序列。包含有20 个种，其中侧孢短芽胞杆菌B.laterosporus23 个，短短芽胞杆菌13 个，波茨坦短芽胞杆菌B.borstelensis7 个，土壤短芽胞杆菌B.agri5 个，有26 个菌株没有种分类信息。短芽胞杆菌基因组大小范围为4.00～7.13 Mb，GC 含量范围为40%～8.5%。从图1 可以看出，不同种之间基因组大小和GC 含量差异明显，如侧孢短芽胞杆菌基因组GC 含量偏低（40%～41.3%），基因组大小也偏小（4.53～5.80 Mb），而短短芽胞杆菌基因组GC含量偏高，为46.8%～47.5%，基因组大小也偏大，为6.02～6.73 Mb，在基因组大小-GC 含量图上能明显分开（图1）。而短短芽胞杆菌，抗生素短芽胞杆菌B.antibioticus, 美丽短芽胞杆菌B.formosus, 强音短芽胞杆菌B.fortis和波塔短芽胞杆菌B.porteri聚在一起，这与平均核苷酸一致性（ANI, Average Nucleotide Identity）和基于全基因组聚类结果一致（结果未公布），表明这几个种亲缘关系较近。

图1 102 个Brevibacillus 菌株的聚类分布图Fig.1 The cluster distribution of 102 Brevibacillus genomes using the genome size as X-axis and GC content as Y-axis

2.2 短短芽胞杆菌基因家族分析

比较了4 个短短芽胞杆菌基因组，用OrthoMCL 对4 个短短芽胞杆菌的编码基因进行了聚类，构建基因家族。首先对4 个菌株的蛋白序列用BLAST 进行两两比对，然后用MCL 进行聚类，共构建了5933 个基因家族，其中4749 个基因家族在4 个菌株都存在，可视为短短芽胞杆菌的核心基因组，其中4718 个基因家族为单拷贝基因家族，每个基因组只含一份（表1，图2）。

表1 4 个短短芽胞杆菌基因家族概况Table 1 The gene families of the 4 B.brevis with complete genomic sequences

图2 短短芽胞杆菌基因家族比较分析Fig.2 OrthoMCL gene families of 4 Brevicillus brevis

2.3 Spo0E 基因预测

从102 个短芽胞杆菌菌株基因组中一共预测得到837 个spo0E基因，平均每个菌株含8.2 个spo0E类似基因，最少的仅有3 个，最多的有23 个。不同种之间spo0E同源基因数量具有一定差异，侧孢短芽胞杆菌平均含有13 个，变化范围为8～23，而短短芽胞杆菌平均含有8 个spo0E基因，变化范围为6～10 个（图3）。Spo0E 蛋白因子的平均长度为68 个氨基酸，最小的24 个氨基酸，最长的为162 个氨基酸，但不同种之间长度差异不显著。菌株B011 预测含有9 个spo0E类似基因，编码肽链的长度范围为52～101个氨基酸。这9 个基因在基因组上为分散分布，多数与旁边基因并不形成操纵子结构。

图3 不同Brevibacillus 种spo0E 类似基因个数Fig.3 The number of spo0E-like genes in different species of Brevibacillus

2.4 系统进化分析

从837 个spo0E类似基因的系统进化树上可以看出，这些基因明显聚成6 个大的分枝，B011 菌株的9个基因分散在不同的5 个分枝上，且每个分枝上均分布有不同的种，说明这些基因在种分化之前就已经存在了。与其他种不同的是，侧孢短芽胞杆菌趋向于分布在单独的分枝上，特别是如图4 所示的分枝4 绝大部分为侧孢短芽胞杆菌成员，说明该分枝上的基因可能为侧孢短芽胞杆菌特异基因。

图4 短杆芽胞杆菌spo0E 类似基因系统进化树Fig.4 The phylogenetic tree of 837 spo0E-like genes from 102 Brevibacillus species

2.5 短杆芽胞杆菌Spo0E 编码蛋白基序特征

经MEME 的基序分析，发现短短杆芽胞杆菌的837 个spo0E基因编码蛋白序列中有825（98.6%）含有基序1，基序1 当中最保守的序列为SQQLD，这与已报道的芽胞杆菌Spo0E 的SQELD 基序高度保守，特别是基序中的第8 位S 和第12 位D 在所有序列中高度保守，推测可能为其活性位点。菌株B011 的9个基因中，有5 个为SQQLD，2 个为SQELD，1 个为SVQLD，1 个为SIQID（图5）。基序2 位于基序1的前面，但缺乏完全保守的位点，在673 个基因中存在该基序。这两个基序均存在于B011 菌株的9 个基因中。

图5 短杆芽胞杆菌Spo0E 蛋白因子基序特征及B011 菌株的9 个Spo0E 类似蛋白序列多序列比对Fig.5 The main two motifs in the Spo0E-like proteins of Brevibacillus and multiple sequences alignment of that of strain B011

2.6 B011 菌株中spo0E 类似基因的表达情况

B011 菌株的9 个spo0E同源基因除了FO446_11425在5 个转录组样本中没有检测到外，其他8 个基因均检测到有表达（图6），其中基因FO446_19435的表达量非常高，在发酵到18 h 达到最高峰，18 h为其对数生长期，其TPM 值达到了5013，随后下降。基因FO446_04050表达模式与FO446_19435类似，均在18 h 达到最大值，表达量也较高，18 h 时TPM 值达到了1003。而FO446_23370、FO446_27395两个基因表达模式相似，但在18 h 时表达量与30、36 h 时的测量值相当，属于中等表达丰度基因，另外4 个基因的表达量较低，其中基因FO446_01695在12 h 表达量相对更高，随后下降。

图6 菌株B011 菌株中9 个spo0E 类似基因表达热图Fig.6 The gene expression heatmap of 9 spo0E-like genes in strain B011

3 讨论

短短芽胞杆菌是能够产生多种植物病原菌拮抗物质的重要生防菌，其活性物质鉴定和基因组测序分析已有多个报道[27,28]。从已报道的具有完成图的短短芽胞杆菌基因组来看，B011 菌株的基因组大小是最小的，蛋白编码基因也是最少的，预测有5663 个蛋白编码基因，而本实验室筛选到的另一个生防菌X23 菌株同为短短芽胞杆菌[29]，基因组大小为6643473，编码6367 个蛋白，比B011 菌株基因组大了约500000，蛋白编码基因也相应多了约700 个基因。B011 菌株是从植物组织内分离得到的植物内生菌，而X23 菌株是从土壤中分离得到的土栖菌，植物内生菌由于所处的生存环境相对于土壤来说要稳定和简单，对于植物内生菌来说，某些基因可能对其生存作用不大，这部分基因在进化的过程中就有可能丢掉或者退化，这也与植物内生菌或者病原微生物的基因组相对于土栖菌来说，其基因组相对来说要小一些的现象是一致的。

在芽胞形成调控基因方面，B011 菌株中含有9 个spo0E家族基因，除了FO446_11425在5 个转录组样本中没有表达外，其余8 个基因均有表达，说明这些基因可能都在行使功能。Spo0E基因在产芽胞细菌中一般是以多基因家族存在的，相互之间存在功能冗余，如在枯草芽胞杆菌中敲除spo0E后对其芽胞形成并没有明显影响，因为在枯草芽胞杆菌中还存在另外两个spo0E同源基因yisI和ynzD，具有类似功能[30]。在所有B011 的9 个spo0E同源基因中均存在芽胞杆菌SQELD 类似基序，说明SQELD 基序可能是芽胞杆菌共有基序，三维结构显示，基序SQELD 位于活性位点中心[31]，其中S 和D 高度保守，中间三个氨基酸QEL 在不同菌株间存在一定变异，如在B011 菌株中有5 个为SQQLD，但表达量最高的两个spo0E同源基因FO446_19435和FO446_04050均为SQELD。B011 菌株中这9 个spo0E同源基因是否具有功能互补？在体内又是如何进行调控的需进一步研究。

鉴于基因FO446_19435表达量非常高，且表达模式是先升高后下降，符合spo0E的一般表达模式，Spo0E 的生物学功能是对Spo0A～P 进行去磷酸化，从而降低细菌体内Spo0A～P 的含量，抑制芽胞的形成。推测基因FO446_19435为B011 菌株的主效spo0E基因，但同时基因FO446_04050可能是spo0E的候选补充基因，它们的表达模式类似，且表达量也较高。FO446_19435与下游基因FO446_19440属于同一个操纵子，与FO446_19435类似，基因FO446_19440的表达量非常高，18 h 的TPM 达到了27181，但基因FO446_19440功能未知，且只在菌株短短芽胞杆菌GZDF3.1 中找到有同源基因，但奇怪的是与FO446_19440同为一个操纵子的FO446_19435并未在短短芽胞杆菌 GZDF3.1 找到同源基因。

枯草芽胞杆菌中spo0E受到多个转录因子调控，受AbrB 和Rok 负调控，同时还受sigma B 转录因子正调控[32]。在B011 菌株中，abrB(FO446_00625)基因表达随着发酵时间是升高的，发酵后期表达量明显上调，总体趋势是与spo0E的表达相反，说明与枯草芽胞杆菌类似，B011 菌株的spo0E可能也受AbrB 负调控，而sigma B (FO446_02910)转录因子的表达量变化不大，sigma B 转录因子主要与胁迫有关[33]，由于本试验转录组样本没有胁迫因素，可能是导致sigmaB转录因子表达量变化不大的原因。B011 菌株中缺少rok基因，rok基因主要与感受态细胞形成相关[34]，而B011 菌株中同时也缺少感受态细胞形成关键基因comK。考虑到B011 菌株的spo0E基因在发酵18 h 时达到最高，而abrB的表达量是一直上升的，所以推测在B011菌株中，spo0E基因除了受AbrB 调控外，可能还受到其他基因调控。