三峡消落带土壤古菌群落组成与代谢功能预测

叶 飞，吴晓冰，姜 毅，吴胜军，王 雨*

（1.中国科学院 重庆绿色智能技术研究院，重庆 400714；2.中国科学院水库水环境重点实验室，重庆 400714）

长期以来，古菌被认为仅存在于强碱、强酸或者高盐等极端环境中。当与极端环境古菌高度相似的“泉古菌”在海洋浮游生物中被发现后[1-2]，一系列后续研究无不证明古菌是地球生物圈中一类丰富且多样的微生物，广泛地分布于地球的各类环境中，改变了人们对古菌的传统认识[3-7]。数量庞大的古菌在生物地球化学过程中的重要作用远超预期。

在碳循环中，产甲烷古菌产生的甲烷占全球甲烷排放量的70%[8]。除此之外，古菌还积极参与硫循环，能将二硫化碳（CS2）快速转化为硫化氢（H2S）和二氧化碳（CO2）[9]。进行化能自养的氨氧化古菌能将氨（NH3）氧化为亚硝酸盐（NO2-）[10]。氨氧化古菌的发现彻底改变了氨氧化过程由细菌驱动的传统认识，它们在自然界中的庞大数量意味着其在自然界物质转化尤其是全球碳氮循环方面扮演着重要角色[4，7]。

三峡水库实行夏落冬涨的调度运行方案，夏季汛期（6—9月）保持145 m 最低水位，汛期末（10月）开始蓄水，冬季（11月）达到175 m最高水位，之后水位逐渐消落，次年5月底降至最低水位，水位变幅大且同自然洪枯规律相反。季节性水位涨落使库区淹没土地周期性露出水面，形成了面积达348.9 km2的消落带。三峡库区消落带由于其所处的特殊空间位置，同时受到人类活动和水位波动的干扰[11-12]，导致消落带土壤环境比较复杂，进而影响消落带微生物群落的分布。然而目前对消落带古菌群落组成的空间分布规律还缺乏深刻的认识。此外，了解古菌群落的潜在功能对于更好地认识其在生物地球化学循环中的角色是至关重要的，有助于理解长江流域氮磷营养盐循环过程，对流域水资源保护有重要意义。近年来越来越多的研究采用宏基因组测序的方法来揭示环境中微生物群落功能[13-14]。作为一种比较经济的生物信息学方法，PICRUSt（phylogenetic investigation of communities by reconstruction of unobserved states）可以通过16S rRNA基因数据推断微生物群落的功能属性，具有较好的预测结果，已广泛应用于人体或哺乳动物[15]、土壤[16]、海水和沉积物[17]等生境研究，是了解古菌群落代谢功能的有效途径。本研究以不同高程消落带土壤为研究对象，分析其古菌群落组成沿高程梯度的分布规律，并采用基于16S rRNA 高通量测序的PICRUSt 功能预测分析消落带土壤古菌群落的代谢功能特征。

1 材料与方法

1.1 样品的采集与保存

本研究选取的白家溪消落带位于长江一级支流澎溪河流域，处于三峡库区腹地的重庆市开州区境内（图1）。该区域属于中亚热带湿润季风气候，年平均气温18.5 ℃，四季分明，年均降雨量大约1 200 mm，雨季长、雨日多、降水丰富。土壤类型以发育自紫色砂岩的紫色土为主。

图1 研究区和采样点示意图Fig.1 Map of the study region and sampling sites

以海拔150 m 为最低点，沿消落带150～175 m海拔高程选择一个垂直样带，每隔5 m 设置采样点，共6个采样点。在消落带水位上升初期的2013年10月10日和2014年9月23日，以及水位下降末期的2014年4月1日和2015年5月6日共进行两年4次土壤样品的采集。水位线以上的高程样点，用不锈钢土钻采集；水位线以下用彼得森采样器采集。每个样点采集3 个重复样，充分混匀成一个混合样，共形成24 个混合样。样品采集后放置于无菌密封袋，置于冰盒中保存运输至实验室。其中一份立即进行理化指标分析，另一份于-20 ℃储存，用于DNA提取和后续分子生物学实验。

