甘南地區(qū)牦牛曲拉中細(xì)菌群落結(jié)構(gòu)

曹 磊，梁春御，曹瑛瑛，文開勇，文鵬程，楊 敏，馮曉蘶，張忠明,*，張衛(wèi)兵

（1.甘肅農(nóng)業(yè)大學(xué)食品科學(xué)與工程學(xué)院，甘肅 蘭州 730070；2.蘭州資源環(huán)境職業(yè)技術(shù)學(xué)院氣象學(xué)院，甘肅 蘭州 730021；3.甘肅食品藥品監(jiān)督管理局，甘肅 蘭州 730070；4.甘肅農(nóng)業(yè)大學(xué)理學(xué)院，甘肅 蘭州 730070）

牦牛曲拉，又稱奶渣，是將牦牛乳脫脂后，在自然條件下進(jìn)行發(fā)酵使酪蛋白凝結(jié)、干燥后制成的一種發(fā)酵乳制品[1]。與新鮮牦牛乳比較，曲拉具有蛋白質(zhì)含量較高、易貯藏、方便運(yùn)輸?shù)奶攸c(diǎn)[1-2]。曲拉不僅是藏區(qū)牧民貯藏食品蛋白的重要手段，還是制作酸奶的發(fā)酵劑。另外，曲拉還可作為生產(chǎn)干酪素的原料[3]。由于曲拉加工環(huán)境相對(duì)開放，因此環(huán)境中的多種微生物參與了發(fā)酵過程[4]。

目前，微生物多樣性的分析方法主要包括傳統(tǒng)純培養(yǎng)方法和基于分子生物學(xué)技術(shù)的非純培養(yǎng)的方法。純培養(yǎng)的方法只能在實(shí)驗(yàn)室條件下對(duì)可培養(yǎng)的微生物進(jìn)行分析，不能全面了解微生物多樣性；非純培養(yǎng)的方法包括變性梯度凝膠電泳[5]、單鏈構(gòu)象多態(tài)性[6]、溫度梯度凝膠電泳[7]、末端限制性片段長(zhǎng)度多態(tài)性分析[8-10]等，其早期也廣泛應(yīng)用于乳制品中微生物多樣性的研究[11-13]。然而，近年來新興的Illumina MiSeq高通量測(cè)序技術(shù)能夠全面準(zhǔn)確地對(duì)研究對(duì)象中微生物種類組成和結(jié)構(gòu)進(jìn)行分析，相較于傳統(tǒng)的純培養(yǎng)方法及以16S rRNA為基礎(chǔ)的非純培養(yǎng)方法，能夠較多的產(chǎn)生測(cè)序覆蓋深度更大的數(shù)據(jù)量，更全面的呈現(xiàn)其微生物多樣性，明確各類微生物的結(jié)構(gòu)比例及優(yōu)勢(shì)菌群。該技術(shù)已在原乳[14]、發(fā)酵酸牛奶[9]、發(fā)酵酸駝奶和發(fā)酵酸馬奶[15]的微生物多樣性研究中得到了成功應(yīng)用。

本研究采用Illumina MiSeq高通量測(cè)序技術(shù)對(duì)采自于甘肅省甘南藏族自治州的曲拉樣品中細(xì)菌多樣性進(jìn)行分析，以期全面地解析曲拉中的細(xì)菌多樣性，系統(tǒng)地了解曲拉中細(xì)菌組成及群落結(jié)構(gòu)。

1 材料與方法

1.1 材料與試劑

曲拉樣品采自于甘肅省甘南藏族自治州合作市那吾鄉(xiāng)瑪崗村（MG1、MG2、MG3）、卡斯合村（KS1、KS2、KS3）、紹瑪村（SM1、SM2、SM3）3 個(gè)村莊的9 個(gè)不同牧民家庭。將所有樣品置于自封袋中編號(hào)，并置于冰盒中運(yùn)輸至實(shí)驗(yàn)室以備實(shí)驗(yàn)。

