养猪发酵床微生物宏基因组基本分析方法

刘　波，王阶平，陈倩倩，刘国红，车建美，陈德局，郑雪芳，葛慈斌

(福建省农业科学院农业生物资源研究所，福建　福州　350003)



养猪发酵床微生物宏基因组基本分析方法

刘波，王阶平，陈倩倩，刘国红，车建美，陈德局，郑雪芳，葛慈斌

(福建省农业科学院农业生物资源研究所，福建福州350003)

养猪微生物发酵床分解猪粪、消除恶臭，将猪粪形成生物菌肥，实现了养猪污染的微生物治理和猪粪的资源化利用，微生物起到关键作用。然而，关于发酵床微生物群落的种类、数量、结构等方面系统研究鲜见报道。课题组采用高通量宏基因组学方法，系统地开展了微生物发酵床的微生物组研究，揭示不同空间、不同深度、不同发酵程度、不同垫料组成、不同季节发酵床中的微生物组成及其区系演替规律。本文就相关的宏基因组测序的样本采集、样本处理、测序原理、分析模型、数据统计、结果表述等分析流程和方法进行了归纳，阐明微生物组操作分类单元(OTU)鉴定，物种累积曲线，核心微生物组，种类丰度主成分分析，种类秩-多度曲线，物种组成丰度柱状图、热图、星图等分析的原理与方法，列举微生物组种类复杂度α-多样性分析、β-多样性分析、种类显著性差异LEfSe分析，冗余分析等实例，为深入研究发酵床微生物群落提供了完整的分析思路和范例。

宏基因组学；微生物群落；微生物发酵床；分析方法

微生物发酵床养猪是利用植物废弃物如谷壳、秸秆、锯糠、椰糠等农业副产品制作发酵床垫料层，添加微生物发酵剂，经发酵后铺垫到猪舍，生猪(母猪、公猪、小猪、育肥猪)生活在垫料上，排泄物混入垫料发酵，利用生猪的拱翻习性和机械翻耕，使猪粪尿和垫料充分混合，通过发酵床的微生物分解猪粪，消除恶臭，发酵床形成的有益微生物群落为限制猪病源的蔓延创造了有利条件，不仅实现了猪场的零排放，而且猪粪和垫料经发酵转化还可以生产优质微生物菌肥[1]。同时，饲喂饲用益生菌，替代抗生素，抑制猪病发生，解决猪肉药残问题。微生物发酵床养猪具有“五省、四提、三无、两增、一少、零污染”优点，即：五省，即省水、省工、省料、省药、省电；四提，提高品质、提高猪抗病力、提早出栏、提高肉料比；三无，无臭味、无蝇蛆、无环境污染；两增，增加经济效益，增加生态效益；一少，减少猪肉药物残留；零污染，猪粪尿由微生物在猪舍内降解消纳，而实现零污染[2]。

微生物发酵床养猪技术综合利用了微生物学、生态学、发酵工程学原理，以活性功能微生物作为物质能量“转换中枢”的一种生态养殖模式。该技术的核心在于利用活性强大的有益功能微生物复合菌群，长期和持续稳定地将动物粪尿废弃物进行降解[3-4]。尽管以往的研究涉及发酵床的微生物群落，如刘波等研究了零排放猪场基质垫层微生物群落脂肪酸生物标记多样性[5]，郑雪芳等利用磷脂脂肪酸生物标记分析猪舍基质垫层微生物亚群落的分化[6]；李娟等进行鸡发酵床不同垫料理化性质及微生物菌群变化规律研究[7]；张学锋等研究了不同深度垫料对养猪土著微生物发酵床稳定期微生物菌群的影响[8]；赵国华等进行了生物发酵床养猪垫料中营养成分和微生物群落研究[9]；王震等研究了发酵床垫料中优势细菌的分离鉴定及生物学特性[10]；但至今，发酵床微生物种类、数量、结构等方面的系统研究未见报道。其原因之一是通过传统的可培养方法分离鉴定发酵床微生物种类的工作量非常大，要完成丰富的微生物种类的调查难度很大，而且，可培养方法不能全面调查微生物的种类。宏基因组测定，为发酵微生物群落的调查提供了方法。

宏基因组学，是研究直接来自环境的微生物基因材料的学科，被认为是微生物学发展中的一个里程碑。通过测序分析微生物组主要包括16Sr RNA和宏基因组两大技术，它不仅使得对未培养或者不可培养的微生物的研究成为可能，也使得研究同一环境中的微生物在自然条件下的相互作用以及微生物和环境因子的相互作用成为可能[11]。宏基因组学应用范围很广，如分析传统食醋发酵过程微生物多样性[12]、牙菌斑微生物群落[13]、牛胃菌群组成[14]、人和动物胃肠道微生物群落[15]、深海样本[16]、温室黄瓜根结线虫发生地土壤微生物[17]、铅锌尾矿酸性废水微生物[18]、普洱茶渥堆发酵过程中微生物[19]、红树林土壤微生物[20]。

宏基因组数据分析主要包括序列处理、分类、注释及统计分析4个环节。随着测序技术的升级,测序成本将逐步降低，而大数据分析将成为核心内容。数据具备标准化、可积累性特点，通过数据建模是未来应用的基础，培养和基于培养的功能验证将是未来的重点之一。微生物组学将阐述并调整环境与微生物组之间的关系，此领域相关研究有巨大的发展空间[21]

宏基因组分析方法是理解环境中微生物组种类、数量、结构的工具。本研究以养猪发酵床垫料微生物宏基因组为材料，阐明宏基因组测序的样本采集、样本处理、数据分析等方法，以期为发酵床微生物群落的分析提供基础和参考。

1　材料与方法

1.1发酵床垫料样本采集及垫料理化性质测定

