鄭慧娟，李偉程，趙飛燕，蔡虹宇，孫志宏，張和平

（內(nèi)蒙古農(nóng)業(yè)大學 乳品生物技術與工程教育部重點實驗室 農(nóng)業(yè)農(nóng)村部乳制品加工重點實驗室內(nèi)蒙古自治區(qū)乳品生物技術與工程重點實驗室 呼和浩特010018）

乳酸乳球菌乳酸亞種（Lactococcus lactis subsp.lactis）是兼性厭氧、不產(chǎn)氣、不產(chǎn)芽孢且無莢膜的革蘭氏陽性菌，是公認安全的食品級微生物[1]。其可將大分子營養(yǎng)物質(zhì)降解成易被人體消化吸收的小分子營養(yǎng)物質(zhì)，其中2，3-丁二醇等物質(zhì)可賦予發(fā)酵制品特殊的滋氣味，同時產(chǎn)生的乳酸等還可抑制有害菌生長，延長發(fā)酵制品的貨架期[2]。此外乳酸乳球菌乳酸亞種具有多種益生功能，被作為發(fā)酵劑廣泛應用于發(fā)酵乳制品、發(fā)酵果蔬制品與發(fā)酵肉制品等[3]。乳酸乳球菌乳酸亞種較高的經(jīng)濟價值使其成為研究熱點。隨著分子生物學的發(fā)展，測序手段與生物信息學分析的成熟，研究者從基因?qū)用媪私馊樗崛榍蚓樗醽喎N愈發(fā)便利。越來越多乳酸乳球菌乳酸亞種的全基因組序列完成圖被鑒定，為生產(chǎn)提供了遺傳學理論指導[4-5]。

隨著研究地深入，部分研究發(fā)現(xiàn)基因型相似的菌株間表型存在較大的差異性，而這些差異可能由菌株的遺傳表觀修飾造成的[6-7]。其中，DNA甲基化是一種典型遺傳表觀修飾，指DNA 堿基在甲基轉移酶的作用下增加甲基基團化學修飾的過程，其具有不改變DNA 序列的特點。Zhu 等[6]發(fā)現(xiàn)平均核苷酸同源性高于95%的伯克霍爾德氏菌（Burkholderia seminalis）菌株適應于不同生態(tài)位，可能是菌株甲基化組的差異導致的。序列相似度高于98%的短雙歧桿菌（Bifidobacterium breve）擁有多種可保護自身免受外源DNA 侵染甲基轉移酶[7]。Cooper 等[8]通過對不同血清型的大腸桿菌（Escherichia coli）進行基因組測定，發(fā)現(xiàn)一部分菌株毒力的產(chǎn)生與其DNA 的甲基化有關。目前許多甲基化相關研究發(fā)現(xiàn)細菌中的甲基化修飾位點以m6A（N6-methyladenosine，N6-甲基腺嘌呤）為主[9]。這類研究的關注對象多集中于致病菌，而優(yōu)良發(fā)酵劑相關的研究鮮有報道。

BL19 與IMAU10978 均分離自乳源的乳酸乳球菌乳酸亞種，表型差異明顯。其中BL19 具有發(fā)酵劑的特點，包括：發(fā)酵速度快，蛋白水解徹底，能產(chǎn)生具有特殊香氣的揮發(fā)性物質(zhì)等[10]。本文旨在從基因組與甲基化組角度比較BL19 與IMAU10978，探究2 株菌發(fā)酵特性差異的原因，表明乳酸乳球菌乳酸亞種功能基因組特征和m6A 基序類型、分布與功能等甲基化組特征，推測發(fā)酵特性差異的原因，為發(fā)酵工業(yè)的菌株選育提供理論依據(jù)。

1 材料與方法

1.1 試驗菌株與試劑

本研究的實驗菌株BL19 與IMAU10978 均由內(nèi)蒙古農(nóng)業(yè)大學乳品生物技術與工程教育部重點實驗室提供。BL19 分離自俄羅斯布里亞特共和國昆庫爾的嚼口，IMAU10978 分離自中國內(nèi)蒙古自治區(qū)赤峰市的酸牛奶。同時從NCBI（National Coalition Building Institute，https：//www.ncbi.nlm.nih.gov/）下載分離自乳源且被廣泛研究的乳酸乳球菌乳酸亞種KLDS 4.0325、Il1403、FM03、S0 和UC06 的序列，以進行比較分析（序列登錄號如表1所示）。

