養(yǎng)豬發(fā)酵床微生物宏基因組基本分析方法

劉　波，王階平，陳倩倩，劉國紅，車建美，陳德局，鄭雪芳，葛慈斌

(福建省農(nóng)業(yè)科學(xué)院農(nóng)業(yè)生物資源研究所，福建　福州　350003)

?

養(yǎng)豬發(fā)酵床微生物宏基因組基本分析方法

劉波，王階平，陳倩倩，劉國紅，車建美，陳德局，鄭雪芳，葛慈斌

(福建省農(nóng)業(yè)科學(xué)院農(nóng)業(yè)生物資源研究所，福建福州350003)

養(yǎng)豬微生物發(fā)酵床分解豬糞、消除惡臭，將豬糞形成生物菌肥，實(shí)現(xiàn)了養(yǎng)豬污染的微生物治理和豬糞的資源化利用，微生物起到關(guān)鍵作用。然而，關(guān)于發(fā)酵床微生物群落的種類、數(shù)量、結(jié)構(gòu)等方面系統(tǒng)研究鮮見報道。課題組采用高通量宏基因組學(xué)方法，系統(tǒng)地開展了微生物發(fā)酵床的微生物組研究，揭示不同空間、不同深度、不同發(fā)酵程度、不同墊料組成、不同季節(jié)發(fā)酵床中的微生物組成及其區(qū)系演替規(guī)律。本文就相關(guān)的宏基因組測序的樣本采集、樣本處理、測序原理、分析模型、數(shù)據(jù)統(tǒng)計、結(jié)果表述等分析流程和方法進(jìn)行了歸納，闡明微生物組操作分類單元(OTU)鑒定，物種累積曲線，核心微生物組，種類豐度主成分分析，種類秩-多度曲線，物種組成豐度柱狀圖、熱圖、星圖等分析的原理與方法，列舉微生物組種類復(fù)雜度α-多樣性分析、β-多樣性分析、種類顯著性差異LEfSe分析，冗余分析等實(shí)例，為深入研究發(fā)酵床微生物群落提供了完整的分析思路和范例。

宏基因組學(xué)；微生物群落；微生物發(fā)酵床；分析方法

微生物發(fā)酵床養(yǎng)豬是利用植物廢棄物如谷殼、秸稈、鋸糠、椰糠等農(nóng)業(yè)副產(chǎn)品制作發(fā)酵床墊料層，添加微生物發(fā)酵劑，經(jīng)發(fā)酵后鋪墊到豬舍，生豬(母豬、公豬、小豬、育肥豬)生活在墊料上，排泄物混入墊料發(fā)酵，利用生豬的拱翻習(xí)性和機(jī)械翻耕，使豬糞尿和墊料充分混合，通過發(fā)酵床的微生物分解豬糞，消除惡臭，發(fā)酵床形成的有益微生物群落為限制豬病源的蔓延創(chuàng)造了有利條件，不僅實(shí)現(xiàn)了豬場的零排放，而且豬糞和墊料經(jīng)發(fā)酵轉(zhuǎn)化還可以生產(chǎn)優(yōu)質(zhì)微生物菌肥[1]。同時，飼喂飼用益生菌，替代抗生素，抑制豬病發(fā)生，解決豬肉藥殘問題。微生物發(fā)酵床養(yǎng)豬具有“五省、四提、三無、兩增、一少、零污染”優(yōu)點(diǎn)，即：五省，即省水、省工、省料、省藥、省電；四提，提高品質(zhì)、提高豬抗病力、提早出欄、提高肉料比；三無，無臭味、無蠅蛆、無環(huán)境污染；兩增，增加經(jīng)濟(jì)效益，增加生態(tài)效益；一少，減少豬肉藥物殘留；零污染，豬糞尿由微生物在豬舍內(nèi)降解消納，而實(shí)現(xiàn)零污染[2]。

微生物發(fā)酵床養(yǎng)豬技術(shù)綜合利用了微生物學(xué)、生態(tài)學(xué)、發(fā)酵工程學(xué)原理，以活性功能微生物作為物質(zhì)能量“轉(zhuǎn)換中樞”的一種生態(tài)養(yǎng)殖模式。該技術(shù)的核心在于利用活性強(qiáng)大的有益功能微生物復(fù)合菌群，長期和持續(xù)穩(wěn)定地將動物糞尿廢棄物進(jìn)行降解[3-4]。盡管以往的研究涉及發(fā)酵床的微生物群落，如劉波等研究了零排放豬場基質(zhì)墊層微生物群落脂肪酸生物標(biāo)記多樣性[5]，鄭雪芳等利用磷脂脂肪酸生物標(biāo)記分析豬舍基質(zhì)墊層微生物亞群落的分化[6]；李娟等進(jìn)行雞發(fā)酵床不同墊料理化性質(zhì)及微生物菌群變化規(guī)律研究[7]；張學(xué)鋒等研究了不同深度墊料對養(yǎng)豬土著微生物發(fā)酵床穩(wěn)定期微生物菌群的影響[8]；趙國華等進(jìn)行了生物發(fā)酵床養(yǎng)豬墊料中營養(yǎng)成分和微生物群落研究[9]；王震等研究了發(fā)酵床墊料中優(yōu)勢細(xì)菌的分離鑒定及生物學(xué)特性[10]；但至今，發(fā)酵床微生物種類、數(shù)量、結(jié)構(gòu)等方面的系統(tǒng)研究未見報道。其原因之一是通過傳統(tǒng)的可培養(yǎng)方法分離鑒定發(fā)酵床微生物種類的工作量非常大，要完成豐富的微生物種類的調(diào)查難度很大，而且，可培養(yǎng)方法不能全面調(diào)查微生物的種類。宏基因組測定，為發(fā)酵微生物群落的調(diào)查提供了方法。

宏基因組學(xué)，是研究直接來自環(huán)境的微生物基因材料的學(xué)科，被認(rèn)為是微生物學(xué)發(fā)展中的一個里程碑。通過測序分析微生物組主要包括16Sr RNA和宏基因組兩大技術(shù)，它不僅使得對未培養(yǎng)或者不可培養(yǎng)的微生物的研究成為可能，也使得研究同一環(huán)境中的微生物在自然條件下的相互作用以及微生物和環(huán)境因子的相互作用成為可能[11]。宏基因組學(xué)應(yīng)用范圍很廣，如分析傳統(tǒng)食醋發(fā)酵過程微生物多樣性[12]、牙菌斑微生物群落[13]、牛胃菌群組成[14]、人和動物胃腸道微生物群落[15]、深海樣本[16]、溫室黃瓜根結(jié)線蟲發(fā)生地土壤微生物[17]、鉛鋅尾礦酸性廢水微生物[18]、普洱茶渥堆發(fā)酵過程中微生物[19]、紅樹林土壤微生物[20]。

宏基因組數(shù)據(jù)分析主要包括序列處理、分類、注釋及統(tǒng)計分析4個環(huán)節(jié)。隨著測序技術(shù)的升級,測序成本將逐步降低，而大數(shù)據(jù)分析將成為核心內(nèi)容。數(shù)據(jù)具備標(biāo)準(zhǔn)化、可積累性特點(diǎn)，通過數(shù)據(jù)建模是未來應(yīng)用的基礎(chǔ)，培養(yǎng)和基于培養(yǎng)的功能驗(yàn)證將是未來的重點(diǎn)之一。微生物組學(xué)將闡述并調(diào)整環(huán)境與微生物組之間的關(guān)系，此領(lǐng)域相關(guān)研究有巨大的發(fā)展空間[21]

宏基因組分析方法是理解環(huán)境中微生物組種類、數(shù)量、結(jié)構(gòu)的工具。本研究以養(yǎng)豬發(fā)酵床墊料微生物宏基因組為材料，闡明宏基因組測序的樣本采集、樣本處理、數(shù)據(jù)分析等方法，以期為發(fā)酵床微生物群落的分析提供基礎(chǔ)和參考。

1　材料與方法

