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

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

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

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

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

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型電泳槽,北京市六一儀器廠；凝膠成像系統(tǒng),美國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月采自甘肅省甘南藏族自治州合作市那吾鄉(xiāng)塔瓦(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區(qū)域進行擴增，ITS區(qū)擴增引物分別為ITS5F(GGAAGTAAAAGTCGTAACAAGG)和ITS1R(GCTGCGTTCTTCATCGATGC)。PCR條件如下：預變性為95 ℃、5 min，然后95 ℃、30 s、56 ℃、30 s、72 ℃、30 s共25個循環(huán)，72 ℃退火10 min，最后保存在4 ℃條件下。通過瓊脂糖凝膠電泳對PCR產(chǎn)物進行檢測，然后用試劑盒進行回收進行質(zhì)量檢測并建庫。最后由上海派森諾生物科技股份有限公司在MiSeq測序平臺進行雙端測序。

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

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

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

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

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

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

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

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

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

2 結(jié)果與分析

2.1 Alpha多樣性分析

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

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

Shannon指數(shù)和Simpson指數(shù)是綜合衡量物種多樣性的指數(shù)，其值越高，物種多樣性越豐富，反之物種多樣性越少。9個樣品中指數(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ù)呈正相關，而群落多樣性與群落豐富度之間不存在相關性。

2.2 稀疏曲線

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

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

2.3 豐度等級曲線

9份不同曲拉樣品的豐度等級曲線(rank-abundance curve)如圖2所示。由圖2可知，9個樣品的曲線趨勢相似。在水平方向，各樣品曲線寬度反映豐富度，在橫軸上的寬度，體現(xiàn)出不同樣品的物種豐度可能有較大的差異，其中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份樣品的共有優(yōu)勢門(相對豐度>1%)，平均相對豐度為96.544%；擔子菌門在S01、S02、S03、S63、S121、S122和S123樣品中為優(yōu)勢屬(相對豐度大于1%)，平均相對豐度為3.025%；接合菌門、壺菌門、球囊菌門在9份樣品中的豐度很低，為樣品中非優(yōu)勢門；另外，還有一些在門水平上未鑒定出。

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

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

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

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

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

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

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

圖4 各樣品屬水平菌群分布相對豐度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維度上有明顯的聚集和分離趨勢，樣品S01、S02和S03聚為第一類(A)，S61、S62和S63聚為第二類(B)，S121、S122和S123聚為第三類(C)。A和C樣本中各點分布比較緊密，表明樣品個體間差異性較小，微生物群落結(jié)構(gòu)相似；B樣本中各點分布疏散，樣品個體間差異性較大。原因可能是曲拉是自然狀態(tài)下發(fā)酵并且制作工藝比較粗放。這與李偉程等對傳統(tǒng)發(fā)酵乳制品中微生物多樣性研究結(jié)果一致[31]。

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

聚類分析主要是以等級樹的形式展示樣品之間的相似度，通過聚類樹的分枝長度衡量聚類效果的好壞。圖6為基于加權(quán)UniFrac距離的曲拉樣品的聚類圖，由圖6可知，不同來源的樣品聚集到不同的類別，說明不同來源的樣品微生物多樣性存在一定的差異性。在S121、S122和S123曲拉樣品之間，其分枝長度最短，樣本之間相似度最高；其次為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)計檢驗思路，通過對距離矩陣進行置換檢驗，從而評價原始樣本組間差異的大小及其統(tǒng)計學顯著性。表2為基于加權(quán)UniFrac距離的Adonis/PERMANOVA分析，由表可知，F(xiàn)=MSt/MSe=16.34,根據(jù)df1=2，df2=6查F檢驗表可得，F(xiàn)0.01(2,6)=10.92，而F>F0.01(2,6)，P<0.01，表明不同來源的樣品之間差異極顯著。而“Pr(>F)”是通過999次置換檢驗獲得的P值，P值越小組間差異性就越強。本研究中不同組間的樣品真菌組差異極顯著(P=0.001)。

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

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

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

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

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

3 結(jié)論

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