屎肠球菌自溶的基因表达差异分析

聂俊辉， 郭 浩， 王 通， 郭建军， 曾 静， 袁 林*

（1.江西省科学院微生物研究所，江西南昌 330096；2.上饶市农林水科学研究中心，江西上饶 334000；3.江西鑫维生物技术有限公司，江西南昌 330096）

屎肠球菌（Enterococcus faecium）为革兰氏阳性菌，属于肠球菌属，是一种兼性厌氧菌，广泛存在于传统乳制品、发酵食品、人类和动物的肠道当中（景宇超等，2018；Kivan，c等，2016；王晓蕊等，2016；Hassanzadazar等，2014）。屎肠球菌能产生诸如细菌素、有机酸、消化酶等多种抗菌物质，能够抑制病原微生物繁殖，平衡肠道的微生态环境（罗强等，2021；欧秀玲等，2021；田丰松等，2016；王永等，2013）。由屎肠球菌制成的饲用微生态制剂能够调节动物的肠道菌群平衡，减少腹泻等多种消化道系统疾病的发生率，提高动物对饲料的利用率，对提高动物的养殖性能有着良好的促进作用（王卫卫等，2021；肖冉等，2020；夏玉等，2014）。

关于屎肠球菌活菌制剂的制备，影响其生产能力的因素有多种。改善外界的环境因素，如调节培养基配方、优化发酵条件等，能够有效提高屎肠球菌的生物量和存活率（陈志娜等，2021；刘斯斯等，2021）。同时，菌体自身的一些生理特点同样能够会对菌剂的制备造成影响（胡卫国等，2016）。屎肠球菌发酵过程中易自溶的特点限制了其产能的提升，制约了其大规模的制备。因此，通过对屎肠球菌自溶这一生理现象的研究，能够为菌株的高性能生产提供理论依据和改进方向。

RNA测序（RNA-Seq）也称为转录组测序，该技术是转录组学研究的重要手段。通过对研究对象某一特定状态下转录出来的几乎所有RNA进行测序，获得该状态下研究对象的基因表达水平，能够从基因的整体水平反映基因功能及生物学过程（宋尚桥等，2020；Stark等，2019）。目前在对细菌生理及功能的研究中也广泛运用了高通量测序技术。郭金凤等（2020）使用转录组技术分析了副干酪乳杆菌在乙醇胁迫下的基因表达情况，得到了6个可能与抵制乙醇对自身损害相关的基因，为阐明乳杆菌的耐乙醇机制提供了一定的理论依据。宋雪飞等（2017）使用转录组技术分析了植物乳杆菌在盐胁迫下的基因表达变化，为提高植物乳杆菌在工业中生产中的耐受性积累了科学依据。

鉴于屎肠球菌自溶机制在屎肠球菌制备中的重要性，本研究对自溶状态下的屎肠球菌进行转录组测序，对比正常状态下的菌体，从基因表达层面探究屎肠球菌自溶状态下基因表达的情况，分析屎肠球菌自溶的分子机制及在该过程中的关键通路，为屎肠球菌的培养优化、基因工程改良等奠定理论基础。

1 材料与方法

1.1 试验材料

1.1.1 菌株 屎肠球菌EXW275由江西省科学院微生物研究所分离保存。

1.1.2 培养基 MRS培养基（郭建军等，2019）：葡萄糖20.0 g、胰化蛋白胨10.0 g、酵母提取物4.0 g、牛肉膏8.0 g、柠檬酸铵2.0 g、七水磷酸氢二钾2.0 g、三水醋酸钠4 g、七水硫酸镁0.6 g、一水硫酸锰0.25 g、Tween 80 1.0 mL，蒸馏水定容至1000 mL，调节pH至6.2，115℃灭菌20 min。需要时另加入15 g琼脂制成固体培养基。

1.2 试验方法