E.Z.N.A.Soil DNA試劑盒 美國Omega公司；Qubit 2.0 DNA檢測(cè)試劑盒 美國Invitrogen公司；Q5高保真DNA聚合酶 美國New England Biolabs公司；凝膠回收試劑盒 美國Axygen公司；TruSeq Nano DNA LT Library Prep Kit 美國Illumina公司。

1.2 儀器與設(shè)備

Pico-21型臺(tái)式離心機(jī) 美國Thermo Fisher科技公司；DYY-6C型電泳儀、DYCZ-21型電泳槽 北京市六一儀器廠；凝膠成像系統(tǒng) 美國UVP公司；Q32866型Qubit 2.0分光光度計(jì) 美國Invitrogen公司；T100TMThermal Cyeler型聚合酶鏈?zhǔn)椒磻?yīng)（polymerase chain reaction，PCR）儀 伯樂生命醫(yī)學(xué)產(chǎn)品（上海）有限公司；MiSeq System SY-410-1003高通量測(cè)序儀 美國Illumina公司。

1.3 方法

1.3.1 曲拉微生物總DNA的提取

取0.2 g曲拉樣品，采用E.Z.N.A.Soil DNA Kit D5625-01提取試劑盒，參照說明書從樣品中提取DNA；然后使用Qubit 2.0分光光度計(jì)對(duì)提取的DNA進(jìn)行定量，并通過0.8%瓊脂糖凝膠電泳確定DNA提取的完整性。

1.3.2 PCR擴(kuò)增及測(cè)序

以稀釋后的基因組總DNA為模板，取30 ng進(jìn)行PCR靶向擴(kuò)增16S rRNA基因。采用V3-V4區(qū)338F（5’-ACTCCTACGGGAGGCAGCAG-3’）和806R（5’-GGACTACHVGGGTWTCTAAT-3’）分別為特異性引物。16S rRNA基因PCR熱循環(huán)譜如下：95 ℃預(yù)變性5 min；95 ℃變性30 s，56 ℃退火30 s，72 ℃延伸30 s，25 個(gè)循環(huán)；72 ℃退火10 min，最終保持在4 ℃條件下。PCR產(chǎn)物使用2%瓊脂糖凝膠電泳進(jìn)行檢測(cè)，切割后使用Axygen公司的凝膠回收試劑盒回收。將擴(kuò)增產(chǎn)物純化后，利用分光光度計(jì)精確定量后制備測(cè)序文庫，最后在MiSeq測(cè)序平臺(tái)進(jìn)行雙端測(cè)序。PCR擴(kuò)增、測(cè)序文庫制備及測(cè)序工作由上海派森諾生物科技股份有限公司完成。

1.3.3 高通量測(cè)序數(shù)據(jù)處理

采用Mothur（version 1.31.2）和QIIME（version 1.7.0）軟件處理及分析高通量測(cè)序數(shù)據(jù)[16-17]。先對(duì)測(cè)序結(jié)果進(jìn)行整合，使每次讀數(shù)的堿基平均質(zhì)量不小于Q20。然后將整合后的序列使用FLASH（version 1.2.7）軟件根據(jù)重疊堿基進(jìn)行配對(duì)連接，并丟棄無法配對(duì)的序列[18]。使用UCHIME（version 4.2）軟件鑒定和去除嵌合體序列[19-20]。將序列分配入對(duì)應(yīng)樣本，使用UCLUST算法進(jìn)行序列聚類，以97%的序列相似度作為劃分閾值，將樣品聚類成可操作分類單位（operational taxonomic units，OTU）[18]。

1.3.4 群落多樣性和統(tǒng)計(jì)分析

利用Mothur（version 1.31.2）軟件進(jìn)行α多樣性分析，其中包括Simpson指數(shù)、Chao1指數(shù)和Shannon指數(shù)，并在不同的分類水平上對(duì)群落結(jié)構(gòu)進(jìn)行了統(tǒng)計(jì)分析。

