分析方法对细菌群落16S rRNA基因扩增测序分析结果的影响

钟辉 刘亚军 王滨花 和梦洁 吴兰

（1. 南昌大学生命科学学院，南昌 330031；2. 南昌大学 鄱阳湖环境与资源利用教育部重点实验室，南昌 330031）

微生物是地球上数量最多、物种和功能多样性最高的生物，广泛分布于自然环境中，维持着生态系统的物质循环和能量流动［1-2］。16S rRNA基因具有高度的保守性及特异性，被广泛应用于细菌群落结构、组成和多样性的研究［3-4］。随着测序技术在环境微生物组学方面的大量应用，人们对高通量测序数据的可重复性和可比性提出了更高的要求。

以往的研究发现，受16S rRNA基因序列高变区的保守程度不同，测定不同高变区会显著影响细菌群落系统发育的多样性，其中V3-V4区特异性好，数据库信息全，是目前细菌多样性分析注释的常用选择［5］。此外，测序通量的大小也是影响细菌群落信息的重要因素，例如，Illumina MiSeq之所以逐渐取代了最初的454 GS Junior（Roche）测序平台，正是由于Illumina MiSeq具有较高的测序通量，测序结果能够更为精准的反应庞杂的微生物群落特征。除了测序片段和平台，测序深度也是影响可获得微生物群落信息的重要影响因素，较低的测序深度会导致一些低丰度物种的信息丢失［6-8］。然而，以往的研究多集中在上机测序过程，而对于海量的下机数据，不同的处理方式也可能会影响最终的分析结果，尤其是最小分类单元的划分越来越受到科研工作者重视［9-10］。

最小分类单元是指将16S rRNA基因扩增序列按特定的相似性阈值分配到操作分类单位中（例如OTU）［11］。最初，由于可比对的数据信息有限，将获得的基因序列按照97%相似性进行聚类分析，确定最小分类单元，这在一定程度上掩盖PCR及测序过程中产生的错误，同时简化了基因序列的分析。97%的相似性阈值仍然是当前最小分类单元最常用的划分方式［9］。然而，随着16S rRNA基因测序的广泛应用、相关数据库的不断完善，有研究认为将序列相似性阈值提升到98%［12］。99%［13］甚至100%［14］，可以更准确反应微生物群落的真实信息。此外，通过DADA2模型算法可以纠正PCR和测序过程中的错误，并将去重复后得到的单核苷酸变体（ASV）作为最小分类单元［15］，该划分方式也越来越多地用于环境微生物组学的分析。

目前，关于环境微生物的研究主要集中在其多样性、组成以及与环境的协同进化上［16］。微生物群落作为一个整体，其多样性与组成不仅能够直接表征环境生物的物种遗传信息，还在一定程度上体现其在自然环境中扮演的重要生态功能差异［17］。考虑到微生物的多功能性和对环境变化的高敏感性，探索影响其时空变化特征的环境变量，成为当前微生物生态学研究的重要方向。然而，目前在细菌基因序列分析过程中，对于最小分类单元划分方法，存在多种方法并行的奇特研究现状［18］。而不同分类单元划分方式对细菌群落16S分析结果的影响，尤其是对细菌群落多样性、组成以及其与环境因子关系的影响仍待进一步探究。基于此，我们选择了5个不同生境下的微生物群落样本，在Illumina MiSeq平台上对16S rRNA基因的V3-V4区进行高通量测序。我们同时采用97%、98%、99%、100%相似性聚类OTU和非聚类的ASV方法对测序结果进行分析，以评估不同分类单元划分方式对微生物群落结构以及多样性的影响。

1 材料与方法

1.1 材料

