微生物是陆地生态系统中物种丰富的生物群落之一，在各类土壤中发挥着重要的生态系统功能，包括许多分解和地球化学循环过程，以及可以影响地表植被物种多样性和生产力，并且能够保持土壤肥力，是维持陆地生态系统稳定性与可持续性的关键驱动者^[1-5]。生物多样性已被认为是影响生态系统功能的主要因素^[6-7]。然而，我们对土壤微生物群落结构和多样性的认识主要来自于表层土壤样本，因为绝大多数研究只集中在土壤表层大约20 cm处，对整体土壤剖面的关注还比较缺乏。

土壤不是单一的环境，而是由一系列物理化学成分组成的连续可变的生境，为不同的微生物提供了独特的生存空间。许多研究表明土壤微生物作为重要的指示因子，在表层土壤中，微生物群落在很大程度上受到土壤微环境因子如pH、温度、营养水平(N、P、K、S)和水分等变化的影响^{[1-2, 8]}。即使土壤环境中仅相隔几微米到几毫米的距离，但它们的非生物特征、微生物丰度、微生物活动速率和微生物群落组成可能有很大的不同^[8]。然而成熟的土壤剖面有时可深达数百米，并且大量的微生物也会在地下较深的土壤中栖息^[9-10]。最近，有研究表明深度是影响微生物组成和结构的重要因素，并且强调了土壤深度作为构建土壤微生物群落环境梯度的重要性^[11]。这是由于土壤质地以及压实过程沿剖面向下变化时，减少了土壤总孔隙空间和土壤氧含量，从而限制了微生物的生长，并且各种土壤理化性质在剖面上的梯度变化，共同使得土壤成为具有许多环境梯度的高度异质的生境，从而导致微生物组合的垂直分层^[12-13]，也造成微生物多样性、群落组成和功能会随土壤深度的不同而变化。Fierer等^[14]通过比较垂直剖面的微生物多样性发现，随着土壤深度(0−200 cm)的增加，微生物种群数量显著减少。Chu等^[15]对西藏高原土壤细菌群落的深度分布研究揭示了表层和次表层土壤微生物群落结构明显不同，与总碳和碳氮比的变化密切相关。Tao等^[16]研究了氨氧化微生物在黄土高原土壤剖面上的垂直分布，其结果表明氨氧化细菌(ammonium-oxidizing bacteria，AOB)随土壤深度的增加而显著降低，改变了AOB的多样性和功能，可能对土壤氨氧化群落及其功能产生深远影响。然而，人们对土壤剖面上微生物总量和多样性垂直分布特征以及群落结构随深度梯度变化的了解仍然有限。

微生物群落构建的决定性因素通常被认为包含两个主要驱动因素：(1)非生物因素，可以由pH、营养元素、水分等多种因子决定，也可以由空间尺度上的扩散(如地理距离和深度)限制决定；(2)生物因素，即微生物之间的相互作用决定了微生物群落的复杂结构。土壤微生物在复杂的生态网络中共存，物种间存在着不同类型的共生、竞争、捕食、偏利或偏害共生等^[17]。大多数生态功能特征，例如凋落物的降解、温室气体的产生以及封闭的营养循环等都不是单一物种的作用结果，而是微生物物种间密切协作和相互竞争的产物。所以，理解群落中不同物种之间的相互作用以及它们对环境变化的反应并揭示其组装背后的过程，是理解生物多样性、群落稳定性和生态系统功能模式的基础^[18-19]。例如，Tian等^[20]研究了原生林和次生林土壤的细菌物种网络，生态网络分析发现，相较于次生林，原生林中微生物种类间具有更大的生态位共享和更多的相互作用，通过这样的网络交互作用，一个生态系统内的所有成员能够完成单个种群无法完成的生态功能(如养分循环、维持生态系统稳定性)。Jiao等^[21]研究了中国东部不同生境和地区相邻的玉米和水稻地里的古菌、细菌和真菌群落的聚集情况，结果显示群落组装过程的平衡调节了中国东部农业土壤微生物群落的物种共存。揭示群落聚集和物种共存之间的联系是了解生态系统多样性和功能的关键，是微生物生态学中一个基本但较少被研究的课题^[21]。

