高通量測序分析不同采集保存溫度下甘南牧區牦牛酸奶真菌多樣性

張文齊，麻和平，劉彩云，王 潔，彭章普，邵建寧,

（1.甘肅省科學院生物研究所 甘肅省微生物資源開發利用重點實驗室，甘肅蘭州 730000；2.甘肅省科學院，甘肅蘭州 730000）

傳統發酵乳制品制備多采用自然發酵，原料的不同、獨特的制備方法和環境氣溫的差異等造就了傳統發酵乳制品微生物群落的豐富性和多樣性[1?3]。近年來，有關傳統發酵乳制品微生物多樣性研究，既有傳統的分離培養研究，也有高通量測序、宏基因組學分析等非培養研究[4?5]。利用傳統的分離技術從樣品中獲得的微生物，代表具體培養條件所能分離到的微生物，無法分離到樣品中所有微生物，不能完整地揭示樣品的菌群組成[6?7]。高通量測序、宏基因組學分析等非培養研究無需分離培養微生物，能直接反映樣品中微生物群落構成和多樣性[8?11]。

目前，對傳統發酵乳制品中參與發酵微生物的分離和菌群多樣性相關研究較多[12?14]，然而，通過高通量測序分析比較不同采集保存溫度下的傳統發酵乳制品微生物多樣性，優化傳統發酵乳制品采集保存溫度卻鮮有報道。Sun等[15]通過高通量焦磷酸測序對17 個自然發酵乳制品中的細菌和真菌群落多樣性進行研究，所有樣品中細菌屬47 個，真菌屬43 個；王遠微等[16]對川西北部分牧區的10 份傳統發酵牦牛酸奶樣品進行酵母菌的分離，通過常規形態特征和26S rRNA基因測序分析鑒定出16 株馬克斯克魯維酵母菌；袁鳳霞等[17]將前期從新疆、西藏和青海等地區傳統自制乳制品中分離的27 株酵母菌株分別與保加利亞乳桿菌發酵酸乳，篩選優良共生酵母菌株。還有研究表明傳統發酵奶制品在發酵過程中乳酸菌和酵母菌的共同作用使得發酵后的奶制品具有獨特的風味和口感[18?20]。

本研究以常溫、冰袋、干冰和液氮不同保存條件下采集的甘南牧區牧民自制牦牛酸奶樣品為研究對象，采用Illumina NovaSeq高通量測序技術，對牦牛酸奶擴增真菌ITS2 區進行測序，通過分析比較不同采集保存溫度下樣品的真菌群落結構和真菌多樣性，探究基于高通量測序不同采集保存溫度對傳統發酵乳制品微生物真菌多樣性分析的影響，為傳統發酵乳制品采樣條件優化提供理論指導，并為科學保護、開發與利用傳統發酵乳制品中優良微生物資源提供參考依據。

1 材料與方法

1.1 材料與儀器

牦牛酸奶樣品 20 份甘南牧區牦牛酸奶樣品的采樣具體信息見表1；QIAamp Fast DNA Stool Mini Kit、QIAquick Gel Extraction Kit德國QIAGEN公司；PowerSoil DNA Isolation Kit美國MOBIO公司；KAPA HiFi Hotstart ReadyMix Kit美國KAPA Biosystems公司；2× Phanta Master Mix南京諾唯贊生物科技股份有限公司；AxyPrep DNA凝膠回收試劑盒 美國AXYGEN公司。

表1 甘南牧區牦牛酸奶采樣信息Table 1 Information of yak yogurt samples collected from Gannan pastoral area

5424R低溫高速離心機 德國Eppendorf公司；VDRTEX-5 渦旋混合器 海門市其林貝爾儀器制造有限公司；Tanon 3500 凝膠成像系統 上海天能科技有限公司；DYY-6C電泳儀 上海天能科技有限公司；ND-2000C微量紫外分光光度計 美國NanoDrop公司；Applied Biosystems VeritiPro PCR儀 美國Life Technologies公司。

1.2 實驗方法