1.1發(fā)酵床墊料樣本采集及墊料理化性質(zhì)測定

試驗(yàn)地點(diǎn)位于福清漁溪現(xiàn)代設(shè)施農(nóng)業(yè)樣本工程示范基地微生物發(fā)酵床大欄養(yǎng)豬舍，大欄發(fā)酵床養(yǎng)殖面積為1 617 m2(33 m ×49 m)，深度80 cm，發(fā)酵床墊料由70%椰糠和30%谷殼組成。發(fā)酵床飼養(yǎng)1 617頭育肥豬，飼養(yǎng)密度為1頭·m-2。按下表將發(fā)酵床劃分為8 ×4個方格，冬季(11月)和春季(3月)取樣兩次，空間取樣在各個區(qū)域內(nèi)隨機(jī)挖取10 cm深度的墊料樣本，深度取樣分別取0、20、40、60 cm深度的發(fā)酵床樣本。冬季采樣空間平面樣本32個點(diǎn)和3個點(diǎn)4個深度樣本12個，春季采樣空間平面樣本32個點(diǎn)和2個點(diǎn)4個深度樣本8個，共84個樣本。分別測定各墊料樣本的溫度、pH值、含水量、含鹽量、硝態(tài)氮和粗纖維含量等理化性質(zhì)。

1.2發(fā)酵床墊料的發(fā)酵程度等級劃分

將墊料鋪平，將色差計的出光口正對墊料，使其不露光，測定墊料跟空白墊料的色差值(△E)。每份墊料重復(fù)3次，取平均值作為該份墊料的最終值。在大量分析數(shù)據(jù)的基礎(chǔ)上，快速而準(zhǔn)確地判定各樣本的發(fā)酵程度，并進(jìn)行發(fā)酵等級劃分。其中0<△E<30 為發(fā)酵等級一級；30<△E<50為發(fā)酵等級二級；50<△E<65為發(fā)酵等級三級；65<△E<80為發(fā)酵等級四級。

1.3宏基因組高通量測序

1.3.1總DNA的提取按土壤DNA提取試劑盒FastDNA SPIN Kit for Soil的操作指南，分別提取各墊料樣本的總DNA，于-80℃冰箱凍存?zhèn)溆谩?/p>

1.3.216S rDNA和ITS測序文庫的構(gòu)建采用擴(kuò)增原核生物16S rDNA的V3～V4區(qū)的通用引物U341F和U785R對各墊料樣本的總DNA進(jìn)行PCR擴(kuò)增，并連接上測序接頭，從而構(gòu)建成各墊料樣本的真細(xì)菌和古菌16S rDNA V3～V4區(qū)測序文庫。采用擴(kuò)增真菌5.8S和28S rDNA之間的轉(zhuǎn)錄間區(qū)的通用引物ITS-F1-12和ITS-R1-12對各墊料樣本的總DNA進(jìn)行PCR擴(kuò)增，并連接上測序接頭，從而構(gòu)建成各墊料樣本的真菌ITS測序文庫。

1.3.3高通量測序使用Illumina MiSeq測序平臺，采用PE300測序策略，每個樣本至少獲得10萬條reads。

1.4宏基因組測序數(shù)據(jù)質(zhì)控

通過PADNAseq軟件[22]利用重疊關(guān)系將雙末端測序得到的成對短序列(reads)拼接成1條序列，得到高變區(qū)的長reads。然后使用撰寫的程序?qū)ζ唇雍蟮膔eads進(jìn)行處理而獲取去雜短序列(clean reads)：去除平均質(zhì)量值低于20的reads；去除reads含N的堿基數(shù)超過3個的reads；reads長度范圍為250～500 nt。

1.5操作分類單元(OTU)聚類分析

為便于下游的物種多樣性分析，將標(biāo)簽(Tags)聚類成操作分類單元(Operational Taxonomic Units，OTU)。首先把Tags中的單例標(biāo)簽(singletons，即對應(yīng)的reads序列只有1條)過濾掉，因?yàn)檫@些singletons可能由于測序錯誤造成，故將這部分序列去除，不加入聚類分析。然后，利用usearch軟件，在0.97相似度下進(jìn)行聚類，對聚類后的序列進(jìn)行嵌合體過濾后，得到用于物種分類的OTU，每個OTU被認(rèn)為可代表1個物種[23]。

1.6抽平處理

為避免因各樣本數(shù)據(jù)大小不同而造成分析時的偏差，在樣本達(dá)到足夠測序深度的情況下，根據(jù)α-多樣性(Alpha diversity)指數(shù)分析結(jié)果對每個樣本進(jìn)行隨機(jī)抽平處理。抽平的參數(shù)必須在保證測序深度足夠的前提下去選取。

1.7核心微生物組(Core microbiome)分析

根據(jù)樣本的共有OTU以及OTU所代表的物種，可以找到核心微生物組(Core microbiome，即覆蓋90%樣本的微生物組)。OTU維恩圖(Venn)分析。Venn圖可以很好地反映組間共有以及組內(nèi)特有的OTU數(shù)目。利用R語言編寫的VennDiagram軟件里的venn.diagram函數(shù)實(shí)現(xiàn)。OTU的主坐標(biāo)分析(PCoA)：PCoA分析可以初步反映出不同處理或不同環(huán)境間的樣本可能表現(xiàn)出分散和聚集的分布情況，從而可以判斷相同條件的樣本組成是否具有相似性。利用R語言編寫的ade4軟件里的dudi.pca函數(shù)實(shí)現(xiàn)。

1.8物種分類和豐度分析