为了评估不同分类分辨率方法对微生物菌群组成和多样性的影响，选择5种生境下的环境样品进行研究，分别是森林土壤（FS）、湿地土壤（WS）、农田土壤（CS）、湖泊沉积物（LS）和湖泊水体（LW）样本。在每种生境下选择4个样地进行样品采集，每个样地采集3个重复，共计60个样本。将采集好的土壤样品低温运回实验室进行分装和预处理。对于土壤样品，一部分土样自然风干供理化性质分析；一部分保存在-80℃用于DNA提取测序分析。对于水体样品，理化指标在样品带回实验室24 h 内进行测定；同时取1 L 水样通过微孔滤膜（孔径0.22 mm）进行抽滤，滤膜-80℃保存用于DNA提取测序分析。

所有样本的pH、TOC（总有机碳）、TN（总氮）、TP（总磷）、NH+-N（氨态氮）、NO3--N（硝态氮）的测定参照刘亚军等［19］描述的方式（表1-2）。使用PowerSoil DNA 提取试剂盒提取了样本中的总微生物 DNA［20］。使用超微量分光光度计检验 DNA 的浓度以及纯度，通过用通用细菌引物338F（5′-ACTCCTACGGGAGGCAGCAG-3′）和806R（5′-GGACTACHVGGGTWTCTAAT-3′）扩增16S rRNA 基因序列的V3-V4区来构建扩增子文库。合格的文库被送往派森诺（上海）有限公司通过Illumina Miseq 2500平台（2×300 bp）进行测序。原始测序数据已上传至NCBI数据库，登录号分别为森林 土壤（PTJNA388530）、湿地土壤（PRJNA559977）、农田土壤（PRJNA329019）、湖泊沉积物（PRJNA632434）和湖泊水体（PRJNA528375）。

表1 土壤和沉积物样本信息表及环境参数Table1 Information and environmental parameters of soil and sediment samples

1.2 方法

1.2.1 基础序列分析 使用QIIME2（v2020.11）［21］平台对测序数据进行处理分析，不同生境的样本单独进行分析。首先使用cutadapt［22］插件切除序列两端的引物片段，引物切除时设置最大碱基差异率为15%，最低序列长度为225 bp（300×75%），引物去除后的序列分别按以下方式分别进行分析得到4个OTU数据集以及一个ASV数据集。

表2 水体样本信息表及环境参数Table2 Water samples information and environmental parameters

1.2.2 序列分析 OTU数据集构建：分别使用Vsearch［23］插 件 的join-pairs、quantity-filter、dereplicate-sequences模块对去除引物的双端序列进行拼接、质控和去冗余。处理过程中除以下参数外均使用默认参数：拼接设置最小重叠区域（minovlen）为15 bp，重叠区域最大碱基差异（maxdiffs）为3 bp；质量控制设置滑动窗口（quality-window）为30 bp，质量阈值（min-quality）为20。再通过Vsearch插件cluster-feature-de-novo模块对非冗余序列按97%、98%、99%以及100%的相似性阈值聚类，得到4个不同聚类阈值的分类操作单元数据集（后文按聚类阈值简称97 OTU、98 OTU、99 OTU以及100 OTU数据集），并过滤掉reads数为1的OTU。

ASV数据集构建：使用dumx summarize插件对去引序列的碱基质量进行统计，并按照序列的质量分布情况，使用DADA2［15］插件对序列进行质控、纠错以及嵌合体的去除（处理过程中除截断长度设置为260/240 bp外均采用默认参数），最后得到单核苷酸序列变体（ASV）数据集。

1.2.3 物种注释 将各数据集的代表序列使用Silva 138数据库［24］的16S序列预先训练的朴素贝叶斯分类器［25］进行物种注释，而后根据注释结果去除线粒体、叶绿体以及非细菌序列。使用MAFFT［26］算法对数据集的代表序列进行对齐，使用FASTTREE［27］软件构建了各数据集的系统发育树。以样地序列数最低样本为基准，分别对各生境下的所有数据集进行抽样，而后使用q2-diversity插件进行了细菌群落的α和β多样性分析。

