謝芳，謝華德，李孟偉，彭開屏，楊承劍

(中國農業(yè)科學院廣西水牛研究所廣西水牛遺傳繁育重點實驗室，南寧530001)

0 引 言

中國是世界第三大水牛養(yǎng)殖大國，其水牛存欄數(shù)僅次于印度和巴基斯坦。近年來，在市場需求和政策扶持下，奶水牛產(chǎn)業(yè)已經(jīng)成為廣西的優(yōu)勢特色產(chǎn)業(yè)[1]，水牛奶產(chǎn)量也位居全國之最。水牛奶中的蛋白質、脂肪、乳糖、維生素等含量均高于黑白花奶，其免疫球蛋白G、溶菌酶、類胰島素生長因子等更是普通牛乳的5～10倍[2]。與荷斯坦奶牛相比，水牛適應能力強，具有較強的抗病性和免疫力，迄今尚未有水牛瘋牛病例報道[3]，因此，水牛奶具有較好的食用安全性和較高的商業(yè)價值。廣西屬亞熱帶季風氣候區(qū)，年平均氣溫在23℃左右，年平均降雨量達1 300 mm，正是這種長年濕熱的氣候特點孕育了口感醇香獨特的水牛奶。水牛奶也因其豐富的營養(yǎng)而成為微生物生長的溫床[4]，而在自然條件下，溫度和濕度則是影響微生物生長的主要因素，季節(jié)起著決定作用[5]。然目前，尚未見有關廣西地區(qū)水牛乳中微生物組成及水牛乳品質受季節(jié)變化影響的報道。

近年來，隨著分子生物學的發(fā)展，高通量測序作為“下一代”測序技術已迅速發(fā)展，它避開了微生物分離培養(yǎng)的過程，能檢測到純培養(yǎng)和非培養(yǎng)不能發(fā)現(xiàn)的低豐度細菌種類，相較傳統(tǒng)測序技術，具有準確率高、耗時短以及信息處理量大等優(yōu)點[6]。目前，高通量測序已應用于各個領域，包括海洋、土壤、植物和腸道環(huán)境，而利用高通量測序對乳品進行研究的報道較少[7]，以往通過傳統(tǒng)測序可以檢測到可培養(yǎng)和含量較高的菌落，但很難從樣品中獲得曾經(jīng)存活或難以分離的微生物，進而不能全面了解微生物的多樣性[8]。因此，本研究擬利用16SrDNA高通量測序技術，對廣西地區(qū)冬、夏兩季水牛乳中的細菌群落進行比較分析，以期全面了解冬、夏兩季水牛乳中菌落結構和多樣性的差異，旨在為后續(xù)水牛乳加工、安全評估等提供依據(jù)。

1 材料與方法

1.1 材料與儀器

E.Z.N.A.?Soil DNA Kit DNA抽提試劑盒，美國Omega Bio-Tek公司；biowest agArose瓊脂糖，西班牙Biowest；NEXTFLEX?Rapid DNA-Seq Kit建庫試劑盒，美國Bioo Scientific；MiSeq Reagent Kit v3測序試劑盒，美國Illumina；Eppendorf 5430 R離心機、Eppendorf 5424R高速臺式冷凍離心機，德國Eppendorf；BioTek ELx800酶標儀，美國Biotek；QuantusTMFluorometer微型熒光計，美國Promega公司；DYY-6C電泳儀，北京市六一儀器廠；ABI GeneAmp?9700型PCR儀，美國ABI；Illumina M iseq測序儀，美國Illumina。

1.2 樣本采集

水牛乳樣本均取自廣西水牛研究所水牛種畜場，隨機選取10頭健康三品種高代雜交水牛（摩拉×尼里-拉菲×本地沼澤型水牛），于7月中旬（夏季）和12月下旬（冬季）采集樣品，每頭?，F(xiàn)場各采集3管平行，運回實驗室后將3管奶樣混勻為1管混合奶樣，冬、夏季各10份混合樣本，共20個樣本，分別定義為夏季組（X）和冬季組（W），編號為F1、F2、F3、F4、F5、F6、F7、F8、F9、F10和T1、T2、T 3、T 4、T 5、T6、T7、T8、T 9、T10。采樣前用溫水清洗乳頭，人工佩戴無菌手套擠奶，樣本經(jīng)50 mL無菌PC離心管密封好，置于采樣冷藏保溫箱迅速帶回實驗室，于30 min內進行生乳中細菌總脫氧核糖核酸（DNA）的粗提取。