1.2.1 牦牛酸奶樣品的采集 實地采集牧民自制牦牛酸奶，每份樣品采集4 管，每管30 mL，密封后標記樣品號，分別置于常溫、冰袋、干冰和液氮容器中保存，運回實驗室后，分別于常溫、4、?20 ℃和?80℃保存備用。樣品采集自甘肅省甘南藏族自治州卓尼縣（Z）、碌曲縣（L），常溫采集實驗室常溫保存樣品分組為Z-0，冰袋采集實驗室4 ℃保存樣品分組為Z-4，干冰采集實驗室?20 ℃保存樣品分組為Z-20，液氮采集實驗室于?80 ℃保存樣品分組為Z-80。

1.2.2 牦牛酸奶樣品中總DNA提取 取2 mL牦牛酸奶樣品，采用試劑盒QIAamp Fast DNA Stool Mini Kit，按照說明書的要求提取牦牛酸奶的總DNA，利用微量紫外分光光度計和1%瓊脂糖凝膠電泳，進行總DNA完整性和濃度質檢。

1.2.3 PCR擴增與測序 使用KAPA HiFi Hotstart ReadyMix PCR kit將提取到的總DNA作為模板進行PCR擴增，采用引物為ITS3 GCATCGATGAAGAACGCAGC（F2045）和ITS4 TCCTCCGCTTATT-GATATGC（R2390）擴增真菌ITS2 區（引物中已帶特定barcode），PCR反應條件為：95 ℃ 3 min；98 ℃20 s，58 ℃ 15 s，72 ℃ 20 s，30 個循環；72 ℃ 5 min。PCR反應體系為：2×KAPA Library Amplification Ready Mix 15 μL，引物（10 μmol/L）各為1 μL，模板DNA為50 ng，加ddH2O補齊至30 μL。

擴增PCR產物用2%瓊脂糖凝膠電泳檢測，使用AxyPrep DNA凝膠提取試劑盒切膠回收PCR產物，文庫質檢合格后，使用Qubit 2.0 進行文庫定量，并根據每個樣品的數據量要求，進行相應比例的混合，采用Illumina NovaSeq PE250 進行上機測序，Illumina NovaSeq高通量測序由上海銳翌生物科技有限公司完成。

1.2.4 數據處理 ITS2 序列被控制在220～500 bp之間，保證每條reads平均質量值不低于20，且含N堿基數不超過3 個。將序列完全一樣的Clean Reads根據其豐度大小進行排序，將其中的Singletons過濾掉。使用UPARSE（http://drive5.com/uparse/）根據97%相似度進行OTU聚類，并使用Userach（版本7.0.1090）鑒定和移除嵌合體序列。每個OTU都有一個代表性的序列，使用RDP數據庫（http://rdp.cme.msu.edu/），置信度閾值設置為0.8，利用RDP Classifer（http://rdp.cme.msu.edu/）將每個代表性序列進行物種注釋。不同樣本對應的Reads數量差距較大，為避免因樣品數據大小不同而造成分析時的偏差，在樣品達到足夠的測序深度的情況下，對每個樣品進行隨機抽平處理。測序深度使用Alpha多樣性中的coverage指數來衡量，OTU profiling table、Alpha多樣性指數和Beta多樣性指數通過QIIME（version 1.9.1）的python腳本實現。

2 結果與分析

2.1 樣品序列信息及OTU分析

對從采集至實驗室8 d后的牦牛酸奶樣品提取總DNA，經質檢合格，通過樣品真菌ITS2 區測序，20 份牦牛酸奶樣品質控過濾后得到的有效序列為728022 條。采用稀釋曲線評估各樣品測序量是否足夠，結果如圖1 所示，由稀釋曲線可知，隨著序列數的增加，chao1 指數趨于平緩，說明本次樣本測序的深度足夠，每個樣品隨機抽取31495 條reads，抽平足以反映各樣品中真菌物種多樣性，測序量滿足后續生物信息學分析的要求。各個樣品及分組信息，質控后序列數，可觀測的物種數量，chao1、shannon、simpson和coverage指數見表2，所有樣本的測序深度均大于0.99，說明樣本測序結果可以反映樣品的真實情況。所有序列按97%的相似度進行OTU聚類，從20 份樣品中共檢測到702 個OTU，隸屬于6 門、24 綱、61 目、126 科、198 屬。