1.2.1 菌体培养 菌种活化及培养（郭建军等，2019）：将保存的屎肠球菌菌种接种至MRS斜面培养基上活化，37℃恒温培养24 h，连续转接活化两次；再用接种环刮取两环接种至液体MRS培养基中，37℃、120 r/min培养10 h作为种子菌液。在培养瓶中装入发酵培养基后灭菌，无菌条件下向培养瓶内接入种子菌液；初始发酵培养条件为：250 mL培养瓶装100 mL培养基，接种量5%，初始pH 6.2，37℃，150 r/min，摇瓶培养20 h。取此时的菌体作为正常对照组。菌液继续培养传代至自溶发生，取此时菌体作为自溶组。

1.2.2 RNA提取 取菌液离心后去除上清液，保留菌体。使用TRIzol Reagent（Invitrogen）对菌体进行裂解并提取总RNA，操作参照产品说明书。提取总RNA后加入适量DEPC水溶解，使用Nanodrop 2000检测总RNA的浓度及纯度。将合格的总RNA委托上海生工生物工程股份有限公司使用Illumina HiseqTM平台进行转录组测序。

1.2.3 测序数据处理及分析

1.2.3.1 数据评估及质控 将测序得到的原始图像数据文件经CASAVA碱基识别（Base Calling）分析转化为原始测序序列。使用FastQC对原始数据进行质量评估。使用Trimmomatic进行质量剪切，去掉含有带接头的、低质量的序列，得到相对准确的Clean数据。

1.2.3.2 RNAseq测序评估 使用Bowtie2将样本有效数据比对到参考基因组上，统计Mapping信息，并通过RSeQC根据比对结果进行冗余序列分析、插入片段分布等分析。采用Qualimap根据比对结果进行均一性分布检查、基因组结构分布等分析。使用BEDTools进行基因覆盖率统计分析和基因组上测序序列分布等分析。

1.2.3.3 表达水平分析 使用featureCounts和已知的基因模型评估基因的表达量；使用WGCNA进行基因共表达分析；基于样本的表达量矩阵进行样本比较分析等多方向的统计分析和探索。

1.2.3.4 表达差异分析 使用DESeq2进行基因表达差异分析，对表达差异分析结果进行可视化。基于差异分析结果，进行聚类分析；使用topGO进行GO富集分析，并绘制显著性GO有向无环图；使用clusterProfiler进行KEGG通路和COG分类富集分析。

2 结果与分析

2.1 样本测序数据质量控制及序列比对 通过上海生工生物工程有限公司的测序平台分别对屎肠球菌正常状态和自溶状态下的样本进行转录组测序后，共获得89266999个raw reads。经过质量剪切，去掉含有带接头的、低质量的序列，得到相对准确的clean reads。总共获得正常状态下44674190个clean reads，平均reads长度为138.3 bp，GC含量为42.51%，Q10、Q20、Q30 Bases Ratio分别为99.91%、98.33%、93.93%；自溶状态下41888670个clean reads，平均reads长度为137.9 bp，GC含量为42.32%，Q10、Q20、Q30 Bases Ratio分别为99.90%、98.15%、93.30%（表1）。以上结果表明当前两组转录组数据量较为丰富，质量较高。

表1 样本测序数据质量控制及序列比对结果

将测序序列进行基因组定位后分析，并通过RSeQC统计比对结果。结果显示屎肠球菌正常状态下序列能定位到基因组上的占比总计（Total mapped）为96.17%，在参考序列上有多个比对位置的测序序列的数量统计（Multiple mapped）为21.45%，在参考序列上有唯一比对位置的测序序列的数量统计（Unique mapped）为74.72%，整段比对到外显子的测序序列（Non-splice reads）统计为74.72%；屎肠球菌自溶状态下序列能定位到基因组上的占比总计为64.42%，在参考序列上有多个比对位置的测序序列的数量统计为21.30%，在参考序列上有唯一比对位置的测序序列的数量统计为43.11%，整段比对到外显子的测序序列统计为43.11%（表2）。

表2 参考序列比对结果

2.2 基因表达量分析及差异基因筛选 基因表达水平的直接体现就是其转录本的丰度情况，丰度越高，则基因表达水平越高。通过定位到基因组外显子区的测序序列计数来评估基因的表达水平，使用TPM（Transcripts per million）来计量在RNA池中的转录比例。通过对基因表达量密度分析发现，屎肠球菌正常和自溶状态下基因表达的密度出现明显区隔，自溶状态下的基因表达出现了改变（图1）。