1.2 土壤理化性质分析

土壤pH采用pH计（FE20，Mettler Toledo，美国）测定，水土比为1∶5。含水率（Moisture）用烘干质量测定法。土壤NH4+和NO3-的含量用2 mol/L KCl 浸提后，采用连续流动注射分析仪（FIA Star 5000，FOSS Tecator，瑞典）测定。土壤有机质（origanic matter，OM）采用烧失量法测定。总碳（total carbon，TC）、总氮（total nitrogen，TN）和总硫（total sulfur，TS） 用元素分析仪 （Vario EL cube，Elementar，德国）测定[18]。

1.3 DNA的提取和Illumina高通量测序

采用Power Soil®DNA Isolation Kit（Mobio，美国）试剂盒提取土壤样品中的总DNA。用超微量紫外分光光度计NanoVue Plus Spectrophotometer（GE Healthcare，英国）确定DNA浓度。用浓度为1%（w/v）的琼脂糖凝胶电泳检验DNA的质量。提取的DNA溶液保存于-20 ℃冰箱中待用。

针对古菌16S rRNA 基因的V3-V5 高变区，使用带标签序列（barcode）的引物Arch344F（5′-ACGGGGYGCAGCAGGCGCGA-3′） 和 Arch915R（5′-GTGCTCCCCCGCCAATTCCT-3′）[19]进行聚合酶链反应（polymerase chain reaction，PCR 扩增。扩增采用TransGen AP221-02（TransGen，中国）的 20 μL 体系，具体为：5 μL 5 × FastPfu buffer，2 μL dNTPs（脱氧核糖核酸，2.5 mmol/L），0.3 μL FastPfu DNA 聚合酶（2.5 U），正、反向引物（5 μmol/L）各0.8 μL，0.2 μL BSA（牛血清白蛋白，20 μg/μL）（Takara，中国），DNA模板（10 ng），双蒸水补足至20 μL。PCR 扩增程序为：95 ℃预变性3 min，然后95 ℃变性30 s，55 ℃退火30 s，72 ℃延伸45 s，共32 个循环，最后再72 ℃延伸10 min。扩增所用的PCR 仪为ABI GeneAmp®9700。PCR 产物使用AxyPrepDNA 凝胶回收试剂盒（AXYGEN，美国）纯化后用2%琼脂糖凝胶电泳进行检测。参照电泳初步定量结果，将PCR 产物用QuantiFluor™-ST 蓝色荧光定量系统（Promega，美国）进行检测定量，之后按照每个样本的测序量要求，进行相应比例的混合后双端（pair-end，PE）测序（2 × 300）。本研究的Illumina 高通量测序委托上海美吉生物医药科技有限公司①http://www.Majorbio.com在Illumina Miseq PE300 平台上完成。在一个淹水周期内，160 m 样点的淹水和出露时间与165 m 样点非常接近，因此，在进行高通量测序及后续相关分析过程中将160 m高程采集的样品排除在外。

1.4 高通量测序序列处理

MiSeq测序得到的是双端序列数据，首先根据双末端片段（PE reads）之间的重叠（overlap）关系，将成对的片段（reads）拼接（merge）成一条序列，同时对reads 的质量和拼接的效果进行质控过滤，根据序列首尾两端的barcode 和引物序列区分样品得到有效序列，并校正序列方向，进行数据优化。优化后的高质量序列使用Usearch 软件（vsesion 7.0，http://drive5.com/uparse/）在97%的相似度下划分独立操作单元（Operational Taxonomic Units，OTUs）。为了得到每个OTU对应的物种分类信息，采用核糖体数据库项目分类器（RDP classifier）贝叶斯算法对OTU 代表序列进行分类学分析，分析比对使用SILVA rRNA 数据库（Release 128，http://www.arb-silva.de）。为了便于比较各样品间的差异性，根据测序样品中的最小序列数将所有样品的序列稀释到同一测序深度（每个样品21 068条序列）。

