基于16S rDNA高通量測(cè)序分析水牛初乳與常乳中細(xì)菌多樣性

謝 芳，謝華德，唐振華，謝耀鋒，郭艷霞，黃鈺涵，楊承劍

（中國(guó)農(nóng)業(yè)科學(xué)院廣西水牛研究所，廣西水牛遺傳繁育重點(diǎn)實(shí)驗(yàn)室，廣西南寧 530001）

我國(guó)水牛資源豐富，廣西水牛數(shù)量更是位居全國(guó)首位。近年來(lái)，在政府扶植和政策驅(qū)動(dòng)下，奶水牛產(chǎn)業(yè)已經(jīng)成為廣西的特色優(yōu)勢(shì)產(chǎn)業(yè)[1]，水牛奶產(chǎn)量也升至全國(guó)之最。水牛奶因其優(yōu)良的品質(zhì)和香醇的口感而被世界公認(rèn)為“奶中珍品”[2]，其蛋白質(zhì)、脂肪、乳糖、維生素等含量均高于黑白花奶[3]與荷斯坦奶牛相比，水牛適應(yīng)能力強(qiáng)，具有較強(qiáng)的抗病性和免疫力，迄今尚未有水牛瘋牛病例報(bào)道[4]，尤其是水牛初乳，因其含有大量的免疫和生長(zhǎng)因子而一度成為市場(chǎng)新寵，更有研究表明水牛初乳中的免疫球蛋白G、溶菌酶、類胰島素生長(zhǎng)因子等是普通牛乳的5～10 倍[5]，因此，水牛乳具有較好的食用安全性和較高的商業(yè)價(jià)值，深加工優(yōu)勢(shì)明顯。

近年來(lái)，隨著分子生物學(xué)的發(fā)展，高通量測(cè)序作為“下一代”測(cè)序技術(shù)已迅速發(fā)展，它避開了微生物分離培養(yǎng)的過(guò)程，能檢測(cè)到純培養(yǎng)和非培養(yǎng)技術(shù)不能發(fā)現(xiàn)的低豐度細(xì)菌種類，相較傳統(tǒng)的測(cè)序技術(shù)，具有準(zhǔn)確率高、耗時(shí)短以及信息處理量大等優(yōu)點(diǎn)[6]。目前，高通量測(cè)序技術(shù)現(xiàn)已應(yīng)用于各個(gè)領(lǐng)域，包括海洋、土壤、植物和腸道環(huán)境，而利用高通量測(cè)序?qū)θ槠奉I(lǐng)域進(jìn)行研究的報(bào)道較少[7]，以往通過(guò)傳統(tǒng)測(cè)序可以檢測(cè)到可培養(yǎng)和含量較高的菌落，但很難從樣品中獲得曾經(jīng)存活或難以分離的微生物，進(jìn)而不能全面了解微生物的多樣性[8]。牛乳是微生物最好的培養(yǎng)基，這些微生物可能來(lái)自患乳房炎奶牛的乳房、擠奶工具亦或是牛棚環(huán)境[9]，而生奶中致病細(xì)菌以及腐敗細(xì)菌的分布可反映奶水牛的健康狀況、水牛奶的污染情況等[10]。作為生產(chǎn)乳品的基礎(chǔ)，其微生物含量則是衡量水牛乳品質(zhì)的關(guān)鍵指標(biāo)，故了解水牛原料乳中微生物多樣性至關(guān)重要。

本研究擬利用16S rDNA高通量測(cè)序技術(shù)，對(duì)廣西水牛研究所種畜場(chǎng)存欄量較大的三品雜水牛（摩拉×尼里-拉菲×本地沼澤型）初乳和常乳組樣本中的細(xì)菌群落結(jié)構(gòu)進(jìn)行比較分析，以期全面了解水牛初乳和常乳中細(xì)菌群落結(jié)構(gòu)和多樣性差異，旨在為后續(xù)水牛乳加工、安全評(píng)估等提供參考。