使用QIIME軟件對(duì)基于Unweighted和Weighted的UniFrac距離矩陣進(jìn)行UPGMA聚類分析，并使用R（version 2.15.3）軟件進(jìn)行可視化。

通過PICRUSt（version 1.0）軟件對(duì)群落基因功能特征進(jìn)行預(yù)測(cè)[21]。

1.4 數(shù)據(jù)分析

利用Excel 2007軟件進(jìn)行數(shù)據(jù)處理分析并作圖。

2 結(jié)果與分析

2.1 PCR擴(kuò)增、測(cè)序序列及OTU基本情況

由圖1可以看出，PCR擴(kuò)增后得到清晰明顯的條帶，可用于后續(xù)的研究分析。通過細(xì)菌的16S rRNA基因V3-V4區(qū)測(cè)序，9 份樣品共產(chǎn)生216 630 條有效序列，將得到的有效序列經(jīng)過質(zhì)量控制，在剔除一些疑問序列之后得到高質(zhì)量序列為185 910 條，高質(zhì)量序列占有效序列的比例為85.82%。將所有序列按97%的相似度進(jìn)行OTU聚類，得到4 754 個(gè)OTU。由序列數(shù)及OTU聚類可以看出，曲拉中細(xì)菌種類繁多，物種豐富，且不同來源的曲拉樣品中存在一定差異。

圖1 PCR擴(kuò)增產(chǎn)物凝膠電泳圖Fig. 1 Electrophoresis profile of PCR amplified product

2.2 稀疏曲線

稀疏曲線用來評(píng)判每個(gè)樣品在當(dāng)前的測(cè)序深度下是否足以反映該樣品中所包含的細(xì)菌群落多樣性。由圖2可以看出，當(dāng)樣品測(cè)序量較低時(shí)，每個(gè)樣品的Shannon指數(shù)呈顯著上升的趨勢(shì)，說明在當(dāng)前測(cè)序量上樣品物種的多樣性較高，樣品中還有較多的物種沒有被檢測(cè)到。但隨著測(cè)序數(shù)量逐漸增加，各曲拉樣品中Shannon指數(shù)呈緩慢增加的趨勢(shì)，稀疏曲線的斜率逐漸降低，并且隨著測(cè)序數(shù)量的繼續(xù)增加，Shannon曲線基本穩(wěn)定。雖然在當(dāng)前的測(cè)序深度下每個(gè)樣品的稀疏曲線未能完全進(jìn)入平臺(tái)期，但與X軸幾乎接近平行，已經(jīng)達(dá)到飽和。因此，盡管隨著測(cè)序量的增加新的物種可能會(huì)被發(fā)現(xiàn)，但在此測(cè)序水平下，樣品中細(xì)菌的群落多樣性已能夠充分的展現(xiàn)。

圖2 Shannon指數(shù)稀疏分析圖Fig. 2 Sparse analysis of Shannon index

2.3 α多樣性分析

曲拉樣品中微生物α多樣性指數(shù)如表1所示。對(duì)于微生物群落而言，可以用Chao1指數(shù)、Shannon指數(shù)、ACE指數(shù)和Simpson指數(shù)反映樣品中物種的豐富度和多樣性。由表1可以看出，曲拉樣本中細(xì)菌的Chao1指數(shù)平均值為307.06±124.62，ACE指數(shù)為473.93±180.14，Simpson指數(shù)為0.74±0.13，Shannon指數(shù)為3.29±1.35。從微生物多樣性指數(shù)結(jié)果分析可以得出，曲拉樣品中細(xì)菌群落的物種多樣性及總體豐度相對(duì)較高，曲拉中的微生物數(shù)量處于較高的水平。同時(shí)，不同村莊來源的樣品微生物多樣性指數(shù)存在差異，這表明樣品制備環(huán)境可能會(huì)影響樣品中的微生物物種豐富度及多樣性。Matsoni[22]和Tarag[23]對(duì)于發(fā)酵乳的研究中也發(fā)現(xiàn)了相同的結(jié)果。在該研究中所有樣品的覆蓋率范圍為99.91%～99.98%，均大于99%，這表明測(cè)序水平能夠達(dá)到鑒定樣品中的大多數(shù)微生物群落的多樣性。