长期以来，土壤生态调查主要集中在单一的生态系统之中，以调查生物群落多样性和构建机制，如森林^[22]、草地^[23]、农田^[24-25]以及被污染的生态系统^[26-27]等。然而，不同栖息地之间的过渡区域，即交错带，特别是干旱地区，很少受到研究人员的关注。已有研究表明各类交错带(如林际交错带、湖泊交错带、农牧交错带等)都对动物、植物、微生物有一定的影响^[28-31]。农牧交错带是重要的陆地生态系统，因此，迫切需要研究生存于其中的微生物群落，以了解决定这些地区物种组成和群落结构的生态过程，为其合理开发和系统保护提供科学依据。分子生物学研究方法的进步使研究人员能够在一个更加精细的水平内确定土壤微生物多样性，特别是基于微生物遗传物质DNA和RNA的土壤微生物群落分析现在比较普遍，并对土壤系统中特定的微生物建立更为全面的了解^[32-33]。因此，我们使用16S rRNA基因高通量测序来检测农牧交错带土壤原核微生物群落随土壤深度(0−80 cm)的变化，以期阐明微生物是如何通过改变群落水平的组成结构和相互作用来适应土壤的垂直梯度。本研究的主要目的是探讨以下问题：(1)原核生物群落的组成与结构如何随深度变化？(2)原核生物群落相互作用的分子生态网络如何受深度的影响？

1 材料与方法 1.1 土壤样品收集

选取的采样地点在内蒙古自治区多伦县(42°02 N，116°17 E)，属于典型农牧交错区域，该地为暖温带季风气候区域^[34]。该生态系统以多年生草本植物为主，包括小针茅(Stipa klemenzii Roshev.)、冷蒿(Artimesia frigida)、星毛委陵菜(Potentilla acaulis)、砂韭(Allium bidentatum)和糙隐子草(Cleistogenes squarrosa)等^[34]。长期年平均降水量(mean annual precipitation，MAP)约为383 mm，其中90%的降水分布在5−10月的6个月内；月平均气温变化范围为−17.5 ℃ (1月)至18.91 ℃ (7月)，年平均气温保持在2.11 ℃左右^[35]。根据联合国粮农组织的分类，该地区的土壤类型被划分为普通钙质土，是分别由62.7%的沙子、20.3%的泥沙和17.0%的粘土组成的一种土壤；土壤pH为中性(6.84±0.07)，平均容重约为1.31 g/cm^3[36]。为了研究微生物在土壤剖面上的变化，于2017年8月使用长度为1 m的采样器在空间尺度东西跨度达6 000 m的区域收集土壤剖面样品，共设置了8个采样点，并将所得土壤按照10 cm等间距区分为可以代表不同深度的样本。根据取样柱剖面形态特征发现，该地部分样本80 cm以下已呈散状砂石结构，为了保持样本重复一致，重点分析其中0−80 cm共64个样本，并且依次命名为0−10、10−20、20−30、30−40、40−50、50−60、60−70、70−80 cm。土壤被挑除植物根系和筛除沙砾后立即送到实验室，其中部分土壤4 ℃保存用于理化参数测定，另一部分储存在−80 ℃冰箱用于DNA提取。

1.2 土壤理化参数测定

对于每个样本，土壤湿度的测量方法是将15 g土壤放入温度为55 ℃的烤箱中烘干，直到其重量恒定，土壤含水率(%)根据烘箱干燥后的失重占原重的百分比来计算。土壤pH值是在土壤悬液中测量的，悬浮液由土壤:水=1:2.5 (质量体积比)形成，根据酸度计操作标准测试。本研究中所需要的其他理化参数由中国科学院南京土壤研究所进行测量。将风干后的土壤样品与重铬酸钾标准液和浓硫酸混合、消煮，硫酸亚铁滴定测定土壤总有机碳(total organic carbon，TOC)。用过硫酸钾溶液在高压锅中121 ℃下溶解土壤样品40 min，冷却后取出，用紫外分光光度计测定消化液浓度，计算确定总氮(total nitrigen，TN)含量。把土样加入2 mol/L氯化钾溶液振荡1 h后800×g离心15 min，静置10 min后通过0.45 μm微孔滤膜过滤，分别在波长420 nm和210 nm处测定上清液的吸光度，计算得到土壤铵态氮(NH₄⁺-N)和硝态氮(NO₃^--N)的含量。

1.3 土壤微生物全基因组提取、PCR扩增以及测序