1 材料與方法

1.1 材料與儀器

生水牛乳樣本 均取自廣西水牛研究所水牛種畜場(chǎng)；E.Z.N.A.? Soil DNA Kit DNA抽提試劑盒，美 國(guó)Omega Bio-Tek公 司；biowest agArose瓊 脂糖 西班牙Biowest；NEXTFLEX ? Rapid DNASeq Kit建庫(kù)試劑盒 美國(guó)Bioo Scientific；MiSeq Reagent Kit v3 測(cè)序試劑盒 美國(guó)Illumina；磷酸緩沖鹽溶液 賽默飛世爾科技；蔗糖緩沖液 北京雷根生物；溶菌酶、變?nèi)芫?上海聯(lián)碩寶為生物；蛋白酶K、RNase、PE Buffer上海生工生物；氯化鈉、冰異丙醇、苯酚、氯仿、異戊醇 上海北諾生物科技；糖原 北京百奧萊博。

Eppendorf 5430 R離心機(jī)、Eppendorf 5424R高速臺(tái)式冷凍離心機(jī) 德國(guó)Eppendorf；NanoDrop2000超微量分光光度計(jì) 美國(guó)Thermo Fisher Scientific；BioTek ELx800 酶標(biāo)儀美國(guó)Biotek；Quantus?Fluorometer微型熒光計(jì) 美國(guó)Promega公司；DYY-6C電泳儀 北京市六一儀器廠；ABI GeneAmp?9700 型PCR儀 美國(guó)ABI；Illumina Miseq測(cè)序儀 美國(guó)Illumina。

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

1.2.1 樣本采集 生水牛乳樣本隨機(jī)選取9 頭處于圍產(chǎn)期的三品種高代雜交水牛（摩拉×尼里-拉菲×本地沼澤型水牛），于每頭牛產(chǎn)犢后的第1 d開始取樣，連續(xù)采集3 d，每頭牛采集2 管平行，將每頭牛采集到的6 管奶樣混勻?yàn)榛旌夏虡樱x為初乳（CR4），共9 個(gè)樣本，編號(hào)為C1、C2、C3、C4、C5、S2、S3、S4、S5。于每頭牛產(chǎn)犢后的第7 d開始采集奶樣，連續(xù)3 d，每頭采集2 管平行，同樣將6 管奶樣混勻成混合奶樣，定義為常乳（PR），共9 個(gè)樣本，編號(hào)為M1、M2、M3、M4、M5、M6、M7、M8、M9。采樣前用溫水清洗乳頭，人工佩戴無(wú)菌手套擠奶，去掉前數(shù)把奶后再留樣，樣本經(jīng)50 mL無(wú)菌PC離心管密封好，置于采樣冷藏保溫箱迅速帶回實(shí)驗(yàn)室，于30 min內(nèi)進(jìn)行生乳中細(xì)菌總脫氧核糖核酸（DNA）的粗提取。

1.2.2 總DNA粗提和PCR擴(kuò)增 取40 mL生水牛乳樣本，添加8 mL乙二胺四乙酸二鈉（EDTA）（500 mmol/L，pH8.0），在6500 r/min下4 ℃冷凍離心30 min，用滅菌脫脂棉去掉乳脂層，棄上清液，加入磷酸緩沖鹽溶液（PBS）至10 mL，12000 r/min高速冷凍離心10 min，用1 mL蔗糖緩沖液沖洗兩次（每次12000 r/min，離心10 min），后懸浮400 μL蔗糖緩沖液及2U變?nèi)芫睾?00 μ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 次以上，然后18000×g、4 ℃離心15 min提取上清，重復(fù)3 次，最后加入700 μL冰異丙醇（?20 ℃預(yù)冷）和20 μg糖原（glycogen,Roche Diagnostic），混勻放置5 min，然后-20 ℃沉淀過(guò)夜。取出過(guò)夜樣本，18000×g離心30 min獲取沉淀，用70%冰乙醇清洗兩遍（500～750 μL，彈起白色沉淀，13000 r/min離心10～15 min），傾倒上清，將離心管倒放，風(fēng)干DNA，加入50 μL PE Buffer,70 ℃孵育5 min，用NanoDrop2000 測(cè)定DNA濃度。

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

