甘南牦牛曲拉中真菌群落結(jié)構(gòu)

馬江，文鵬程，羅俏俏，曹磊，朱艷，楊敏，張衛(wèi)兵*，張忠明*

1(甘肅農(nóng)業(yè)大學(xué) 食品科學(xué)與工程學(xué)院，甘肅 蘭州，730070)2(甘肅農(nóng)業(yè)大學(xué) 理學(xué)院，甘肅 蘭州，730070)

牦牛曲拉是藏族牧民將牦牛乳脫脂、自然發(fā)酵、脫水、干燥等工藝制備得到的一種發(fā)酵乳制品[1-2]。有研究報(bào)道，我國(guó)曲拉總產(chǎn)量約3萬t，其中1/3的留作牧民食用[3-4]。與新鮮牛乳相比，曲拉中蛋白質(zhì)含量高達(dá)75.0%以上，脂肪含量為4.0%～7.0%，營(yíng)養(yǎng)價(jià)值極高[5]，同時(shí)也具有調(diào)節(jié)人體腸道微生物和改善人體腸道菌群的作用[6]。甘南地區(qū)的曲拉是以牧民家庭自制為主，其制作環(huán)境相對(duì)開放，其中含有多種真菌資源。

目前，有關(guān)乳制品中真菌資源的報(bào)道較多。張曉旭[7]通過傳統(tǒng)培養(yǎng)法從新疆和內(nèi)蒙古曲拉中分離出91株酵母菌，對(duì)其生物學(xué)特性和發(fā)酵特性進(jìn)行比較，結(jié)果表明，地域不同的同種酵母亦表現(xiàn)不同的特性；楊俊俊[8]從牦牛曲拉中鑒定出畢赤酵母、釀酒酵母、乳酸克魯維酵母等為主要優(yōu)勢(shì)菌群；李先勝等[9]對(duì)西藏地區(qū)11份曲拉進(jìn)行分離篩選，Saccharomycescerevisiae、Kluyveromycesmarxianus、Debaryomyceshansenii、Candidazeylanoides和Torulasporadelbrueckii為曲拉樣品中的優(yōu)勢(shì)屬，且不同樣品之間菌落差異較大；次頓等[10]通過變性梯度凝膠電泳法(polymerase chain reaction-denaturing gradient gel electrophoresis，PCR-DGGE)從拉薩酥油中鑒定出優(yōu)勢(shì)真菌菌群為假絲酵母屬、亞羅酵母屬和畢赤酵母；烏仁圖雅[11]和張冬蕾[12]通過焦磷酸測(cè)序技術(shù)對(duì)傳統(tǒng)牦牛酸奶中真菌多樣性進(jìn)行研究，結(jié)果表明，牦牛酸奶中優(yōu)勢(shì)門均為子囊菌門。以上研究采用的方法主要為傳統(tǒng)培養(yǎng)法、變性梯度凝膠電泳、454焦磷酸測(cè)序等，具有通量低、操作復(fù)雜和準(zhǔn)確率低等缺陷。新興的Illumina高通量測(cè)序技術(shù)[13-14]具有操作簡(jiǎn)單、成本較低的優(yōu)勢(shì)，并且采用邊合成邊測(cè)序原理，結(jié)果可信度高。因此，通過Illumina高通量測(cè)序技術(shù)能夠全面而準(zhǔn)確地了解研究對(duì)象中微生物種類組成和結(jié)構(gòu)。

本研究采用Illumina高通量測(cè)序技術(shù)對(duì)采自甘肅省甘南藏族自治州的曲拉樣品中真菌多樣性進(jìn)行分析，以期全面解析曲拉中的真菌組成及群落結(jié)構(gòu)，為曲拉的安全生產(chǎn)提供理論指導(dǎo)，同時(shí)為適合曲拉發(fā)酵微生物的篩選奠定基礎(chǔ)。

1 材料與方法

1.1 材料與試劑

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

1.2 儀器與設(shè)備

Pico-21型臺(tái)式離心機(jī),Thermo Fisher；DYY-6C型電泳儀、DYCZ-21型電泳槽,北京市六一儀器廠；凝膠成像系統(tǒng),美國(guó)UVP公司；Q32866型Qubit 2.0分光光度計(jì),Invitrogen公司；T100TM Thermal Cyeler型PCR儀,BIO-RAD公司；MiSeq System SY-410-1003高通量測(cè)序儀,美國(guó)Illumina公司。