1.5 PICRUSt功能预测与统计分析

基于古菌群落16S rRNA 基因代谢功能的预测通过PICRUSt v1.0.0 软件实现[16]。参照Greengenes数据库（13_5_release）将用于PICRUSt 预测的序列划分为OTUs（97%相似性）。在比对过程中剔除与参照序列不符合的序列。根据OTU 对应基因组中16S rRNA 基因拷贝信息，将每个OTU 对应序列数除以其16S拷贝数进行标准化。最后，将标准化的数据乘以其对应的基因组中基因含量，从而实现代谢功能预测，产生KEEG（Kyoto Encyclopedia of Genes and Genomes） 同源 序 列 （KEEG orthologs，KOs）表。预测结果的准确度可以通过参数NSTI（Nearest Sequenced Taxon Index）值的大小来评估。NSTI 表征某样品中所有微生物OTU 与其亲缘关系最近的已测序基因组间系统进化距离的评价值。因此，该值越小表示预测结果越准确。

不同高程土壤理化性质差异利用IBM SPSS Statistics 20.0 for Windows 的单因素方差分析（One-way ANOVA）进行比较。使用R 语言（version 3.0.1）的“vegan”包计算基于布雷-柯蒂斯距离（Bray-Curtis Distance）的不同高程土壤古菌群落的差异性（Analysis of Similarity，ANOSIM；置换检验等于999），并进行古菌群落组成变化与土壤性质关系的Mantel检验。主坐标轴分析（Princi-pal Coordinates Analysis，PCoA）和冗余分析（Redundancy Analysis，RDA） 通过 CANOCO 5 软件实现。图形绘制使用ORIGIN 9.0软件完成。

2 结果

2.1 土壤性质沿高程的分布

消落带土壤NH4+和NO3-浓度、pH 值以及含水率均表现出沿采样高程的梯度分布规律（图2）。其中，NH4+，pH 值和含水率表现出随着采样高程降低而增加的趋势。但150 m高程样点的NH4+浓度和含水率要明显低于155 m 高程（P＜0.05），这种情况很可能是因为150 m 高程样点位于河床附近，采集的土壤样品中含有大量的砂石。此外，消落带土壤中的NO3-浓度总体上随着采样高程的升高而增加。而土壤中OM，TC，TN 和TS 的变化都与采样高程没有关系。

图2 消落带4次采样不同高程土壤性质的均值Fig.2 The mean values of the soil properties at different elevations in four sampling surveys.

2.2 古菌群落组成变化及驱动因子

通过基于OTU 水平的ANOSIM 分析，比较不同高程样点间古菌群落组成的差异（表1），可以发现各高程土壤古菌群落组成与相邻高程没有显著性差异，但是与相隔较远高程的土壤差异显著。如175 m高程样点的古菌群落组成与170 m没有显著差异（P＞0.05），而与165 m，155 m 和150 m 高程的土壤古菌群落都有显著差异（P＜0.05）。170 m高程样点与165 m 高程样点的古菌群落无显著差异 （P＞0.05），而与155 m 和150 m 高程呈现显著差异（P＜0.05）。说明消落带土壤古菌群落组成沿高程梯度呈现出明显的空间分异特征。Mantel检验结果表明，消落带古菌群落组成的变化主要与土壤NH4+，TN和含水率显著正相关（P＜0.05，表2）。

表1 基于ANOSIM分析的不同高程间土壤古菌群落组成比较Tab.1 Comparison of soil archaeal community compositions between different sampling elevations based on ANOSIM analysis

表2 基于Mantel检验的土壤性质与古菌群落组成变化的相关关系Tab.2 The correlations of soil properties with the variation of archaeal community composition based on Mantel test

2.3 代谢功能预测