1.2.4 统计分析 所有统计分析都是在R（v.4.0.3）中使用vegan［28］，stats以及multcomp［29］软件包完成，可视化则是通过ggplot2［30］软件包完成。所有统计检验的显著性均在α = 0.05的水平上确定。

为了评估不同分类分辨率是否会显著影响群落的α多样性，使用方差分析（ANOVA）以及Tukey多重检验鉴定了不同分类分辨率数据集的系统发育多样性（Faith’s PD指数）和物种多样性指数（Chao1及Shannon）的差异。基于群落Bray-curtis距离进行的主坐标分析（PCoA）以及多元方差分析（PERMANOVA）评估样本间群落的差异性和相似性。此外，还通过冗余分析（RDA）量化了环境因子对不同分类分辨率下群落结构变化的解释程度［9］。

2 结果

2.1 不同分类分辨率对细菌群落α多样性的影响

微生物群落的α多样性体现了群落内物种多样性以及丰度。我们选择最为常用的Chao1、Shannon以及PD指数用来表征细菌群落α多样性，所有细菌群落α多样性均是在最小分类单元（OTU或者ASV）上进行（表3）。研究发现，所有生境下OTU数据集的Shannon和Chao1指数均受最小分类单元划分方式的显著影响（P ＜ 0.05），随着分类分辨率［9］（分类单元的精细程度）的上升而上升，表现为100 OTU ＞ 99 OTU ＞ 98 OTU ＞ 97 OTU。而PD指 数 受OTU数据集分类分辨率的影响较小，仅在湿地土壤中有显著差异，表现为98 OTU（135.4 ± 20.4）最高，其次是97 OTU（131.8 ± 20）和99 OTU（131.7 ± 21），而100 OTU（112.2 ± 17.6）最低。

表3 细菌群落的α多样性Table 3 Alpha diversity of bacterial community

值得注意的是，所有生境下基于ASV分类水平所得的Chao1指数均显著低于基于OTU分类方法得到的Chao1指数（P ＜ 0.05）。而对于Shannon指数，单从数值上看，所有生境下ASV数据集的Shannon指数大小均介于97 OTU和100 OTU之间。ASV数据集的PD指数和Chao1指数类似，在数值上均低于OTU数据集，且在森林和湿地土壤中差异显著（P ＜ 0.05）。总之，提高OTU数据集的分类分辨率，可以获得更高的Shannon和Chao1指数，而ASV划分方法会在一定程度上降低Chao1和PD（系统发育多样性）指数。

2.2 不同分类分辨率对群落组成的影响

为了研究不同分类单元划分方式对细菌群落组成的影响，同时在门和属水平对物种的相对丰度进行了单因素方差分析（ANOVA）。结果表明，在门水平（图1-A），不同分类分辨率对湿地土壤、农田土壤、湖泊沉积物和湖泊水体样本的物种组成（门水平）无显著影响（P ＞ 0.05），仅对森林土壤中Armatimonadota和Bdellovibrionota的相对丰度有显著影响（P ＜ 0.05）。其中Bdellovibrionota的相对丰度随着分类分辨率（OTU数据集）的上升而下降，同时基于ASV分析方法注释得到的丰度最低。随着分类分辨率的上升，Armatimonadota并未表现出规律性变化，其在99 OTU中相对丰度最高，而基于ASV分析方法并未注释出相应的门。