按照操作说明，使用专门针对分离土壤有机体基因组DNA的MP Fast DNA Spin Kit for Soil从土壤中提取微生物总DNA。所有样品浓度达到10 ng/μL，且OD₂₆₀/OD₂₈₀的值在1.8−2.0范围内，保证DNA质量合格可以进行16S rRNA基因的PCR扩增实验。

用16S rRNA基因扩增通用引物515F (5′-GTGCCAGCMGCCGCGGTAA-3′)和806R (5′-GGACTACHVGGGTWTCTAAT-3′)扩增提取后合格的DNA。引物序列中包含简并碱基，M代表C或A；H代表A、C、T；V代表A、C、G；W代表A或T。该引物集针对原核生物16S rRNA基因V3−V4高变区为靶点，并具有识别样品的特异性标签。PCR反应体系(50 μL)：Taq DNA聚合酶(5 U/μL) 0.5 μL，10×PCR缓冲液5 μL，核苷酸混合物(2.5 nmol/μL) 1.5 μL，正、反向引物(10 μmol/L)各1.5 μL，模板DNA (20−30 ng/μL) 1 μL，超纯水补至50 μL。PCR反应条件：94 ℃ 1 min；94 ℃ 20 s，57 ℃ 25 s，68 ℃ 45 s，30个循环；68℃ 10 min，4 ℃保存。使用琼脂糖凝胶电泳的方法分离扩增后的DNA片段，对产物进行鉴定。

根据不同的样品标签，将不同样品回收后的DNA以相同的质量进行混合(每个文库约50个样品)，定量到20 ng/μL以上。采用KAPA HyperPrep Kit DNA文库构建试剂盒制备测序文库。之后，根据MiSeq基因测序试剂盒说明对样品进行上机前准备。最后，将准备好的测序文库在Illumina MiSeq平台进行正反向双端测序。

1.4 序列处理