1.2.3 Illumina Miseq測(cè)序 使用2%瓊脂糖凝膠回收PCR產(chǎn)物，用AxyPrep DNA Gel Extraction Kit進(jìn)行純化，Tris-HCl洗脫，2%瓊脂糖電泳檢測(cè)。利用QuantiFluor?-ST進(jìn)行檢測(cè)定量。根據(jù)Illumina MiSeq平臺(tái)標(biāo)準(zhǔn)操作規(guī)程，將純化后的擴(kuò)增片段構(gòu)建PE 2*300 文庫(kù)。文庫(kù)構(gòu)建步驟：a.接“Y”字形接頭；b.使用磁珠篩選去除接頭自連片段；c.利用PCR擴(kuò)增進(jìn)行文庫(kù)模板的富集；d.磁珠回收PCR產(chǎn)物得到最終的文庫(kù)。利用Illumina公司的Miseq PE300 平臺(tái)進(jìn)行測(cè)序（上海美吉生物醫(yī)藥科技有限公司）。原始數(shù)據(jù)上傳至美國(guó)國(guó)家生物信息中心（NCBI）數(shù)據(jù)庫(kù)進(jìn)行比對(duì)分析。

1.2.4 細(xì)菌組成與豐度差異分析

1.2.4.1 稀釋曲線分析 構(gòu)建方法是統(tǒng)計(jì)每個(gè)樣本中，每個(gè)OTU所含的序列數(shù)，將OTUs按豐度（所含有的序列條數(shù)）由大到小等級(jí)排序，再以O(shè)TU的排序等級(jí)為橫坐標(biāo)，以每個(gè)OTU中所含的序列數(shù)（也可用OTU中序列數(shù)的相對(duì)百分含量）為縱坐標(biāo)作圖。

1.2.4.2 OUT、韋恩圖分析 選用相似水平為97%的OTU對(duì)序列進(jìn)行分類操作單元（OUT）聚類。軟件:采用R語(yǔ)言工具統(tǒng)計(jì)和作圖。

1.2.4.3 群落柱形圖（Bar圖）直觀呈現(xiàn)各樣本在（如域、界、門、綱、目、科、屬、種、OTU等）分類學(xué)水平上含有哪些優(yōu)勢(shì)物種以及樣本中各優(yōu)勢(shì)物種的相對(duì)豐度（所占比重）。軟件：基于tax_summary_a文件夾中的數(shù)據(jù)表，利用R語(yǔ)言工具作圖。

1.2.5 細(xì)菌多樣性比較分析

1.2.5.1 Heatmap圖 根據(jù)物種或樣本間豐度的相似性進(jìn)行聚類，并將結(jié)果呈現(xiàn)在群落Heatmap圖上，可使高豐度和低豐度的物種分塊聚集，通過(guò)顏色變化來(lái)反映不同分組（或樣本）在各分類學(xué)水平上群落組成的相似性和差異性。軟件及算法：R語(yǔ)言vegan包。

1.2.5.2 主成分分析（PCA分析）通過(guò)分析不同樣本群落組成可以反映樣本間的差異和距離，PCA運(yùn)用方差分解，將多組數(shù)據(jù)的差異反映在二維坐標(biāo)圖上，坐標(biāo)軸取能夠最大反映樣品間差異的兩個(gè)特征值。如樣本物種組成越相似，反映在PCA圖中的距離越近。軟件：R語(yǔ)言PCA統(tǒng)計(jì)分析和作圖。