M17 培養(yǎng)基，英國Oxoid 公司；RNA 酶，天根生化科技有限公司；Wizard?Promega Genomic DNA Purification Kit，美國Promega 公司；核酸染料，北京全式金生物技術有限公司。

1.2 儀器、設備

Qubit2.0 熒光計，美國Life Technologies 公司；LRH-250 電熱恒溫培養(yǎng)箱、HWS28 電熱恒溫水浴鍋，上海一恒科技有限公司；5810R 臺式高速冷凍離心機，德國Eppendorf 公司；ZHJH-C1214C型智能安全型超凈工作臺，上海智城分析儀器制造有限公司；PacBio SMRT RSⅡ測序平臺，Pacific Biosciences 公司。

1.3 菌株活化與菌泥制備

將BL19 與IMAU10978 通過M17 肉湯培養(yǎng)基活化至3 代，培養(yǎng)條件為37℃，24 h。對第3 代培養(yǎng)液離心以收獲菌泥，4 000×g，5 min。之后用PBS 洗滌菌泥2 次，7 000×g 離心2 min。

1.4 DNA 提取與純化

按照Wizard?Genomic DNA Purification Kit說明書對潔凈菌泥進行DNA 提取。

通過0.6×瓊脂糖凝膠對提取的DNA 進行電泳，以判定所提取DNA 條帶的完整性與純凈性。用0.46×磁珠對DNA 進行吸附。棄廢液，并用75%乙醇對吸附DNA 的磁珠進行快速洗滌以純化DNA。最后用回溶液將純化后的DNA 從磁珠上回溶。并用Qubit 檢測純化后DNA 的質(zhì)量濃度，通過所得DNA 體積與質(zhì)量濃度相乘獲得所純化DNA 的最終質(zhì)量。

1.5 全基因組測序及組裝

上機要求每株菌的DNA 體積在150 μL 以內(nèi)，質(zhì)量大于7 ng。通過PacBio RSⅡ測序平臺對高質(zhì)量的乳酸乳球菌乳酸亞種DNA 進行全基因組測序。通過SMRT?Portal（V2.7）軟件的方案RS_HGAP_Assembly.3 進行質(zhì)控和菌株基因組組裝。使用circlator（V1.5.5）對下機數(shù)據(jù)進行環(huán)化已得到基因組完成圖。

1.6 基因組比較分析

1.6.1 基因組分析 計算試驗菌株的平均核苷酸一致性（Average Nucleotide Identity，ANI）與總核苷酸一致性（Total Nucleotide Identity，TNI）[11-12]。

1.6.2 功能基因組分析 將上述組裝好的核酸序列文件上傳至RAST（Repaid Annotion using Subsystem Technology http：//rast.nmpdr.org/rast.cgi）進行注釋，并下載對應的氨基酸注釋文件。

將BL19 與IMAU10978，以及從NCBI 下載的FM03、Il1403、KLDS 4.0325、S0、UC06 菌株的氨基酸注釋序列與蛋白相鄰類的聚簇（Cluster of Orthologous Groups of proteins，COG）數(shù)據(jù)庫比對，設置參數(shù)：E 值小于10-3。采納問詢序列與數(shù)據(jù)庫序列相似度大于40%且問詢序列比對上長度與序列總長度的比值大于80%的比對結果[13]。

將各菌株的蛋白注釋序列與碳水化合物活性酶（Carbohydrate-Active enZYmes，CAZy）數(shù)據(jù)庫比對。采納E 值小于10-3、相似度大于40%且問詢序列比對上長度與序列總長度的比值大于80%的比對結果[13]。

1.7 甲基化組分析

使用SMRT?Portal（V2.7）軟件中RS_Motification_and_motify_Analysis.1 方案對BL19 與IMAU10978 進行甲基化分析，參考序列選擇對應組裝效果好的序列。對基因組的甲基化位點進行預測分析，并分析特定基序的分布及其功能。

1.8 數(shù)據(jù)作圖

使用R（V3.5.0）軟件中的“pheatmap”程序包繪制熱圖。使用OriginPro 2015 繪制折線圖、柱狀圖與3D 色帶圖（3D Ribbons）。