物種注釋分析：首先，分別從各個OTU中挑選出一條序列作為該OTU的代表序列，將該代表序列與已知物種的16S rRNA數(shù)據(jù)庫(網(wǎng)站Green Genes，http://greengenes.lbl.gov)[24]進(jìn)行比對，從而對每個OTU進(jìn)行物種鑒定；然后，根據(jù)每個OTU中序列的條數(shù)，得到各個OTU的豐度值。物種豐度分析：在門、綱、目、科、屬水平，將每個注釋上的物種或OTU在不同樣本中的序列數(shù)整理在一張表格，形成物種豐度的柱狀圖、星圖及統(tǒng)計表等。物種熱圖分析：物種熱圖利用顏色梯度可以很好反映出樣本在不同物種下的豐度大小以及物種聚類、樣本聚類信息，可利用R語言的gplots包的heatmap.2函數(shù)實(shí)現(xiàn)。

1.9樣本復(fù)雜度分析

1.9.1單個樣本復(fù)雜度分析α-多樣性也稱為生境內(nèi)的多樣性(within-habitat diversity)，是對單個樣本中物種多樣性的分析，常用的度量指標(biāo)有：測定物種指數(shù)(observed species)、chao1指數(shù)、香農(nóng)指數(shù)(Shannon)、辛普森指數(shù)(Simpson)以及譜系多樣性指數(shù)(phylogenetic diversity，PD_whole_tree)等。利用QIIME軟件計算樣本的α-多樣性指數(shù)值，并做出相應(yīng)的稀釋曲線[25]。稀釋曲線是利用已測得16S rDNA序列中已知的各種OTU的相對比例，來計算抽取n個(n小于測得reads序列總數(shù))reads時各α-多樣性指數(shù)的期望值，然后根據(jù)一組n值(一般為1組小于總序列數(shù)的等差數(shù)列)與其相對應(yīng)的α-多樣性指數(shù)的期望值做出曲線來。并作出α-多樣性指數(shù)的統(tǒng)計表格。

1.9.2樣本間復(fù)雜度比較分析β-多樣性也稱為生境間的多樣性(between-habitat diversity)，反映了不同樣本在物種多樣性方面存在的差異大小。分析各類群在樣本中的含量，進(jìn)而計算出不同樣本間的β-多樣性值。本分析中通過QIIME軟件，采用迭代算法，分別在加權(quán)物種分類豐度信息和非加權(quán)物種分類豐度信息的情況下進(jìn)行差異計算，得到最終的統(tǒng)計分析結(jié)果表并做出組間的距離箱線圖(Box Plot)及主成分分析(Principal Component Analysis，PCA)展示圖。

1.10顯著性差異分析

1.10.1判別效應(yīng)分析(LEfSe)LEfSe分析即LDA Effect Size分析，LEfSe采用線性判別分析(LDA，Linear discriminative analysis)來估算每個組分(物種)豐度對組分間差異效果影響的大小，可以實(shí)現(xiàn)多個分組之間的比較，還進(jìn)行分組比較的內(nèi)部進(jìn)行亞組比較分析，從而找到組間在豐度上有顯著差異的物種標(biāo)記(biomaker)，找出對組分劃分產(chǎn)生顯著性差異性影響的群落或物種[26]。本研究采用LEfSe Tools進(jìn)行分析[27]?；痉治鲞^程如下：首先，在多組樣本中采用Kruskal-Wallis非參數(shù)多組檢驗(yàn)法來檢測不同分組間豐度差異顯著的物種；然后，在上一步中獲得的顯著差異物種，用成組的威爾科克森(Wilcoxon)秩和檢驗(yàn)來進(jìn)行組間差異分析；最后，用線性判別分析(LDA)對數(shù)據(jù)進(jìn)行降維和評估差異顯著的物種的影響力(即LDA score)。LEfse分析可以進(jìn)行本地分析也可以在線分析，在線分析的網(wǎng)址是：https://huttenhower.sph.harvard. edu/galaxy/。

1.10.2組間差異分析使用秩和檢驗(yàn)的方法對不同分組之間的進(jìn)行顯著性差異分析，以找出對組間劃分產(chǎn)生顯著性差異影響的物種。本分析對于兩組間的差異分析采用R語言stats包的wilcox.test函數(shù)，對于兩組以上的組間差異分析采用R語言stats包的kruskal.test函數(shù)。

1.10.3物種典型分析(CCA)/冗余分析(RDA)即典型分析(canonical correspondence analusis，CCA)/冗余分析(redundancy analysis RDA)。RDA是基于線性模型，CCA是基于單峰模型。CCA/RDA主要用來檢測環(huán)境因子、樣本、菌群三者之間的關(guān)系或者兩兩之間的關(guān)系。可以基于樣本的所有OTU，也可基于樣本的優(yōu)勢物種或者差異物種。CCA采用R語言編寫的vegan軟件的cca函數(shù)，RDA采用vegan軟件的rda函數(shù)。

1.11發(fā)酵床微生物高通量測序與統(tǒng)計工作流程

養(yǎng)豬發(fā)酵床微生物組高通量測序(圖1)和生物信息學(xué)分析工作流程(圖2)分別歸納為：通過Illumina平臺(Hiseq或者M(jìn)iseq)進(jìn)行配對末端(Paired-end)測序，通過reads之間的覆蓋部分(overlap)關(guān)系拼接成長reads，并對拼接后的reads進(jìn)行質(zhì)控，除去平均質(zhì)量值低于20的reads；除去reads含N的堿基數(shù)超過3個的reads；reads長度范圍為250～500 nt，得到clean reads(圖1)。

測序處理后得到的clean reads(去雜短序列)，進(jìn)入生物信息學(xué)分析。通過拼接→去雜統(tǒng)計→OTU聚類→抽平處理→序列選擇→分類比對→豐度統(tǒng)計等，得到發(fā)酵床微生物種類的OTU，為微生物群落分析提供基礎(chǔ)(圖2)。

2　結(jié)果與分析

2.1微生物發(fā)酵床細(xì)菌OTU提取

2.1.1發(fā)酵床細(xì)菌OTU提取將序列完全一樣的去低值后reads(clean reads)歸為一種標(biāo)簽(tag)，統(tǒng)計每條tag對應(yīng)的豐度(即reads數(shù)目)，然后將tags根據(jù)其豐度大小進(jìn)行排序，將其中的單例標(biāo)簽(singletons，即對應(yīng)reads只有一條的序列)過濾掉，因?yàn)閟ingletons可能由于測序錯誤造成。利用usearch7.0軟件，在0.97相似度下進(jìn)行聚類，對聚類后的序列進(jìn)行嵌合體過濾后，最后將所有去低值reads比對到OTU序列上，將能比對上OTU的reads提取出來，得到最終的比對reads(mapped reads)(表1)。從表1可知，發(fā)酵床不同深度層次的墊料微生物tags、singleton、對比序列數(shù)(Mapped reads)、clean tags、OTUs存在顯著差異。用微生物發(fā)酵床夏季和冬季2個季節(jié)墊料進(jìn)行微生物宏基因組測序，2個季節(jié)共采樣84個樣本，分析結(jié)果表明，總體趨勢是發(fā)酵床表層和底層的相應(yīng)數(shù)值較小，中間層的相應(yīng)數(shù)值較大。統(tǒng)計各個樣本每個16S rRNA OTU中的豐度信息，一共產(chǎn)生154 556個16S rRNA-OTU(表1)，平均每個樣本有(1862.12±505.96)個OTU。分析結(jié)果表明養(yǎng)豬發(fā)酵床墊料中的微生物種類非常豐富。

表1　養(yǎng)豬發(fā)酵床樣本數(shù)據(jù)和16S rRNA部分OTU統(tǒng)計

注：第一列sample name是樣本名稱，S1.1.1代表樣點(diǎn)S1,1的表層，順序下去為20 cm、40 cm、60 cm層的樣本。

2.1.2發(fā)酵床細(xì)菌OTU對原始數(shù)據(jù)進(jìn)行質(zhì)量檢驗(yàn)(QC)后，用usearch7.0軟件對數(shù)據(jù)進(jìn)行去嵌合體和聚類的操作，usearch7.0軟件聚類時，先將reads按照豐度從大到小排序，通過97%相似度的標(biāo)準(zhǔn)聚類，得到16S rRNA OTU，每個OTU被認(rèn)為可代表一個細(xì)菌物種。在聚類過程中利用從頭測序(de novo)方法去除嵌合體(chimeras)。接下來對每個樣本的標(biāo)簽(tags)進(jìn)行隨機(jī)抽平處理以保證數(shù)據(jù)質(zhì)量，提取對應(yīng)的OTU序列。然后使用微生物生態(tài)學(xué)數(shù)量評估軟件(Qiime，Quantitative Insights Into Microbial Ecology，QIIME v 1.9.0)，做α-多樣性指數(shù)稀釋曲線。根據(jù)稀釋曲線選擇合理的抽平參數(shù)，利用Qiime軟件對得到的抽平后的OTU進(jìn)行分析，首先從OTU中分別提取一條reads作為代表序列，將該代表序列與核苷酸數(shù)據(jù)庫(RDP，Ribosomal Database Project)的16S rRNA比對，從而對每個OTU進(jìn)行細(xì)菌物種分類。歸類后，根據(jù)每個OTU中序列的條數(shù)，從而得到OTU物種含量豐度，最后根據(jù)該OTU豐度表進(jìn)行后續(xù)分析(表1)。

2.1.3發(fā)酵床細(xì)菌OTU抽平處理由于不同樣本對應(yīng)的reads數(shù)量差距較大，為了保證后期分析結(jié)果合理，對每個樣本的數(shù)據(jù)進(jìn)行隨機(jī)抽平處理，抽平參數(shù)根據(jù)α-多樣性指數(shù)的稀釋曲線來確定。α-多樣性反映的是單個樣本內(nèi)部的物種多樣性，包括測定物種指數(shù)(observed species)、豐富度Chao1指數(shù)，香農(nóng)多樣性指數(shù)(Shannon)以及辛普森優(yōu)勢度指數(shù)(Simpson)，譜系多樣性指數(shù)(PD_whole_tree)等。測定物種指數(shù)(observed species)和chao1指數(shù)反映樣本中細(xì)菌群落的豐富度(species richness)，即簡單指細(xì)菌群落物種總的數(shù)量，而不考慮細(xì)菌群落每個物種的豐度情況。這2個指數(shù)對應(yīng)的稀釋曲線還可以反映樣本測序量是否足以覆蓋所有類群。為了對各樣品進(jìn)行定量比較分析，根據(jù)發(fā)酵床細(xì)菌群落α-多樣性分析結(jié)果，結(jié)合測序飽和度和樣本信息完整性，對各樣本的測序數(shù)據(jù)進(jìn)行隨機(jī)抽平處理。根據(jù)分析結(jié)果，我們從發(fā)酵床部分樣本的測序結(jié)果中隨機(jī)抽取60 623條reads，進(jìn)行各樣本的各OTU豐度計算，分析部分結(jié)果見表2?？梢钥闯霭l(fā)酵床S(1，1)點(diǎn)表層墊料的細(xì)菌種類(OTU)為1 643，20 cm深度墊料細(xì)菌種類(OTUs)為2 055，等等，說明發(fā)酵床不同深度墊料的細(xì)菌種類(OTU)差異顯著。