1.2.5.3 偏最小二乘法判別分析（PLS-DA）是多變量數(shù)據(jù)分析技術(shù)中的判別分析法。通過(guò)對(duì)主成分適當(dāng)?shù)男D(zhuǎn)，PLS-DA可以有效地對(duì)組間觀察值進(jìn)行區(qū)分，并且能夠找到導(dǎo)致組間區(qū)別的影響變量。這種分析方法有利于發(fā)現(xiàn)組間的異同點(diǎn)。軟件：R語(yǔ)言mixOmics包中plsda分析和作圖。

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

使用Trimmomatic軟件原始測(cè)序序列進(jìn)行質(zhì)控，使用FLASH軟件進(jìn)行拼接：過(guò)濾reads尾部質(zhì)量值20 以下的堿基，設(shè)置50 bp的窗口，如果窗口內(nèi)的平均質(zhì)量值低于20，從窗口開始截去后端堿基，過(guò)濾質(zhì)控后50 bp以下的reads，去除含N堿基的reads；根據(jù)PE reads之間的overlap關(guān)系，將成對(duì)reads拼接（merge）成一條序列，最小overlap長(zhǎng)度為10 bp；拼接序列的overlap區(qū)允許的最大錯(cuò)配比率為0.2，篩選不符合序列；根據(jù)序列首尾兩端的barcode和引物區(qū)分樣本，并調(diào)整序列方向，barcode允許的錯(cuò)配數(shù)為0，最大引物錯(cuò)配數(shù)為2。使用的 UPARSE軟件（version 7.1 http://drive5.com/uparse/)，根據(jù)97%的相似度對(duì)序列進(jìn)行分類操作單元（OUT）聚類并剔除嵌合體。利用RDP classifier（http://rdp.cme.msu.edu/）對(duì)每條序列進(jìn)行物種分類注釋，比對(duì) Silva數(shù)據(jù)庫(kù)（SSU128），設(shè)置比對(duì)閾值為70%。

菌群多樣性分析及繪圖：均采用R語(yǔ)言進(jìn)行統(tǒng)計(jì)分析和作圖[11]。

2 結(jié)果與分析

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

通過(guò)對(duì)初乳和常乳組共計(jì)18 個(gè)樣本進(jìn)行16S rDNA測(cè)序、數(shù)據(jù)質(zhì)控后，共得到有效序列783372條，初乳組和常乳組分別得到優(yōu)化序列條數(shù)為392724 條和390648 條。如圖1 所示，兩組樣本的優(yōu)化序列長(zhǎng)度均主要分布在401～440 bp，集中在421～440 bp，與設(shè)計(jì)引物擴(kuò)增長(zhǎng)度300 bp左右接近，可滿足序列分析要求。

圖1 樣本有效序列長(zhǎng)度分布Fig.1 Distribution of effective column length of sample order

在97%相似性水平下對(duì)所有樣本序列進(jìn)行聚類分析[12]，結(jié)果見圖2，初乳組和常乳組共注釋出1176 個(gè)OTU，兩組共有的OTU為669，各自獨(dú)有的OTU分別為79 和428。初乳組共產(chǎn)生748 個(gè)OTU，常乳組共產(chǎn)生1097 個(gè)OTU，由此可知，常乳組OTU數(shù)遠(yuǎn)高于初乳組，兩組水牛乳樣本間存在大量共有的細(xì)菌菌群，且常乳組中細(xì)菌多樣性高于初乳組。

圖2 樣本OTU分布韋恩圖Fig.2 Venn diagram of sample OTUs distribution

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

圖3 Sobs（a）、Shannon（b）指數(shù)稀釋性曲線Fig.3 Rarefaction curves for Sobs (a) and Shannon (b) index

2.2 細(xì)菌組成與豐度差異分析