表1 微生物多樣性指數(shù)Table 1 Microbial diversity indexes

2.4 樣品細(xì)菌群落在門水平的分類及比較

表2 各樣品門水平菌群分布相對(duì)豐度Table 2 Relative abundance of bacterial flora distribution in Qula samples at phylum level

對(duì)9 份曲拉樣品從門的分類水平進(jìn)行鑒定，結(jié)果如表2所示。在曲拉樣品中共檢測(cè)出8 個(gè)門，分屬于厚壁菌門（Firmicutes）、變形菌門（Proteobacteria）、放線菌門（Actinobacteria）、酸桿菌門（Acidobacteria）、藍(lán)藻門（Cyanobacteria）、擬桿菌門（Bacteroidetes）、綠彎菌門（Chlorof l exi）和假單胞菌門（Gemmatimonadetes）。由表2可知，厚壁菌門為曲拉樣品中的第1優(yōu)勢(shì)門，其相對(duì)豐度在44.515%～81.010%的范圍內(nèi)，平均相對(duì)豐度為60.751%；其次為變形菌門，為樣品中的第2優(yōu)勢(shì)門，其相對(duì)豐度介于18.973%～55.363%之間，平均相對(duì)豐度為39.138%；放線菌門和酸桿菌門豐度較低，平均相對(duì)豐度僅為0.064%和0.026%，為樣品中的非優(yōu)勢(shì)門；藍(lán)藻門（平均相對(duì)豐度為0.012%）、擬桿菌門（平均相對(duì)豐度為0.003%）、綠彎菌門和假單胞菌門（平均相對(duì)豐度為0.001%）的豐度都很低，為曲拉樣品中的稀有門。厚壁菌門和變形菌門也是牛乳、馬乳、駝乳及其發(fā)酵酸乳中的優(yōu)勢(shì)門[15]，與本實(shí)驗(yàn)的研究結(jié)果一致。另外，在發(fā)酵牛乳[24]、發(fā)酵牦牛乳[25]、發(fā)酵馬乳[15,26]中均檢測(cè)到綠彎菌門、酸桿菌門、藍(lán)藻門，且均為非優(yōu)勢(shì)門，而在先前的研究中沒有檢測(cè)到假單胞菌門。

2.5 樣品細(xì)菌群落在屬水平的分類及比較

表3 各樣品屬水平菌群分布相對(duì)豐度Table 3 Relative abundance of bacterial flora distribution in Qula sample at genus level

對(duì)9 份曲拉樣品從屬的分類水平進(jìn)行鑒定，共檢測(cè)出55 個(gè)屬，結(jié)果如表3所示。從屬分類水平來看，在曲拉樣品中檢測(cè)出的乳桿菌屬（Lactobacillus）、醋酸桿菌屬（Acetobacter）、乳球菌屬（Lactococcus）、假單胞菌屬（Pseudomonas）和明串珠菌屬（Leuconostoc）為豐度相對(duì)較高的屬，其相對(duì)豐度均大于5%。厚壁菌門的乳桿菌屬的相對(duì)豐度在11.628%～80.946%的范圍內(nèi)，平均相對(duì)豐度為43.027%，為曲拉樣品中的第1優(yōu)勢(shì)屬；變形菌門的醋酸桿菌屬的相對(duì)豐度在5.057%～52.408%的范圍內(nèi)，平均相對(duì)豐度為24.739%，為第2優(yōu)勢(shì)屬；厚壁菌門的乳球菌屬的相對(duì)豐度相對(duì)較高，平均相對(duì)豐度為11.402%，為第3優(yōu)勢(shì)屬。假單胞菌屬和明串珠菌屬的平均相對(duì)豐度分別為7.971%和5.088%。不同來源的曲拉樣品在屬水平上細(xì)菌組成存在差異。另外，樣品中還包含許多在屬水平上未鑒定的屬（unclassified），未鑒定出屬也具有較高的豐度，其相對(duì)豐度在1.942%～10.928%的范圍內(nèi)，平均相對(duì)豐度為4.922%。