圖1 所有樣品稀釋曲線圖Fig.1 Rarefaction curve of all samples

表2 20 份牦牛酸奶樣品的序列信息及α 多樣性指數Table 2 Sequence information and α diversity of 20 yak yoghurt samples

2.2 不同采集保存溫度下樣品真菌群落結構的比較

2.2.1 組間真菌群落結構比較 Z2-20、Z4-20 樣品可觀測的物種數量分別為432、259，明顯高于其同一牦牛酸奶不同采集保存溫度下的其它樣品，但在干冰采集實驗室?20 ℃保存的其它樣品Z1-20、Z3-20、L2-20 可觀測的物種數量未出現明顯高于同一牦牛酸奶不同采集保存溫度下的其它樣品，故排除Z2-20、Z4-20 異常樣本，對不同采集保存溫度樣品組間α多樣性指數進行了kw檢驗分析，結果見表3，組間α多樣性各指數均無顯著差異（P>0.05）。對不同采集保存溫度樣品的真菌群落進行了主成分（PCA）分析，PCA可以初步反映出不同樣品可能表現出分散和聚集的分布情況，結果見圖2，常溫采集組（Z-0）、冰袋采集組（Z-4）、干冰采集組（Z-20）、液氮采集組（Z-80）的樣品分布均較為聚集，從而可以判斷出不同采集保存溫度樣品真菌群落組成差異不大，群落組成較為接近。

表3 組間樣品α 多樣性指數比較Table 3 Comparison of α diversity indices among groups

圖2 各組樣品的PCA分析Fig.2 PCA analysis of the different groups

2.2.2 不同采集保存溫度樣品門水平物種豐度分析 在真菌門水平，樣品中一共鑒定出6 個門水平物種，圖3 顯示不同采集保存溫度樣品組在門水平真菌相對豐度物種，主要5 個門，其余相對豐度太低的合并為other。常溫采集組的牦牛酸奶真菌ITS2 區測序，共檢測出4 個真菌菌門，其中子囊菌門（Ascomycota）平均豐度值為99.25%，擔子菌門（Basidiomycota）平均豐度值為0.73%，接合菌門（Zygomycota）平均豐度值為0.0052%，未分類菌門平均豐度值為0.010%。冰袋采集組的牦牛酸奶經高通量測序共檢測出5 個真菌菌門，其中子囊菌門（Ascomycota）平均豐度值為 98.39%，擔子菌門（Basidiomycota）平均豐度值為1.26%，接合菌門（Zygomycota）平均豐度值為0.043%，壺菌門（Chytridiomycota）平均豐度值為0.022%，未分類菌門平均豐度值為0.28%。干冰采集組的牦牛酸奶經高通量測序共檢測出5 個真菌菌門，其中子囊菌門（Ascomycota）平均豐度值為 96.04%，擔子菌門（Basidiomycota）平均豐度值為3.28%，接合菌門（Zygomycota）平均豐度值為 0.19%，壺菌門（Chytridiomycota）平均豐度值為0.057%，未分類菌門平均豐度值為0.44%。液氮采集組的牦牛酸奶經高通量測序共檢測出5 個真菌菌門，其中子囊菌門（Ascomycota）平均豐度值為 82.09%，擔子菌門（Basidiomycota）平均豐度值為16.90%，接合菌門（Zygomycota）平均豐度值為 0.36%，壺菌門（Chytridiomycota）平均豐度值為0.013%，未分類菌門平均豐度值為0.64%。ITS2 區測序分析的甘南牧區牦牛酸奶樣品中子囊菌門（Ascomycota）豐度占總數的82.09%～99.25%，為樣品中真菌優勢菌門。通過門水平真菌群落結構分析比較，不同采集保存溫度樣品組間真菌群落組成基本相同，僅組成比例不同。

圖3 門分類水平真菌相對豐度Fig.3 Relative abundance of fungi at phylum level