基于OTU注釋結(jié)果，分析樣本在不同分類水平上的群落組成。本次測(cè)序，Coverage（覆蓋率）為0.999，可代表99.9%的樣本信息。兩組樣本共注釋到1176 個(gè)OUT，26 個(gè)細(xì)菌門，97 個(gè)目，197 個(gè)科，495 個(gè)屬和767 個(gè)種。初乳組樣本共得到19 個(gè)門，52 個(gè)目，113 個(gè)科，292 個(gè)屬和437 個(gè)種。而常乳兩組樣本共得到7 個(gè)細(xì)菌門，45 個(gè)目，84 個(gè)科，203 個(gè)細(xì)菌屬和330 個(gè)菌種，相對(duì)而言，初乳組的細(xì)菌多樣性高于常乳組。為可視化呈現(xiàn)各樣本的物種豐度信息，本研究采用柱形圖來(lái)展示各樣本在門水平和屬水平的菌落分布及差異性[15]，結(jié)果如圖4、圖5 所示。由圖4 可知，在門分類水平上，以相對(duì)豐度≥0.01%為閾值，將<0.01%的菌門歸為others，初乳和常乳組共得到5 個(gè)細(xì)菌門，從高到低依次為變形菌門（Proteobacteria）、擬桿菌門（Bacteroidetes）、厚壁菌門（Firmicutes）、放線菌門（Actinobacteria ）、Patescibacteri。其中變形菌門（Proteobacteria）、厚壁菌門（Firmicutes）、放線菌門（Actinobacteria ）在常乳組中占比較多，而擬桿菌門（Bacteroidetes）、厚壁菌門（Firmicutes）則在初乳組中較占優(yōu)勢(shì)，同一菌門在不同樣本中的相對(duì)豐度差異較大。變形菌門、擬桿菌門、厚壁菌門為兩組樣本中的絕對(duì)優(yōu)勢(shì)菌門，相對(duì)豐度分別為35.4%、27.6%、24.8%，合計(jì)達(dá)88%。

圖4 門水平菌落豐度圖Fig.4 Abundance map of horizontal colony of phylum

在屬水平上，豐度≥0.1%的細(xì)菌屬有15 個(gè)（圖5），其中含量較高的有9 個(gè)屬：金黃桿菌屬（Chryseobacterium）、不動(dòng)桿菌屬（Acinetobacter）、假單胞菌屬（Pseusomonas）、巨型球菌屬（Macrococcus）、棲水菌屬（Enhydrobacter）、乳球菌屬（Lactococcus）、腸球菌屬（Enterococcus）、蒼白桿菌屬（Ochrobactrum）、鏈球菌屬（Streptococcus），各屬的相對(duì)豐度依次為24.7%、9.1%、6.6%、5.2%、5.2%、5.0%、4.7%、3.4%、2.9%，豐度和為61.6%，另外6 個(gè)屬的豐度和為25.2%，其它<0.1%的菌屬歸為others，豐度為13.2%。

圖5 屬水平菌落豐度圖Fig.5 Abundance map of horizontal colony of genus

由圖5 可知，不同樣本中主要菌屬占比差異較大，相對(duì)而言，常乳組樣本中菌群多樣性要高于初乳組。在種水平上，將初乳和常乳兩組樣本中各自的同名菌種豐度求和，其占比如表1 所示，兩組樣本間優(yōu)勢(shì)細(xì)菌豐度差異較大。初乳組中的炭疽金桿菌（Chryseobacterium anthropi）、溶酪大球菌（Macrococcus caseolyticus）、奧斯陸莫拉菌（Moraxella osloensis）占比較大，而這幾種菌在常乳組中的占比相對(duì)較少，尤其是炭疽金桿菌和溶酪大球菌。常乳組中占優(yōu)勢(shì)的菌為意大利腸球菌（Enterococcus italicus）、烏爾辛不動(dòng)桿菌（Acinetobacter ursingii）、炭疽金桿菌（Chryseobacterium anthropi）。而別雷斯不動(dòng)桿菌、假單胞菌、喬斯特金桿菌這三株菌在初乳組樣本中的相對(duì)豐度較高，但在常乳組樣本中卻未檢出或相對(duì)豐度小于<0.1%。