腸球菌屬（Enterococcus）和葡萄球菌屬（Staphylococcus）為樣品中的低豐度屬，其平均相對(duì)豐度分別為0.054%和0.001%。腸球菌屬及葡萄球菌屬屬于食源性病原菌屬，說明曲拉相對(duì)粗放的生產(chǎn)環(huán)境有一定的影響。

不同乳制品的微生物組成有著豐富的多樣性，如傳統(tǒng)酸奶[23]、發(fā)酵牛乳[24,26-27]、發(fā)酵馬乳[26]、發(fā)酵牦牛乳[25]和開非爾[28]的優(yōu)勢(shì)屬多為乳桿菌屬，與本實(shí)驗(yàn)的研究結(jié)果一致。奶酪中的優(yōu)勢(shì)屬各不相同，愛爾蘭奶酪[29]為乳球菌屬，丹麥原奶酪[28]中優(yōu)勢(shì)菌株歸屬于乳桿菌屬和鏈球菌屬。有些發(fā)酵乳制品樣品中也檢測(cè)到埃希氏桿菌屬（Escherichia）和沙門氏菌屬（Salmonella）[23]，葡萄球菌屬和志賀氏菌屬（Shellogella）[14]等也有檢測(cè)到。

2.6 物種豐度熱圖

豐度熱圖顏色越綠則說明所含的菌屬較少。由圖3可以看出，9 個(gè)樣本的前5 種菌屬的顯示顏色比較紅，含量較多，分別是乳桿菌屬、醋酸桿菌屬、乳球菌屬、假單胞菌屬、明串珠菌屬。另外，一些沒有與數(shù)據(jù)庫匹配的未鑒定出屬也具有較高的豐度。不同來源的曲拉樣品菌屬中有的樣品顏色比較紅，有的樣品顏色比較綠，說明每個(gè)樣品中菌群的豐度存在差異，因此，不同來源的樣品中菌群多樣性存在差異性。

圖3 屬水平物種豐度熱圖Fig. 3 Heatmap of species abundance at genus level

2.7 基于UniFrac距離的樣本聚類分析

圖4 曲拉中菌群結(jié)構(gòu)UniFrac非加權(quán)主坐標(biāo)分析Fig. 4 Unweighted UniFrac principal coordinate analysis of bacterial communities in Qula samples

采用基于UniFrac距離的Weighted和Unweighted主坐標(biāo)分析對(duì)樣品中細(xì)菌群落結(jié)構(gòu)進(jìn)行比較，結(jié)果分別見圖4、5。圖4為基于非加權(quán)UniFrac距離的主坐標(biāo)分析，其第1主成分、第2主成分和第3主成分的貢獻(xiàn)率分別為47.68%、28.60%和4.76%，不同來源的樣品在PC1和PC2維度上明顯的呈聚集和分離趨勢(shì)，個(gè)別樣品之間存在部分交疊，能夠?qū)⒉煌瑯悠份^好的區(qū)分，而在PC3維度上所有樣品可以達(dá)到很好區(qū)分效果。

圖5為基于加權(quán)UniFrac距離的主坐標(biāo)分析，其第1主成分、第2主成分和第3主成分的貢獻(xiàn)率分別是71.75%、27.93%和0.15%，不同樣品在PC1和PC2維度上有樣品之間存在交疊現(xiàn)象，而在PC3維度上所有樣品可以很好地區(qū)分。