16S rRNA基因的原始序列以FASTQ格式被提交到我们构建在Galaxy平台上一个公开的集成了多种生物信息学分析工具的序列分析平台(http://mem.rcees.ac.cn:8080)进行进一步的分析^[37]。序列Reads根据标签被分配到不同的样品中并允许一个错配发生，然后去掉标签和引物序列。使用FLASH步骤合并相同序列的正向和反向Reads^[38]。使用Btrim对Reads进行过滤，保持平均质量分数大于20，最小长度为140 bp^[39]。任何带有简并碱基的序列都被删除，只保留长度在245−260 bp之间的Reads。之后，使用Greengenes参考数据集作为嵌合体检测的参考^[40]，使用UPARSE^[41]将序列聚类为可操作分类单元(operational taxonomic units，OTU)，序列相似度阈值为97%。单条序列将被保留，因为稀有物种会导致β多样性发生重大变化^[42]。利用核糖体数据库项目16S分类器将16S rRNA基因代表性序列分类注释^[43]，其置信度为50%。基于每个OTU的代表性序列，使用FastTree构建16S rRNA基因的系统发育树^[44-45]。

1.5 荧光定量qPCR

根据Zhang等^[46]对土壤总DNA定量的描述，利用CFX96实时定量PCR系统测定土壤原核微生物16S rRNA基因的丰度。qPCR引物序列为：Forward：5′-CCTACGGGAGGCAGCAG-3′；Reverse：5′-TTACCGCGGCTGCTGGCAC-3′。qPCR反应体系(20 μL)：模板DNA 10 ng，2×SuperReal PreMix Plus 10 μL，正、反向引物(10 μmol/L) 1.5 μL，超纯水补至20 μL。标准曲线由连续10倍稀释包含目标扩增子的标准质粒组成。qPCR反应条件：95 ℃ 15 min；95 ℃ 15 s，55 ℃ 30 s，72 ℃ 32 s，共40个循环。PCR结束时，在72-98 ℃下以0.5 ℃为增量进行熔解曲线分析。16S rRNA基因的扩增效率在90%−110%之间，并且R²高于0.99。所有样品都有技术重复一式三份，并且无菌水阴性对照显示无特异性扩增。

1.6 主要试剂和仪器

土壤基因组DNA提取试剂盒、FastPrep匀浆仪，MP生物医学公司；ExTaq酶，TaKaRa公司；DNA纯化试剂盒，Omega公司；SuperReal PreMix Plus荧光定量试剂盒，天根生化科技(北京)有限公司；KAPA HyperPrep Kit DNA文库构建试剂盒，KAPA Biosystems公司。PCR仪，Bio-Rad公司；NanoDrop 2000超微量分光光度计、Qubit荧光仪，Thermo Fisher公司；MiSeq高通量基因测序系统，Illumina公司。

1.7 分子生态网络

为了评估物种在不同土壤深度的共存情况，基于16S rRNA基因的焦磷酸测序数据，将其应用于确定不同深度水平下的微生物群落网络结构，通过在线平台(http://ieg2.ou.edu/MENA)构建了由原核微生物群落成员组成的分子生态网络(molecular ecological networks，MENs)^[47]。首先，通过SparCC相关性计算成对OTU的相似矩阵。SparCC相关性与物种多样性无关，并对所有参数值进行了改进，计算的相关性更加精确^[48]。其次，基于随机矩阵理论(Random Matrix Theory，RMT)的方法自动计算分子网络。与其他网络构建方法相比，这种方法的显著特点是网络是自动定义的，并且对噪声具有很强的鲁棒性^[47]。最后，进行网络拓扑结构分析，包括网络拓扑特征、模块检测、基于模块的特征基因分析和模块角色识别。为了描述网络的拓扑结构，计算了一系列基本的指标，包括网络节点(node)、边(edge)、平均度(avgK)、平均聚类系数(avgCC)、模块化(modularity)等。根据节点的模块内连接度(Zi)和模块间连接度(Pi)将网络中所有节点的模块角色分为4种类型，包括模块中心点(module hubs)、连接节点(connectors)、网络中心点(network hubs)以及外围节点(peripherals)。通常模块中心点、连接节点和网络中心点被当作分子生态网络中的关键物种。分子生态网络图可视化通过Gephi 0.9.2实现。

1.8 统计分析

采用ANOVA方差分析和Tukey事后比较对各土壤理化参数和原核微生物群落的多样性指数进行了检验，并且在P < 0.05的显著性水平判断是否存在显著差异。使用Bray-Curtis距离矩阵来表示原核生物群落之间的差异，并且Bray-Curtis矩阵被用于非度量多维尺度分析(non-metric multidimensional scaling，NMDS)来可视化群落组成的差异。利用置换多元方差分析(permutational multivariate analysis of variance，PERMANOVA)测试了不同土壤深度间群落组成的差异^[49]。使用Mantel和Partial Mantel检验^[50]以及典范对应分析(canonical correspondence analysis，CCA)来确定地理距离和土壤变量对原核生物群落差异的相对贡献。在进行计算之前，将对所有变量进行标准化，使每个变量的权重相等。所有统计分析均使用R软件(V 3.6)进行^[51]。

2 结果与分析 2.1 土壤理化性质的变化

不同土壤深度的各因子参数比较如表 1所示。随着土壤深度的增加，pH值大小有略微增加(7.26−7.81)，但是其他各参数含量都是逐渐减少。所有参数在不同深度的方差分析结果显示，含水率在整个土壤剖面没有显著的差异(P > 0.05)，而pH与各类有机物在表层与深层之间表现出明显的不同(P < 0.05)。特别是土壤总氮和总有机碳含量，在浅层土壤内(0−30 cm)的变化也极其显著。总体而言，土壤理化性质具有一定的垂直梯度变化，从表层到深层呈现营养物质递减的趋势。

2.2 微生物数量、组成与结构

64个土壤样本共得到473余万条高质量原核微生物16S rRNA基因序列，这些序列被划分为37 720个OTU。由于各样本序列差异较大，将所有样本统一到相同的总Reads数量(20 291个)对每个样本进行随机重抽，重抽后的OTU表格进行后续的群落分析。使用原核生物OTU数量、微生物总量表示群落的α多样性。每个样本中的微生物OTU数量在2 135−4 108之间变化，随土壤剖面向下发展，每个层次的平均OTU数目明显减少(图 1A)。荧光定量qPCR获取的16S rRNA基因拷贝数结果显示(图 1B)，0−10 cm土壤深度的基因平均拷贝数最大，70−80 cm最少，并随土壤深度的增加而减少，总体呈指数级的变化(0−10 cm的10¹⁰到70−80 cm的10⁸)，说明微生物总量在深度上的差异比较大。总体而言，微生物OTU数量和总量随深度的增加均梯度减少(图 1)。但是通过ANOVA Turky事后比较检验发现接近地表的三层土壤中的微生物α多样性没有明显的差异，然而表层与深层的微生物α多样性有极显著的差异，即较深层的微生物α多样性显著小于表层土壤。10个微生物丰度最高的类群中，9个属于细菌，只有一个类群属于古菌(图 1C)。在所有样本中，放线菌(18%−23%)、奇古菌(13%−23%)、变形菌(12%−16%)和酸杆菌(11%−18%)是群落中优势菌群，占整个土壤剖面原核生物总数的60%以上。与OTU和总量的变化相反，微生物分类学组成在土壤剖面中似乎受深度变化的影响有限。虽然奇古菌、未分类、拟杆菌和芽单胞菌在部分深度间有显著性差异(P < 0.05)，但是没有呈现出梯度变化。以上结果说明，在较小的深度差异上，微生物多样性也可以表现出一定的差异，即表层的物种数量和总量与深层显著不同，相邻土壤中的微生物变化较小，但一致表现出多样性随深度增加而下降的趋势。然而深度对原核生物的组成的影响较小，表现为多数微生物相对丰度比较均匀。

为了进一步探究不同土壤深度下原核微生物群落之间的差异，使用基于Bray-Curtis距离矩阵的非度量多维尺度分析描述土壤剖面群落组成的总体格局(图 2)。结果显示在平面空间相距6 000 m不同取样点获取的八层土壤剖面样本，其微生物类群之间依然呈现出垂直梯度的分异。同时此结果也显示在相邻的两层土壤中，部分样本微生物群落的位置发生重叠。通过对菌群组成的差异分析证实了这一结果，表明菌群结构存在一定程度的明显差异(表 2)。检验结果与NMDS图一致，不同深度之间的微生物群落有显著差异，但是相邻两层之间没有显著的差异性。在任何深度的土壤中，微生物的聚类距离越近，表明群落拥有更加相似的结构组成。其中，表层的微生物聚集紧密，随着深度的增加，下层的微生物位置越来越分散(图 2)，说明在深度增加的过程中，微生物的群落在每一层土壤的变异是增大的。以上结果说明在土壤内部，虽然深度差较小，但还是可以在表层土壤和深层之间发现显著的多样性变化和结构组成差异。表层群落组成基本一致，而深层存在较强的成分变异性，并且这种变异随着深度逐渐增大。

2.3 驱动微生物群落变异的环境因子

为了研究土壤理化性质对原核微生物群落的影响，首先应用Mantel和Partial mantel检验确定不同的环境因子对土壤原核生物群落结构(Bray-Curtis距离)的影响强度(表 3)。如果不排除环境因子之间的相互影响，统计结果表明，包括含水率、pH、总氮、氨态氮、硝态氮和总有机碳在内的所有理化参数都对原核生物群落结构施加显著的影响，共同改变了微生物群落结构。然而排除其他因子对特定某个土壤参数的影响，Partial mantel结果显示氨态氮和硝态氮的影响极小，而含水率、pH、总氮与总有机碳的作用依然很显著。进一步通过典范对应分析方法对所有的原核生物群落结构与土壤理化性质的关系进行分析(图 3)。原核微生物群落可以被CCA1和CCA2轴较好地解释(F=0.29，P=0.001)，其群落结构与土壤中总氮(TN)、土壤总有机碳(TOC)和土壤pH等理化性质存在较强的相关关系，说明在剖面上影响最为显著的因素可能是总氮、总有机碳和pH值。从图 3中可以看出，在各种环境因子的影响下，原核生物群落大致沿CCA1轴分布，表明CCA1可能是这些因素叠加作用下所代表的土壤深度影响。所以深度的影响是包括了各种土壤剖面环境参数变化在内的综合作用，是影响原核生物群落分布的重要因素。

2.4 微生物生态网络

为了了解不同深度水平的原核微生物相互作用及其变化，通过筛选并保留了在50%的样本中同时出现的OTU构建了不同深度下微生物物种间的互作网络(MENs) (图 4)。MENs在不同土壤深度条件下均以相同的相似阈值0.913构建，使不同网络的拓扑结构系数可以直接进行比较(表 4)。一个节点(node)代表一个OTU，每条边(edge)分别表示由其连接的两个节点之间的相关性。整体拓扑指数显示，所有网络连通性分布曲线均与幂律模型(R²值为0.785−0.888)吻合良好，表明网络的无标度(scale-free)属性。这意味着MENs中较少的微生物与其他微生物之间的联系较多，而大多数微生物之间的关联较少^[47]。这些网络拓扑指数与具有相同节点数和边的随机网络有显著差异(P < 0.01，表 4)，表明网络结构不是随机无序构建的，其相互连接具有显著的确定性意义。

土壤不同深度的原核微生物互作网络的连接紧密程度存在差异，并且这种差异随深度的增加而增大(图 4)。然而各种网络拓扑指数也一致显示，表层土壤中原核微生物的相互关联与深层土壤中有很大差别(表 4)。更多的节点和边意味着微生物间的网络更大且其连接可能更复杂。由于网络节点和边都与深度具有负相关关系(图 5)，表层明显拥有较多的节点个数与强相关性关系，增加了OTU之间连接的密度，所以使得表层土壤微生物形成了比深层更大与复杂的分子网络模式(表 4和图 4)。此外，由于深层原核生物间相互关联较少，容易造成网络的相互分离(图 4)。网络复杂度的降低还反映在从表层到深层逐渐降低的平均度(average degree)、聚类系数(average clustering coefficient)以及增大的测量距离(harmonic geodesic distance)，这些网络拓扑参数意味着网络连接程度降低，反映了微生物MENs是逐渐变为松散的(图 5)。综合以上结果可以看出，微生物在土壤剖面的垂直分布过程中，深处的生物群落相互联系比表层相对较弱。

根据微生物物种之间的互作关系，它们会形成相互独立或联系的模块(module)出现在群落中，所以对复杂网络中模块的识别可以揭示这些网络的组织结构，并且判断其如何进一步成为紧密联系的群落。任何一层原核生物群落的MENs模块化值(modularity)均在0.779−0.932之间，显著高于相应的随机化网络的模块化值，因此所构建的所有MENs都表现为模块化属性(表 4)。每个不同深度的MENs被划分为一些包含紧密相关原核生物类群的离散模块，分别得到49−74个相互分离的模块(表 4和图 4)。模块内共存的类群被认为具有共同的环境偏好或是以相互促进的方式发生互作。为了判断网络模块之间的变化以及如何受土壤环境影响，在每一层土壤原核生物MENs中挑选前5个最大的网络模块进行分析。模块中的物种组成是变化的，它们包括不同的OTU数目和分类(图 6)。0−10 cm最大的模块7包含8个不同的门62个OTU，10−20 cm最大的模块4包含8个不同的门53个OTU，20−30 cm最大的模块3包含7个不同的门38个OTU，30−40 cm最大的模块4包含9个不同的门43个OTU，40−50 cm最大的模块5包含9个不同的门36个OTU，50−60 cm最大的模块1包含7个不同的门23个OTU，60−70 cm最大的模块1包含6个不同的门17个OTU，70−80 cm最大的模块1包含8个不同的门17个OTU。Mantel检验用来揭示所有挑选的模块与各土壤环境变量之间的关系(图 6)。模块丰度与土壤因子的比较表明，10−20 cm模块4丰度的增加可能与土壤含水率和氨态氮的含量降低有关，然而在氨态氮、硝态氮和总有机碳较低的40−50 cm，模块5的丰度却相对较高。此外，聚类分析表明不同深度的模块可以被划分到相同的分类簇(如60−70 cm模块1、30−40 cm模块4和40−50 cm模块9)，而且这些模块所受环境条件影响相当，表明它们可能具有相同或类似的生态位。研究结果表明，网络的模块化结构主要由环境因素所决定。

每个节点都会在网络中扮演不同的拓扑角色，分析模块中节点作用对识别关键种群具有重要意义。为了识别各节点在网络中可能扮演的拓扑角色，根据节点的模块内连通性和模块间连通性值将节点区分为外围节点、连接节点、模块中心点和网络中心点四类(图 7A)。每个网络中的大多数节点都是外围节点，因为它们的大部分连接都在各自的模块中(图 7A)。虽然网络中心点被当作是超级通才(super generalists)，但是没有在任何网络中发现网络中心点。此外，连接节点是与其他模块高度连接的节点，模块中心点是与自己模块中的许多其他节点高度连接的节点，它们被认为是网络中的通才(generalists)^[47]。我们在6/8的不同MENs发现了不同的连接节点(1−11个)，然而没有在最深处的两个网络中找到连接节点。同时，在好几个MENs中观察到模块中心点：0−10 cm包含4个、10−20 cm包含3个、20−30 cm包含1个、30−40 cm包含3个、40−50 cm包含3个，其他层的网络中的没有发现模块中心点(图 7B)。这些模块中心点都属于细菌，其中，7/14是酸杆菌门，3/14是放线菌门，1/14拟杆菌门，1/14未分类，2/14疣微菌门。由于连接节点和模块中心点是关键的分类单元，并且在网络拓扑结构中具有重要作用。在本研究中属于酸杆菌门的连接节点和模块中心点为16/44，并且几乎每一层都可以发现该分类群的功能节点(0−10、10−20、20−30、30−40和40−50 cm)，鉴于此，酸杆菌门类微生物可能是该农牧交错带草地生态系统中最为重要的关键物种。

3 讨论

微生物群落多样性的变化、物种间相互作用及其参与的生物地球化学循环过程被证实在深层土壤中也持续存在。农牧交错带是一个脆弱的生态系统，作为两种不同生境之间的过渡区域，生存于土壤不同深度的微生物群体对农牧交错带的生态功能发挥着巨大的作用。最近的研究表明，交错带对环境变化具有高度的响应性，尤其是在干旱和半干旱地区^[30]。然而土壤微生物对环境变化极为敏感，而且可以同时对各种变化作出响应，因此，研究其地下微生物群落的构成机制和对环境因子的响应显得越来越重要^[52]。确定不同土壤深度微生物群落的组成和特征，将有助于了解重要的地下微生物活动。

3.1 土壤深度影响原核微生物的总量、群落多样性与物种组成结构

本研究表明，在农牧交错带草地土壤中，原核微生物的总量和多样性沿土壤垂直梯度上呈现极为规律的梯度变化。在这一区域母质层上方的土层范围内(0-80 cm)，原核微生物的物种丰富度(OTU richness)与微生物总量都随土壤深度的增加(每10 cm)而显著减少(图 1)。微生物总量和α多样性的这种减少趋势与以往的许多其他研究结果相一致^{[14, 16, 53-54]}。之前的研究均表明，土壤微生物群落的丰富度、多样性和组成具有层位特异性，在不同的地层中存在差异^{[11, 55]}，但这些研究中土壤分层都较宽，因此梯度变化趋势不连续或者规律不清晰。本研究以每10 cm为间隔，其结果显示相邻的深度间微生物群落表现出的差异不显著(表 2)，但整体上这种分化差异是与深度成正相关关系并逐渐扩大，即深层土壤的原核生物群落β多样性差异大于浅层的微生物结构变异(图 2)。这些连续变化趋势对认识土壤分层也有着重要的参考价值。土壤深度对微生物群落施加的影响包括土壤质地、物理特征(压力、空隙度等)和各种理化属性(pH、含氧量、含水率、营养物质等)共同的作用。研究的土壤剖面具有较强的环境梯度，多种土壤因子随深度的增加发生变化(表 1)。这些土壤环境参数与微生物群落的相关性表明(表 3和图 3)，所有因子都会直接或间接地影响群落的结构，特别是微生物数量随深度的增加成指数级下降，且与总氮、硝态氮和总有机碳含量随深度的增加而下降之间具有一致的趋势。因此，这些微生物群落的垂直差异在很大程度上可以归因于各类资源可利用性随土壤深度的下降。以上结果共同表明，土壤深度是影响微生物多样性与结构的重要因素，而土壤深度造成土壤环境的梯度变化是微生物物种垂直分布特征的最重要因素。

3.2 微生物互作网络的垂直变化特征

生态网络相关关系可以用来代表物种构建过程中的生物因素，由于微生物群落的复杂性，评估微生物互作关系一直都是一个不小的挑战^[56]。本研究采用针对测序数据优化的SparCC算法，结合分子生态网络分析(MENA)构建不同土壤剖面下的微生物物种互作网络。结果显示，与表层的土壤原核生物群落相比，深层的群落形成了相对更小、更简单且分散的网络结构，而且这种变化是随着土壤深度的增加呈现明显的梯度变化(图 4和表 4)，说明原核生物之间的相互作用关系随深度增加逐渐减弱。微生物数量在土壤深处相对较少，并且许多研究也证明各类微生物随土壤深度增加而迅速降低^{[13-14, 16]}。因此，生活于深处的微生物数量少且迁移速率低，它们接触产生联系的机会变少，这可能导致了其网络大小及复杂程度的不同变化。虽然之前的研究发现农田的土壤剖面中微生物生态网络的结构较为稳定，没有随深度体现出规律性的变化^[57]。这种观测结果的区别可能来自于生态系统的差异和不同系统中微生物对环境因子的敏感程度。农田的土壤受长期耕种和人为管理的影响，其基质分布更均匀，而本研究针对的是半干旱区的农牧交错带，其土壤分层中微生物群落不容易受到人为耕作的影响，进而保留了物种间互作关系上的差别。

模块化是生态网络中一个重要的概念，其量化了网络在多大程度上可以分解成更小的组件，可能指示控制群落结构的生态过程^[58]。本研究中不同水平的MENs模块化指数较高、相互隔离的模块数目众多，反映了原核生物群落特殊的相互关系、差异的环境生态位或重叠等^[59]。通过观察网络模块与土壤环境因子的相关性，发现部分网络可以被聚类到一起，表明这些模块内的微生物可能具有相近的生态位，而且受影响的程度相似(图 6)。部分模块的丰度变化可能与营养、资源的变化相关(图 6)，因为资源可利用性是塑造网络结构的重要驱动因素，揭示了环境属性的变化会改变生态网络的结构组成^[60]。本研究中总有机碳和总氮随深度的衰减一定会造成下层微生物营养欠缺，刺激了下层微生物间的竞争，同时也改变了网络的结构与组成。

虽然微生物网络以细菌OTU为主，但同时也有古菌的出现(图 6)，这突出了在土壤微生物网络中功能菌之间联系的重要性。例如，Tu等^[61]对二氧化碳升高情况下草地生态系统中固氮微生物群落变化的研究，揭示了固氮生物类群与其他微生物种之间共存的状况(固氮螺菌/放线菌、中间根瘤菌/束缚菌、慢生根瘤菌/酸杆菌)，说明在生态系统中，部分微生物对其他类群微生物具有潜在的吸引力，通过加强相互之间的联系发挥功能作用。本研究中，通过识别网络中的关键物种，发现酸杆菌门是多个土壤剖面中的优势菌属，虽然酸杆菌门类群的丰度不是最高且随深度降低(图 1)，但是依然可以在不同深度的互作网络中占据重要地位(图 7)。酸杆菌是土壤中分布最广的细菌类群，在不同的土壤环境中大约有20%的微生物群落^[62]。酸杆菌种群具有极其广泛的多样性，目前共分为26个亚类(subdivisions，GP)^[63]。这些酸杆菌亚类具有碳降解、氮代谢、表达多种活性转运体和产生胞外多糖的能力，并且生活方式存在差异^[64]。把这些关键的酸杆菌细分到亚类，我们得到Gp4、Gp6、Gp7和Gp16亚群。这些微生物会根据不同的环境选择压力改变其生存策略，并且完成在营养匮乏环境下的自身养分吸收，进一步协同其他微生物物种发挥作用，这对农牧交错带中酸杆菌微生物功能发挥和系统稳定保持可能至关重要。

4 结论

在本研究中，我们调查了典型农牧交错带中土壤剖面上的原核微生物总量、多样性和互作网络结构与环境因素的交互作用，揭示了随土壤剖面发展变化的微生物群落模式。原核生物的多样性与深度成负相关，这与许多其他研究结果一致。原核生物群落结构变异会逐渐扩大，这与地下的环境差异关系显著。然而分子生态网络方法揭示了原核生物之间的相关关系也随深度变化，并且相关性会减弱，相较于浅层土壤，表现为不再紧密聚集的原核微生物分子网络结构。这些结果为我们提供了农牧交错带微生物动态变化的重要见解，有助于理解支撑群落结构和生态过程的关键驱动因素。我们推测一旦当地受到较大的环境扰动，底部的微生物群落可能难以快速调整适应，将使生态系统功能发生显著变化。因此，未来可以进一步增加对微生物时间序列的宏基因组学研究，以及其他微生物功能的研究，有望加深我们对土壤微生物群落复杂性及其与生物地球化学因子关系的理解，为生态保护策略的制定提供科学依据。