2 試驗結果

2.1 基因組概況

菌株BL19 的基因組大小為2.72 Mb，包含5個質(zhì)粒，GC 含量為35.19%，菌株IMAU10978 的基因組大小為2.82 Mb，包含6 個質(zhì)粒，GC 含量為35.20%。從NCBI 下載的菌株UC06、S0、KLDS 4.0325、Il1403、FM03 的基因組大小為（2.57 ±0.17）Mb，其GC 含量為（35.26 ±0.07）%。除了Il1403 與S0，所有菌株的基因組均含有多個質(zhì)粒。

表1 試驗菌株的基本信息Table 1 Basic information of experimental strains

ANI 反映基因組間的相似性，比DNA 雜交（DNA-DNA hybridization，DDH）更簡單便利，常被用于進化距離的計算[14]。由圖1可知，試驗菌株之間的ANI 相似度均高達98.24%，TNI 均高達77.63%。這2 個最低值均高于乳酸乳球菌乳酸亞種的平均值，說明試驗選取菌株具有較高的基因組相似性[15]。

圖1 試驗菌株的核苷酸一致性Fig.1 Nucleotide identity of test strains

2.2 功能預測

2.2.1 RAST 注釋RAST 注釋結果（圖2a）顯示乳酸乳球菌乳酸亞種基因組以碳水化合物、氨基酸及其衍生物與蛋白代謝相關基因為主。其中BL19 與IMAU10978 蛋白質(zhì)折疊、蛋白質(zhì)生物合成、蛋白質(zhì)加工和改性相關基因數(shù)量基本一致。同時2 株菌均包含低聚果糖和棉子糖利用與乳糖和半乳糖的攝取和利用等碳水化合物相關基因，其數(shù)量不存在明顯差別。

2.2.2 COG 注釋 試驗選取菌株的基因組約83.60%得到注釋。圖2b 發(fā)現(xiàn)試驗菌株具有類似的特點，即E（氨基酸運輸和代謝）、J（翻譯，核糖體結構和生物合成）與G（碳水化合物轉移和代謝）3 種已知功能大類在基因組中占最高比例。其中，氨基酸運輸和代謝功能相關的基因數(shù)量相同。

圖2 試驗菌株的功能注釋Fig.2 Functional annotations of test strains

2.2.3 CAZy 注釋對試驗菌株進行CAZy 注釋。結果（圖2c）表明，7 株乳源乳酸乳球菌乳酸亞種均注釋到碳水化合物酯酶家族（Carbohydrate esterase family，CE）、糖苷水解酶家族（Glycoside hydrolase family，GH）、糖基轉移酶家族（Glycosyltransferase family，GT）和輔助酶家族（Auxiliary activities family，AA）的基因，包括CE3、CE10、GH2、GT3、GH20、GH125、CE1、GH1、GH77、CE9與GT5 等酶的基因，其中以GH1（包括β-葡萄糖苷酶與β-半乳糖苷酶等）的基因居多。

3 種功能注釋結果均顯示試驗菌株在功能基因組層面表現(xiàn)出極高的相似性。與菌株生長代謝相關的基因在基因組中占更高的比例，包括碳水化合物與氨基酸代謝等。Kelleher 等[16]認為乳酸乳球菌乳酸亞種區(qū)別于乳酸乳球菌乳脂亞種的特點為富含豐富的碳水化合物相關基因。其COG 注釋顯示乳酸乳球菌乳酸亞種的碳水化合物代謝和氨基酸代謝相關基因明顯高于乳脂亞種，分別占基因組14%和15%。其結果略高于本研究的結果，這可能是不同的篩選條件導致的，總之2 種結果的整體規(guī)律是一致的。Sun 等[17]發(fā)現(xiàn)乳酸菌含有豐富多樣的GH 與GT 基因，本研究的注釋結果相似也發(fā)現(xiàn)了類似規(guī)律。

2.3 甲基化比較結果

由于試驗菌株基因組較為相似，無法很好地解釋BL19 與IMAU10978 的發(fā)酵特性差異，因而從甲基化組角度對2 株菌進行更深入地分析。由圖3可知，2 株菌堿基A 的甲基化質(zhì)量值（Modification quality value，Mod QV）均明顯高于其它堿基，同時堿基A 與C 均高于甲基化檢測閾值。甲基化質(zhì)量值可用來預示位點發(fā)生甲基化的概率，甲基化檢測閾值通常為Mod QV 30（P=0.001）[18]。同時，Murray 等[19]也認為SMRT 測序技術檢測甲基化是方便可靠。