表2　發(fā)酵床部分樣本微生物群落OTU抽平統(tǒng)計

2.1.4發(fā)酵床OTU稀釋曲線微生物多樣性分析中需要驗(yàn)證測序數(shù)據(jù)量是否足以反映樣本中的物種多樣性，稀釋曲線(豐富度曲線)可以用來檢驗(yàn)這一指標(biāo)，評價測序量是否足以覆蓋所有類群，并間接反映樣本中物種的豐富程度。稀釋曲線是利用已測得16S rDNA序列中已知的各種OTU的相對比例，來計算抽取n個(n小于測得reads序列總數(shù))reads時出現(xiàn)OTU數(shù)量的期望值，然后根據(jù)一組n值(一般為一組小于總序列數(shù)的等差數(shù)列，如0 ×104、2×104、4 ×104、6×104、8 ×104等)與其相對應(yīng)的OTU數(shù)量的期望值做出曲線來。當(dāng)曲線趨于平緩或者達(dá)到平臺期時，也就可以認(rèn)為測序深度已經(jīng)基本覆蓋到樣本中所有的物種；反之，則表示樣本中物種多樣性較高，還存在較多未被測序檢測到的物種。圖3展示了發(fā)酵床墊料樣本細(xì)菌chao指數(shù)和測定物種指數(shù)(observed species)的稀釋曲線圖。結(jié)果顯示，各樣本均具有較好的物種豐富度，且不同樣本的物種豐富度存在明顯差異；在測序達(dá)到一定深度后，各樣本的chao指數(shù)和測定物種指數(shù)(observed species)的稀釋曲線趨于平緩或者達(dá)到平臺期，表明測序深度已經(jīng)基本覆蓋到樣本中所有的物種。同時，與其他樣品相比，圖3中最下方的4條曲線表示這4個樣品中的物種多樣性明顯偏低。

2.1.5發(fā)酵床細(xì)菌OTU群落α-多樣性香農(nóng)多樣性指數(shù)(Shannon)以及辛普森優(yōu)勢度指數(shù)(Simpson)反映群落的物種多樣性(species diversity)，受樣本細(xì)菌群落中物種豐富度(species richness)和物種均勻度(species evenness)的影響。相同物種豐富度的情況下，群落中各物種均勻度越大，則認(rèn)為群落多樣性越高，同樣，Shannon指數(shù)和Simpson指數(shù)值越大，說明個體分分布越均勻。如果每一個體都屬于不同的種，Shannon指數(shù)和Simpson指數(shù)就大，如果每一個體都屬于同一種，則Shannon指數(shù)和Simpson指數(shù)就小。圖4展示了發(fā)酵床墊料部分樣本的Shannon指數(shù)和Simpson指數(shù)稀釋曲線圖。結(jié)果表明，絕大多數(shù)樣本的Shannon指數(shù)和Simpson指數(shù)稀釋曲線很快就達(dá)到了平臺期，說明各個樣本均具有較好的物種多樣性，而且不同個體屬于不同種的可能性大，即物種均勻度大。

2.1.6發(fā)酵床部分樣本微生物(OTU)物種譜系演化α-多樣性指數(shù)(PD_whole_tree and- Goods_coverage)稀釋曲線該指數(shù)反映了樣本中細(xì)菌物種譜系演化的差異。譜系多樣性指數(shù)(PD_whole_tree)指數(shù)越大說明物種譜系演化差異越大。測序深度指數(shù)(Goods_coverage)反映了測序的深度，測序深度指數(shù)(Goods_coverage)指數(shù)越接近于1，說明測序深度已經(jīng)基本覆蓋到樣本中所有的物種。圖5展示了樣本譜系多樣性指數(shù)(PD_whole_tree)和測序深度指數(shù)(Goods_coverage)的稀釋曲線。圖中一條曲線代表一個樣本，橫軸表示從某個樣本中隨機(jī)抽取的clean reads數(shù)目，縱軸表示該reads數(shù)目對應(yīng)的alpha多樣性指數(shù)的大小。在測序深度指數(shù)(Goods_coverage)的稀釋曲線中，隨著測序深度的增加，當(dāng)曲線趨于平緩時，表示此時測序深度已經(jīng)基本覆蓋到樣本中所有的物種，測序數(shù)據(jù)量比較合理。從圖5可以看出，每一個樣品均獲得了高質(zhì)量的測序深度，能覆蓋樣品中的所有物種信息。

2.2微生物發(fā)酵床核心微生物組(Core microbiome)分析

2.2.1發(fā)酵床微生物組共有OTU數(shù)與覆蓋樣本數(shù)關(guān)系基于16S rRNA基因信息的系統(tǒng)分類結(jié)果與基于全基因組信息的分類結(jié)果很相似，并且消除了克隆問題，可以綜合研究各個可變區(qū)，分析方法也相對成熟，廣泛應(yīng)用于核心微生物組的研究。以微生物發(fā)酵床墊料樣本為例，對采集的夏季和冬季的84個樣本分別進(jìn)行宏基因組分析，列出各樣本的細(xì)菌種類OTUs，對各樣本共有的OTUs進(jìn)行統(tǒng)計。分析結(jié)果(圖6)表明微生物發(fā)酵床墊料各樣本共有細(xì)菌種類(OTUs)數(shù)與樣本數(shù)的關(guān)系，圖中橫坐標(biāo)表示的是覆蓋一定比例以上(如>50%、>60、>70%、>80%、>90%等)樣本的共有細(xì)菌種類(OTUs)數(shù)目的比例，縱坐標(biāo)是統(tǒng)計的覆蓋大于此比例樣本的共有細(xì)菌種類(OTUs)數(shù)目。例如，若1個細(xì)菌種類(OTUs)覆蓋50%以上的樣本，則在>0.5的橫坐標(biāo)所對應(yīng)的縱坐標(biāo)的細(xì)菌種類(OTUs)數(shù)目為1200種(OTUs)，覆蓋60%以上的樣本細(xì)菌種類數(shù)目(OTUs)為880種，以此類推，覆蓋樣本越高的細(xì)菌種類(OTUs)其數(shù)量越少，覆蓋100%樣本的細(xì)菌種類(OTUs)則非常少(圖6)。

2.2.2微生物發(fā)酵床冬季和春季共有和特有微生物組(OTUs)分析在0.97的相似度下，得到了發(fā)酵床墊料每個樣本的OTUs個數(shù)，利用韋恩圖(Venn)展示多樣本共有和各自特有細(xì)菌物種(OTUs)數(shù)目，直觀展示樣本間物種重疊情況。結(jié)果表明，微生物發(fā)酵床冬季樣本(S1)(40個)和春季樣本(S2)(44個)鑒定到的細(xì)菌種類(OTUs)數(shù)目分別為5 472和6 858個，冬季和春季特有的OTUs數(shù)目分別為598和1 984個，春季的OTUs比冬季略豐富，而共有的OTU數(shù)目為4 874個，說明墊料的微生物群落具有一定的穩(wěn)定性(圖7)。

