水旱轮作对土壤微生物群落构建过程的影响机制①

刘文静，张建伟，邱崇文，包远远，冯有智*，林先贵

(1 土壤与农业可持续发展国家重点实验室(中国科学院南京土壤研究所)，南京 210008；2 中国科学院大学，北京 100049；3 广东海纳农业研究院，广东惠州 510000)

土壤微生物在不同土壤类型、不同土壤发生发育时期呈现特征性的群落结构差异，从而适应生存环境变化并发挥不同的生态功能[1]。当前研究工作已逐步从探究不同生境下的微生物群落分布格局深入到其内在群落构建机制。群落构建过程决定微生物群落组成和多样性，进一步影响其生态功能。明确群落生态构建过程和微生物群落结构以及功能之间的关系是近年来微生物生态学和群落生态学的研究热点[2]。Vellend[3]在微生物群落构建机制理论研究领域做出了先驱性的贡献，其理论进一步被Stegen 等[4]发展为群落构建的4 个生态过程：同质性选择和异质性选择过程(确定性过程)以及扩散限制和同质性扩散过程(随机性过程)。传统的生态位理论认为确定性因素，例如物种特性、种间关系(如竞争、共存和寄生)以及环境因素(如pH、温度和水分等)控制群落结构的形成。确定性过程在某种特定环境条件下可以直接决定此生物群落结构。相反，随机性(中性)理论认为所有的物种都是等生态位的，群落结构由出生、死亡、扩散以及定殖等随机过程控制，环境过滤作用与群落构建生态过程无直接关联。目前普遍认为随机性和确定性过程共同控制微生物群落构建，都起到不可或缺的作用。此生态过程理论被广泛用于探究时空尺度下的微生物群落演替机制的相关研究[5-9]。

生态过程理论现已初步应用到评价农田生态系统可持续性的研究中。可持续性表征生态系统在较长时期内维持某种功能稳定的能力[10-11]。对于农田生态系统，土壤肥力是评价农业可持续性的重要指标之一，高肥力土壤通常能够维持稳定的粮食产量以及较少的肥料施用量[12-14]。相对于旱作生态系统，水旱轮作的稻田生态系统因其较高的粮食生产能力以及土壤涵养能力被认为是典型的可持续农田生态系统[15-16]。除气候、土壤母质、水文灌溉条件等非生物因素外，土壤生物尤其微生物是维系农业土壤可持续性的关键因素。微生物对土壤结构改善、有机物质降解和养分循环等发挥至关重要的作用[17-18]。经典微生物生态功能理论认为，数量越多的微生物将提供更多功能特性，丰富的生态功能有助于可持续生态系统的形成与稳定[19-20]。稻田生态系统被认为是一类可持续农田生态系统，其微生物多样性被广泛研究，但群落构建机制鲜有报道。本研究旨在通过探究同一土壤母质的旱作生态系统、稻田水旱轮作生态系统以及长期淹水植藕的水生生态系统微生物群落构建过程的差异，以明确稻田可持续性生态系统的微生物学机制。

在极端且剧烈变化的生境中，某个特定环境因子将成为影响微生物群落演替的限制性因素，在此条件下，确定性过程主导微生物群落构建过程；随机性过程则在相对温和且扰动较小的环境中发挥主导作用，微生物群落变化更多来自于群落内部的基因漂变等随机因素[4,6]。农业土壤由于受到不断的扰动(例如翻耕、施肥等)，只有特定的微生物能够耐受[21]，所以相对于水生生态系统，旱地和稻田的微生物群落构建以确定性为主，这是本研究的第一个科学假设。稻田以及水生生态系统微生物由于淹水环境的保护受外部扰动以及环境气候变化的扰动较小，加之水的流动性会促进微生物随机性扩散，所以本研究第二个科学假设是相对于旱作生态系统，水旱轮作稻田微生物群落构建过程中随机性过程的影响将会增大，而由于长期淹水状态，水生生态系统中随机性过程可能占主导性作用。通过比较水旱生态系统的群落构建差异将有助于明确可持续性农田生态系统的形成和维系机制，以便更好地利用和调控微生物以应对未来环境干扰并服务于农业生产。

1 材料与方法

1.1 样品采集

试验样品于2018 年8 月在广东省惠州市惠城区水口街道下源村采集。其中，水生生态系统用于莲藕种植(长期淹水)；水旱轮作生态系统开展双季稻种植，冬闲时为旱地(长期水旱轮作)；旱作生态系统用于果树种植(长期旱作)。每个生态系统采集8 个样品，共计24 个样品用于研究。样品带回实验室后，于阴凉通风处风干以除去过多的水分。当样品表面未见明显水分后，用冻存管取样分装并保存于 -20℃ 冰箱中用于提取土壤细菌基因组DNA，剩余样品完全风干并过筛用于测定土壤理化性质。