2.2.3 不同采集保存溫度樣品屬水平物種豐度分析 在真菌屬水平，圖4 顯示了不同采集保存溫度牦牛酸奶樣品在屬水平真菌相對豐度前20 的物種，其它物種合并為other。常溫采集組的牦牛酸奶相對豐度大于1%的菌屬包括酵母菌屬（Saccharomyces，63.49%）、克魯維酵母菌屬（Kluyveromyces，31.89%）、畢赤酵母屬（Pichia，3.39%）。冰袋采集組的牦牛酸奶相對豐度大于1%的菌屬包括酵母菌屬（Saccharomyces，61.59%）、克魯維酵母菌屬（Kluyveromyces，16.08%）、畢赤酵母屬（Pichia，10.01%）、假絲酵母屬（Candida，9.48%）。干冰采集組的牦牛酸奶相對豐度大于1%的菌屬包括酵母菌屬（Saccharomyces，23.98%）、克魯維酵母菌屬（Kluyveromyces，27.35%）、畢赤酵母屬（Pichia，35.21%）、Davidiella屬（3.12%）、馬拉色氏霉菌屬（Malassezia，1.76%）。液氮采集組的牦牛酸奶相對豐度大于1%的菌屬包括酵母菌屬（Saccharomyces，37.58%）、克魯維酵母菌屬（Kluyveromyces，22.68%）、畢赤酵母屬（Pichia，13.91%）、Davidiella屬（2.53%），馬拉色氏霉菌屬（Malassezia，12.28%）。測序分析的甘南牧區牦牛酸奶樣品中酵母菌屬（Saccharomyces）豐度占總數的23.98%～63.49%，克魯維酵母菌屬（Kluyveromyces）豐度占總數的16.08%～31.89%，畢赤酵母屬（Pichia）豐度占總數的3.39%～35.21%，是相對豐度較高的真菌屬，為樣品中的優勢菌屬。通過屬水平真菌群落結構分析比較，不同采集保存溫度下牦牛酸奶樣品組間真菌群落組成相似，僅組成比例不同。

圖4 屬水平前20 名真菌的相對豐度Fig.4 Relative abundance of the top 20 fungi at genus level

2.2.4 組間物種差異分析 屬水平的組間物種差異分析，采用秩和檢驗得到不同采集保存溫度下牦牛酸奶樣品組間的差異細菌屬共2 個，其中液氮采集組，Davidiella、Hypocreales unidentified2 個菌屬顯著高于其它組別的真菌群落菌屬（P<0.05），結果見表4。為了直觀展示各組及各樣品間真菌屬差異，對差異真菌屬進行熱圖（Heatmap）分析和主成分分析（PCA），Heatmap圖是以顏色梯度來表征二維矩陣或表格中的數據大小，反映不同分組或樣本在各分類水平上群落組成的相似性與差異性，由圖5 可知，液氮采集組樣品的群落組成和其它樣品組差異較大。PCA圖能夠直觀地把樣品之間相似或差異的程度呈現在二維坐標系上，在圖6 中常溫采集組和冰袋采集組樣品表現出聚集的分布，表明常溫、冰袋采集組樣品的群落組成更為接近；干冰采集組和液氮采集組樣品表現出聚集的分布，表明干冰、液氮采集組樣品的群落更為接近。

圖5 屬水平差異物種的HeatmapFig.5 Heatmap of different genus among groups

圖6 屬水平差異物種的PCA分析Fig.6 PCA analysis of different genus among groups

表4 差異顯著物種Table 4 Significantly different genus among groups

3 討論與結論