1.3 總DNA粗提和PCR擴增

取40 mL生水牛乳樣本，添加8 mL乙二胺四乙酸二鈉（EDTA）（500 mmol/L，p H 8.0），在6 500 r/min下4℃冷凍離心30 min，用滅菌脫脂棉去掉乳脂層，棄上清液，加入磷酸緩沖鹽溶液（PBS）至10 mL，12 000 r/min高速冷凍離心10 min，用1 mL蔗糖緩沖液沖洗兩次（每次12 000 r/min，離心10 min），后懸浮400μL蔗糖緩沖液及2U變溶菌素和800μg溶菌酶，37℃培養(yǎng)1 h；取上述培養(yǎng)后的樣液，分別加入5μL SDS，12μL EDTA，5μL蛋白酶K及78μL蔗糖緩沖液，55℃孵育1 h，再加入65μL蔗糖緩沖液和135μL氯化鈉（5 mmol/L），最后加入2μL RNase于37℃孵育30 min，用700μL苯酚：氯仿：異戊醇（25∶24∶1）消除蛋白，上下?lián)u勻10次以上，然后18 000×g、4℃離心15 min提取上清，重復3次，最后加入700μL冰異丙醇（(-20℃預冷）和20μg糖原（glycogen,Roche Diagnostic），混勻放置5 min，然后-20℃沉淀過夜。取出過夜樣本，18 000×g離心30 min獲取沉淀，用70%冰乙醇清洗兩遍（500～750μL，彈起白色沉淀，13 000 r/min離心10～15 min），傾倒上清，將離心管倒放，風干DNA，加入50μL PE buffer,70℃孵育5 min，用Nano-Drop2000測定DNA濃度。

用NanoDrop2000檢測DNA純度和濃度，為避免序列太長影響測序，同時兼顧擴增特異性，選用338F（5’-ACTCCTACGGGAGGCAGCAG-3’）和806R（5’-GGACTACHVGGGTWTCTAAT-3’）對16SrDNA基因V3-V4可變區(qū)進行PCR擴增，反應條件如下：95℃預變性3 min，27個循環(huán)（95℃變性30 s，55℃退火30 s,72℃延伸30 s)，然后72℃穩(wěn)定延伸10 min，最后在4℃進行保存[9]。聚合酶鏈式反應（PCR）反應體系為：5×TransStart FastPfu緩沖液4μL，2.5mmol/L d NTPs 2μL，上游引物（5μmol/L）0.8μL，下游引物（5μmol/L）0.8μL,TransStart FastPfu DNA聚合酶0.4μL，模板DNA 10 ng，補足至20μL。每個樣本3個PCR重復，將3個重復的PCR產(chǎn)物混合，使用2%瓊脂糖凝膠電泳檢測產(chǎn)物。

1.4 Illumina Miseq測序

使用2%瓊脂糖凝膠回收PCR產(chǎn)物，用AxyPrep DNA Gel Extraction Kit(AxygenBiosciences,Union City,CA,USA)進行純化，Tris-HCl洗脫，2%瓊脂糖電泳檢測。利用QuantiFluorTM-ST(Promega,USA)進行檢測定量。根據(jù)Illumina MiSeq平臺(Illumina,SanDiego,USA)標準操作規(guī)程，將純化后的擴增片段構建PE 2*300文庫。利用Illumina公司的Miseq PE300平臺進行測序。原始數(shù)據(jù)上傳至美國國家生物信息中心(NCBI)數(shù)據(jù)庫進行比對分析。

1.5 數(shù)據(jù)處理