1.3 實(shí)驗(yàn)方法

1.3.1 曲拉樣品的采集

曲拉樣品于2017年9月采自甘肅省甘南藏族自治州合作市那吾鄉(xiāng)塔瓦(S01、S02、S03)、加拉(S61、S62、S63)和瑪崗村(S121、S122、S123)3個(gè)村莊的9個(gè)不同牧民家庭，將所有樣品裝進(jìn)自封袋中，置于冷藏箱中運(yùn)輸至實(shí)驗(yàn)室以備試驗(yàn)。

1.3.2 曲拉微生物總DNA的提取

曲拉樣品中微生物總DNA提取采用E.Z.N.A.Soil DNA Kit D5625-01試劑盒，按照使用說明從樣品中提取DNA。

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

使用真菌特異性引物對(duì)曲拉樣品所提取的DNA的ITS區(qū)域進(jìn)行擴(kuò)增，ITS區(qū)擴(kuò)增引物分別為ITS5F(GGAAGTAAAAGTCGTAACAAGG)和ITS1R(GCTGCGTTCTTCATCGATGC)。PCR條件如下：預(yù)變性為95 ℃、5 min，然后95 ℃、30 s、56 ℃、30 s、72 ℃、30 s共25個(gè)循環(huán)，72 ℃退火10 min，最后保存在4 ℃條件下。通過瓊脂糖凝膠電泳對(duì)PCR產(chǎn)物進(jìn)行檢測(cè)，然后用試劑盒進(jìn)行回收進(jìn)行質(zhì)量檢測(cè)并建庫。最后由上海派森諾生物科技股份有限公司在MiSeq測(cè)序平臺(tái)進(jìn)行雙端測(cè)序。

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