牦牛酸奶是由牦牛奶通過傳統發酵方法制得，是經過長期自然馴化優良食品微生物的天然載體，含有極具開發利用價值的微生物資源[21]。本研究采用Illumina NovaSeq高通量測序技術測序甘南牧區牦牛酸奶樣品真菌ITS2 區，分析牦牛酸奶中真菌菌群組成，比較不同采集保存溫度樣品真菌多樣性，對20 份甘南牧區牦牛酸奶樣品高通量測序結果分析表明，子囊菌門（Ascomycota）為樣品中真菌優勢菌群門，酵母菌屬（Saccharomyces）、克魯維酵母菌屬（Kluyveromyces）、畢赤酵母屬（Pichia）為優勢菌屬。Zhang等[22]報道青海地區牦牛酸奶與其它發酵酸奶相比含有更多的乳酸菌和酵母菌，其中青海西北部的高原牦牛酸奶和青海東部的環湖牦牛酸奶中乳酸菌是優勢菌，而青海南部的環湖牦牛酸奶中酵母菌是優勢菌。Wu等[23]研究發現，3 種乳酸菌以及5 種酵母菌是西藏不同海拔地區牦牛酸奶中主要發酵菌群。孫志宏等[24]運用宏基因組學技術對西藏地區自然發酵乳中微生物多樣性進行分析，結果表明在真菌門水平上子囊菌門（Ascomycota）為最優勢菌門，在屬水平上優勢菌屬為畢赤酵母屬（Pichia）、耐堿酵母屬（Galactomyces）、酵母菌屬（Saccharomyces）。烏仁圖雅[25]、張冬蕾[26]、楊彥榮[27]通過焦磷酸測序技術對傳統發酵乳制品酸牦牛奶、酸牛奶進行了真菌多樣性研究，結果表明，酸牦牛奶、酸牛奶中真菌優勢菌門均為子囊菌門。從以上研究結果中可以看出牦牛乳制品中的真菌菌群盡管比較復雜，但是都以酵母菌為主，只是優勢的酵母菌的種類有所不同，本次甘南牧區牦牛酸奶樣品真菌菌群組成分析結果與上述學者研究結論較為一致。

采用Illumina高通量測序技術分析比較了4 種不同采集保存溫度下牦牛酸奶樣品真菌多樣性。目前，尚未見不同采集保存溫度對傳統發酵乳制品微生物多樣性分析影響的報道，呼斯楞等[28]在內蒙古呼倫貝爾地區傳統發酵乳中乳酸菌的多樣性分析中，樣品采集采用低溫保存，帶回實驗室于?80 ℃超低溫保藏；陳芝蘭等[29]在西藏地區傳統發酵乳中乳酸菌多樣性及微生物數量分析中，將牧民自然發酵制作的酸奶放入冰盒內保存；劉振東等[30]在西藏不同產區曲拉細菌群落結構的比較分析中將采集樣品放入冰盒保存并立即運送到實驗室，置于?80 ℃冰箱中備用；張敏等[31]在基于16S rDNA高通量測序方法比較新疆西北部地區乳品中微生物的多樣性中，將采集樣品放入車載冰箱運回實驗室，在?80 ℃超低溫冰箱保存；馬江等[32]在甘南牦牛曲拉中真菌群落結構研究中，將樣品置于冷藏箱中，運輸至實驗室。以上乳制品微生物多樣性研究采樣均以低溫采集為主。在未來研究中，為了更大可能的獲得牦牛酸奶中不同種類的微生物多樣性信息和定向的分離出目標微生物，可在借鑒已有最新測序技術在傳統發酵乳制品微生物基因組和生物多樣性中的研究結果，進一步優化不同的采樣方法與分離策略，更好地開展傳統發酵乳制品中優良微生物資源的保護與發掘研究工作。

采用Illumina NovaSeq高通量測序技術，對4 種不同保存溫度條件下采集自甘南牧區20 份牦牛酸奶樣品真菌ITS2 區進行測序，經OTU聚類分析，采集的甘南牧區牦牛酸奶樣品中的真菌分屬于6 門、24 綱、61 目、126 科、198 屬。子囊菌門（Ascomycota）在采集牦牛酸奶樣本真菌中占比最高，其豐度值變化范圍為82.09%～99.25%，是20 份樣品中的真菌優勢菌門。酵母菌屬（Saccharomyces）豐度值變化范圍為23.98%～63.49%，克魯維酵母菌屬（Kluyveromyces）豐度值變化范圍為16.08%～31.89%，畢赤酵母屬（Pichia）豐度值變化范圍為3.39%～35.21%，是樣品中的真菌優勢菌屬。基于Illumina NovaSeq高通量測序分析比較，常溫、冰袋、干冰和液氮4 種不同保存溫度條件下采集自甘南牧區牦牛酸奶真菌多樣性之間無顯著差異，真菌群落組成基本相同，僅組成比例不同。