进一步并将有显著差异（P ＜ 0.05）的门/属的丰度（不同划分方式下该门/属相对丰度的均值）进行统计后发现，相对于门水平，更多的注释属受到分类单元划分方法的影响（表4和图1）。具体而言，在不同的分类单元划分方式下，森林、湿地、农田、湖泊沉积物和水体样品中分别有75（2.9%）、30（0.35%）、38（3.21%）、15（14.9%）和18（1.32%）个属的相对丰度存在显著差异（P ＜ 0.05）。进一步对丰度较高（top10）差异属的分析发现，湖泊沉积物的Vibrionimonas（11.4%）是相对丰度最高的差异属，其相对丰度随分类分辨率上升而下降，而ASV方法注释的丰度与100 OTU近似。除Vibrionimonas外，其他差异属的相对丰度均较低（小于0.2%），它们在不同样地中表现出不一样的规律，例如，农田土壤和湖泊水体OTU数据集中的WWE3并不会随分类分辨率的变化而发生改变，而ASV方法注释得到的相对丰度显著高于OTU数据集；森林和湿地土壤中g_0319_6G20的相对丰度随分类分辨率上升而下降，而ASV注释的丰度最低；Bdellovibrio在湿地土壤中的丰度会随分类分辨率的上升而下降，ASV则近似于100 OTU，但在森林土壤中表现为98 OTU ＞ 97 OUT ＞ 99 OTU ＞ 100 OTU ＞ ASV。总而言之，分类方法的改变对于门水平的物种丰度影响较小，更多的影响体现在属水平，大多数受影响的属均不属于优势属（top100），且在不同生境下受影响的模式可能并不一致。

表4 最小分类单元划分方法对细菌群落门和属水平物种丰度的影响Table 4 Effects of the minimum taxonomy unit division method on the abundances of bacterial community phylum and genus

图1 最小分类单元划分方法对细菌群落门（A）和属（B）相对丰度的影响Fig.1 Effects of the minimum taxonomy unit division method on the relative abundance of bacteria community phylum（A）and genus（B）

2.3 不同分类方法对微生物群落β多样性的影响

基于不同分类分辨率间可比较的最低分类组成（属水平）进行聚类分析（图2-A）和主坐标分析（PCoA，图2-B）。聚类分析表明，在同一生境类型下，同一样地不同分类分辨率所得的细菌群落结构更为接近（图2-A）。其中98 OTU和99 OTU聚在一起，ASV和100 OTU注释的分类群更为接近。类似的，主坐标分析也发现，同一生境下不同样地（S1、S2、S3和S4）的细菌群落沿第一轴分开（图2-B）。值得注意的是，最小分类单元的划分对细菌群落也造成了一定的影响，使得同一细菌群落在不同的分类分辨率下沿第二轴依次排列。为进一步分析最小分类单元的划分方式对细菌群落β多样性的影响，我们计算了生境下各样本间的Bray- Curtis距离指数，并通过ANOVA（单因素方差分析）和Tukey（事后多重检验）检验了各分辨率间β多样性的差异（图3-C）。分析发现所有生境下OTU数据集的β多样性均随着分类分辨率的上升而上升，表现为100 OTU ＞ 99 OTU ＞ 98 OTU ＞ 97 OTU。值得注意的是，ASV分类方法得到了较高的β多样性（图3-C），与100 OTU无显著差异（P ＞ 0.05）。总而言之，提高分类分辨率会提升群落的β多样性，但对整体群落间的差异影响较小。

图2 基于细菌群落（属水平）bray-curtis距离的β多样性分析Fig.2 Analysis of β diversity based on the bacterial community（genus level）bray-curtis dissimility

2.4 不同分辨率下的微生物群落与环境因子的 联系

