甘南牦牛曲拉中真菌群落結構

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

1(甘肅農業大學 食品科學與工程學院，甘肅 蘭州，730070)2(甘肅農業大學 理學院，甘肅 蘭州，730070)

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

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

本研究采用Illumina高通量測序技術對采自甘肅省甘南藏族自治州的曲拉樣品中真菌多樣性進行分析，以期全面解析曲拉中的真菌組成及群落結構，為曲拉的安全生產提供理論指導，同時為適合曲拉發酵微生物的篩選奠定基礎。

1 材料與方法

1.1 材料與試劑

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

1.2 儀器與設備

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

1.3 實驗方法

1.3.1 曲拉樣品的采集

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

1.3.2 曲拉微生物總DNA的提取

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

1.3.3 PCR擴增及測序

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

1.3.4 高通量測序數據處理

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

1.3.5 群落多樣性和統計分析

利用Mothur(V.1.31.2)軟件進行Alpha多樣性分析，并在不同的分類水平上對群落結構進行了統計分析;

使用R軟件對Weighted的UniFrac距離矩陣分別進行NMDS分析，通過二維排序圖描述群落樣本的結構分布;

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

使用Mothur軟件，計算優勢屬之間的Spearman等級相關系數，對其中rho>0.6且P<0.01的相關優勢屬構建關聯網絡，并導入軟件進行可視化。

1.4 數據分析

利用Excel 2010和Origin 2018軟件進行數據處理分析并作圖。

2 結果與分析

2.1 Alpha多樣性分析

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

表1 測序結果及真菌Alpha多樣性指數表Table 1 Sequencing results and fungal Alpha diversity index

Shannon指數和Simpson指數是綜合衡量物種多樣性的指數，其值越高，物種多樣性越豐富，反之物種多樣性越少。9個樣品中指數最高的分別為S63(3.53)和S01(0.86)，表明樣品S63和S01中真菌OTU的多樣性較高；S62的指數值最低，分別為1.20和0.37，反映了S62中真菌多樣性較低；Chao1指數和ACE指數主要側重于體現稀有群落的豐富度，指數越大，表明群落的豐富度越高。而曲拉樣品中Chao1指數和ACE指數最高的均為S63，分別為254.22和255.77，并且OTU也是最多的(391)。可以看出，樣品的Chao1指數和ACE指數大小與OTU數呈正相關，而群落多樣性與群落豐富度之間不存在相關性。

2.2 稀疏曲線

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

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

2.3 豐度等級曲線

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

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

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

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

子囊菌門為優勢菌門這一結果與新疆地區傳統發酵酸牛乳、對韓國酒精飲料和中國白酒中的真菌多樣性研究結果一致[23-25]。另外，在新疆阿圖什和烏什傳統發酵酸奶[26]中均檢出擔子菌門、接合菌門、壺菌門、球囊菌門和羅茲菌門，并且壺菌門、球囊菌門和羅茲菌門均為非優勢門。

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

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

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

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

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

莖點霉屬在我國屬于檢疫性病菌[29]，這是由于制作的環境存在安全隱患，并且在發酵豆醬中也檢測到該菌[30]，目前未有危害報道。

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

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

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

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

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

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

2.7 置換多元方差分析

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

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

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

2.8 優勢物種互作Spearman關聯網絡分析

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

圖7 優勢屬的關聯網絡圖Fig.7 Diagram of the associated network of dominant genus注：節點代表各優勢屬，以不同的顏色標識，節點之間的連接表明兩個屬之間存在相關性，紅線表明正相關，綠線表明負相關。通過某節點的連接越多，表明該屬與菌群中其他成員的關聯越多

3 結論

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