图1 基因表达量密度曲线图

通过对两组样本进行主成分分析（PCA）和主坐标分析（PCoA）发现，PC1、PC2、PC3分别为72.53%、10.79%、8.43%，PCo1、PCo2、PCo3分别为99.47%、0.35%、0.12%（图2）。以上结果说明两组样本间组内重复性良好，组间差异明显，组间差异显著大于组内样本间差异，屎肠球菌自溶状态下基因表达出现明显改变。

图2 PCA分析与PCoA分析

经过对两组样本差异分析后，共筛选出差异表达基因783个，以屎肠球菌正常状态下的转录基因作为对照组，自溶状态下的屎肠球菌基因表达上调个数为689，下调个数为94，上调基因数远远大于下调基因数（图3）。

图3 屎肠球菌自溶状态下差异表达基因展示图

2.3 差异基因的功能注释 GO（Gene Ontology）是一个国际标准化的基因功能分类体系，提供了一套动态更新的标准词汇表来全面描述生物体中基因和基因产物的属性（Carbon等，2019）。GO总共有三个ontology，分别描述基因的分子功能（MF）、细胞组分（CC）、参与的生物过程（BP）（Blake等，2015）。将差异表达的基因进行GO功能统计分析，在三大功能分类中，分别包含了22（BP）、17（CC）、12（MF）个功能分类（图4）。在生物过程中，差异基因分类在代谢过程、细胞过程、定位、生物调节的数量最多，分别为403、400、93、93个，分别占分类注释基因的34.8%、35.75%、26.65%、37.80%；在细胞组成中，差异基因分类在细胞、细胞部分、胞膜、胞膜部分的数量最多，分别为434、431、213、187个，分别占分类注释基因的39.10%、39.18%、26.89%、25.69%；在分子功能中，差异基因分类在催化活性、结合、转运活性、结构分子活性的数量最多，分别为424、401、76、56个，分别占分类注释基因的32.07%、35.02%、28.57%、96.43%（图4）。

图4 差异基因功能注释

对屎肠球菌自溶下的差异表达基因进行显著性分析，富集程度最高的30个GO分类如图5所示，富集程度最高的几个分类为细胞、细胞部分、生物合成、有机物质生物合成几个方面。这说明屎肠球菌自溶状态下的胞内生物合成途径出现显著改变，与自溶的发生具有相关性。

图5 差异基因功能富集分析

2.4 差异基因的KEGG分析 京都基因与基因组百科全书（KEGG）是基因组破译方面的数据库，在给出染色体中一套完整基因的情况下，其可以对蛋白质交互（互动）网络在各种细胞活动起的作用作出预测（Kanehisa等，2016；Kanehisa等，2000）。通过对屎肠球菌自溶下的差异表达基因进行KEGG富集分析发现，屎肠球菌自溶状态下在核糖体、氨基酰基-tRNA生物合成、嘌呤代谢、肽聚糖生物合成、嘧啶代谢、同源重组、氨基酸合成等7个途径具有显著性（P<0.05）（表3）。在对自溶下表达下调的基因进行聚类后则发现，下调的基因主要富集在碳代谢、糖酵解/糖异生、丙酮酸代谢等能量代谢相关通路（图6）。

表3 差异表达基因KEGG显著性富集分析

图6 下调差异表达基因KEGG显著性富集分析

3 讨论