2.2.3發(fā)酵床微生物組(OTU)豐度主成分分析主成分分析(Principal Component Analysis，PCA)是一種分析和簡化數(shù)據(jù)集的技術(shù)。主成分分析經(jīng)常用于減少數(shù)據(jù)集的維數(shù)，保持?jǐn)?shù)據(jù)集中的對方差貢獻(xiàn)最大的特征。通過保留低階主成分，忽略高階主成分而做到的。保留的低階成分往往能夠最大程度的保留數(shù)據(jù)的重要特征。PCA運(yùn)用降維的思想，通過分析不同樣本OTUs(97%相似性)組成可以反映樣本的差異和距離，將多維數(shù)據(jù)的差異反映在二維坐標(biāo)圖上，坐標(biāo)軸取值采用對方差貢獻(xiàn)最大的前2個特征值。如果2個樣本距離越近，則表示這2個樣本的組成越相似。不同處理或不同環(huán)境間的樣本，可表現(xiàn)出分散和聚集的分布情況，從而可以判斷相同條件的樣本組成是否具有相似性。以微生物發(fā)酵床2個季節(jié)的部分樣本進(jìn)行分析說明，結(jié)果表明，微生物發(fā)酵床2個季節(jié)的部分樣本的微生物組成具有一定的相似性(圖8)。分析結(jié)果表明，2個季節(jié)共有的種類主要集中在第1相限，種類較多，表明微生物群落具有穩(wěn)定性；春季S2獨(dú)有的種類主要分布在第3相限，分布較為分散，表明中間差異較大；冬季S1獨(dú)有的種類主要分布在第4相限，分布較為集中，表明中間差異較小。

2.2.4發(fā)酵床微生物組(OTU)秩-多度曲線秩-多度曲線(Rank_Abundance Curve)能同時解釋樣本所含物種的豐富度和均勻度。物種的豐富程度由曲線橫軸上長度來反映，物種的相對豐度由縱坐標(biāo)長度來反映；物種組成的均勻程度由曲線整體斜率來反映，曲線斜率越大則表示樣本中各物種所占比例差異較大，分布越均勻。如圖9所示，微生物發(fā)酵床冬季和春季84個樣本的曲線在橫軸上均在103左右，說明各樣本具有較高的豐富度；各曲線的斜率相對較大，表明各樣本的均勻度較好，且樣本間的均勻度存在一定差異。

2.2.5發(fā)酵床微生物組(OTU)物種累積曲線物種累積曲線(species accumulation curves)用于描述隨抽樣量增加物種數(shù)量增加的狀況，是理解和預(yù)測物種豐富度(species richness)的有效工具，利用物種累積曲線不僅可以判斷抽樣量充分性，在抽樣量充分的前提下，運(yùn)用物種累積曲線可以對物種豐富度進(jìn)行預(yù)測。對微生物發(fā)酵床冬季和春季84個墊料樣本細(xì)菌物種(OTUs)測序結(jié)果進(jìn)行抽樣，分析結(jié)果如圖10所示。隨著抽樣reads的增加，墊料樣本的物種累積曲線迅速趨于平緩，最后達(dá)到平臺期，說明樣本測序深度符合要求，同時說明各樣品具有較好的物種豐富度。

2.2.6發(fā)酵床微生物組(OTU)分類階元數(shù)量分布統(tǒng)計結(jié)果(表3)表明從各個OTU中挑選出豐度最高的1條序列，作為該OTU的代表序列。使用核苷酸數(shù)據(jù)庫(RDP)比對方法，將該代表序列與已知物種的16S rRNA數(shù)據(jù)庫進(jìn)行比對，從而對每個OTU進(jìn)行物種歸類。微生物發(fā)酵床樣本共鑒定出細(xì)菌OTUs總數(shù)為7 456，每個樣本可鑒定的OTUs范圍為1 643～2 467，其中可鑒定到科水平OTUs有4 886，可鑒定到屬水平OTUs有2 351，可鑒定到種水平OTUs有246。

表3　微生物發(fā)酵床選擇樣本細(xì)菌OTU統(tǒng)計

2.2.7發(fā)酵床微生物組(OTU)豐度柱狀圖分析微生物發(fā)酵床所有84個樣本門分類階元物種(OTUs)含量前20位的物種的組成與豐度(圖11)結(jié)果表明，豐度在前5位的分別是放線菌門(Actinobacteria)>厚壁菌門(Firmicutes)>變形菌門(Proteobacteria)>擬桿菌門(Bacteroidetes)>綠彎菌門(Chloroflexi)，與一般的土壤樣本(變形菌門含量最高)存在明顯差異，說明微生物發(fā)酵床具有獨(dú)特的生境類型。

2.2.8發(fā)酵床冬季和春季微生物組(OTU)豐度分析發(fā)酵床冬季和春季微生物組(OTU)豐度如圖12所示，將發(fā)酵床冬季和春季分別作為樣本，兩個樣本細(xì)菌物種門分類階元(OTU)的相對豐度在前5位的也是放線菌門Actinobacteria>厚壁菌門Firmicutes>變形菌門Proteobacteria>擬桿菌門Bacteroidetes>綠彎菌門Chloroflexi。

2.2.9發(fā)酵床微生物組(OTU)物種熱圖發(fā)酵床微生物組(OTU)物種熱圖分析結(jié)果如圖13所示，物種熱圖是以顏色梯度來代表數(shù)據(jù)矩陣中數(shù)值的大小，并能根據(jù)物種或樣本豐度相似性進(jìn)行聚類的一種圖形展示方式。聚類結(jié)果加上樣本的處理或取樣環(huán)境分組信息，可以直觀觀察到相同處理或相似環(huán)境樣本的聚類情況，并直接反映了樣本的群落組成的相似性和差異性。微生物發(fā)酵床冬季和春季樣本中的細(xì)菌種類分別在門，綱，目，科，屬，種分類等級進(jìn)行heatmap聚類分析?？v向聚類表示所有物種在不同樣本間的相似情況，距離越近，枝長越短，說明樣本的物種組成及豐度越相似。橫向聚類表示該物種在各樣本豐度相似情況，與縱向聚類一樣，距離越近，枝長越短，說明兩物種在各樣本間的組成越相似。

2.2.10組(OTU)物種星圖(star)分析圖14展示了發(fā)酵床冬季和春季84份樣本屬水平豐度在前10位的細(xì)菌物種(OTU)整理的星圖(star)，每個星形圖中的扇形代表一個物種，用不同顏色區(qū)分，用扇形的半徑來代表物種相對豐度的大小，扇形半徑越長代表此扇形所對應(yīng)的物種的相對豐度越高。結(jié)果表明不同樣本在屬水平的物種組成存在明顯差異性，而且，2個季節(jié)樣本中的棒桿菌屬Corynebacterium、葡萄球菌屬Staphylococcus、梭菌屬Clostridium是多數(shù)樣本的優(yōu)勢屬，這3個屬主要來源于人和動物，凸顯了養(yǎng)殖微生物發(fā)酵床中的微生物組成與區(qū)別于一般土壤樣本微生物組成的差異。

2.3發(fā)酵床微生物組(OTU)α-多樣性種類復(fù)雜度分析