试验地点位于福清渔溪现代设施农业样本工程示范基地微生物发酵床大栏养猪舍，大栏发酵床养殖面积为1 617 m2(33 m ×49 m)，深度80 cm，发酵床垫料由70%椰糠和30%谷壳组成。发酵床饲养1 617头育肥猪，饲养密度为1头·m-2。按下表将发酵床划分为8 ×4个方格，冬季(11月)和春季(3月)取样两次，空间取样在各个区域内随机挖取10 cm深度的垫料样本，深度取样分别取0、20、40、60 cm深度的发酵床样本。冬季采样空间平面样本32个点和3个点4个深度样本12个，春季采样空间平面样本32个点和2个点4个深度样本8个，共84个样本。分别测定各垫料样本的温度、pH值、含水量、含盐量、硝态氮和粗纤维含量等理化性质。

1.2发酵床垫料的发酵程度等级划分

将垫料铺平，将色差计的出光口正对垫料，使其不露光，测定垫料跟空白垫料的色差值(△E)。每份垫料重复3次，取平均值作为该份垫料的最终值。在大量分析数据的基础上，快速而准确地判定各样本的发酵程度，并进行发酵等级划分。其中0<△E<30 为发酵等级一级；30<△E<50为发酵等级二级；50<△E<65为发酵等级三级；65<△E<80为发酵等级四级。

1.3宏基因组高通量测序

1.3.1总DNA的提取按土壤DNA提取试剂盒FastDNA SPIN Kit for Soil的操作指南，分别提取各垫料样本的总DNA，于-80℃冰箱冻存备用。

1.3.216S rDNA和ITS测序文库的构建采用扩增原核生物16S rDNA的V3～V4区的通用引物U341F和U785R对各垫料样本的总DNA进行PCR扩增，并连接上测序接头，从而构建成各垫料样本的真细菌和古菌16S rDNA V3～V4区测序文库。采用扩增真菌5.8S和28S rDNA之间的转录间区的通用引物ITS-F1-12和ITS-R1-12对各垫料样本的总DNA进行PCR扩增，并连接上测序接头，从而构建成各垫料样本的真菌ITS测序文库。

1.3.3高通量测序使用Illumina MiSeq测序平台，采用PE300测序策略，每个样本至少获得10万条reads。

1.4宏基因组测序数据质控

通过PADNAseq软件[22]利用重叠关系将双末端测序得到的成对短序列(reads)拼接成1条序列，得到高变区的长reads。然后使用撰写的程序对拼接后的reads进行处理而获取去杂短序列(clean reads)：去除平均质量值低于20的reads；去除reads含N的碱基数超过3个的reads；reads长度范围为250～500 nt。

1.5操作分类单元(OTU)聚类分析

为便于下游的物种多样性分析，将标签(Tags)聚类成操作分类单元(Operational Taxonomic Units，OTU)。首先把Tags中的单例标签(singletons，即对应的reads序列只有1条)过滤掉，因为这些singletons可能由于测序错误造成，故将这部分序列去除，不加入聚类分析。然后，利用usearch软件，在0.97相似度下进行聚类，对聚类后的序列进行嵌合体过滤后，得到用于物种分类的OTU，每个OTU被认为可代表1个物种[23]。

1.6抽平处理

为避免因各样本数据大小不同而造成分析时的偏差，在样本达到足够测序深度的情况下，根据α-多样性(Alpha diversity)指数分析结果对每个样本进行随机抽平处理。抽平的参数必须在保证测序深度足够的前提下去选取。

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

根据样本的共有OTU以及OTU所代表的物种，可以找到核心微生物组(Core microbiome，即覆盖90%样本的微生物组)。OTU维恩图(Venn)分析。Venn图可以很好地反映组间共有以及组内特有的OTU数目。利用R语言编写的VennDiagram软件里的venn.diagram函数实现。OTU的主坐标分析(PCoA)：PCoA分析可以初步反映出不同处理或不同环境间的样本可能表现出分散和聚集的分布情况，从而可以判断相同条件的样本组成是否具有相似性。利用R语言编写的ade4软件里的dudi.pca函数实现。

1.8物种分类和丰度分析

