分析方法對細菌群落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劃分方式。