圖3 IMAU10978（a）與BL19（b）的甲基化位點動力學檢測散點圖Fig.3 Scatter plot of kinetic detection of modification sites IMAU10978（a）and BL19（b）

2.3.1m6A 修飾位點及基序鑒定結果IMAU10978與BL19 的基因組中m6A 修飾位點分別有9 377與4 541 個。2 株菌中m6A 修飾的特有基序類型的統(tǒng)計情況如表2。BL19 基因組中甲基化率最高的m6A 甲基化基序為ACAAYC，其甲基化率為99.02%。IMAU10978 基因組中甲基化率最高的m6A 甲基化基序為ATGNNNNNGTTA 和GAAHNNN NNNTCC，其甲基化率分別為98.04%和97.50%。

表2 IMAU10978 與BL19 的甲基化基序信息Table 2 Methylation motif information of IMAU10978 and BL19

2.3.2m6A 修飾基序的分布對上述m6A 修飾基序進行分布統(tǒng)計，發(fā)現(xiàn)3 種m6A 修飾基序在基因間區(qū)與基因區(qū)分布如圖4a所示，3 種基序均大量分布于基因間區(qū)，推測m6A 修飾位點呈一種隨機分布。與分布于基因間區(qū)的m6A 修飾位點相比，分布于功能基因組的m6A 修飾位點可能對菌株產(chǎn)生更嚴重的影響，從而導致發(fā)生這類變化的菌被淘汰的更多[20]。

對3 種m6A 修飾基序的分布情況進行近一步的分析。由圖4b～d所示，BL19 的ACAAYC 與IMAU10978 的ATGNNNNNGTTA 均在染色體基因組中隨機分布，而GAAHNNNNNNTCC 則在染色體基因組的起始處分布頻次出現(xiàn)極短暫的小高峰，隨后迅速達到平穩(wěn)。3 種m6A 修飾基序的分布趨勢整體表現(xiàn)為隨機分布，在一定程度證實了基序隨機分布的推論。Zhu 等[21]在結核桿菌的甲基化基序分布趨勢中發(fā)現(xiàn)相似規(guī)律。

圖4 3 種基序的分布Fig.4 Distribution of three motifs

2.3.3m6A 修飾基因的注釋對分別含有3 種m6A修飾基序的基因進行功能注釋。由圖5可知，3 種基序中m6A 修飾位點的分布規(guī)律存在差異。BL19特有基序ACAAYC 主要分布于可移動元件蛋白（Mobile element protein）、復制蛋白（Replication protein）、噬菌體蛋白（Phage protein）、DNA 原酶（DNA primase，EC：2.7.7.-）等功能基因上。而IMAU11978 的基序ATGNNNNNGTTA 和GAAHNNNNNNTCC 則主要分布于流動元素蛋白（Mobile element protein）、位點特異性重組酶/DNA 轉化酶相關蛋白（Site-specific recombinase DNA invertase Pin related protein）、精氨酸/鳥氨酸反轉運子ArcD（Arginine/ornithine antiporter ArcD）、PTS 系統(tǒng)甘露醇特異性IIC 組分（PTS system mannitol-specific IIC component）與非特征化MFS 型轉運體（Uncharacterized MFS-type transporter）等功能基因上。

圖5 3 種特有基序的功能注釋結果Fig.5 Functional annotations of three specific motifs

與BL19 相比，IMAU11978 的基因組上含有較多m6A 修飾的肽轉運載體基因，這可能是2 株菌蛋白水解特性存在差異的一個原因。同時，2 株菌m6A 修飾基序在PTS 系統(tǒng)相關基因上的分布存在差異，IMAU11978 中m6A 修飾基因包括：甘露醇特異性IIC 組分（10 個）、細胞二糖特異性IIB 成分（Cellobiose-specific IIB component，EC：2.7.1.205，4 個）、海藻糖特異性IIA 成分（Trehalose-specific IIA component，EC：2.7.1.201，2 個）、纖維二糖特異性IIC 組分（Cellobiose-specific IIC component，1 個）與甘露醇特異性IIA 成分（Mannitolspecific IIA component，EC：2.7.1.197，1 個），共計18 個。BL19 中m6A 修飾基因有PTS 系統(tǒng)的甘露醇特異性IIC 組分（5 個）與乳糖特異性IIA 成分（Lactose-specific IIA component，EC：2.7.1.207，1個）。推測m6A 修飾位點在PTS 系統(tǒng)相關基因的差異分布可能導致了菌株能量代謝的差異[22]。