表1 豐度較高的細(xì)菌及其在門水平上的歸類Table 1 Bacteria with high abundance and their classification at phylum level

在屬分類水平，取樣本平均總豐度前30 的菌屬，按其物種和樣本層級(jí)聚類方式繪制群落豐度聚類熱圖（圖6），使高豐度和低豐度菌落分區(qū)聚集[16]，通過(guò)顏色梯度變化來(lái)反映不同分組樣本在屬水平群落組成的相似性和差異性。由圖6 可知，初乳組和常乳組整體并沒有分區(qū)聚集，但兩組樣本中優(yōu)勢(shì)菌群的豐度明顯不同，金黃桿菌（Chryseobacterium）、巨型球菌屬（Macrococcus）、棲水菌屬（Enhydrobacter）在初乳組樣本中的含量顯著高于常乳組，而假單胞菌屬（Pseudomonas）、蒼白桿菌屬（Ochrobactrum）、腸球菌屬（Enterococcus）在常乳組中的豐度均高于初乳組。金黃桿菌（Chryseobacterium）和不動(dòng)桿菌（Acinetobacter）在兩組樣本中均占比較高，葡萄球菌屬（Staphylococcus）在兩組樣本中均有低豐度出現(xiàn)。

圖6 樣本基于屬水平的菌落聚類熱圖Fig.6 Community clustering heatmaps of samples based on genus level

2.3 細(xì)菌多樣性比較分析

采用樣本層級(jí)聚類樹來(lái)分析初乳和常乳組間群落組成的相似性。18 個(gè)樣本基于OTU的樣本層級(jí)聚類樹如圖7 所示，不同顏色代表不同分組，樹枝長(zhǎng)度和寬度分別表示樣本間的距離和遠(yuǎn)近，相似度越高，距離越小[17]。通過(guò)樣本聚類樹可以清楚地看出樣本間分支距離遠(yuǎn)近，直觀反映出各樣本間菌群的相似度和差異性[18]。由圖7 可知，除去個(gè)別異常樣本，初乳組樣本（紅色）主要聚成一類，相似度較高，其中C1 和C2 是最相似的兩個(gè)樣本。常乳組（藍(lán)色）則分成兩類，其中M2、M3、M6 聚成一類，而M1、M5、M7、M8、M9 聚成另一類，其中M7、M9 相似度最高，這兩聚類樣本間相似度差異較大。整體而言，初乳組與常乳組組間樣本差異不明顯，初乳組組內(nèi)樣本間菌落組成相似度較高，而常乳組組內(nèi)樣本間菌落組成差異較大。

圖7 樣本層級(jí)聚類分析Fig.7 Sample level cluster analysis

對(duì)初乳（CR4）與常乳（PR）樣本中細(xì)菌屬（相對(duì)豐度閾值>0.1%）進(jìn)行主成分（PCA）分析，以便找出兩組樣本間菌群組成的相似性和差異性。樣本通過(guò)降維后，在主成分上會(huì)顯示相對(duì)坐標(biāo)點(diǎn)，樣本點(diǎn)越接近，表明兩樣本物種組成越相似[19]，同理，各點(diǎn)之間距離越大，表明它們之間的菌群差異越大[20]，由圖8 可知，初乳（CR4）與常乳（PR）兩組樣本代表的點(diǎn)都按各自分組相對(duì)聚集，但兩組樣本點(diǎn)之間的距離較近，在空間上沒有完全分開，這說(shuō)明初乳組與常乳組組間樣本個(gè)體之間菌群組成有一定的差異，導(dǎo)致少部分樣本菌群結(jié)構(gòu)相似，大部分樣本菌群存在差異性。主成分1（PC1）和主成分2（PC2）對(duì)樣本差異性的貢獻(xiàn)分別達(dá)到了41.21%和17.66%，共計(jì)58.87%，可以代表大部分變量信息。由圖8 可知，初乳組與常乳組在第一主成分PC1 所在的X軸方向，兩組樣本的微生物組成有明顯的分離趨勢(shì)，這說(shuō)明PC1 因素有較高的可能性影響了樣本的組成，表明結(jié)果數(shù)據(jù)組內(nèi)重復(fù)性好，組內(nèi)樣本差異性大，且PC1 貢獻(xiàn)值為41.21%，可以作為有效的主成分，來(lái)解釋兩組樣本的組間差異。