1.2 土壤理化性质测定

土壤有机碳(TOC)采用重铬酸钾氧化法测定；土壤全氮(TN)采用凯氏定氮法测定；土壤全磷(TP)和全钾(TK)采用氢氟酸-次氯酸消煮，分别由钼蓝比色法和火焰光度计法测定；土壤有效氮(AN)采用2 mol/L KCl 提取，流动分析仪测定；土壤有效磷(AP)采用0.5 mol/L NaHCO3提取，钼蓝比色法测定；土壤有效钾(AK)采用10% HNO3+1% HF 消煮，电感耦合等离子体发射光谱法测定[22]；土壤pH 采用土水质量比1∶2.5 浸提pH 计(Thermo，USA)测定[23]。

1.3 土壤DNA 提取

土壤样品DNA 依照土壤基因组DNA 提取试剂盒(MP Biochemicals, Solon, OH, USA)使用步骤提取。提取的基因组DNA 溶解于50 μl TE(tris-EDTA)缓冲液中，并用1% 琼脂糖凝胶检测是否降解。之后，采用核酸定量仪 NanoDrop ND-1000(Thermo Secientific, Wilmington, DE, USA)对DNA 进行定量检测。

1.4 16S rRNA 基因PCR 扩增及高通量测序

本研究使用引物519F(3′-CAGCMGCCGCGGT AATWC-5′)和907R(3′-CCGTCAATTCMTTTRAGTT T-5′)对细菌16S rRNA 基因V4 ～ V5 可变区进行扩增。正向引物3′端加入5 bp 已知碱基序列的分类标签(Barcode)以区分不同样品。50 μl PCR 体系包括：1.25 μmol/L 的dNTP，2 μl(15 μmol/L)正向和反向引物，2 μmol/L TaqDNA 聚合酶(TaKaRa, Japan)和1 μl(50 ng)的DNA 模板。所有样品扩增3 份后混合检测以消除PCR 过程中的偏好性和随机性误差，阴性对照以无菌水代替DNA 模板。PCR 扩增条件如下：94℃ 5 min，30 个循环(94℃ 30 s，55℃ 30 s，72℃45 s)，最后72℃延长10 min。所有PCR 产物均通过1% 琼脂糖凝胶进行条带特异性检测。使用QIAquick PCR Purification Kit(Qiagen)纯化试剂盒纯化PCR 扩增产物，并用NanoDrop ND-2000(Thermo Scientific，USA)测定产物浓度。将不同样品的PCR 产物按照等摩尔原则进行混合，使用IlluminaMiSeq 测序平台进行高通量测序。