使用Trimmomatic軟件原始測序序列進行質控，使用FLASH軟件進行拼接：過濾reads尾部質量值20以下的堿基，設置50 bp的窗口，如果窗口內的平均質量值低于20，從窗口開始截去后端堿基，過濾質控后50 bp以下的reads，去除含N堿基的reads。使用的UPARSE軟件（version 7.1 http://drive5.com/uparse/)，根據(jù)97%的相似度對序列進行分類操作單元(OUT)聚類并剔除嵌合體。利用RDP classifier(http://rdp.cme.msu.edu/)對每條序列進行物種分類注釋，比對Silva數(shù)據(jù)庫（SSU 128），設置比對閾值為70%。菌群多樣性分析及繪圖：均采用R語言進行統(tǒng)計分析和作圖。

2 結果與分析

2.1 測序數(shù)據(jù)、OTU及稀釋曲線分析

通過對冬季組和常夏季組共計20個樣本進行16SrDNA測序、數(shù)據(jù)質控后，共得到有效序列690880條，冬季組和夏季組分別得到優(yōu)化序列條數(shù)為345440條。兩組樣本的優(yōu)化序列長度均主要分布在401～440 bp，集中在421～440 bp，與設計引物擴增長度300 bp左右接近，可滿足序列分析要求。

圖1 樣本有效序列長度分布

在97%相似性水平下對兩組所有樣本序列進行聚類分析[10]，結果見圖2，冬季組和夏季組共注釋出1126個OTU，兩組共有的OTU為502，冬季組共產(chǎn)生1105個OUT，夏季組共產(chǎn)生523個OTU，各自獨有的OTU分別為603和21。由此可見，冬季組OTU數(shù)遠高于夏季組，兩組水牛乳樣本間存在大量共有的細菌菌群，且冬季組中細菌多樣性高于夏季組。

圖2 樣本OTU分布韋恩圖

稀釋性曲線主要反映測序深度，根據(jù)曲線是否達到平緩來判斷本次測序數(shù)量是否足夠，也間接反映出樣本中菌群的豐富度和多樣性[11]。Sobs和Shannon指數(shù)分別用于評估樣本群落的豐富度和多樣性[12]。兩組樣本基于OTU的Sobs和Shannon稀釋曲線結果如圖3所示。各樣本的群落豐富度Sobs指數(shù)稀釋曲線大部分還在上升期，但群落多樣性Shannon指數(shù)稀釋曲線在2000條時就已經(jīng)平緩，說明增加測序數(shù)量可能會發(fā)現(xiàn)少量新的OTU，但樣本群落多樣性不會增加[13]，說明本次測序數(shù)量可滿足后續(xù)分析要求。

圖3 Sobs（a）、Shannon（b）指數(shù)稀釋性曲線

2.2 細菌組成與豐度差異比較

基于OTU注釋結果，分析樣本在不同分類水平上的群落組成。本次測序，Coverage（覆蓋率）為0.999，可代表99.9%的樣本信息。冬季組樣本共得到19個細菌門，92個目，185個科，456個細菌屬和733個細菌種。而夏季組樣本共得到14個細菌個門，63個目，127個科，274個細菌屬和424個菌種，相較而言，冬季組的細菌多樣性顯著高于夏季組。為可視化呈現(xiàn)各樣本的物種豐度信息，本研究采用柱形圖來展示各樣本在門水平和屬水平的菌落分布及差異性[14]，結果如圖4、圖5所示。由圖4可知，在門分類水平上，以相對豐度≥0.01%為閾值，將＜0.01%的菌門歸為others，初乳和常乳組共得到6個細菌門，從高到低依次為變形菌門（Proteobacteria）、厚壁菌門（Firmicutes）、擬桿菌門（Bacteroidetes）、放線菌門（Actinobacteria）、Patescibacteria、Deinococcota。其 中 厚 壁 菌 門（Firmicutes）、放線菌門（Actinobacteria）、（Deinococcota）在冬季組中占比較多，而變形菌門（Proteobacteria）、擬桿菌門（Bacteroidetes）、（Patescibacteria）則在夏季組中較占優(yōu)勢，同一菌門在不同樣本中的相對豐度差異較大。變形菌門、厚壁菌門為兩組樣本中的絕對優(yōu)勢菌門，相對豐度分別為55.4%和22.9%，合計達78.3%。