2.3.1發(fā)酵床單個樣本微生物組(OTU)種類復(fù)雜度分析為了分析單個微生物發(fā)酵床樣本的細(xì)菌種類復(fù)雜度，分別進(jìn)行了衡量各個樣品細(xì)菌種類組成和豐度復(fù)雜度的α-多樣性指數(shù)常用度量指標(biāo)[包括測定物種指數(shù)(observed species)、chao1指數(shù)、香農(nóng)指數(shù)(Shannon)、辛普森指數(shù)(Simpson)、測序深度指數(shù)(Goods_coverage)、譜系多樣性指數(shù)(PD_whole_tree]的計算。表4展示了部分冬季和春季部分樣本的α-多樣性指數(shù)的統(tǒng)計結(jié)果，可以看出，測序深度指數(shù)(Goods_coverage)和Simpson指數(shù)均接近于100%或1，表明有較好的α-多樣性。

2.3.2發(fā)酵床冬季和春季微生物組(OTU)箱線圖分析圖15展示了發(fā)酵床冬季(S1)和春季(S2)樣本α-多樣性指數(shù)箱線圖(Box-plot)，能更直觀顯示各樣品的α-多樣性指數(shù)差異。其中，中間的黑色橫線是數(shù)據(jù)的中位數(shù)(median)，即數(shù)據(jù)中占據(jù)中間位子的數(shù)，即數(shù)據(jù)中有一半大于中位數(shù)(在其之上)，另一半小于中位數(shù)(在其之下)。盒子的上下兩邊稱為上下四分位數(shù)(Hinges)，其意義為：數(shù)據(jù)中有四分之一的數(shù)目大于上四分位數(shù)(即盒子的上邊)，即在紅色盒子之上；另外有四分之一的數(shù)目小于下四分位數(shù)(即盒子的下邊)，也就是在盒子之下。也就是說有一半的數(shù)目在中間封閉盒子的范圍內(nèi)，有一半分布在盒子上下兩邊。在盒子上下兩邊分別有一條縱向的線段，叫觸須線。上截止橫線是變量值本體最大值，下截止橫線是變量值本體最小值。本體指的是除異常值和極值以外的變量值，稱為本體值。異常值標(biāo)記為o，極值標(biāo)記為*。高于觸須線上截止橫線的值的取值范圍為：①異常值：x>上四分位數(shù)+1.5IQR；②極值：x>上四分位數(shù)+3.0IQR；低于觸須線下截止橫線的值的取值范圍為：①異常值：x<下四分位數(shù)-1.5IQR；②極值：x<下四分位數(shù)-3.0IQR；從而表明盒子外面數(shù)值點(diǎn)的分布。IQR(interquartile range)=上四分位數(shù)-下四分位數(shù)。

從圖15可以看出：①春季的度量比冬季的分散得多；②冬季和春季的Simpson指數(shù)幾乎沒有差異，而且數(shù)值均較高，春季的測序深度指數(shù)(Goods_coverage)指數(shù)略優(yōu)于冬季，其他4個指數(shù)均為冬季優(yōu)于春季(主要根據(jù)中位數(shù)和盒子的位置來判斷)。

2.3.3發(fā)酵床微生物組(OTU)α-多樣性指數(shù)秩和檢驗(yàn)對細(xì)菌物種(OTU)α-多樣性指數(shù)進(jìn)行秩和檢驗(yàn)(rank sum test)，分析不同條件下顯著差異的α-多樣性指數(shù)。分析結(jié)果見表5，表明冬季和春季的Shannon和Simpson指數(shù)均無顯著差異，其P-value值分別為0.1269和0.1036(遠(yuǎn)大于0.05或0.01)，說明在多樣性上冬季和春季微生物組無顯著差異；冬季和春季微生物組的測定物種(observed_species)、測序深度指數(shù)(Goods_coverage)、chao1指數(shù)差異顯著(P<0.05)，譜系多樣性指數(shù)(PD_whole_tree)差異極顯著(P<0.01)。

表5　微生物發(fā)酵床細(xì)菌物種(OTU)α-多樣性指數(shù)秩和檢驗(yàn)

2.4發(fā)酵床樣本間微生物組(OTU)β-多樣性種類復(fù)雜度分析

2.4.1發(fā)酵床樣本間微生物組(OTU)UniFrac距離通過利用系統(tǒng)進(jìn)化的信息來比較樣本間的細(xì)菌群落差異，其計算結(jié)果可以作為一種衡量β-多樣性指數(shù)，它考慮了物種間的進(jìn)化距離，該指數(shù)越大表示樣本間的差異越大。微生物發(fā)酵床部分樣品β-多樣性UniFra距離數(shù)據(jù)矩陣如表6(weighted_unifrac)和表7(unweighted_unifrac)所示。

2.4.2發(fā)酵床樣本間微生物組(OTU)UniFrac距離聚類分析分析結(jié)果見圖16。展示了發(fā)酵床樣品的UniFrac距離分布熱圖(heatmap)和聚類圖，通過對UniFrac距離的聚類，具有相似的β-多樣性的樣本聚類在一起，反映了樣本間的相似性。

表6　加權(quán)物種豐度信息計算得到的樣本β-多樣性數(shù)據(jù)矩陣(weighted_unifrac)

表7　非加權(quán)物種豐度信息計算得到的樣本β-多樣性數(shù)據(jù)矩陣(unweighted_unifrac)

2.4.3發(fā)酵床樣本間微生物組(OTU)主成分分析(PCA)為了進(jìn)一步展示樣本間物種多樣性差異，使用PCA分析的方法展示各個樣本間的差異大小。圖17給出了PCA對樣本間物種多樣性的分析結(jié)果，如果兩個樣本距離較近，則表示這2個樣本的物種組成較相近。從結(jié)果來看，與加權(quán)的相比，非加權(quán)的冬季和春季樣本間物種多樣性的差異更加顯著。

2.5發(fā)酵床微生物組(OTU)差異性分析

2.5.1發(fā)酵床微生物組(OTU)差異效應(yīng)判別分析(LEfSe)分析結(jié)果見圖17。LEfSe(LDA Effect Size)分析結(jié)果可以用LDA值(Linear Discriminant Analysis)作分布柱狀圖和進(jìn)化分支圖(圖18)。在LDA值分布柱狀圖中，展現(xiàn)了發(fā)酵床不同組中(即紅色代表春季；綠色代表冬季)微生物組(OTU)豐度有顯著差異的物種，即具有統(tǒng)計學(xué)差異的微生物標(biāo)記(Biomaker)，柱狀圖的長度代表差異物種的影響值大小。在進(jìn)化分支圖中的著色原則：無顯著差異的物種統(tǒng)一著色為黃色，差異物種(Biomarker)跟隨組別進(jìn)行著色，紅色節(jié)點(diǎn)表示在春季組別起到重要作用的微生物類群，綠色節(jié)點(diǎn)表示在冬季組別中起到重要作用的微生物類群，其他圈顏色意義類同。圖中英文字母表示的物種名稱在右側(cè)圖例中進(jìn)行展示。由內(nèi)至外輻射的圓圈代表了由門至屬(或種)的分類級別。在不同分類級別上的每一個小圓圈代表該水平下的一個分類，小圓圈直徑大小與相對豐度大小呈正比。從結(jié)果可以看出，不同季節(jié)各自對應(yīng)不同的差異物種(Biomarker)，而且，春季的差異物種(Biomarker)明顯多于冬季。

2.5.2發(fā)酵床不同分類階元微生物組(OTU)差異分析分析結(jié)果見表8。分別從門、綱、目、科、屬階元，通過秩和檢驗(yàn)對冬季和春季樣本篩選顯著差異(P<0.05，Kruskal-Wallis test)的物種，表8是在冬季和春季組間20個差異顯著的不同分類水平OTUs(P<0.05)，共找到317個，其中g(shù)代表屬(Genus)，f代表科(Family)，o代表目(Order)，c代表綱(Class)。

表8　發(fā)酵床冬季和春季不同分類階元微生物組(OTU)差異分析

2.5.3發(fā)酵床屬分類階元微生物組(OTU)差異分析通過Kruskal-Wallis test分析可以找出在冬季和春季組間有明顯差異(P<0.05)的屬共有136個，表9展示了部分結(jié)果，可看到各屬的基本信息。

表9　發(fā)酵床冬季和春季屬分類階元微生物組差異分析

注：第一列是物種名稱；第二列至第四列分別為兩組樣本的均值；最后一列是秩和檢驗(yàn)的P值。

2.5.4發(fā)酵床冬季和春季微生物組(OTU)熱圖和主成分分析為了直觀地展示這些具有顯著差異的屬，分別對它們進(jìn)行PCA分析(圖19)和熱圖(圖20)。結(jié)果可以看出，冬季和春季的差異屬可以分別聚類在一起(圖19)。

2.6發(fā)酵床微生物組(OTU)種類冗余(RDA)分析

梯度分析理論，如基于單峰模型的典型分析(CCA，Canonical correspondence analysis)和基于線性模型的冗余分析(RDA，Redundancy analysis)，可用于多指標(biāo)之間的對比，進(jìn)行描述環(huán)境變遷的可信度(顯著性)及解釋的能力(重要性)等問題的統(tǒng)計學(xué)定量刻畫。其中RDA是一種直接梯度分析方法，能從統(tǒng)計學(xué)的角度來評價一個或一組變量與另一組多變量數(shù)據(jù)之間的關(guān)系。圖21展示了微生物發(fā)酵床墊料中的微生物組與季節(jié)因子(即氣溫)RDA分析(即關(guān)聯(lián)性分析)結(jié)果，為了保證圖中字符不重疊，物種名取了前8個字符來代表，圖中標(biāo)出了豐度為前20位的物種。藍(lán)色線條與紅色線條成銳角，表示這些物種與季節(jié)因子(氣溫)呈正相關(guān)；與紅色線條成鈍角，表示這些物種與季節(jié)因子(氣溫)呈負(fù)相關(guān)。微生物組(OTU)落在第二象限的種類與季節(jié)呈正相關(guān)，如短芽胞桿菌屬Brevibacillus、短桿菌屬Brachybacterium、周氏桿菌屬Zhouia、糖多孢菌屬Saccharopolyspora、海桿菌屬M(fèi)arinobacter，碳酸噬胞菌屬Aequorivita等，它們與季節(jié)因子(氣溫)均具有顯著的正相關(guān)性，其他種類為負(fù)相關(guān)性。

3　討　論

我國畜禽糞便是農(nóng)業(yè)面源污染的主要來源，每年僅豬糞的產(chǎn)生量就突破50億t，由于處理不當(dāng)對環(huán)境造成嚴(yán)重污染，已經(jīng)成為經(jīng)濟(jì)發(fā)達(dá)地區(qū)或水環(huán)境敏感地區(qū)優(yōu)先控制的污染源。微生物發(fā)酵床能徹底地處理畜禽糞便污染，其作用機(jī)理在于發(fā)酵床墊料中的微生物[1]。盡管以往的研究有涉及發(fā)酵床微生物群落，如發(fā)酵床微生物群落脂肪酸標(biāo)記多樣性[5]、發(fā)酵床微生物脂肪酸亞群落分化[6]、雞發(fā)酵床微生物菌群變化[7]、穩(wěn)定期發(fā)酵床微生物菌群特性[8]、發(fā)酵床微生物群落研究[9]、發(fā)酵床優(yōu)勢細(xì)菌鑒定[10]等，但至今發(fā)酵床與微生物的密切相關(guān)性知之甚少。有鑒于此，本研究利用高通量的宏基因組學(xué)分析方法，全面而系統(tǒng)地開展了微生物發(fā)酵床的微生物組研究，揭示不同空間、不同深度、不同發(fā)酵程度、不同墊料組成、不同季節(jié)、不同管理的發(fā)酵床中的微生物組成、區(qū)系演替規(guī)律，以期探明豬糞降解機(jī)理、資源化利用機(jī)制以及污染治理原理。宏基因組分析方法的出現(xiàn)，為發(fā)酵床微生物組分析提供了強(qiáng)有力的工具，為使這一技術(shù)盡快地應(yīng)用到發(fā)酵床微生物組的分析上，作者從樣本采集、樣本處理、測序原理、分析模型、數(shù)據(jù)統(tǒng)計、結(jié)果表述等方面系統(tǒng)地進(jìn)行了介紹，以推進(jìn)發(fā)酵床微生物組的研究。

發(fā)酵床微生物組(OTU)α-多樣性種類復(fù)雜度分析，也稱為生境內(nèi)的多樣性(within-habitat diversity)，是對單個樣本中物種多樣性的分析，已經(jīng)有很多度量指標(biāo)(包括測定物種指數(shù)observed species、chao1指數(shù)、香農(nóng)指數(shù)Shannon、辛普森指數(shù)Simpson以及譜系多樣性指數(shù)phylogenetic diversity，PD_whole_tree等)來進(jìn)行評價[28]。β-多樣性也稱為生境間的多樣性(between-habitat diversity)，反映了不同樣本在物種多樣性方面存在的差異大小，可衡量群落之間的差別，在許多研究領(lǐng)域，尤其是生態(tài)學(xué)研究中，具有重要的意義[29]。然而，目前β-多樣性的評價方法不多，在這些方法中，加權(quán)UniFrac(weighted UniFrac)和非加權(quán)UniFrac(unweighted UniFrac)的應(yīng)用非常廣泛[30]。其中加權(quán)UniFrac考慮了物種的豐度，非加權(quán)UniFrac則不考慮物種豐度。與P-檢驗(yàn)類似，UniFrac分析的先決條件是一個包含所有物種的有根、各枝長已知的系統(tǒng)發(fā)生樹。2個群落之間的UniFrac距離定義為：對于系統(tǒng)發(fā)育樹所有分枝，考查其指向的葉節(jié)點(diǎn)是否只存于同一個群落，那些葉節(jié)點(diǎn)只存在于同一群落的枝的枝長和占整個樹的枝長總和的比例。直觀來講，就是計算了僅被一個群落占據(jù)的進(jìn)化歷史的相對大小，這個值越大，說明兩個群落中獨(dú)立的進(jìn)化過程越多，也就說明這兩個群落的差別越大。若兩個群落完全相同，那么它們沒有各自獨(dú)立的進(jìn)化過程，UniFrac距離為0；若2個群落在進(jìn)化樹中完全分開，即它們是完全獨(dú)立的兩個進(jìn)化過程，那么UniFrac距離為1[31]。

本研究就相關(guān)的宏基因組測序的樣本采集、樣本處理、測序原理、分析模型、數(shù)據(jù)統(tǒng)計、結(jié)果表述等分析流程和方法進(jìn)行了歸納，闡明微生物組操作分類單元(OTU)鑒定，物種累積曲線，核心微生物組，種類豐度主成分分析，種類秩-多度曲線，物種組成豐度柱狀圖、熱圖、星圖等分析的原理與方法，列舉微生物組種類復(fù)雜度α-多樣性分析、β-多樣性分析、種類顯著性差異LEfSe分析，冗余分析等實(shí)例，為深入研究發(fā)酵床微生物群落提供完整的分析思路和范例。

[1]劉波，朱昌雄. 微生物發(fā)酵床零污染養(yǎng)豬技術(shù)的研究與應(yīng)用[J]. 北京：中國農(nóng)業(yè)科學(xué)技術(shù)出版社，2009.

[2] 劉波，藍(lán)江林，唐建陽，等. 微生物發(fā)酵床菜豬大欄養(yǎng)殖豬舍結(jié)構(gòu)設(shè)計[J]. 農(nóng)業(yè)科學(xué)與技術(shù):英文版，2014，9(5): 1521-1525.

[3] 王遠(yuǎn)孝，錢輝，王恬. 微生物發(fā)酵床養(yǎng)豬技術(shù)的研究與應(yīng)用[J].中國畜牧獸醫(yī)，2011，38(5): 206-209.

[4] 藍(lán)江林，劉波，宋澤瓊，等. 微生物發(fā)酵床養(yǎng)豬技術(shù)研究進(jìn)展[J]. 生物技術(shù)進(jìn)展，2012，(6): 411-416.

[5] 劉波，鄭雪芳，林營志，等. 零排放豬場基質(zhì)墊層微生物群落脂肪酸生物標(biāo)記多樣性分析[J]. 生態(tài)學(xué)報，2008，28(12)：5488-5498.

[6] 鄭雪芳，劉波，林營志，等. 利用磷脂脂肪酸生物標(biāo)記分析豬舍基質(zhì)墊層微生物亞群落的分化[J]. 環(huán)境科學(xué)學(xué)報，2009，29(11):2306-2317.

[7] 李娟，石緒根，李吉進(jìn)，等. 雞發(fā)酵床不同墊料理化性質(zhì)及微生物菌群變化規(guī)律的研究[J]. 中國畜牧獸醫(yī)，2014，41(2): 139-143.

[8]張學(xué)峰，周賢文，陳群，等. 不同深度墊料對養(yǎng)豬土著微生物發(fā)酵床穩(wěn)定期微生物菌群的影響[J]. 中國獸醫(yī)學(xué)報，2013，33(9): 1458-1462.

[9] 趙國華，方雅恒，陳貴. 生物發(fā)酵床養(yǎng)豬墊料中營養(yǎng)成分和微生物群落研究[J]. 安徽農(nóng)業(yè)科學(xué)，2015，48(8)：98-101.

[10]王震，尹紅梅，劉標(biāo)，等. 發(fā)酵床墊料中優(yōu)勢細(xì)菌的分離鑒定及生物學(xué)特性研究[J]. 浙江農(nóng)業(yè)學(xué)報，2015，27(1)：87-91.

[11]常秦. 宏基因組數(shù)據(jù)分析中的統(tǒng)計方法研究[D]. 濟(jì)南：山東大學(xué)，2012.

[12]聶志強(qiáng)，韓玥，鄭宇，等. 宏基因組學(xué)技術(shù)分析傳統(tǒng)食醋發(fā)酵過程微生物多樣性[J]. 食品科學(xué)，2013，34(15):198-203.

[13]陳林. 老年人根面齲患者和健康人牙菌斑微生物群落的宏基因組學(xué)研究[D]. 武漢：武漢大學(xué)，2014.

[14]彭帥. 應(yīng)用宏基因組方法檢測豬致病微生物及分析牛胃菌群組成[D]. 長春：吉林大學(xué)，2015.

[15]許波，楊云娟，李俊俊，等. 宏基因組學(xué)在人和動物胃腸道微生物研究中的應(yīng)用進(jìn)展[J]. 生物工程學(xué)報，2013，29(12):1721-1735.

[16]江夏薇. 基于嗜耐鹽菌基因組分析與深海宏基因組文庫的酯酶研究[D]. 杭州：浙江大學(xué)，2013.

[17]趙志祥，蘆小飛，陳國華，等. 溫室黃瓜根結(jié)線蟲發(fā)生地土壤微生物宏基因組文庫的構(gòu)建及其一個殺線蟲蛋白酶基因的篩選[J]. 微生物學(xué)報，2010，50(8):1072-1079.

[18]韓玉姣. 凡口鉛鋅尾礦酸性廢水微生物宏基因組及宏轉(zhuǎn)錄組研究[D]. 廣州：中山大學(xué)，2013.

[19]呂昌勇. 普洱茶渥堆發(fā)酵過程中微生物宏基因組學(xué)的測定與分析[D]. 昆明：昆明理工大學(xué)，2013.

[20]蔣云霞. 基于紅樹林土壤微生物資源研發(fā)的宏基因組學(xué)平臺技術(shù)的建立與應(yīng)用初探[D]. 廈門：廈門大學(xué)，2007.

[21]盛華芳，周宏偉. 微生物組學(xué)大數(shù)據(jù)分析方法、挑戰(zhàn)與機(jī)遇[J]. 南方醫(yī)科大學(xué)學(xué)報，2015，35(7)：931-934.

[22]MASELLA A P， BARTRAM A K，TRUSZKOWSKI J M，et al. PANDAseq: paired-end assembler for illumina sequences[J]. BMC Bioinformatics，2012，13:31.

[23]EDGAR R C. UPARSE: Highly accurate OTU sequences from microbial amplicon reads[J]. Nat Methods， 2013，10(10):996-998.

[24]MCDONALD D，PRICE M N，GOODRICH J，et al. An improved Green genes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaea[J]. ISMEJ，2012，(6): 610-618.

[25]PAUL F K，JOSEPHINE Y A. Bacterial diversity in aquatic and other environments: what 16S rDNA libraries can tell us[J]. FEMS Microbiol Ecol，2004，47:161-177.

[26]ZHANG C，LI S，YANG L，et al. Structural modulation of gut microbiota in life-long calorie-restricted mice[J]. Nat Commun，2013，(4):2163.

[27]SEGATA N，IZARD J，WALDRON L，et al. Metagenomic biomarker discovery and explanation[J]. Genome Biol，2011，12(6): R60.

[28]BASSIOUNI A，CLELAND E J，PSALTIS A J，et al. Sinonasal microbiome sampling: a comparison of techniques[J]. PLoS One，2015，10(4):e0123216.

[29]LAROCHE F，JARNE P，PERROT T，et al. The evolution of the competition-dispersal trade-off affects α- and β-diversity in a heterogeneous metacommunity[J]. Proc Biol Sci，2016，283(1829). pii: 20160548.

[30]CHANG Q，LUAN Y, SUN F. Variance adjusted weighted UniFrac: a powerful beta diversity measure for comparing communities based on phylogeny[J]. BMC Bioinformatics，2011，12:118.

[31]LOZUPONE CA，HAMADY M，KELLEY S T，et al. Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities[J]. Appl Environ Microbiol，2007，73(5):1576-1585.

(責(zé)任編輯：林海清)

Metagenomic Analysis of Microbial Community in a Microbial Fermentation-bed for Pig Raising

LIU Bo, WANG Jie-ping, CHEN Qian-qian, LIU Guo-hong, CHE Jian-mei, CHEN De-ju,ZHENG Xue-fang, GE Ci-bin

(AgriculturalBio-resourcesResearchInstitute,Fujianacademyofagriculturalsciences,Fujzhou,Fujian350003,China)

The microbial community can decompose pig feces and eliminate funk odor in the microbial fermentation-bed for pig raising. Thus, the pig raising technology in the microbial fermentation-bed can not only achieve zero release of the pig raising pollutants, but also carry out resourceful utilization of the pig feces resources. Moreover, some key techniques in the microbial fermentation-bed for pig raising are directly or indirectly related to microorganisms. However, very few researches have been done on the microbial species, quantity, structure, and other aspects of the microbial fermentation-bed. Our group had comprehensively and systematically performed the microbiomics analysis of the microbial fermentation-bed to reveal its microbial community and the ecological succession in different space and depth with different fermenting degrees and litter composition in different seasons. In this paper, we summarized the analytical methods and procedures for sample collection, sample processing, and data analysis of the microbial fermentation-bed metagenomic analysis, of which included bacterial operational taxonomic units (OTUs) identification; species accumulation curve; core microbiome; principal component analysis of bacterial abundance; rank-abundance curve; histogram, heatmap, star analyses of bacterial species composition and abundance; α- and β-diversity analyses of bacterial species complexity; linear discriminative analysis (LDA) coupled with effect size measurement for significant difference, and bacterial species redundancy analysis (RDA). Our work could provide the complete analytical methods and ideas for more deep and broad researches on the microbial community in microbial fermentation-bed.

metagenomics; microbial community; microbial fermentation-bed; analytical procedure

劉波，王階平，陳倩倩，等.養(yǎng)豬發(fā)酵床微生物宏基因組基本分析方法[J].福建農(nóng)業(yè)學(xué)報，2016，31(6)：630-648.

LIU B，WANG J-P，CHEN Q-Q，et al.Metagenomic Analysis of Microbial Community in a Microbial Fermentation-bed for Pig Raising[J].FujianJournalofAgriculturalSciences，2016，31(6)：630-648.

2016-04-12初稿；2016-05-22修改稿

劉波(1957-),博士，研究員, 研究方向：微生物生物技術(shù)與農(nóng)業(yè)生物藥物(E-mail:fzliubo@163.com)

國家自然科學(xué)基金(31370059)；福建省科技重大專項(xiàng)(2015NZ0003-1)；福建省種業(yè)創(chuàng)新與產(chǎn)業(yè)化工程項(xiàng)目(FJZZZY-1544)；福建省科技計劃項(xiàng)目——省屬科研院所基本科研專項(xiàng)(2015R1018-1)

Q 933

A

1008-0384(2016)06-630-19