16S rRNA 测序基因序列使用 Quantitative Insights Into Microbial Ecology(QIIME)平台进行分析。质量低于20、长度小于200 bp 的序列被去除，并根据Barcode 序列与样品进行匹配。使用cutadapt和trimmomatic 分别去除扩增引物和接头序列，使用Usearch 算法按照97% 的相似性对序列进行聚类并进行嵌合体去除，将样品划分为不同的OTU，其中每个 OTU 中丰度最高的序列被选作代表序列。物种注释选择SILVA 119 database(http://www.arb-silva.de/download/archive/qiime/)数据库。使用FastTree 建立系统发育树用以计算系统发育水平多样性。

测序共得到601 514 条高质量16S rRNA 基因序列，样品序列在8 742 ～ 46 254 条。为了进行样品间多样性比较，将所有样品的序列统一抽平到8 742 条计算多样性和物种组成。

1.5 生态过程计算

由于系统发育保守性，本研究通过平均最近物种距 离 (between-community mean nearest taxon distance ， βMNTD) 和 最 近 物 种 指 数(betweencommunity nearest taxon index，βNTI)计算样品间的系统发育组成变化(系统发育多样性)。βNTI 表示βMNTD 的观测值和基于零模型(系统发育树结构不变，物种位置随机置换)计算的随机值之间的差异。该分析使用R 软件“picante”工具包的“comdistnt”功能来完成。βNTI 在 -2 ～ 2 表示随机构建过程，而 ＜-2或 ＞2 表示确定性选择过程。

此外，本研究利用表征Bray-Curtis 距离观测值和基于零模型的随机值之间差异的Bray-Curtis-based Raup-Crick(RCbray) 来计算基于扩散影响的微生物群落构建过程[4,24]。试验观测到的Bray-Curtis 值和零模型的偏差被标准化为 -1 ～ 1。RCbray值＜-0.95 或 ＞0.95表示观测值与零模型预期存在显著差异。

确定性构建过程由同质性选择(βNTI＜-2)和异质性选择(βNTI＞2)组成。随机性过程包括扩散限制(|βNTI|＜2，RCbray＞0.95)和 同 质 性 扩 散(|βNTI|＜2，RCbray＜-0.95)。如果|βNTI|＜2，|RCbray|＜0.95，则群落构建过程由未知因素控制。

1.6 数据分析

基于物种水平的Observed OTUs 以及基于系统发育水平的phylogenetic diversity (PD 多样性指数)[25]用于表征群落alpha 多样性。基于Bray-Curtis 距离(OTUs 水平)和βMNTD 距离(系统发育水平)的非度量多维尺度分析(NMDS)用以展示群落结构差异。置换多元方差分析(permutational multivariate analysis of variance，PERMANOVA)用以检验群落结构差异的显著性。以上分析使用R 软件“vegan”包。方差分析(ANOVA)用以多样性指数的组间差异分析，使用SPSS16.0 计算。冗余分析 (redundancy analysis，RDA)用以进行群落结构与环境因子的相关性分析，揭示驱动群落结构变化的主要环境因子。本文所有统计分析，显著和极显著水平分别为P＜0.05 和P＜0.01。使用MENA 在线平台(http://ieg4.rccc.ou.edu/mena)进行共发生网络分析。

2 结果

2.1 细菌群落结构

如图1 所示，水生、旱作和水旱轮作生态系统样品的细菌群落结构均呈现显著分异，且水旱轮作系统样品细菌群落处于水生和旱作生态系统样品之间。PERMANOVA 分析结果表明，3 种样地两两之间细菌群落结构差异均极显著(表1)。

2.2 细菌群落物种组成

根据OTU 物种注释的结果，本研究在门水平上挑选出相对丰度最高的优势细菌用于展示不用样点间的细菌物种组成差异。如图2 所示，在门分类水平上，3 种生态系统中占比最高的物种分别为：变形菌门(Proteobacteria)、酸杆菌门(Acidobacteria)、绿弯菌门(Chloroflexi)、厚壁菌门(Firmicutes)、拟杆菌门(Bacteroidetes)、放线菌门(Actinobacteria)、硝化螺旋菌门(Nitrospirae)。其中，变形菌门相对丰度在水生生态系统中所占比例最大(39.75%)，水旱轮作(29.5%)、旱作生态系统中较低(27.13%)。同样，拟杆菌门的相对丰度呈现相同的规律，在水生、水旱轮作、旱作生态系统中的占比分别为 6.35%，5.35%和4.65%。而酸杆菌门相对丰度呈现相反的规律，在水生(15.05%)、水旱轮作(29.5%)、旱作(35.6%)生态系统呈增加趋势。

图1 基于Bray-Curtis 和βMNTD 距离的NMDS 分析Fig. 1 Nonmetric multidimensional scaling (NMDS) analysis based on Bray-Curtis and βMNTD distances

表1 不同生态系统细菌群落结构差异Table 1 Dissimilarity tests of bacterial communities among different ecosystems using permutational multivariate analysis

图2 不同生态系统下土壤细菌门水平的物种组成Fig. 2 Relative abundances of bacterial community composition at phylum level under different ecosystems

2.3 细菌alpha 多样性指数

本文分别计算了不同生态系统的 Observed OTUs 值和PD 值，分别用来表征物种水平和系统发育水平的alpha 多样性。结果表明，两种水平下细菌alpha 多样性均呈现水生＞水旱轮作＞旱作生态系统的趋势(P＜0.05)(图3)。

图3 不同生态系统细菌群落alpha 多样性指数Fig. 3 Alpha diversity indexes under different ecosystems

2.4 土壤化学性质对细菌群落结构的影响

本文将表征细菌群落结构的各指标分别与土壤化学性质进行相关性分析，以探究土壤化学性质对细菌群落结构的影响。土壤化学性质包括TOC、TN、TP、TK、AN、AP、AK 和pH，其中Observed OTUs和PD 值分别与理化因子进行Pearson 相关性分析，Bray-Curtis和βMNTD距离矩阵分别与各理化因子进行Mantel 相关性分析。结果显示，Observed OTUs和 PD 值均与 AN 极显著相关(P＜0.01)。同样，Bray-Curtis 和βMNTD 距离矩阵与土壤肥力(TOC、TN、TP 和AN)以及pH 呈显著正相关(P＜0.05)，其中与AN 的相关性最强(表2)。

表2 表征细菌群落结构的各指标与土壤化学属性间的相关性分析Table 2 Correlations between bacterial community structure indices and soil chemical properties

土壤化学性质与水生、水旱轮作、旱作生态系统细菌群落的RDA 结果如图4 所示。不同生态系统的细菌群落各自聚集在一起。主坐标 1(RDA1)解释了总变异的56.21%，主坐标2(RDA2)解释了总变异的 28.5%，二者总共解释了总变异的84.71%。不同生态系统间细菌群落结构变异与土壤肥力(TOC、AN、TN 和AK)的相关性最强：TOC、AN、TN 与水生和水旱轮作生态系统细菌群落呈现较好的正相关关系，AK 与旱作生态系统群落结构呈现较好的正相关关系。与Mantel 相关性结果一致，AN 对群落结构变异影响最大。土壤化学性质指标在水生和水旱轮作生态系统下均高于旱作生态系统(表3)。

图4 不同生态系统细菌群落与土壤化学性质的RDA 排序图Fig. 4 Redundancy analysis (RDA) of bacterial community and soil chemical properties

表3 不同生态系统土壤化学性质差异Table 3 Variations of soil chemical properties among different ecosystems

2.5 分子生态网络分析

本文分析了水生、水旱轮作、旱作生态系统下共发生网络的各拓扑结构参数，包括节点数(OTU数)(total nodes)、连接数(total links)、平均路径长度(the average path)、平均聚类系数(average clustering coefficient) 、 传递系数(transitivity) 、 连接度(connectedness)以及模块化指数(modularity)，各参数具体数值见表4。其中，模块化指数均大于0.4，表明群落具有模块结构；节点数在水生生态系统中最高，水旱轮作生态系统中次之，旱作生态系统中最低；连接度指数与节点数呈现相同的趋势；而传递系数在旱作生态系统下最高，说明旱作生态系统下网络稳定性最低。这些结果说明，细菌在淹水条件的水生和水旱轮作生态系统中的种间交互性要高于旱作生态系统，这可能是由于厌氧环境下电子受体的多样化导致了细菌彼此间互作性增强和稳定性提高。

表4 不同生态系统细菌群落共生网络拓扑参数Table 4 Topological parameters of bacterial network of each ecosystem

2.6 驱动群落结构变化的群落构建过程

为了探寻驱动不同生态系统下细菌群落结构差异的原因，本研究耦合βNTI 和RCbray数据分析群落构建的生态过程，结果如图5 所示，驱动水生生态系统微生物群落构建的生态过程主要由扩散限制过程(32.1%)控制，同质性选择和同质性扩散过程各占17.9%；水旱轮作下，同质性选择过程占50%，同质性扩散过程占25%，扩散限制过程占14.3%；旱作生态系统的群落构建主要由同质性选择过程(78.6%)驱动，扩散限制过程占10.7%。此外本研究中3 个生态系统均未发现异质性选择过程。从水生到旱作生态系统，同质性选择过程明显增加(从 17.9% 增加到78.6%)，扩散性限制过程则降低(从32.1% 降低到10.7%)。进一步分析群落构建的确定性过程(同质性选择 + 异质性选择)、随机性过程(同质性扩散 + 扩散限制)、同质性过程(同质性选择 + 同质性扩散)和异质性过程(异质性选择 + 扩散限制)(表5)，结果显示旱作和水旱轮作生态系统由确定性过程主导，而水生生态系统由随机性过程主导；确定性过程(从17.9%到78.6%)和同质性过程(从35.7% 到78.6%)由水生至旱作生态系统增加，随机性过程(从50% 到10.7%)和异质性过程(从32.1% 到10.7%)则降低。

图5 不同生态系统下微生物群落构建生态过程Fig. 5 Relative contributions of ecological processes under different ecosystems

表5 不同生态系统下微生物群落构建的确定性、随机性、同质性以及异质性过程Table 5 Deterministic, stochastic, homogenizing and differentiating processes influencing community assembly under different ecosystems

3 讨论

与本研究的第一个科学假设一致，确定性过程主导旱作与水旱轮作稻田生态系统下的微生物群落构建。其中旱作生态系统的确定性过程比例(78.6%)最高，而稻田生态系统中的确定性过程比例(50.0%)降低，而随机性过程比例(39.3%)升高；在长期淹水的水生生态系统下随机性比例持续升高，并占据主导(50.0%)。在旱作及水旱轮作生态系统下，频繁的农业耕作活动，例如翻耕、施肥等管理措施，使得土壤物理和化学因素变化显著，剧烈的外部环境干扰形成一定的生态位格局，环境过滤作用选择某些微生物生存下来，因此群落构建由确定性过程主导，而该确定性过程因为灌溉而降低。水稻土由于近半年的灌溉，淹水厌氧环境限制了有机质分解消耗，使土壤肥力得以积累[26]，增加的可利用微生物能量以及养分将大量微生物从限制性的生存压力下释放出来，因此随机性过程增加[8,21]。有研究表明，随机性过程在环境压力较小的生境下起主导作用，较低的有机物质分解率、水体连通性等特性导致其内部物质周转缓和，进而缓解了微生物对外部环境扰动的影响[4]。该推论在长期淹水的水生生态系统中也得到了印证：其微生物群落构建由随机性过程主导(表5)。因此，本研究的第二个科学问题也得到了证实。具体到4 个生态过程，水生和水旱轮作生态系统由于淹水耕作而形成独特的氧化还原电位梯度，使得微生物只能局限在特定的生态位中定殖，而旱作生态系统不存在该氧化还原电位梯度，所以扩散性限制在旱作生态系统中比例最小，水旱轮作生态系统次之，水生生态系统中最高。旱作生态系统下，频繁的翻耕活动以及特定施肥方式导致了特定微生物才能在此环境过滤作用下存活，因此具有最高的同质性选择作用，而受外界环境扰动相对较小的水生生态系统下，则同质性选择比例最低(图5)。同理，异质性过程由旱作到水生生态系统呈现了增加趋势(表5)。水生、水旱轮作生态系统中较高的异质性过程有助于微生物群落形成更多生态功能。

无论在水生、水旱轮作生态系统还是在旱作生态系统，随机性过程和确定性过程都共同驱动群落结构的形成(图5 和表5)。确定性过程可以强化特定生境下的微生物功能，而随机性过程能够产生更多样的生态功能以维持生态系统功能的稳定性和持续性[2]。同时，随机性微生物群落构建过程可以对于剧烈环境变化所带来扰动起到一定的缓冲作用。因此，随机性过程对于可持续生态系统功能有正向反馈机制[27-30]。水生及水旱轮作生态系统下增加的随机性过程可提供更丰富的物种多样性以及生态系统功能。这一结论同样从微生物群落多样性的结果得以证实：无论在物种水平还是系统发育水平，细菌群落alpha 多样性均呈现由旱作到水生生态系统显著增加(图3)。可见，有水分参与形成的生态系统利于形成物种更加丰富的微生物群落结构。同理，从网络分析的拓扑参数中也能明显发现，水生和水旱轮作生态系统下的平均路径长度、平均聚类系数以及连接度指数参数均显著高于旱作生态系统 (表4)。该结果说明在厌氧环境由于氧化还原电位的梯度，产生了更多的生态位[31]，可以容纳多种微生物共存[32]。同时，水生和水旱轮作生态系统下较低的传递系数表明了更稳定的群落物种间交互关系[33]。微生物间较高和较稳定的交互性有利于可持续的生态系统构建。

土壤微生物是土壤养分的储存库以及植物生长可利用养分的重要驱动者。土壤微生物群落是反映土壤健康与土壤生产力的一个重要指标，并快速响应于土壤理化性质的演变[34-35]。在较高的随机性过程影响下，微生物多样性高、群落互作网络紧密、生态功能丰富等均有利于养分物质的循环和积累，因此水生和水旱轮作生态系统下的土壤肥力状况显著优于旱作生态系统(表3)。因此，与前人研究结论相吻合[36]，土壤肥力(TOC、TN、TP、AN)与群落结构呈现较强的相关性(表2 和图4)。但是本研究从群落构建的层面更好地揭示了微生物群落和土壤肥力的耦合性。同时，也应该看到，可能正是这些群落构建机制，稻田生态系统成为了可持续性生态系统，而水生生态系统更是对全球物种贮存、气候调节、保护生态环境具有不可替代的重要作用[37]。下一步工作，应探究驱动群落构建过程的环境因子，明确可持续农田生态系统微生物群落结构、功能和构建机制之间的相互关系，为可持续农田生态系统的构建提供理论支撑。

4 结论

本研究利用生态零模型计算3 种生态系统下的土壤微生物群落构建机制，发现水旱轮作形成的缓和内部环境及土壤肥力积累，使稻田相比于旱作生态系统，随机性过程增加，确定性过程下降。随机性过程维持微生物物种多样性以及物种生态网络的交互度和稳定性，增强了微生物抵抗外部环境扰动能力，该过程有助于维持稻田生态系统功能的稳定性和可持续性。