圖4 門水平菌落豐度圖

在屬水平上，兩組樣本豐度≥0.01的細菌屬共有33個（圖5）。其中，冬季組檢測到豐度≥2%的有12個細菌屬，從高到低依次為：不動桿菌屬（Acinetobacter）、芽孢桿菌屬（Bacillus）、棒桿菌屬（Corynebacterium）、庫特氏菌屬（Kurthia）、土地桿菌屬（Pedobacter）、腸桿菌屬（Enterobacter）、金黃桿菌屬（Chryseobacterium）、醋桿菌屬（Acetobacter）、棲水菌屬（Enhydrobacter）、微桿菌屬（Microbacterium）、假單胞菌屬（Pseudomonas）、克雷伯氏菌屬（Klebsiella）；而夏季組豐度≥2%的細菌屬有7個，即：不動桿菌屬（Acinetobacter）、乳球菌屬（Lactococcus）、金黃桿菌屬（Chryseobacterium）、芽孢桿菌屬（Bacillus）、伊麗莎白菌屬（Elizabethkingia）、TM 7a菌屬、棲水菌屬（Enhydrobacter）。

圖5 屬水平菌落豐度圖

將冬季和夏季兩組樣本中豐度≥0.01%的相同菌屬豐度占比求和后進行比較，結果如表1所示，夏季組中的優(yōu)勢菌屬為：不動桿菌屬（Acinetobacter）、乳球菌屬（Lactococcus）、金黃桿菌屬（Chryseobacterium），而冬季組的優(yōu)勢菌屬則為：不動桿菌屬（Acinetobacter）、芽孢桿菌屬（Bacillus）、棒桿菌屬（Corynebacterium）、庫特氏菌屬（Kurthia）、土地桿菌屬（Pedobacter）。由此可知，冬、夏兩季水牛乳中的優(yōu)勢菌群差異顯著，相同菌屬在不同樣本中的占比均不同，冬季組樣本的菌群多樣性顯著高于夏季組。

表1 冬夏兩季水牛乳中相同菌屬及其相對豐度

在屬分類水平，取兩組樣本平均總豐度前30的菌屬，按其物種和樣本層級聚類方式繪制群落豐度聚類熱圖（圖6），使高豐度和低豐度菌落分區(qū)聚集[15]，通過顏色梯度變化來反映不同分組樣本在屬水平群落組成的相似性和差異性。由圖6可知，夏季組和冬季組樣本中優(yōu)勢菌群的豐度呈現(xiàn)明顯不同，顏色梯度交替出現(xiàn)，整體沒有分區(qū)聚集趨勢，不動桿菌屬（Acinetobacter）、乳球菌屬（Lactococcus）、金黃桿菌屬（Chryseobacterium）在夏季組樣本中的含量顯著高于冬季組，而芽孢桿菌屬（Bacillus）、棒桿菌屬（Corynebacterium）、庫特氏菌屬（Kurthia）、土地桿菌屬（Pedobacter）在冬季組中的豐度均遠高于夏季組。不動桿菌（Acinetobacter）在兩組樣本中均占比較高，葡萄球菌屬（Staphylococcus）在兩組樣本中均有低豐度出現(xiàn)。

圖6 樣本基于屬水平的菌落聚類熱圖

2.3 細菌多樣性比較分析