圖8 樣本基于屬水平的主成分分析Fig.8 The sample principal component analysis on the genus level

PLS-DA（偏最小二乘法判別分析Partial Least Squares Discriminant Analysis）也是一種常見的檢驗(yàn)樣本組間差異的方法，此方法是根據(jù)分組方案，弱化組內(nèi)樣本差異，突出組間樣本差異[21]，更便于發(fā)現(xiàn)組間樣本的差異性。由PLS-DA（圖9）顯示，除去M8、S3 兩個(gè)離群嚴(yán)重的異常樣本，其余樣本點(diǎn)均按初乳與常乳所代表的分組相對(duì)聚集，兩組樣本點(diǎn)在COMP1（X軸）和COMP2(Y軸)方位能夠很好地區(qū)分開來(lái)，說(shuō)明初乳與常乳兩組樣本在組間的微生物組成存在明顯差異。這種一目了然的分離效果在之前的PCA分析圖中是沒有得到體現(xiàn)的。初乳組（CR4）所代表的樣本點(diǎn)都比較集中，說(shuō)明初乳組組內(nèi)樣本之間菌群組成較相似，而常乳組（PR）所代表的樣本點(diǎn)都比較分散，這說(shuō)明常乳組組內(nèi)樣本間菌群組成有一定的差異性。結(jié)合PCA和PLS-DA可知，初乳與常乳組組間差異顯著，組內(nèi)差異不明顯。

圖9 樣本基于屬水平的PLS-DA分析Fig.9 Sample PLS-DA analysis based on the genus level

3 討論

相較傳統(tǒng)的涂片、劃線等細(xì)菌培養(yǎng)和API試劑條等其它檢測(cè)方法，16S rDNA高通量測(cè)序更全面和準(zhǔn)確地反映了樣本中菌群的多樣性和豐度信息，解決了以往不可培養(yǎng)及低豐度菌在傳統(tǒng)方法中無(wú)法分離鑒別的局限[22]。依據(jù)測(cè)序結(jié)果，本次試驗(yàn)Coverage（覆蓋率）為0.999，顯示樣本細(xì)菌99.9%的信息得到了反映，且樣本的Shannon稀釋性曲線都趨于平緩，說(shuō)明本次實(shí)驗(yàn)的測(cè)序深度合理，測(cè)序量足夠，可代表樣本中絕大數(shù)的菌群多樣性信息。

本次測(cè)序，在門水平上，初乳和常乳組主要分屬于變形菌門（Proteobacteria）、擬桿菌門（Bacteroidetes）和厚壁菌門（Firmicutes），在屬水平上，初乳和常乳組主要菌群種類不盡相同，初乳組中的主要菌群有金黃桿菌屬（Chryseobacterium）、不動(dòng)桿菌屬（Acinetobacter）和假單胞菌屬（Pseudomonas），而常乳組中的優(yōu)勢(shì)菌屬主要為金黃桿菌屬（Chryseobacterium）、乳球菌屬（Lactococcus ）和腸球菌屬（Enterococcus）。黃衛(wèi)強(qiáng)[23]在研究母乳微生物群落構(gòu)成中，鑒定出的優(yōu)勢(shì)菌門與本試驗(yàn)結(jié)果相同，而在屬水平上的差異較大。這可能是因?yàn)槿巳楹退Ｈ橐蚱贩N和乳源不同所致。Mazmanian等[24]的研究表明，Bacteroides（擬桿菌）在人類初乳中含量非常豐富，它可能在新生兒腸道早期定群中發(fā)揮主要作用。同樣，本次試驗(yàn)水牛初乳樣本中，擬桿菌的占比也較大。