通过OTU/ASV水平的微生物群落组成与6个环境因子（pH、TOC、TN、TP、NH4+-N和NO3--N）的RDA（冗余分析）分析发现，随着分类分辨率的上升，环境因子对OTU数据集的解释度逐渐增加，表现为100 OTU ＞ 99 OTU ＞ 98 OTU ＞ 97 OTU。此外，基于ASV分析方法得到的数据结果也与环境因子有较高的匹配度，在所有生境下均高于99 OTU（图3-A）。进一步对环境因子的贡献度进行排序发现（图3-B），在所有生境下97、98和99 OTU各因子的贡献排序一致，不同于100 OTU和ASV的分析结果。比如，在湿地土壤中，97、98和99 OTU数据集中，影响细菌群落结构最主要的环境因子是TN，而100 OTU和ASV的分析结果是pH。对于农田土壤，97、98和99 OTU的分析结果认为NO3--N对细菌群落的影响要大于TP和TOC，而100 OTU和ASV的分析结果则相反。类似的，对于湖泊沉积物和水体样品，NH4+-N对细菌群落的重要性也受到不同分析方法的影响。然而，在森林土壤中，分类方法的变化不会影响环境因子的相对次序。总而言之，随着分类分辨率的上升，环境因子对微生物群落组成变化的解释力也逐渐提高。同时，受最小分类单元划分方式的影响，在探究细菌群落驱动因素时，不同环境因子的相对重要性可能会发生改变。

图3 环境因子对不同分类分辨率数据集的群落组成变异的解释程度 Fig.3 Variation of community composition explained by environment factors in different taxonomy resolution dataset

3 讨论

α多样性常用于表征群落中的个体总数、物种数目、物种均匀性以及系统发育的进化广度等生物学信息［31］。基于物种数目和物种丰度获得Chao1和Shannon指数，是目前最常用的微生物群落α多样性指数。本研究发现随着分类分辨率提升，Shannon和Chao1指数会显著上升，但对PD多样性指数的影响较小（表3）。在Chao1和Shannon指数的计算 中，物种的数目有较大权重，而OTU分类分辨率的上升会增加最小分类单元的数目，因而较高的分类分辨率往往有更高的Shannon和Chao1指数 值［10，32］。此外，有研究认为，提高分类分辨率有利于我们发现一些难以检测到的低丰度物种，然而较高的分类分辨率（如100 OTU）可能无法排除PCR和测序错误，进而高估微生物多样性［33-34］。与Shannon和Chao1指数不同，PD多样性指数基于序列的相似性［35］，在本研究中，OTU数据集分辨率提升产生的OTU代表序列与更低分辨率的OTU代表序列的差异较小（小于3%），分类分辨率的变化对系统发育多样性影响并不显著。ASV数据集的Chao1和PD指数显著低于OTU，但Shannon指数与97 OTU没有显著差异，这表明Shannon指数在不同的分类方法之间更趋于一致。与OTU聚类不同，ASV的划分并不是不按照序列相似性进行聚类，而是通过错误概率模型对可能存在的测序错误进行纠正，这在一定程度上减少了假阳性序列［18］。此外，两种方法使用不同的嵌合体过滤方法，这也可能是导致ASV与OTU数据集多样性指数差异的原因［36］。ASV数据集的Chao1以及PD指数显著低于OTU，但Shannon指数介于97 OTU以及100 OTU之间，表明从OTU方法转到ASV方法，稀有物种所受的影响更大［37-38］。此外，我们还发现在多样性较低的（LS、LW）样地ASV数据集与97 OTU、98 OTU和99 OTU数据集中Shannon和PD指数并没有显著差异，这表明在多样性较低的样地中ASV方法和OTU方法都能较好地估计群落的丰富度［18］。总而言之，对下机数据（碱基序列）最小分类单元的划分方式会影响细菌群落α多样性的最终分析结果。