物种注释分析：首先，分别从各个OTU中挑选出一条序列作为该OTU的代表序列，将该代表序列与已知物种的16S rRNA数据库(网站Green Genes，http://greengenes.lbl.gov)[24]进行比对，从而对每个OTU进行物种鉴定；然后，根据每个OTU中序列的条数，得到各个OTU的丰度值。物种丰度分析：在门、纲、目、科、属水平，将每个注释上的物种或OTU在不同样本中的序列数整理在一张表格，形成物种丰度的柱状图、星图及统计表等。物种热图分析：物种热图利用颜色梯度可以很好反映出样本在不同物种下的丰度大小以及物种聚类、样本聚类信息，可利用R语言的gplots包的heatmap.2函数实现。

1.9样本复杂度分析

1.9.1单个样本复杂度分析α-多样性也称为生境内的多样性(within-habitat diversity)，是对单个样本中物种多样性的分析，常用的度量指标有：测定物种指数(observed species)、chao1指数、香农指数(Shannon)、辛普森指数(Simpson)以及谱系多样性指数(phylogenetic diversity，PD_whole_tree)等。利用QIIME软件计算样本的α-多样性指数值，并做出相应的稀释曲线[25]。稀释曲线是利用已测得16S rDNA序列中已知的各种OTU的相对比例，来计算抽取n个(n小于测得reads序列总数)reads时各α-多样性指数的期望值，然后根据一组n值(一般为1组小于总序列数的等差数列)与其相对应的α-多样性指数的期望值做出曲线来。并作出α-多样性指数的统计表格。

1.9.2样本间复杂度比较分析β-多样性也称为生境间的多样性(between-habitat diversity)，反映了不同样本在物种多样性方面存在的差异大小。分析各类群在样本中的含量，进而计算出不同样本间的β-多样性值。本分析中通过QIIME软件，采用迭代算法，分别在加权物种分类丰度信息和非加权物种分类丰度信息的情况下进行差异计算，得到最终的统计分析结果表并做出组间的距离箱线图(Box Plot)及主成分分析(Principal Component Analysis，PCA)展示图。

1.10显著性差异分析

1.10.1判别效应分析(LEfSe)LEfSe分析即LDA Effect Size分析，LEfSe采用线性判别分析(LDA，Linear discriminative analysis)来估算每个组分(物种)丰度对组分间差异效果影响的大小，可以实现多个分组之间的比较，还进行分组比较的内部进行亚组比较分析，从而找到组间在丰度上有显著差异的物种标记(biomaker)，找出对组分划分产生显著性差异性影响的群落或物种[26]。本研究采用LEfSe Tools进行分析[27]。基本分析过程如下：首先，在多组样本中采用Kruskal-Wallis非参数多组检验法来检测不同分组间丰度差异显著的物种；然后，在上一步中获得的显著差异物种，用成组的威尔科克森(Wilcoxon)秩和检验来进行组间差异分析；最后，用线性判别分析(LDA)对数据进行降维和评估差异显著的物种的影响力(即LDA score)。LEfse分析可以进行本地分析也可以在线分析，在线分析的网址是：https://huttenhower.sph.harvard. edu/galaxy/。

1.10.2组间差异分析使用秩和检验的方法对不同分组之间的进行显著性差异分析，以找出对组间划分产生显著性差异影响的物种。本分析对于两组间的差异分析采用R语言stats包的wilcox.test函数，对于两组以上的组间差异分析采用R语言stats包的kruskal.test函数。

1.10.3物种典型分析(CCA)/冗余分析(RDA)即典型分析(canonical correspondence analusis，CCA)/冗余分析(redundancy analysis RDA)。RDA是基于线性模型，CCA是基于单峰模型。CCA/RDA主要用来检测环境因子、样本、菌群三者之间的关系或者两两之间的关系。可以基于样本的所有OTU，也可基于样本的优势物种或者差异物种。CCA采用R语言编写的vegan软件的cca函数，RDA采用vegan软件的rda函数。

1.11发酵床微生物高通量测序与统计工作流程

养猪发酵床微生物组高通量测序(图1)和生物信息学分析工作流程(图2)分别归纳为：通过Illumina平台(Hiseq或者Miseq)进行配对末端(Paired-end)测序，通过reads之间的覆盖部分(overlap)关系拼接成长reads，并对拼接后的reads进行质控，除去平均质量值低于20的reads；除去reads含N的碱基数超过3个的reads；reads长度范围为250～500 nt，得到clean reads(图1)。

测序处理后得到的clean reads(去杂短序列)，进入生物信息学分析。通过拼接→去杂统计→OTU聚类→抽平处理→序列选择→分类比对→丰度统计等，得到发酵床微生物种类的OTU，为微生物群落分析提供基础(图2)。

2　结果与分析

2.1微生物发酵床细菌OTU提取

2.1.1发酵床细菌OTU提取将序列完全一样的去低值后reads(clean reads)归为一种标签(tag)，统计每条tag对应的丰度(即reads数目)，然后将tags根据其丰度大小进行排序，将其中的单例标签(singletons，即对应reads只有一条的序列)过滤掉，因为singletons可能由于测序错误造成。利用usearch7.0软件，在0.97相似度下进行聚类，对聚类后的序列进行嵌合体过滤后，最后将所有去低值reads比对到OTU序列上，将能比对上OTU的reads提取出来，得到最终的比对reads(mapped reads)(表1)。从表1可知，发酵床不同深度层次的垫料微生物tags、singleton、对比序列数(Mapped reads)、clean tags、OTUs存在显著差异。用微生物发酵床夏季和冬季2个季节垫料进行微生物宏基因组测序，2个季节共采样84个样本，分析结果表明，总体趋势是发酵床表层和底层的相应数值较小，中间层的相应数值较大。统计各个样本每个16S rRNA OTU中的丰度信息，一共产生154 556个16S rRNA-OTU(表1)，平均每个样本有(1862.12±505.96)个OTU。分析结果表明养猪发酵床垫料中的微生物种类非常丰富。

表1　养猪发酵床样本数据和16S rRNA部分OTU统计

注：第一列sample name是样本名称，S1.1.1代表样点S1,1的表层，顺序下去为20 cm、40 cm、60 cm层的样本。

2.1.2发酵床细菌OTU对原始数据进行质量检验(QC)后，用usearch7.0软件对数据进行去嵌合体和聚类的操作，usearch7.0软件聚类时，先将reads按照丰度从大到小排序，通过97%相似度的标准聚类，得到16S rRNA OTU，每个OTU被认为可代表一个细菌物种。在聚类过程中利用从头测序(de novo)方法去除嵌合体(chimeras)。接下来对每个样本的标签(tags)进行随机抽平处理以保证数据质量，提取对应的OTU序列。然后使用微生物生态学数量评估软件(Qiime，Quantitative Insights Into Microbial Ecology，QIIME v 1.9.0)，做α-多样性指数稀释曲线。根据稀释曲线选择合理的抽平参数，利用Qiime软件对得到的抽平后的OTU进行分析，首先从OTU中分别提取一条reads作为代表序列，将该代表序列与核苷酸数据库(RDP，Ribosomal Database Project)的16S rRNA比对，从而对每个OTU进行细菌物种分类。归类后，根据每个OTU中序列的条数，从而得到OTU物种含量丰度，最后根据该OTU丰度表进行后续分析(表1)。

2.1.3发酵床细菌OTU抽平处理由于不同样本对应的reads数量差距较大，为了保证后期分析结果合理，对每个样本的数据进行随机抽平处理，抽平参数根据α-多样性指数的稀释曲线来确定。α-多样性反映的是单个样本内部的物种多样性，包括测定物种指数(observed species)、丰富度Chao1指数，香农多样性指数(Shannon)以及辛普森优势度指数(Simpson)，谱系多样性指数(PD_whole_tree)等。测定物种指数(observed species)和chao1指数反映样本中细菌群落的丰富度(species richness)，即简单指细菌群落物种总的数量，而不考虑细菌群落每个物种的丰度情况。这2个指数对应的稀释曲线还可以反映样本测序量是否足以覆盖所有类群。为了对各样品进行定量比较分析，根据发酵床细菌群落α-多样性分析结果，结合测序饱和度和样本信息完整性，对各样本的测序数据进行随机抽平处理。根据分析结果，我们从发酵床部分样本的测序结果中随机抽取60 623条reads，进行各样本的各OTU丰度计算，分析部分结果见表2。可以看出发酵床S(1，1)点表层垫料的细菌种类(OTU)为1 643，20 cm深度垫料细菌种类(OTUs)为2 055，等等，说明发酵床不同深度垫料的细菌种类(OTU)差异显著。

表2　发酵床部分样本微生物群落OTU抽平统计

2.1.4发酵床OTU稀释曲线微生物多样性分析中需要验证测序数据量是否足以反映样本中的物种多样性，稀释曲线(丰富度曲线)可以用来检验这一指标，评价测序量是否足以覆盖所有类群，并间接反映样本中物种的丰富程度。稀释曲线是利用已测得16S rDNA序列中已知的各种OTU的相对比例，来计算抽取n个(n小于测得reads序列总数)reads时出现OTU数量的期望值，然后根据一组n值(一般为一组小于总序列数的等差数列，如0 ×104、2×104、4 ×104、6×104、8 ×104等)与其相对应的OTU数量的期望值做出曲线来。当曲线趋于平缓或者达到平台期时，也就可以认为测序深度已经基本覆盖到样本中所有的物种；反之，则表示样本中物种多样性较高，还存在较多未被测序检测到的物种。图3展示了发酵床垫料样本细菌chao指数和测定物种指数(observed species)的稀释曲线图。结果显示，各样本均具有较好的物种丰富度，且不同样本的物种丰富度存在明显差异；在测序达到一定深度后，各样本的chao指数和测定物种指数(observed species)的稀释曲线趋于平缓或者达到平台期，表明测序深度已经基本覆盖到样本中所有的物种。同时，与其他样品相比，图3中最下方的4条曲线表示这4个样品中的物种多样性明显偏低。

2.1.5发酵床细菌OTU群落α-多样性香农多样性指数(Shannon)以及辛普森优势度指数(Simpson)反映群落的物种多样性(species diversity)，受样本细菌群落中物种丰富度(species richness)和物种均匀度(species evenness)的影响。相同物种丰富度的情况下，群落中各物种均匀度越大，则认为群落多样性越高，同样，Shannon指数和Simpson指数值越大，说明个体分分布越均匀。如果每一个体都属于不同的种，Shannon指数和Simpson指数就大，如果每一个体都属于同一种，则Shannon指数和Simpson指数就小。图4展示了发酵床垫料部分样本的Shannon指数和Simpson指数稀释曲线图。结果表明，绝大多数样本的Shannon指数和Simpson指数稀释曲线很快就达到了平台期，说明各个样本均具有较好的物种多样性，而且不同个体属于不同种的可能性大，即物种均匀度大。

2.1.6发酵床部分样本微生物(OTU)物种谱系演化α-多样性指数(PD_whole_tree and- Goods_coverage)稀释曲线该指数反映了样本中细菌物种谱系演化的差异。谱系多样性指数(PD_whole_tree)指数越大说明物种谱系演化差异越大。测序深度指数(Goods_coverage)反映了测序的深度，测序深度指数(Goods_coverage)指数越接近于1，说明测序深度已经基本覆盖到样本中所有的物种。图5展示了样本谱系多样性指数(PD_whole_tree)和测序深度指数(Goods_coverage)的稀释曲线。图中一条曲线代表一个样本，横轴表示从某个样本中随机抽取的clean reads数目，纵轴表示该reads数目对应的alpha多样性指数的大小。在测序深度指数(Goods_coverage)的稀释曲线中，随着测序深度的增加，当曲线趋于平缓时，表示此时测序深度已经基本覆盖到样本中所有的物种，测序数据量比较合理。从图5可以看出，每一个样品均获得了高质量的测序深度，能覆盖样品中的所有物种信息。

2.2微生物发酵床核心微生物组(Core microbiome)分析

2.2.1发酵床微生物组共有OTU数与覆盖样本数关系基于16S rRNA基因信息的系统分类结果与基于全基因组信息的分类结果很相似，并且消除了克隆问题，可以综合研究各个可变区，分析方法也相对成熟，广泛应用于核心微生物组的研究。以微生物发酵床垫料样本为例，对采集的夏季和冬季的84个样本分别进行宏基因组分析，列出各样本的细菌种类OTUs，对各样本共有的OTUs进行统计。分析结果(图6)表明微生物发酵床垫料各样本共有细菌种类(OTUs)数与样本数的关系，图中横坐标表示的是覆盖一定比例以上(如>50%、>60、>70%、>80%、>90%等)样本的共有细菌种类(OTUs)数目的比例，纵坐标是统计的覆盖大于此比例样本的共有细菌种类(OTUs)数目。例如，若1个细菌种类(OTUs)覆盖50%以上的样本，则在>0.5的横坐标所对应的纵坐标的细菌种类(OTUs)数目为1200种(OTUs)，覆盖60%以上的样本细菌种类数目(OTUs)为880种，以此类推，覆盖样本越高的细菌种类(OTUs)其数量越少，覆盖100%样本的细菌种类(OTUs)则非常少(图6)。

2.2.2微生物发酵床冬季和春季共有和特有微生物组(OTUs)分析在0.97的相似度下，得到了发酵床垫料每个样本的OTUs个数，利用韦恩图(Venn)展示多样本共有和各自特有细菌物种(OTUs)数目，直观展示样本间物种重叠情况。结果表明，微生物发酵床冬季样本(S1)(40个)和春季样本(S2)(44个)鉴定到的细菌种类(OTUs)数目分别为5 472和6 858个，冬季和春季特有的OTUs数目分别为598和1 984个，春季的OTUs比冬季略丰富，而共有的OTU数目为4 874个，说明垫料的微生物群落具有一定的稳定性(图7)。

2.2.3发酵床微生物组(OTU)丰度主成分分析主成分分析(Principal Component Analysis，PCA)是一种分析和简化数据集的技术。主成分分析经常用于减少数据集的维数，保持数据集中的对方差贡献最大的特征。通过保留低阶主成分，忽略高阶主成分而做到的。保留的低阶成分往往能够最大程度的保留数据的重要特征。PCA运用降维的思想，通过分析不同样本OTUs(97%相似性)组成可以反映样本的差异和距离，将多维数据的差异反映在二维坐标图上，坐标轴取值采用对方差贡献最大的前2个特征值。如果2个样本距离越近，则表示这2个样本的组成越相似。不同处理或不同环境间的样本，可表现出分散和聚集的分布情况，从而可以判断相同条件的样本组成是否具有相似性。以微生物发酵床2个季节的部分样本进行分析说明，结果表明，微生物发酵床2个季节的部分样本的微生物组成具有一定的相似性(图8)。分析结果表明，2个季节共有的种类主要集中在第1相限，种类较多，表明微生物群落具有稳定性；春季S2独有的种类主要分布在第3相限，分布较为分散，表明中间差异较大；冬季S1独有的种类主要分布在第4相限，分布较为集中，表明中间差异较小。

2.2.4发酵床微生物组(OTU)秩-多度曲线秩-多度曲线(Rank_Abundance Curve)能同时解释样本所含物种的丰富度和均匀度。物种的丰富程度由曲线横轴上长度来反映，物种的相对丰度由纵坐标长度来反映；物种组成的均匀程度由曲线整体斜率来反映，曲线斜率越大则表示样本中各物种所占比例差异较大，分布越均匀。如图9所示，微生物发酵床冬季和春季84个样本的曲线在横轴上均在103左右，说明各样本具有较高的丰富度；各曲线的斜率相对较大，表明各样本的均匀度较好，且样本间的均匀度存在一定差异。

2.2.5发酵床微生物组(OTU)物种累积曲线物种累积曲线(species accumulation curves)用于描述随抽样量增加物种数量增加的状况，是理解和预测物种丰富度(species richness)的有效工具，利用物种累积曲线不仅可以判断抽样量充分性，在抽样量充分的前提下，运用物种累积曲线可以对物种丰富度进行预测。对微生物发酵床冬季和春季84个垫料样本细菌物种(OTUs)测序结果进行抽样，分析结果如图10所示。随着抽样reads的增加，垫料样本的物种累积曲线迅速趋于平缓，最后达到平台期，说明样本测序深度符合要求，同时说明各样品具有较好的物种丰富度。

2.2.6发酵床微生物组(OTU)分类阶元数量分布统计结果(表3)表明从各个OTU中挑选出丰度最高的1条序列，作为该OTU的代表序列。使用核苷酸数据库(RDP)比对方法，将该代表序列与已知物种的16S rRNA数据库进行比对，从而对每个OTU进行物种归类。微生物发酵床样本共鉴定出细菌OTUs总数为7 456，每个样本可鉴定的OTUs范围为1 643～2 467，其中可鉴定到科水平OTUs有4 886，可鉴定到属水平OTUs有2 351，可鉴定到种水平OTUs有246。

表3　微生物发酵床选择样本细菌OTU统计

2.2.7发酵床微生物组(OTU)丰度柱状图分析微生物发酵床所有84个样本门分类阶元物种(OTUs)含量前20位的物种的组成与丰度(图11)结果表明，丰度在前5位的分别是放线菌门(Actinobacteria)>厚壁菌门(Firmicutes)>变形菌门(Proteobacteria)>拟杆菌门(Bacteroidetes)>绿弯菌门(Chloroflexi)，与一般的土壤样本(变形菌门含量最高)存在明显差异，说明微生物发酵床具有独特的生境类型。

2.2.8发酵床冬季和春季微生物组(OTU)丰度分析发酵床冬季和春季微生物组(OTU)丰度如图12所示，将发酵床冬季和春季分别作为样本，两个样本细菌物种门分类阶元(OTU)的相对丰度在前5位的也是放线菌门Actinobacteria>厚壁菌门Firmicutes>变形菌门Proteobacteria>拟杆菌门Bacteroidetes>绿弯菌门Chloroflexi。

2.2.9发酵床微生物组(OTU)物种热图发酵床微生物组(OTU)物种热图分析结果如图13所示，物种热图是以颜色梯度来代表数据矩阵中数值的大小，并能根据物种或样本丰度相似性进行聚类的一种图形展示方式。聚类结果加上样本的处理或取样环境分组信息，可以直观观察到相同处理或相似环境样本的聚类情况，并直接反映了样本的群落组成的相似性和差异性。微生物发酵床冬季和春季样本中的细菌种类分别在门，纲，目，科，属，种分类等级进行heatmap聚类分析。纵向聚类表示所有物种在不同样本间的相似情况，距离越近，枝长越短，说明样本的物种组成及丰度越相似。横向聚类表示该物种在各样本丰度相似情况，与纵向聚类一样，距离越近，枝长越短，说明两物种在各样本间的组成越相似。

2.2.10组(OTU)物种星图(star)分析图14展示了发酵床冬季和春季84份样本属水平丰度在前10位的细菌物种(OTU)整理的星图(star)，每个星形图中的扇形代表一个物种，用不同颜色区分，用扇形的半径来代表物种相对丰度的大小，扇形半径越长代表此扇形所对应的物种的相对丰度越高。结果表明不同样本在属水平的物种组成存在明显差异性，而且，2个季节样本中的棒杆菌属Corynebacterium、葡萄球菌属Staphylococcus、梭菌属Clostridium是多数样本的优势属，这3个属主要来源于人和动物，凸显了养殖微生物发酵床中的微生物组成与区别于一般土壤样本微生物组成的差异。

2.3发酵床微生物组(OTU)α-多样性种类复杂度分析

2.3.1发酵床单个样本微生物组(OTU)种类复杂度分析为了分析单个微生物发酵床样本的细菌种类复杂度，分别进行了衡量各个样品细菌种类组成和丰度复杂度的α-多样性指数常用度量指标[包括测定物种指数(observed species)、chao1指数、香农指数(Shannon)、辛普森指数(Simpson)、测序深度指数(Goods_coverage)、谱系多样性指数(PD_whole_tree]的计算。表4展示了部分冬季和春季部分样本的α-多样性指数的统计结果，可以看出，测序深度指数(Goods_coverage)和Simpson指数均接近于100%或1，表明有较好的α-多样性。

2.3.2发酵床冬季和春季微生物组(OTU)箱线图分析图15展示了发酵床冬季(S1)和春季(S2)样本α-多样性指数箱线图(Box-plot)，能更直观显示各样品的α-多样性指数差异。其中，中间的黑色横线是数据的中位数(median)，即数据中占据中间位子的数，即数据中有一半大于中位数(在其之上)，另一半小于中位数(在其之下)。盒子的上下两边称为上下四分位数(Hinges)，其意义为：数据中有四分之一的数目大于上四分位数(即盒子的上边)，即在红色盒子之上；另外有四分之一的数目小于下四分位数(即盒子的下边)，也就是在盒子之下。也就是说有一半的数目在中间封闭盒子的范围内，有一半分布在盒子上下两边。在盒子上下两边分别有一条纵向的线段，叫触须线。上截止横线是变量值本体最大值，下截止横线是变量值本体最小值。本体指的是除异常值和极值以外的变量值，称为本体值。异常值标记为o，极值标记为*。高于触须线上截止横线的值的取值范围为：①异常值：x>上四分位数+1.5IQR；②极值：x>上四分位数+3.0IQR；低于触须线下截止横线的值的取值范围为：①异常值：x<下四分位数-1.5IQR；②极值：x<下四分位数-3.0IQR；从而表明盒子外面数值点的分布。IQR(interquartile range)=上四分位数-下四分位数。

从图15可以看出：①春季的度量比冬季的分散得多；②冬季和春季的Simpson指数几乎没有差异，而且数值均较高，春季的测序深度指数(Goods_coverage)指数略优于冬季，其他4个指数均为冬季优于春季(主要根据中位数和盒子的位置来判断)。

2.3.3发酵床微生物组(OTU)α-多样性指数秩和检验对细菌物种(OTU)α-多样性指数进行秩和检验(rank sum test)，分析不同条件下显著差异的α-多样性指数。分析结果见表5，表明冬季和春季的Shannon和Simpson指数均无显著差异，其P-value值分别为0.1269和0.1036(远大于0.05或0.01)，说明在多样性上冬季和春季微生物组无显著差异；冬季和春季微生物组的测定物种(observed_species)、测序深度指数(Goods_coverage)、chao1指数差异显著(P<0.05)，谱系多样性指数(PD_whole_tree)差异极显著(P<0.01)。

表5　微生物发酵床细菌物种(OTU)α-多样性指数秩和检验

2.4发酵床样本间微生物组(OTU)β-多样性种类复杂度分析

2.4.1发酵床样本间微生物组(OTU)UniFrac距离通过利用系统进化的信息来比较样本间的细菌群落差异，其计算结果可以作为一种衡量β-多样性指数，它考虑了物种间的进化距离，该指数越大表示样本间的差异越大。微生物发酵床部分样品β-多样性UniFra距离数据矩阵如表6(weighted_unifrac)和表7(unweighted_unifrac)所示。

2.4.2发酵床样本间微生物组(OTU)UniFrac距离聚类分析分析结果见图16。展示了发酵床样品的UniFrac距离分布热图(heatmap)和聚类图，通过对UniFrac距离的聚类，具有相似的β-多样性的样本聚类在一起，反映了样本间的相似性。

表6　加权物种丰度信息计算得到的样本β-多样性数据矩阵(weighted_unifrac)

表7　非加权物种丰度信息计算得到的样本β-多样性数据矩阵(unweighted_unifrac)

2.4.3发酵床样本间微生物组(OTU)主成分分析(PCA)为了进一步展示样本间物种多样性差异，使用PCA分析的方法展示各个样本间的差异大小。图17给出了PCA对样本间物种多样性的分析结果，如果两个样本距离较近，则表示这2个样本的物种组成较相近。从结果来看，与加权的相比，非加权的冬季和春季样本间物种多样性的差异更加显著。

2.5发酵床微生物组(OTU)差异性分析

2.5.1发酵床微生物组(OTU)差异效应判别分析(LEfSe)分析结果见图17。LEfSe(LDA Effect Size)分析结果可以用LDA值(Linear Discriminant Analysis)作分布柱状图和进化分支图(图18)。在LDA值分布柱状图中，展现了发酵床不同组中(即红色代表春季；绿色代表冬季)微生物组(OTU)丰度有显著差异的物种，即具有统计学差异的微生物标记(Biomaker)，柱状图的长度代表差异物种的影响值大小。在进化分支图中的着色原则：无显著差异的物种统一着色为黄色，差异物种(Biomarker)跟随组别进行着色，红色节点表示在春季组别起到重要作用的微生物类群，绿色节点表示在冬季组别中起到重要作用的微生物类群，其他圈颜色意义类同。图中英文字母表示的物种名称在右侧图例中进行展示。由内至外辐射的圆圈代表了由门至属(或种)的分类级别。在不同分类级别上的每一个小圆圈代表该水平下的一个分类，小圆圈直径大小与相对丰度大小呈正比。从结果可以看出，不同季节各自对应不同的差异物种(Biomarker)，而且，春季的差异物种(Biomarker)明显多于冬季。

2.5.2发酵床不同分类阶元微生物组(OTU)差异分析分析结果见表8。分别从门、纲、目、科、属阶元，通过秩和检验对冬季和春季样本筛选显著差异(P<0.05，Kruskal-Wallis test)的物种，表8是在冬季和春季组间20个差异显著的不同分类水平OTUs(P<0.05)，共找到317个，其中g代表属(Genus)，f代表科(Family)，o代表目(Order)，c代表纲(Class)。

表8　发酵床冬季和春季不同分类阶元微生物组(OTU)差异分析

2.5.3发酵床属分类阶元微生物组(OTU)差异分析通过Kruskal-Wallis test分析可以找出在冬季和春季组间有明显差异(P<0.05)的属共有136个，表9展示了部分结果，可看到各属的基本信息。

表9　发酵床冬季和春季属分类阶元微生物组差异分析

注：第一列是物种名称；第二列至第四列分别为两组样本的均值；最后一列是秩和检验的P值。

2.5.4发酵床冬季和春季微生物组(OTU)热图和主成分分析为了直观地展示这些具有显著差异的属，分别对它们进行PCA分析(图19)和热图(图20)。结果可以看出，冬季和春季的差异属可以分别聚类在一起(图19)。

2.6发酵床微生物组(OTU)种类冗余(RDA)分析

梯度分析理论，如基于单峰模型的典型分析(CCA，Canonical correspondence analysis)和基于线性模型的冗余分析(RDA，Redundancy analysis)，可用于多指标之间的对比，进行描述环境变迁的可信度(显著性)及解释的能力(重要性)等问题的统计学定量刻画。其中RDA是一种直接梯度分析方法，能从统计学的角度来评价一个或一组变量与另一组多变量数据之间的关系。图21展示了微生物发酵床垫料中的微生物组与季节因子(即气温)RDA分析(即关联性分析)结果，为了保证图中字符不重叠，物种名取了前8个字符来代表，图中标出了丰度为前20位的物种。蓝色线条与红色线条成锐角，表示这些物种与季节因子(气温)呈正相关；与红色线条成钝角，表示这些物种与季节因子(气温)呈负相关。微生物组(OTU)落在第二象限的种类与季节呈正相关，如短芽胞杆菌属Brevibacillus、短杆菌属Brachybacterium、周氏杆菌属Zhouia、糖多孢菌属Saccharopolyspora、海杆菌属Marinobacter，碳酸噬胞菌属Aequorivita等，它们与季节因子(气温)均具有显著的正相关性，其他种类为负相关性。

3　讨　论

我国畜禽粪便是农业面源污染的主要来源，每年仅猪粪的产生量就突破50亿t，由于处理不当对环境造成严重污染，已经成为经济发达地区或水环境敏感地区优先控制的污染源。微生物发酵床能彻底地处理畜禽粪便污染，其作用机理在于发酵床垫料中的微生物[1]。尽管以往的研究有涉及发酵床微生物群落，如发酵床微生物群落脂肪酸标记多样性[5]、发酵床微生物脂肪酸亚群落分化[6]、鸡发酵床微生物菌群变化[7]、稳定期发酵床微生物菌群特性[8]、发酵床微生物群落研究[9]、发酵床优势细菌鉴定[10]等，但至今发酵床与微生物的密切相关性知之甚少。有鉴于此，本研究利用高通量的宏基因组学分析方法，全面而系统地开展了微生物发酵床的微生物组研究，揭示不同空间、不同深度、不同发酵程度、不同垫料组成、不同季节、不同管理的发酵床中的微生物组成、区系演替规律，以期探明猪粪降解机理、资源化利用机制以及污染治理原理。宏基因组分析方法的出现，为发酵床微生物组分析提供了强有力的工具，为使这一技术尽快地应用到发酵床微生物组的分析上，作者从样本采集、样本处理、测序原理、分析模型、数据统计、结果表述等方面系统地进行了介绍，以推进发酵床微生物组的研究。

发酵床微生物组(OTU)α-多样性种类复杂度分析，也称为生境内的多样性(within-habitat diversity)，是对单个样本中物种多样性的分析，已经有很多度量指标(包括测定物种指数observed species、chao1指数、香农指数Shannon、辛普森指数Simpson以及谱系多样性指数phylogenetic diversity，PD_whole_tree等)来进行评价[28]。β-多样性也称为生境间的多样性(between-habitat diversity)，反映了不同样本在物种多样性方面存在的差异大小，可衡量群落之间的差别，在许多研究领域，尤其是生态学研究中，具有重要的意义[29]。然而，目前β-多样性的评价方法不多，在这些方法中，加权UniFrac(weighted UniFrac)和非加权UniFrac(unweighted UniFrac)的应用非常广泛[30]。其中加权UniFrac考虑了物种的丰度，非加权UniFrac则不考虑物种丰度。与P-检验类似，UniFrac分析的先决条件是一个包含所有物种的有根、各枝长已知的系统发生树。2个群落之间的UniFrac距离定义为：对于系统发育树所有分枝，考查其指向的叶节点是否只存于同一个群落，那些叶节点只存在于同一群落的枝的枝长和占整个树的枝长总和的比例。直观来讲，就是计算了仅被一个群落占据的进化历史的相对大小，这个值越大，说明两个群落中独立的进化过程越多，也就说明这两个群落的差别越大。若两个群落完全相同，那么它们没有各自独立的进化过程，UniFrac距离为0；若2个群落在进化树中完全分开，即它们是完全独立的两个进化过程，那么UniFrac距离为1[31]。

本研究就相关的宏基因组测序的样本采集、样本处理、测序原理、分析模型、数据统计、结果表述等分析流程和方法进行了归纳，阐明微生物组操作分类单元(OTU)鉴定，物种累积曲线，核心微生物组，种类丰度主成分分析，种类秩-多度曲线，物种组成丰度柱状图、热图、星图等分析的原理与方法，列举微生物组种类复杂度α-多样性分析、β-多样性分析、种类显著性差异LEfSe分析，冗余分析等实例，为深入研究发酵床微生物群落提供完整的分析思路和范例。

[1]刘波，朱昌雄. 微生物发酵床零污染养猪技术的研究与应用[J]. 北京：中国农业科学技术出版社，2009.

[2] 刘波，蓝江林，唐建阳，等. 微生物发酵床菜猪大栏养殖猪舍结构设计[J]. 农业科学与技术:英文版，2014，9(5): 1521-1525.

[3] 王远孝，钱辉，王恬. 微生物发酵床养猪技术的研究与应用[J].中国畜牧兽医，2011，38(5): 206-209.

[4] 蓝江林，刘波，宋泽琼，等. 微生物发酵床养猪技术研究进展[J]. 生物技术进展，2012，(6): 411-416.

[5] 刘波，郑雪芳，林营志，等. 零排放猪场基质垫层微生物群落脂肪酸生物标记多样性分析[J]. 生态学报，2008，28(12)：5488-5498.

[6] 郑雪芳，刘波，林营志，等. 利用磷脂脂肪酸生物标记分析猪舍基质垫层微生物亚群落的分化[J]. 环境科学学报，2009，29(11):2306-2317.

[7] 李娟，石绪根，李吉进，等. 鸡发酵床不同垫料理化性质及微生物菌群变化规律的研究[J]. 中国畜牧兽医，2014，41(2): 139-143.

[8]张学峰，周贤文，陈群，等. 不同深度垫料对养猪土著微生物发酵床稳定期微生物菌群的影响[J]. 中国兽医学报，2013，33(9): 1458-1462.

[9] 赵国华，方雅恒，陈贵. 生物发酵床养猪垫料中营养成分和微生物群落研究[J]. 安徽农业科学，2015，48(8)：98-101.

[10]王震，尹红梅，刘标，等. 发酵床垫料中优势细菌的分离鉴定及生物学特性研究[J]. 浙江农业学报，2015，27(1)：87-91.

[11]常秦. 宏基因组数据分析中的统计方法研究[D]. 济南：山东大学，2012.

[12]聂志强，韩玥，郑宇，等. 宏基因组学技术分析传统食醋发酵过程微生物多样性[J]. 食品科学，2013，34(15):198-203.

[13]陈林. 老年人根面龋患者和健康人牙菌斑微生物群落的宏基因组学研究[D]. 武汉：武汉大学，2014.

[14]彭帅. 应用宏基因组方法检测猪致病微生物及分析牛胃菌群组成[D]. 长春：吉林大学，2015.

[15]许波，杨云娟，李俊俊，等. 宏基因组学在人和动物胃肠道微生物研究中的应用进展[J]. 生物工程学报，2013，29(12):1721-1735.

[16]江夏薇. 基于嗜耐盐菌基因组分析与深海宏基因组文库的酯酶研究[D]. 杭州：浙江大学，2013.

[17]赵志祥，芦小飞，陈国华，等. 温室黄瓜根结线虫发生地土壤微生物宏基因组文库的构建及其一个杀线虫蛋白酶基因的筛选[J]. 微生物学报，2010，50(8):1072-1079.

[18]韩玉姣. 凡口铅锌尾矿酸性废水微生物宏基因组及宏转录组研究[D]. 广州：中山大学，2013.

[19]吕昌勇. 普洱茶渥堆发酵过程中微生物宏基因组学的测定与分析[D]. 昆明：昆明理工大学，2013.

[20]蒋云霞. 基于红树林土壤微生物资源研发的宏基因组学平台技术的建立与应用初探[D]. 厦门：厦门大学，2007.

[21]盛华芳，周宏伟. 微生物组学大数据分析方法、挑战与机遇[J]. 南方医科大学学报，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.

(责任编辑：林海清)

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

刘波，王阶平，陈倩倩，等.养猪发酵床微生物宏基因组基本分析方法[J].福建农业学报，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-),博士，研究员, 研究方向：微生物生物技术与农业生物药物(E-mail:fzliubo@163.com)

国家自然科学基金(31370059)；福建省科技重大专项(2015NZ0003-1)；福建省种业创新与产业化工程项目(FJZZZY-1544)；福建省科技计划项目——省属科研院所基本科研专项(2015R1018-1)

Q 933

A

1008-0384(2016)06-630-19