研究[25?26]表明，金黃桿菌、不動(dòng)桿菌、假單胞菌屬均具有一定的抗感染能力，且假單胞菌屬和腸球菌屬與糖類和酶類代謝存在重要聯(lián)系。綜合本實(shí)驗(yàn)結(jié)果，水牛初乳和常乳組所擁有的這些核心菌屬及其所具備的抗感染、修復(fù)能力以及糖和酶代謝功能，不論是對(duì)奶水牛產(chǎn)后乳房炎的修復(fù)還是犢牛早期腸道菌群定殖、消化等，都起到重要作用。

不論是在初乳還是常乳組中，炭疽金桿菌Chryseobacterium anthropi均被檢測(cè)到占有較高比例（見表1），而炭疽桿菌屬于需氧芽孢桿菌屬[27]，能引起人畜共患的炭疽病，牛、羊、駱駝等草食動(dòng)物是其主工傳染源[28]。炭疽桿菌對(duì)社會(huì)公共衛(wèi)生和經(jīng)濟(jì)發(fā)展的危害，迄今仍占相當(dāng)大的比重[29]。通過(guò)此次實(shí)驗(yàn)，說(shuō)明水牛所種畜牧場(chǎng)大部分奶水牛可能或已經(jīng)感染了炭疽病，應(yīng)馬上采取針對(duì)性的防治措施應(yīng)對(duì)。

同樣是水牛乳樣本，在本實(shí)驗(yàn)檢出的部分低豐度菌群，在某些研究中卻是優(yōu)勢(shì)菌群[30]，有研究[31?32]顯示，水牛乳中Staphylococcus（葡萄球菌）、Streptococci（鏈球菌）在陽(yáng)性組中所占比例較大，其組合感染是引起奶水牛乳腺病變的主要原因。本試驗(yàn)在兩組樣本中均檢測(cè)到低豐度的Staphylococcus（葡萄球菌）和Streptococci（鏈球菌），盡管這些奶水牛在樣本采集時(shí)并沒表現(xiàn)出乳腺炎的病癥[33]。這將為本所種畜廠奶水牛乳腺炎針對(duì)性的預(yù)防治療提供依據(jù)。

4 結(jié)論

采用16S rDNA高通量測(cè)序技術(shù)，對(duì)廣西水牛研究所種畜廠內(nèi)的三品雜奶水牛初乳與常乳組樣本中的細(xì)菌組成及多樣性進(jìn)行了比較分析，結(jié)果顯示，兩組樣本共注釋到1176 個(gè)OUT，26 個(gè)細(xì)菌門，97 個(gè)目，197 個(gè)科，495 個(gè)屬和767 個(gè)種。其中初乳組樣本共得到292 個(gè)細(xì)菌屬和437 個(gè)細(xì)菌種。初乳組較常乳組細(xì)菌多樣性豐富。通過(guò)PCA和PLSDA分析顯示，初乳與常乳兩組樣本，組內(nèi)差異較小，組間差異明顯，各樣本優(yōu)勢(shì)菌群和豐度差異較大。通過(guò)本次測(cè)序，得知了水牛研究所種畜場(chǎng)三品雜水牛乳中的菌群分布狀況，對(duì)后續(xù)改善該種畜場(chǎng)奶源衛(wèi)生條件、防止原料乳污染、乳腺炎的前期預(yù)警等均有指導(dǎo)意義。