自溶是指在一定条件下菌体自身释放被称为自溶素（Autolysin）的肽聚糖水解酶（PGHs），水解细胞壁肽聚糖的网络结构导致其裂解并向周围环境释放胞内物质的过程，是细菌程序性死亡的一种常见方式（Lewis 等，2000；Shockman 等 ，1996；Hltje等，1995）。影响细菌自溶的因素复杂多样，培养基的成分、温度、补料速率、接种量等均可能对自溶的发生产生影响（Vegarud等，1983），如培养基中碳源耗尽时会导致S.thermophilus和L.lactis的自溶，环境中NaCl浓度、乙醇含量或热休克等可以导致S.thermophilus的自溶（Husson-Kao等，1999；Riefe等，1997），过高的补料速率和接种量导致大肠杆菌自溶的发生（饶犇等，2020）。尽管过去对细菌自溶的研究有很多，但是关于屎肠球菌自溶的原因及基因表达规律尚不明确。深入了解屎肠球菌自溶下基因的表达情况，对于挖掘探索自溶有关基因，分析导致屎肠球菌自溶的环境因素具有重要的意义。

细菌的自溶受到多种因素的影响，不同的菌种对环境变化的敏感性不同，涉及的机制复杂而多变。屎肠球菌自溶的分子机制尚不完全清楚，遗传基础了解相对较少。不同于过往从自溶酶等直接因素入手进行乳酸菌自溶的研究，本研究选用屎肠球菌正常状态下的菌体以及自溶发生时的菌体进行转录组测序，通过基因表达水平及聚类分析比较，从基因转录表达的角度探索自溶发生的机制。基于RNA-Seq的转录组测序技术结合生物信息学分析为全面了解屎肠球菌自溶发生的基因表达改变提供了技术支撑，为屎肠球菌自溶的分子机制提供了研究基础。

本研究将所有的差异基因进行GO功能聚类分析，探索自溶下屎肠球菌的功能及代谢改变。差异基因富集程度最高的几个分类为细胞、细胞组分、生物合成、有机物质生物合成，揭示了屎肠球菌自溶状态下菌体的细胞状态改变以及菌体自身合成代谢的改变。KEGG通路分析则发现，大量差异基因富集在核糖体、氨基酰基-tRNA生物合成、嘌呤代谢、肽聚糖生物合成、嘧啶代谢、同源重组等途径当中，说明自溶发生时菌体内的生物合成途径发生显著改变。氨基酰基-tRNA主要功能是通过将氨基酸与含有相应反密码子的tRNA精确匹配，把氨基酸传递到核糖体当中进行肽链的合成（Gomez等，2020；Ibba等，2000）。差异基因富集在核糖体、氨基酰基-tRNA生物合成途径，说明自溶状态下菌体内蛋白质的合成发生了变化；而差异基因同时富集在嘌呤代谢、嘧啶代谢以及同源重组途径则说明了菌体内核酸代谢途径同样发生了改变。蛋白质与核酸合成的改变可能是屎肠球菌自溶发生的原因，也可能是自溶发生导致的结果，关于这一点，还需后续更加深入的研究来加以判断。

肽聚糖是革兰氏阳性菌细胞壁的重要组成，革兰氏阳性菌的自溶酶多是一些肽聚糖水解酶类（Chapot-Chartier等，2014；Vollmer等，2008）。大量差异基因富集在肽聚糖生物合成途径中也证明了肽聚糖合成与自溶发生的强相关性。此外，将差异基因中下调的基因进行KEGG富集分析发现，下调的差异基因主要富集在碳代谢、糖酵解/糖异生、丙酮酸代谢等能量代谢相关通路中，说明屎肠球菌的自溶与对碳源的利用能力、能量的利用能力相关。以往有研究显示碳源的不足会导致细菌自溶的发生，如培养基中碳源耗尽时会导致S.thermophilus和L.lactis的自溶 （Kawasaki等，1998；Douglas等，1974）。结合屎肠球菌下调差异基因的分析结果，推测有可能是长期培养下屎肠球菌对碳源利用能力的下降导致能量供应不足进而促进了菌体自溶的发生。

4 结论

本研究采用正常状态下和自溶状态下的屎肠球菌进行转录组测序，初步明确了自溶状态下屎肠球菌差异基因表达情况、GO功能分类及KEGG富集通路，分析发现自溶状态下屎肠球菌代谢发生了改变，能量利用出现下降。研究结果为进一步分析屎肠球菌自溶相关机制，改善屎肠球菌自溶发生情况提供了理论研究基础。