本試驗采用樣本層級聚類樹來分析水牛乳冬季和夏季組間群落組成的相似性。20個樣本基于OTU的樣本層級聚類樹如圖7所示，通過樣本聚類樹可以清楚地看出樣本間分支距離遠近，較直觀地反映出各樣本間菌群的相似度和差異性[16-17]。由圖7可知，夏季組樣本（紅色）主要聚在一大類分支上，相似度較高，其中F8和F9是最相似的兩個樣本。冬季組（藍色）樣本則較分散，主要分成三大類，其中T1、T9、T 10、T7這4個樣本聚成一類，與大部分夏季組樣本同聚在一大分支上，而T 8、T 5、T 6則與T3、T4分別單獨聚成兩類分支，且這兩層級樣本完全不在一個分支上，樣本間相似度差異較大。整體而言，冬季組與夏季組樣本組間差異不顯著，夏季組樣本組內間菌落組成相似度較高，而冬季組樣本組內間菌落組成差異較大。

圖7 樣本層級聚類分析

對夏季與冬季組樣本中細菌屬進行主成分（PCA）分析，以顯示兩組樣本間菌群組成的相似性和差異性。樣本通過降維后，在主成分顯示的相對坐標點越接近，表明兩樣本物種組成越相似[18]，同理，各點之間距離越大，表明它們之間的菌群差異越大[19]，由圖8可知，夏季組與冬季組兩組樣本點之間的距離較近，在空間上沒有完全分開，這說明冬季組與夏季組組間樣本菌群結構相似，差異不顯著。主成分1（PC1）和主成分2（PC2）對樣本差異性的貢獻分別達到了51.79%和19.43%，共計71.22%，可以代表大部分變量信息。由圖8可知，夏季組與冬季組在第一主成分PC1所在的X軸方向，兩組樣本的菌落組成沒有明顯的分離趨勢，這說明PC1因素有較高的可能性影響了樣本的組成，則表示結果數(shù)據(jù)組內重復性好，組內樣本差異性大，且PC1貢獻值為51.79%，可以作為有效的主成分，來解釋兩組樣本的組間差異。

圖8 樣本基于屬水平的主成分分析

PLS-DA（偏最小二乘法判別分析Partial Least Squares Discriminant Analysis）是一種常見的檢驗樣本組間差異的方法，此方法主要根據(jù)分組方案，弱化組內樣本差異，突出組間差異[20]。由PLS-DA（圖9）顯示，除去離群嚴重的異常樣本T 3，其余樣本點均按夏季與冬季所代表的分組相對聚集，兩組樣本點在COMP1和COMP2方位很好的區(qū)分開來，這說明夏季與冬季兩組樣本在組間的細菌組成存在明顯差異，且夏季組代表的樣本點都比較集中，說明夏季組組內樣本間菌群組成較相似，而冬季組所代表的樣本點都比較分散，這說明冬季組組內樣本菌群組成有一定的差異性。結合PCA和PLS-DA分析可知，冬季組與夏季組組間差異顯著，組內差異不顯著。

圖9 樣本基于屬水平的PLS-DA分析

3 結 論

本試驗采用16SrDNA高通量測序技術，對廣西水牛研究所種畜廠存欄最多的三品雜水牛在冬、夏兩季水牛乳樣本中的細菌多樣性進行了比較分析，結果顯示：冬季組樣本共得到19個細菌門，456個屬和733個種。而夏季組樣本共得到14個門，274個屬和424個種，相較而言，冬季組的細菌多樣性顯著高于夏季組。兩組樣本優(yōu)勢菌群組成及豐度差異顯著，組內差異小，組間差異大。隨著氣溫升高，夏季組中的乳球菌屬（Lactococcus）、金黃桿菌屬（Chryseobacterium）占優(yōu)勢，而氣溫下降后，冬季組中的芽孢桿菌屬（Bacillus）、庫特氏菌屬（Kurthia）和土地桿菌屬（Pedobacter）則成為優(yōu)勢菌群，這說明水牛乳中細菌多樣性與季節(jié)氣候因素密切相關，這與大部分國外研究結果不同，究其原因，可能是因為廣西地區(qū)冬少夏多，長年濕熱，季候區(qū)分不明顯所致。通過本次測序，全面掌握了水牛研究所種畜場三品雜水牛乳中的益生菌與致病群分布狀況，這對后續(xù)改善該種畜場飼養(yǎng)條件、奶源衛(wèi)生以及奶水牛乳腺炎的前期防治等均有指導意義。