由樣品主坐標(biāo)分析可知，所有樣品沒有緊密地組合在一起，各樣品所代表的點(diǎn)在空間分布有一定的距離，說明不同位置來源的曲拉樣品的細(xì)菌群落組成具有差異性。本研究中使用的曲拉樣品是通過傳統(tǒng)方法生產(chǎn)的，因此，樣品間細(xì)菌群落結(jié)構(gòu)的差異不同很可能歸因于原料、地理位置和其他環(huán)境因素。從圖6可以看出，不同采樣地的樣品聚集到不同的類別，說明不同來源樣品的微生物多樣性存在一定的差異性，這與主坐標(biāo)分析的結(jié)果一致。

圖5 曲拉中菌群結(jié)構(gòu)UniFrac加權(quán)主坐標(biāo)分析Fig. 5 Weighted UniFrac principal coordinate analysis of bacterial communities in Qula samples

圖6 基于Unweighted UniFrac距離矩陣的UPGMA聚類分析圖Fig. 6 Dendrogram from UPGMA cluster analysis based on unweighted UniFrac distance matrix

2.8 樣品中與細(xì)菌的基因代謝有關(guān)的功能特征

表4 與樣品中細(xì)菌的基因代謝有關(guān)的功能特征Table 4 Functional features related to genes involved in metabolism from bacteria in Qula samples

表5 樣本中細(xì)菌的其他功能基因Table 5 Other functional genes from bacteria in Qula samples

PICRUSt軟件用于預(yù)測(cè)樣品中存在的細(xì)菌類群的功能基因及其代謝途徑[30]。對(duì)曲拉樣本中的功能基因及代謝途徑預(yù)測(cè)情況如表4所示。在樣本中的功能基因中，44.04%～50.37%的微生物基因與以下代謝途徑有關(guān)：氨基酸代謝、碳水化合物代謝、能量代謝、核苷酸代謝、輔助因子及維生素代謝、外源物質(zhì)的生物降解與代謝、脂質(zhì)代謝、酶家族、糖生物合成與代謝、其他氨基酸的代謝、萜類化合物與多酮類化合物的代謝以及其他次生代謝產(chǎn)物的生物合成。其中，碳水化合物代謝和氨基酸代謝的基因在曲拉樣品中占主導(dǎo)地位，分別占所有基因的10.40%和8.68%。

樣品中存在的細(xì)菌的其他功能基因與細(xì)胞處理、環(huán)境信息處理、遺傳信息處理、人類疾病、機(jī)體系統(tǒng)等有關(guān)，具體情況如表5所示。其中，涉及膜轉(zhuǎn)運(yùn)的基因占主導(dǎo)地位，占所有基因的12.85%。

以不同樣品中前50 種功能基因繪制熱圖，結(jié)果如圖7所示。相同來源的曲拉樣品之間的功能基因相似程度較高，而不同來源的樣品的差異性較大，這一結(jié)果與前面細(xì)菌群落的聚類分析結(jié)果一致。

圖7 基于樣本前50 個(gè)功能基因的聚類分析熱圖Fig. 7 Heatmap from clustering analysis based on top 50 functional genes in Qula samples

3 結(jié) 論

本研究基于Illumina MiSeq高通量測(cè)序平臺(tái)分析曲拉樣品中細(xì)菌群落結(jié)構(gòu)及多樣性。結(jié)果表明，不同來源的曲拉樣品中的微生物多樣性存在差異性。與傳統(tǒng)方法相比較，高通量測(cè)序的研究結(jié)果能夠更全面地反映曲拉中菌群分布多樣性。曲拉中細(xì)菌群落組成分析表明，曲拉樣品中優(yōu)勢(shì)門為厚壁菌門和變形菌門；優(yōu)勢(shì)屬為乳桿菌屬、醋酸桿菌屬、乳球菌屬，細(xì)菌群落組成表明在不同來源的樣品中群落組成存在差異。樣品中存在的細(xì)菌的功能基因預(yù)測(cè)表明，不同來源樣品的細(xì)菌群落之間存在差異。