MiSeq測(cè)序得到的數(shù)據(jù)采用Mothur(V.1.31.2)和QIIME(V.1.7.0)軟件進(jìn)行處理及分析[15-16]。首先采用滑動(dòng)窗口法對(duì)FASTQ格式的雙端序列逐一進(jìn)行質(zhì)量篩選，然后利用FLASH軟件(v1.2.7)對(duì)質(zhì)量初篩的雙端序列進(jìn)行配對(duì)連接。將連接后的序列識(shí)別分配對(duì)應(yīng)樣本，從而獲得有效序列。在測(cè)序過程中會(huì)產(chǎn)生一些錯(cuò)誤或疑問序列，因此采用QIIME軟件(v1.8.0)[17]識(shí)別疑問序列。通過QIIME軟件(v1.8.0)調(diào)用USEARCH(v5.2.236)檢查并剔除嵌合體序列[18-21]。使用QIIME軟件，調(diào)用UCLUST算法進(jìn)行序列聚類，以97%的序列相似度進(jìn)行歸并和OTU[22]劃分。測(cè)序數(shù)據(jù)在NCBI數(shù)據(jù)庫中的收錄編號(hào)為PRJNA431342(https://www.ncbi.nlm.nih.gov/bioproject/431342)。

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

利用Mothur(V.1.31.2)軟件進(jìn)行Alpha多樣性分析，并在不同的分類水平上對(duì)群落結(jié)構(gòu)進(jìn)行了統(tǒng)計(jì)分析;

使用R軟件對(duì)Weighted的UniFrac距離矩陣分別進(jìn)行NMDS分析，通過二維排序圖描述群落樣本的結(jié)構(gòu)分布;

使用QIIME軟件進(jìn)行UPGMA聚類分析和Adonis/PERMANOVA多元方差分析;

使用Mothur軟件，計(jì)算優(yōu)勢(shì)屬之間的Spearman等級(jí)相關(guān)系數(shù)，對(duì)其中rho>0.6且P<0.01的相關(guān)優(yōu)勢(shì)屬構(gòu)建關(guān)聯(lián)網(wǎng)絡(luò)，并導(dǎo)入軟件進(jìn)行可視化。

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

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

2 結(jié)果與分析

2.1 Alpha多樣性分析

9份不同曲拉樣品真菌Alpha多樣性指數(shù)如表1所示，通過真菌的ITS區(qū)測(cè)序，9份樣品共產(chǎn)生高質(zhì)量序列663 037條，將所有序列按97%的相似度進(jìn)行OTU聚類，得到1 667個(gè)OTU。由序列數(shù)及OTU聚類可以看出，曲拉中真菌種類繁多，物種豐富，且不同家庭手工制作的曲拉樣品中存在一定差異。

表1 測(cè)序結(jié)果及真菌Alpha多樣性指數(shù)表Table 1 Sequencing results and fungal Alpha diversity index

Shannon指數(shù)和Simpson指數(shù)是綜合衡量物種多樣性的指數(shù)，其值越高，物種多樣性越豐富，反之物種多樣性越少。9個(gè)樣品中指數(shù)最高的分別為S63(3.53)和S01(0.86)，表明樣品S63和S01中真菌OTU的多樣性較高；S62的指數(shù)值最低，分別為1.20和0.37，反映了S62中真菌多樣性較低；Chao1指數(shù)和ACE指數(shù)主要側(cè)重于體現(xiàn)稀有群落的豐富度，指數(shù)越大，表明群落的豐富度越高。而曲拉樣品中Chao1指數(shù)和ACE指數(shù)最高的均為S63，分別為254.22和255.77，并且OTU也是最多的(391)?？梢钥闯觯瑯悠返腃hao1指數(shù)和ACE指數(shù)大小與OTU數(shù)呈正相關(guān)，而群落多樣性與群落豐富度之間不存在相關(guān)性。

2.2 稀疏曲線

9份不同曲拉樣品的香濃指數(shù)稀疏曲線如圖1所示。由圖1可知，不同樣品的香濃指數(shù)隨著測(cè)序量的增加而呈現(xiàn)上升趨勢(shì)，說明在此測(cè)序水平下，樣品中真菌微生物的多樣性較高，并且樣品中還有較多的物種還沒有被檢測(cè)到；當(dāng)測(cè)序量較高時(shí)，香濃曲線逐漸與X軸接近平行，說明在此測(cè)序水平下，樣品中真菌的群落多樣性已能夠充分的展現(xiàn)。

圖1 香濃指數(shù)稀疏分析圖Fig.1 The sparse analysis diagram of Shannon index

2.3 豐度等級(jí)曲線

9份不同曲拉樣品的豐度等級(jí)曲線(rank-abundance curve)如圖2所示。由圖2可知，9個(gè)樣品的曲線趨勢(shì)相似。在水平方向，各樣品曲線寬度反映豐富度，在橫軸上的寬度，體現(xiàn)出不同樣品的物種豐度可能有較大的差異，其中S63豐富度最高，S62、S123豐富度最低。另外，曲線的形狀反映樣品的均勻度，曲線越平緩，群落組成的均勻度越高，曲線越陡峭，則群落中各OTU間的豐度差異越大，均勻度越低。圖2中S63均勻度最高，群落中各OTU間的差異最小。

圖2 豐度等級(jí)曲線圖Fig.2 The diagram of rank abundance curve

2.4 樣品中真菌群落在門水平的比較

圖3為9份曲拉樣品從門的分類水平進(jìn)行鑒定。在曲拉樣品中共檢測(cè)出6個(gè)門，分屬于子囊菌門(Ascomycota)、擔(dān)子菌門(Basidiomycota)、接合菌門(Zygomycota)、壺菌門(Chytridiomycota)、球囊菌門(Glomeromycota)和羅茲菌門(Rozellomycota)。由圖2可知，子囊菌門為9份樣品的共有優(yōu)勢(shì)門(相對(duì)豐度>1%)，平均相對(duì)豐度為96.544%；擔(dān)子菌門在S01、S02、S03、S63、S121、S122和S123樣品中為優(yōu)勢(shì)屬(相對(duì)豐度大于1%)，平均相對(duì)豐度為3.025%；接合菌門、壺菌門、球囊菌門在9份樣品中的豐度很低，為樣品中非優(yōu)勢(shì)門；另外，還有一些在門水平上未鑒定出。

子囊菌門為優(yōu)勢(shì)菌門這一結(jié)果與新疆地區(qū)傳統(tǒng)發(fā)酵酸牛乳、對(duì)韓國(guó)酒精飲料和中國(guó)白酒中的真菌多樣性研究結(jié)果一致[23-25]。另外，在新疆阿圖什和烏什傳統(tǒng)發(fā)酵酸奶[26]中均檢出擔(dān)子菌門、接合菌門、壺菌門、球囊菌門和羅茲菌門，并且壺菌門、球囊菌門和羅茲菌門均為非優(yōu)勢(shì)門。

圖3 各樣品門水平菌群分布相對(duì)豐度Fig.3 The relative abundance of horizontal flora distributionin phylum level of each sample

2.5 樣品中真菌群落在屬水平的比較

圖4為9份不同曲拉樣品從屬的分類水平進(jìn)行鑒定，共檢測(cè)出123個(gè)屬。從屬分類水平來看，9份曲拉樣品中畢赤酵母屬(Pichia)和雙足囊菌屬(Dipodascus)為共有優(yōu)勢(shì)屬，平均相對(duì)豐度分別為26.635%和11.393%；解脂耶式酵母屬(Yarrowia)在樣品S02、S03、S121、S122和S123中為優(yōu)勢(shì)屬，平均相對(duì)豐度為24.634%；念珠菌屬(Candida)和曲霉屬(Aspergillus)在樣品S01、S02、S03和S63中為優(yōu)勢(shì)屬，平均相對(duì)豐度分別為4.179%和3.640%；毛孢子菌屬(Trichosporon)在S121、S122和S123中為優(yōu)勢(shì)屬，平均相對(duì)豐度為1.711%；Archaeorhizomyces在樣品S01和S02中占優(yōu)勢(shì)，馬拉色霉菌屬(Malassezia)、絲孢畢赤氏酵母屬(Hyphopichia)和莖點(diǎn)霉屬(Phoma)僅在S63中占優(yōu)勢(shì)：Phaeoacremonium只在S01中為優(yōu)勢(shì)屬。

這是由于曲拉主要以家庭作坊式生產(chǎn)為主，發(fā)酵過程易受奶源、制作方法、海拔、地理環(huán)境、氣候環(huán)境、發(fā)酵溫度、發(fā)酵時(shí)間等影響[27]。另外，樣品中還包含許多在屬水平上未鑒定的屬，未鑒定出屬也具有較高的豐度，其豐度值在7.363%～55.276%，平均豐度值為23.818%。

畢赤酵母、解耶氏酵母和曲霉屬也在西藏曲拉檢測(cè)到，這是由于這些酵母是酸凝乳[28]干酪中最常見的菌，能夠利用發(fā)酵乳糖產(chǎn)生乳酸，使凝乳的pH有所升高，有助于發(fā)酵乳制品的成熟。

莖點(diǎn)霉屬在我國(guó)屬于檢疫性病菌[29]，這是由于制作的環(huán)境存在安全隱患，并且在發(fā)酵豆醬中也檢測(cè)到該菌[30]，目前未有危害報(bào)道。

圖4 各樣品屬水平菌群分布相對(duì)豐度Fig.4 The relative abundance of horizontal flora distribution in genus level of each sample

2.6 基于UniFrac距離的NMDS非度量多維尺度分析

圖5是基于加權(quán)UniFrac距離的NMDS分析，由圖5可知，9份不同來源的曲拉樣品在NMDS1和NMDS2維度上有明顯的聚集和分離趨勢(shì)，樣品S01、S02和S03聚為第一類(A)，S61、S62和S63聚為第二類(B)，S121、S122和S123聚為第三類(C)。A和C樣本中各點(diǎn)分布比較緊密，表明樣品個(gè)體間差異性較小，微生物群落結(jié)構(gòu)相似；B樣本中各點(diǎn)分布疏散，樣品個(gè)體間差異性較大。原因可能是曲拉是自然狀態(tài)下發(fā)酵并且制作工藝比較粗放。這與李偉程等對(duì)傳統(tǒng)發(fā)酵乳制品中微生物多樣性研究結(jié)果一致[31]。

圖5 加權(quán)UniFrac NMDS分析的樣本二維排序圖Fig.5 UniFrac NMDS analysis of two-dimensional sorting in samples

聚類分析主要是以等級(jí)樹的形式展示樣品之間的相似度，通過聚類樹的分枝長(zhǎng)度衡量聚類效果的好壞。圖6為基于加權(quán)UniFrac距離的曲拉樣品的聚類圖，由圖6可知，不同來源的樣品聚集到不同的類別，說明不同來源的樣品微生物多樣性存在一定的差異性。在S121、S122和S123曲拉樣品之間，其分枝長(zhǎng)度最短，樣本之間相似度最高；其次為S01、S02和S03，最后為S61、S62和S63曲拉樣品。

圖6 基于加權(quán)UniFrac距離矩陣的飛加權(quán)組平均法聚類分析圖Fig.6 UPGMA cluster analysis graph based on nweightedUniFrac distance matrix

2.7 置換多元方差分析

基于置換的PERMANOVA(permutational multivariate analysis of variance)[32]分析借鑒了ANOVA方差分析多組間差異的統(tǒng)計(jì)檢驗(yàn)思路，通過對(duì)距離矩陣進(jìn)行置換檢驗(yàn)，從而評(píng)價(jià)原始樣本組間差異的大小及其統(tǒng)計(jì)學(xué)顯著性。表2為基于加權(quán)UniFrac距離的Adonis/PERMANOVA分析，由表可知，F(xiàn)=MSt/MSe=16.34,根據(jù)df1=2，df2=6查F檢驗(yàn)表可得，F(xiàn)0.01(2,6)=10.92，而F>F0.01(2,6)，P<0.01，表明不同來源的樣品之間差異極顯著。而“Pr(>F)”是通過999次置換檢驗(yàn)獲得的P值，P值越小組間差異性就越強(qiáng)。本研究中不同組間的樣品真菌組差異極顯著(P=0.001)。

表2 加權(quán)UniFrac距離的置換多元方差分析Table 2 PERMANOVA analysis of Weighted UniFrac distance

注：“***”代表差異極顯著(P<0.001)

2.8 優(yōu)勢(shì)物種互作Spearman關(guān)聯(lián)網(wǎng)絡(luò)分析

關(guān)聯(lián)網(wǎng)絡(luò)基于微生物成員之間相互關(guān)系，對(duì)不同群落成員之間進(jìn)行分析，推斷不同微生物類群之間的的相互作用[33]。本研究使用Spearman等級(jí)相關(guān)系數(shù)計(jì)算牦牛曲拉樣品中屬之間的關(guān)系，并通過Cytoscape[34]軟件可視化。圖7為豐度在前50位的優(yōu)勢(shì)屬關(guān)聯(lián)網(wǎng)絡(luò)圖，由圖7可知，網(wǎng)絡(luò)圖由46個(gè)節(jié)點(diǎn)和106個(gè)邊組成。正相關(guān)和負(fù)相關(guān)之比為105∶1，Aspergillus、Penicillium、Neurospora、Rasamsonia、Apodospora、Placopsis、Chaetosphaeria、Leccinum、Mycosphaerella、Talaromyces、Exophiala、Wardomyces、Gyrophanopsis、Coriolopsis和Russula是網(wǎng)絡(luò)的中樞屬(每個(gè)節(jié)點(diǎn)≥6個(gè)邊)。第1優(yōu)勢(shì)屬畢赤酵母屬與優(yōu)勢(shì)屬毛孢子菌屬呈負(fù)相關(guān)，曲霉屬與Archaeorhizomyces、Guehomyces、Wallemia、Penicillum、Phaeoacremonium、Simplicillium呈正相關(guān)。而解脂耶式酵母、雙足囊菌屬和念珠菌與其他屬之間沒有相關(guān)性。

圖7 優(yōu)勢(shì)屬的關(guān)聯(lián)網(wǎng)絡(luò)圖Fig.7 Diagram of the associated network of dominant genus注：節(jié)點(diǎn)代表各優(yōu)勢(shì)屬，以不同的顏色標(biāo)識(shí)，節(jié)點(diǎn)之間的連接表明兩個(gè)屬之間存在相關(guān)性，紅線表明正相關(guān)，綠線表明負(fù)相關(guān)。通過某節(jié)點(diǎn)的連接越多，表明該屬與菌群中其他成員的關(guān)聯(lián)越多

3 結(jié)論

本研究基于Illumina MiSeq高通量測(cè)序平臺(tái)分析甘南牦牛乳曲拉樣品中真菌群落結(jié)構(gòu)及多樣性。結(jié)果表明，不同來源的曲拉樣品中的微生物多樣性存在差異性。曲拉中真菌群落組成分析表明，曲拉樣品中共有優(yōu)勢(shì)菌門為子囊菌門(Ascomycota)；共有優(yōu)勢(shì)屬為畢赤酵母屬(Pichia)和雙足囊菌屬(Dipodascus)；Beta多樣性結(jié)果表明真菌群落組成在不同來源的樣品中存在差異。Adonis/PERMANOVA多元方差分析表明，不同組間的樣品真菌組差異極顯著(P=0.001)。Spearman關(guān)聯(lián)網(wǎng)絡(luò)圖表明，真菌群落之間正相關(guān)占主導(dǎo)地位。