值得注意的是，BL19 中發(fā)生m6A 修飾基因有乙偶姻特異性的2，3-丁二醇脫氫酶（EC：1.1.1.4），推測這種酶可能在NAD+協(xié)助下，將2，3-丁二醇轉化為乙偶姻。GC-MS（Gas Chromatography-Mass Spectrometer）檢測顯示BL19 含有豐富的2，3-丁二醇[10]。因此推測該基因的m6A 修飾可能導致其表達活性降低，致使BL19 的發(fā)酵乳中2，3-丁二醇得以積累。Ferrocino 等[23]也發(fā)現(xiàn)香腸微生物基因組中的2，3-丁二醇脫氫酶主導的通路與乙偶姻含量呈正相關。

此外，在2 株菌中發(fā)生m6A 修飾的基因都注釋到了豐富的微量元素相關基因與限制修飾系統(tǒng)相關基因。IMAU11978 中HsdS 以及I 型限制修飾系統(tǒng)的限制子單元R（Restriction subunit R，EC：3.1.21.3）與特異性亞單位S（Specificity subunit S）的基因發(fā)生了m6A 修飾，BL19 含有I 型限制性修飾系統(tǒng)特異性亞單位S 的m6A 修飾基因。限制修飾系統(tǒng)（Restriction modification system）由甲基轉移酶與限制性內(nèi)切酶組成，其作用是保護細菌免受外源DNA 的侵染[24]。有研究發(fā)現(xiàn)結核桿菌與雙歧桿菌等原核生物均包含豐富的限制修飾系統(tǒng)[7，21]。其中HsdS 是I 型DNA 甲基轉移酶中負責識別特定DNA 序列的特異性亞基[25]。

綜上所述，BL19 的ACAAYC 和IMAU11978的ATGNNNNNGTTA 和GAAHNNNNNNTCC 基序?qū)δ芑虻挠绊?，包括：肽轉運載體、PTS 系統(tǒng)與2，3-丁二醇脫氫酶等相關基因，可能是導致2株菌表型差異的重要原因。

3 結論

本研究對發(fā)酵特性存在明顯差異的乳酸乳球菌乳脂亞種IMAU10978 與BL19 進行了基因組與甲基化組分析，基因組分析發(fā)現(xiàn)2 株菌的基因組相似性高，ANI 值高達98.24%。甲基化組分析結果顯示IMAU10978 與BL19 的m6A 甲基化差異巨大。BL19 的基序以m6A 修飾的ACAAYC 為主，其甲基化比率高達99.06%，主要分布于可移動元件蛋白、復制蛋白和噬菌體蛋白等功能基因上。I MAU10978 的基序以m6A 修飾的AT GNNNNNGTTA 與GAAHNNNNNNTCC 為主，其甲基化比率均高于97.02%，主要分布于流動元素蛋白、位點特異性重組酶/DNA 轉化酶相關蛋白、精氨酸/鳥氨酸反轉運子ArcD 與PTS 系統(tǒng)甘露醇特異性IIC 組分等功能基因上。2 株菌的肽轉運載體、PTS 系統(tǒng)與2，3-丁二醇脫氫酶相關基因的甲基化差異巨大。這可能是導致BL19 與IMAU10978 在發(fā)酵速度、蛋白水解與產(chǎn)生特殊香氣物質(zhì)等方面存在差異的原因。本研究展現(xiàn)了乳酸乳球菌乳脂亞種的基因組、m6A 甲基化特征及其功能，從遺傳學和表觀遺傳學2 個角度探究了菌株表型產(chǎn)生差異可能的原因，為理解遺傳表觀修飾對乳酸乳球菌乳脂亞種的影響及優(yōu)良發(fā)酵劑的篩選提供堅實的理論基礎和數(shù)據(jù)支持。