屎腸球菌自溶的基因表達差異分析

聶俊輝， 郭 浩， 王 通， 郭建軍， 曾 靜， 袁 林*

（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富集通路，分析發現自溶狀態下屎腸球菌代謝發生了改變，能量利用出現下降。研究結果為進一步分析屎腸球菌自溶相關機制，改善屎腸球菌自溶發生情況提供了理論研究基礎。