本研究中每个样品的NSTI值范围为0.126～0.296，平均为0.238。采用PICRUSt 预测的功能基因组成的变化沿高程可以分为2 个主要的组，分别是170～175 m 和150～155 m 样点组。组内各自对应的两个高程样品的功能组成比较类似。然而，165 m高程样点分布比较零散，此高程功能基因组成表现出一定的随机性[图3（a）]。此外，本研究将具有代谢功能的基因组从KOs 中筛选出来，对土壤中代谢功能组成进行PCoA 分析，发现与全功能组成分布类似，具有代谢功能的基因组成同样沿高程可以分为170～175 m和150～155 m样点2个组[图3（b）]。

图3 采用PICRUSt 预测的土壤样品（a）KEGG orthologies（KOs）和（b）在KEGG功能层次结构中被归为“新陈代谢”的KOs的主坐标轴分析（PCoA）Fig.3 Principal coordinate analysis (PCoA) of soils using PICRUSt predicted (a) KEGG orthologies (KOs) and (b) KOs only grouped as “Metabolism” in KEGG functional hierarchies

为了解消落带土壤中代谢功能变化的驱动因子，本研究采用RDA 方法分析代谢功能组成与土壤性质的相关关系。RDA1 和RDA2 两轴分别解释了代谢功能变化的32.6%和3.6%（图4）。土壤样品预测的代谢功能变化与土壤的含水率、OM和TC显著相关（图4）。但是没有任何一个测得的土壤性质与某一高程样点预测的新陈代谢功能的分布模式有明显的相关关系。

图4 预测的代谢功能群落与土壤性质的冗余分析（RDA）Fig.4 Redundancy analysis (RDA) of predicted metabolism functional communites and soil properties

3 讨论

三峡消落带处于湿地生态系统次生演替初期阶段，大量的本地物种消失、多样性降低[11]。微生物群落作为消落带生态系统的关键组成部分，处于适应变化环境的过程中，最终会形成一个更稳定的群落来应对这种干扰[20-21]。因此微生物群落对于指示消落带生态系统功能的恢复[22]以及生物地球化学循环[23]等方面具有重要作用。

古菌群落组成的PCoA 分析表明，消落带土壤中古菌群落具有明显的空间异质性。这可能跟淹水导致不同高程土壤的氧环境存在较大差异有关。绝大多数的广古菌门（如产甲烷菌）与土壤含水量显著正相关[24]。此外，稳定的氧化还原条件更有利于广古菌门微生物的生长[25]。在三峡库区消落带，较低高程土壤的淹水时间较长，形成的相对稳定的厌氧环境更适合广古菌。而奇古菌门主要由好氧氨氧化古菌组成，短期的淹水导致较高高程的消落带土壤形成较好的微氧生境，为好氧氨氧化古菌提供了良好的生境[26]。除此之外，好氧氨氧化古菌的分布还与氨氮浓度相关[27]。Mantel检验同样证实消落带古菌群落的分布与土壤氨氮浓度和含水率显著正相关。

本研究中各样品的NSTI均值为0.238，预测精度要低于已有研究，这可能是消落带土壤环境包含很多未知的多样性[16]。消落带土壤古菌功能基因组成变化表现出与古菌群落组成变化类似的模式，即沿高程的空间分布。175 m 和170 m 高程土壤古菌功能基因组成类似，而155 m和150 m高程的功能基因组成类似，但较高高程和较低高程之间存在较明显的差异。然而，不同高程功能基因组成的整体变化要弱于古菌群落组成的变化，古菌群落这种预测功能结构可能与其功能性冗余有关[28]。总的来说，与古菌群落结构组成的驱动因子一样，含水率同样是影响预测的代谢功能组成的主要环境因子。古菌中的产甲烷古菌和好氧氨氧化古菌是碳氮循环中的重要组成部分[8，10]，此外某些古菌还是硫循环的主要参与者[9]。其中，产甲烷古菌和氨氧化古菌的分布与土壤含水率密切相关，这在一定程度上反映了古菌代谢功能沿高程的空间分异，较高高程（175～170 m 高程组）主要行使氨氧化功能，而较低高程（155～150 m 高程组）的古菌主要参与碳循环。因此，长时间的淹水可能会导致三峡消落带土壤氮循环和碳循环在空间上的分异。