群落的物种组成是环境微生物群落研究的另一个重要方面，不同组成的微生物群落往往行使着不同的生态功能［39］。目前，环境中微生物群落组成的研究主要是基于门（较高的分类学水平）和属（较低的分类学水平）。我们的研究发现不同的分类方法对较高的分类学水平（门水平）的物种组成影响较小，更多的差异体现在较低的分类学水平（属水平），并且在不同的生境中的差异模式并不一致。这表明分类方法的变化会影响基于属水平的分析结果，尤其是一些低丰度的属（＜ 0.2%）。而最近的研究发现稀有种属在生态功能发挥中起着重要的作用，越来越多的人关注于稀有物种在群落中的作用［40］。例如，在本研究的森林土壤的差异属中，Paenibacillus具有促进植物生长、固氮以及磷酸盐溶解等功能，是一种具有广泛运用前景的农作物促生菌［41-42］，而 Conexibacter可以将硝酸盐还原为亚硝酸盐，在生态修复等方面具有重要的作用［43-44］。湿地土壤与湖泊水体中的差异属Methylorosula是广泛分布于湿地中的一种甲基营养菌，是湿地甲烷排放的重要物种［45］。因此，未来对于低丰度物种的分析，要更加注重最小分类单元的划分方式带来的差异。此外，无论是不同聚类阈值的OTU还是ASV都是希望从种水平重建微生物的分类，当我们在较高的分类学水平分析样地间细菌群落组成差异时，分类方法并不会造成太大的影响。

β多样性常用来量化群落间物种组成差异，对于理解物种多样性的时空变化和群落组成的机制至关重要［46］。在本研究中，提高分类分辨率水平更有利于体现不同细菌群落间的差异。Li等［33］认为，分类分辨率的提升（从97%转到99%）会使得群落的组成进一步细化，从而增大了群落间的差异，这在一定程度上可以帮助我们更好的区分不同的微生物群落。然而，值得注意的是不同样地间的差异大小对于群落β多样性对分类分辨率变化的响应有一定的影响。例如，在β多样性变化幅度较大的湖泊沉积物生境下，分类方法的变化并不会显著影响β多样性，但在群落组成更为接近的农田土壤中，分类分辨率会显著影响其β多样性。这意味着在一些变量梯度均匀的实验（例如精准的模拟控制实验）中，更易于受到最小分类单元划分方式的影响。

微生物对环境条件的变化十分敏感，能够更加迅速的对环境变化做出响应，其群落动态可以在一定程度上表征生态系统功能的变化［20，47］。RDA（冗余分析）常被用来分析环境对微生物群落驱动影响。然而，我们的研究发现最小分类单元的划分方式会影响环境因子对细菌群落变化的解释度，表现为更高分类分辨率具有更好解释度。进一步分析发现在多样性更高的样地（FS、WS和CS）中，ASV划分方式相比于OTU的划分方式往往与环境因子有较好的拟合效果，表明ASV方法可以捕捉到更为准确的群落信息。基于此，我们推测对多样性高的样本（如土壤）更宜采用ASV的划分方式对细菌群落进行探索。而在多样性较低的样地（LS和LW）OTU方法也能够捕捉群落的复杂性。Li等［33］的研究也表明，分类分辨率的提升使得群落组成进一步细分，能够捕捉到更多的群落变异，从而使得它们能够在排序分析中可以更好地被区分，这与我们的结论相印证。除此之外，本研究还发现不同最小分类单元的划分方式，往往会改变特定环境因子的相对重要性。例如，基于100 OTU和ASV 划分方式的结果表明pH是影响湿地土壤细菌群落组成最重要的环境因子，与An等［48］的研究结论一致。而基于较低分类分辨率 （97%，98%和99%）的分析结果认为造成细菌群落组成差异最主要的因素是土壤TN，表明最小分类单元的划分方式会影响环境因子对群落组成相对贡献的结果。因而，为了能够得出更为稳健的生物学结论，未来可能需要同时对下机数据（基因序列）进行多种最小分类单元的划分，采用更多种（如RDA和Mantel等）关联性分析手段分析环境因子与微生物群落之间的联系。

4 结论

通过对比分析不同最小分类单元划分方式下16S rRNA基因扩增测序结果发现，提高分类分辨率能够提高群落的α多样性（Chao1和Shannon指数）和β多样性，同时影响一些低丰度类群（属水平）的组成。此外，还发现分辨率的提升会提升环境因子对群落组成变化的解释度，因而对于物种信息更为复杂的环境样品（如土壤）建议采用ASV划分方式。