根据《中国海洋环境状况公报》(1997–2015)，石油类污染物是中国近海海域主要污染物之一^[1]。生物降解被认为是污染环境中石油烃降解的重要途径，微生物是最重要的降解者之一^[2]。石油组分比较复杂，主要包括饱和烃、芳香烃以及胶质、沥青等非烃类物质。相对于饱和烃，更难被降解的芳香烃具有更高的“三致效应”毒性，属优先检测的有机污染物^[3]。四环及四环以上的多环芳烃(polycyclic aromatic hydrocarbons，PAHs)因其结构特殊并具有较高的热力学稳定性而可在环境中持久存在，其“三致效应”更加明显。因此，研究微生物降解多环芳烃的机理机制是当前石油烃生物修复研究领域的重要方向。目前实验室内已驯化分离到数量可观的多环芳烃降解菌，多属于Proteobacteria、Actinobacteria、Cyanobacteria以及Firmicutes^[4]。这些微生物在纯培养条件下厌氧降解^[5]或好氧降解^[6]多环芳烃的机制也已被很好地阐释。值得注意的是，近年的研究认为在复杂微生物群落中广泛存在由不同微生物之间复杂共生关系(互利共生^[7]、竞争^[8]等)构成的生态网络^[9]，并且生态网络可以影响微生物群落的变迁^[10]。一般认为，存在协同作用的微生物菌群具有相对更高的降解效率^[11]。我们的前期研究也表明微生物之间的共生关系是它们共适应降解石油烃及重建稳定微生物生态系统的重要途径^[12]。因此，研究石油烃修复过程中微生物群落的变化尤其是微生物之间的互作网络对开发更高效的微生物修复技术具有重要的理论指导意义，同时也有助于更好地理解微生物生态系统如何响应环境污染。

由于存在共生关系的微生物在丰度上表现出一定的相关性，例如互利共生微生物的丰度会表现出正相关^[7]，而处于竞争关系的微生物的丰度会表现出负相关^[8]，因此借助微生物的丰度数据我们可以构建不同微生物之间的生态网络进而推测微生物之间的相互作用^{[9, 13-14]}。还可以根据生态网络的拓扑性质了解“关键种”^[15]及互作关系。目前关于污染环境中微生物群落结构变化的研究已有很多，但多数仅仅关注了群落组成及相关功能基因的丰度变化。一方面我们认为优势菌是影响群落功能的主要因素，因此研究焦点便会更多地集中在优势菌上。例如在前期研究中，我们通过建立实验室驯化系统研究了芘和苯甲酸钠胁迫后渤海沉积物中微生物群落结构的变化^[16]。实验结果表明芘和苯甲酸钠显著改变了沉积物中微生物群落的组成，极大降低了海洋沉积物中菌群的丰度和多样性，其中海杆菌属(Marinobacter)、假海源菌属(Pseudidiomarina)及碳酸噬胞菌属(Aequorivita)得到了极大富集。另一方面，低丰度菌在近几年也被认为对群落功能有巨大贡献。例如Adam等人^[17]利用¹³C标记的芘研究土壤中微生物群落对芘胁迫的响应，结果表明土壤中的降解微生物种类繁多，起到关键作用的并不仅仅是丰度较高的优势菌，还包括数目庞大的低丰度菌。因此如何系统地理解群落功能的执行便成为重要的前沿科学问题，基于系统生物学理论的微生物生态网络是其中一种重要手段：微生物生态网络可用于揭示不同分类的微生物之间的生态互作^{[9, 14]}，而网络的拓扑属性被用于研究微生物群落应对环境改变的稳定性^{[10, 12]}。得益于高通量测序技术和生物信息学的迅猛发展，基于16S rDNA的微生物群落数据和各种算法(如eLSA^[18]，CCLasso^[19])极大促进了微生物生态网络研究的发展。同时海量的相关基因组信息及代谢模型数据，也为更深入研究微生物之间的相互作用提供了必要基础^[20]。Jeanbille等人^[21]通过构建基于Spearman相关性的微生物生态网络研究了法国毗邻地中海及大西洋的近海海域表层沉积物中微生物群落对持续性多环芳烃污染的响应，结果表明环境扰动显著影响了微生物生态网络的结构。然而目前关于多环芳烃降解过程中微生物生态网络的研究仍然不足：已报道的研究总数较少，且其中陆地环境研究较多(如赵慧慧等人的研究^[22]表明长期石油污染显著影响了胜利油田周围盐碱化土壤中的微生物分子生态网络的拓扑结构)，海洋环境较少。

芘是由4个苯环构成的PAHs，是检测PAHs污染的指示物及研究PAHs生物降解的模型分子^[23]。因此本文拟构建室内微宇宙体系，以来自渤海表层沉积物的微生物群落为研究对象，以芘为添加污染物，研究芘降解过程中微生物群落结构变迁，重点研究芘降解过程中的微生物生态网络。本研究将有助于更好地理解功能菌群在生态系统内的相互作用及共适应机制这一前沿科学问题，也将有助于探索如何调节功能菌群使其发挥最大降解潜力。

1 材料和方法 1.1 沉积物采集及处理

沉积物于2016年1月17日使用箱式采泥器从渤海(38.95 N 118.95 E)采集，取表层1–2 cm沉积物约500 g置于灭菌烧杯中，于冰盒中搅拌30 min以充分混匀后转移至无菌密封袋中并贮存于4 ℃，48 h内带回实验室进行样品前处理。处理方法如下：取约450 g沉积物再次混匀(方法同上)，利用真空冷冻干燥机(岛津，日本)对其中约30 g沉积物进行冷冻干燥，干燥后的样品置于–80 ℃保存用于总DNA提取。

1.2 试剂及培养基

二氯甲烷、正己烷及无水硫酸钠均为国药分析纯，甲醇为阿拉丁产分析纯，芘(99.9%)、蒽(99.9%)，购自Sigma公司。无机盐培养基(MSM)组成为：(NH₄)₂SO₄ 1.0 g，Na₂HPO₄ 0.8 g，KH₂PO₄ 0.2 g，MgSO₄·7H₂O 0.2 g，FeCl₃·3H₂O 0.005 g，(NH₄)₆Mo₇O₂₄·4H₂O 0.001 g，CaCl₂·2H₂O 0.1 g，NaCl 20 g，溶于1 L去离子水中，pH为7.2，121 ℃灭菌20 min。芘母液：用二氯甲烷配制芘母液，使其浓度为0.4 g/mL。蒽母液：二氯甲烷配制，浓度为0.01 g/mL。

1.3 微宇宙实验设计

实验设计见表 1。取63个250 mL锥形瓶(已灭菌)，设对照组(21个)和实验组(42个)。向各瓶中加入100 mL灭菌MSM培养基、100 µL芘母液，使芘的终浓度为0.4 g/L，待二氯甲烷完全挥发后使用。向实验组各瓶加入10 g沉积物，并将对照组及实验组所有瓶置于摇床中培养(30 ℃，60 r/min)。分别在第0、3、6、12、18、24、30天取样，每次取3个对照组，6个实验组。其中3个对照组及3个实验组用于萃取分析芘降解率，3个实验组冷冻干燥后保存于–80 ℃，用于总DNA提取。

1.4 芘降解率测定

以蒽为内标参照我们的前期研究^[12]进行降解率的测定。主要包括超声萃取后过滤收集液相、分液漏斗收集有机相并用无水硫酸钠除水、旋蒸后用甲醇进行溶剂置换并定容至15 mL、0.45 µm滤膜过滤并–20 ℃保存、HPLC检测。检测条件为：流动相(乙腈:水)=80:20；流速1.0 mL/min；温度(T) 30 ℃；检测波长254 nm；进样量10 µL；色谱柱为C₁₈反相柱(4.6 mm×250 mm，SunFire)。降解率计算公式为：

A：对照组中芘的含量(A₀，0天芘的含量；A_t，t天芘的含量)；

B：实验组中芘的含量(B₀，0天芘的含量；B_t，t天芘的含量)。

1.5 DNA提取及测序

DNA提取使用MoBio PowerSoil^TM kit (MoBio，Carlsbad，CA，USA)试剂盒并参照使用说明进行。本研究共7个取样时间(T0 d、T3 d、T6 d、T12 d、T18 d、T24 d、T30 d)，每个时间点3个平行，共计21个DNA样品。样品被送至北京诺禾致源公司(北京，中国)，使用Illumina Hiseq 2500对V4-V5高变区进行测序。原始序列已上传至SRA数据库，样品编号为SAMN06216651– SAMN06216671。

1.6 生物信息学分析

使用USEARCH 9.0^[24]对原始序列进行聚类分析，主要步骤包括双末端序列装配、去除低质量序列、去重复、去嵌合体、在97%相似性水平上进行操纵分类单元(operational taxonomic unit，OTU)聚类分析、全局比对后生成含有OTU在各样品中序列数的表格，以上步骤的参数采用默认推荐参数值。使用Greengenes (version 13.8) 数据库对获得的OTUs代表序列进行分类水平的鉴定。由于引物是面向细菌的V4-V5区进行设计，且只有6个OTUs被鉴定为古菌(远低于古菌的实际种类数)，因此在界水平上被分类为古菌及不能被分类的OTUs将被移除。利用Mothur软件^[25]使用重抽样的方式将所有样品的测序深度标准化到49000条序列并计算Good’s coverage、Shannon-Wiener、Chao1指数，运算执行1000次后取平均值。计算各OTUs在各样品中的相对丰度并分别在纲及属水平上进行汇总。基于非参数Kruskal-Wallis秩和检验统计在T0 d和T30 d中相对丰度存在显著差异的OTUs (P＜0.05)，将这些OTUs称为差异OTUs。根据它们的相对丰度变化情况，将差异OTUs分为升高组(increased OTUs，i-OTUs)与降低组(decreased OTUs，d-OTUs)，汇总T0 d及T30 d的i-OTUs及d-OTUs数目及总相对丰度，并利用气泡图可视化二者在分类水平上的汇总信息。鉴于我们的研究对象为可能影响芘生物降解的细菌，因此只有在T30 d的样品中平均相对丰度大于0.05%的OTUs被纳入分析。基于OTUs的相对丰度及CCLasso算法^[19]计算微生物之间的相关性(r＞0，P＜0.05) 并构建微生物生态网络(Microbial ecological network，MEnet)；统计在MEnet中i-OTUs与d-OTUs的数目。上述分析中如无特别指出，计算及绘图均借助R软件^[26]进行。

对未能在属水平上获得分类信息的OTUs (属未知OTUs)，借助本地BLAST 2.2.22程序、OTUs的代表序列及ModelSeed数据库(http://modelseed.org/)获得近源基因组(16S序列相似性>97%)的物种信息并以此补充属未知OTUs的属水平分类信息^[27]。构建属水平MEnet (g-MEnet)，则网络中的点代表属，边代表连接的2个属之间存在相关的OTUs的对数(定义为NC值)。根据网络的拓扑属性探究重要菌群及他们之间的关系。

2 结果和分析 2.1 芘的降解率变化趋势

对照组的芘含量在30 d内损失不足1%，表明芘的降解主要由微生物降解驱动。如图 1-A所示，实验组芘的降解率在30 d内持续增加且在30 d时达到(67.07±2.37)%，但降解速率(即降解率的斜率)以3 d为界呈现先增加后降低的趋势。

2.2 微生物群落结构及变迁

经过质控及筛选，21个样品共获得1330269条序列；在97%相似性水平上共获得3712个OTUs。在所有样品中，代表覆盖度水平的Good’s coverage指数均大于98%，表明本研究的测序深度可以很好地反映微生物群落的多样性及群落结构。代表群落物种丰度水平的OTUs获得数(图 1-B)及Chao1指数(图 1-C)在0 d分别为2475.12±94.18、3170.59±67.53，而在30 d则为972.89±180.39、1547.04±217.84，总体趋势表现为持续降低。可同时反映多样性及均匀度的Shannon-Wiener指数(图 1-D)在0 d的值为5.52±0.21，在12 d达到最低值2.98±0.87，随后升高并于第24天趋于稳定(24 d：3.95±0.16，30 d：3.83±0.34，P＞0.05)。

在纲水平上，沉积物中的优势菌(0 d)主要为Alphaproteobacteria [(8.86 ± 1.57) %]、Deltaproteobacteria [(24.72±4.12) %]以及Gammaproteobacteria [(36.73±6.10)%]。30 d后，优势菌主要为Alphaproteobacteria [(45.05±14.83)%]、Gammaproteobacteria [(30.68±7.04)%]及Flavobacteriia [(7.25±6.62)%] (图 2-A)。在属的水平上，沉积物中的优势菌主要包括Desulfococcus [(5.16 ± 1.03) %]、Pseudidiomarina [(2.63±0.79)%]及Marinobacter [(2.54±0.91)%]等，经过30 d培养后的优势属为Flavobacterium [(4.84 ±1.91) %]、Bacillus [(4.46±1.83)%]及Marinobacter [(3.45±1.28)%]等(图 2-B)。通过统计i-OTUs及d-OTUs的分类信息表明：(1) 共计102个i-OTUs与1551个d-OTUs，与前文alpha多样性降低的结果一致；(2) 对i-OTUs而言，在纲水平上主要属于Alphaproteobacteria、Flavobacteriia、Planctomycetia及Bacilli (图 2-C)，在属水平上主要属于Flavobacterium、Bacillus、Planctomyces及Erythrobacter等(图 2-D)；(3) 对d-OTUs而言，在纲水平上主要属于Deltaproteobacteria、Gammaproteobacteria、Anaerolineae及Spirochaetes等(图 2-E)，在属水平上主要属于Desulfococcus、Spirochaeta及Pseudoalteromonas等(图 2-F)。

2.3 微生物生态网络分析

筛选后共计127个OTUs被用于预测MEnet，结果表明MEnet由53个点及317条边组成，其中属于i-OTUs与d-OTUs的分别为37个与0个(表 2)。在g-MEnet中，共有29个点143条边(图 3)。点度最高的点为Unclassified Flavobacteriaceae (degree=27)，其次为Unclassified alpha proteobacterium Q-1 (degree=19) 及Planctomyces (degree=19) 等(表 3)。NC值最高的边为Erythrobacter-Unclassified Flavobacteriaceae (NC=17)，其次为Erythrobacter-Unclassified Alphaproteobacteria (NC=15)、Erythrobacter-Unclassified alpha proteobacterium Q-1 (NC=12) 等(表 3)。

3 讨论

一般来说，包括芘在内的石油烃对微生物群落的影响主要体现在两个方面：一是它们可以作为部分细菌生长所需的唯一碳源，如Bacillus、Cyanobacteria、Cycloclasticus、Flavobacterium、Marinobacter、Mycobacterium、Planctomyces、Pseudomonas、Pseudoalteromonas、Pseudidiomarina、Rhodococcus和Vibrio等^[28]；另一方面，它们对某些微生物具有毒性。因此在受到石油烃胁迫后，微生物群落结构将发生明显改变^[12]。我们前期的研究表明，添加芘胁迫后，渤海表层沉积物中的微生物群落结构发生了明显改变^{[12, 16]}：一部分可能具有降解功能的微生物的相对丰度明显增加，另一部分微生物的相对丰度明显降低甚至变为0进而引起相关群落生态功能的缺失。因此通过分析微生物相对丰度的变化趋势，我们可以推测芘降解过程中的潜在功能菌，也可以获知胁迫对群落生态功能的影响。本研究中，属于i-OTUs的Flavobacterium、Bacillus、Planctomyces等的多环芳烃降解能力均已有报道^[28]。值得注意的是，58.8%的i-OTUs属于Alphaproteobacteria，但它们在属水平上的分类信息多为未知，提示我们对渤海沉积物中多环芳烃降解菌的认识仍非常缺乏。沉积物中的Deltaproteobacteria (如本研究中明显降低的Desulfococcus属)多参与硫酸盐还原^[29]，同时在沉积物环境中，石油烃的厌氧降解多与硫循环过程耦合^[30]。因此本研究中Deltaproteobacteria的相对丰度显著降低，一方面可能与Deltaproteobacteria多厌氧而本研究的培养条件存在氧气输入有关，另外一方面也暗示多环芳烃污染对渤海沉积物中的硫循环过程存在一定的不利影响。

由于绝大多数微生物的不可培养性，使得我们无法采用纯培养的方式获取微生物之间的共生关系。但由于共生关系在丰度相关性上的体现，使得我们可以基于微生物群落组成及它们的相对丰度变化趋势预测微生物之间的共生关系^{[9–10, 14–15]}。微生物之间的共生关系表现在丰度相关性上可以是正相关也可以是负相关。在环境中，非生物因素(如营养盐水平)的生境选择和生物因素(微生物之间的相互作用)是影响丰度水平的重要因素^{[10, 31–32]}。因此，不同细菌之间存在丰度正相关可能是由于他们之间存在功能互补如互利共生(直接相关)，也可能是由于它们有相似的环境要求(间接相关，如资源丰富的有利环境引起的丰度同步增加或者不利环境引起的丰度同步降低均可在相关性上表现出正相关)。值得注意的是，传统计算相关性的算法(如Pearson相关，Spearman相关等)并不能区分直接相关还是间接相关。这样不仅会增加分析的复杂性，还会在分析生态网络的拓扑属性时带来更大的误差。本研究使用的CCLasso算法^[19]借助最小二乘法和增广拉格朗日乘子法移除第三方效应，使得构建的生态网络不包含间接相关。在MEnet中不存在属于d-OTUs的点，表明由于芘胁迫引起的相对丰度同步降低而表现出正相关的OTUs并未包含在我们构建的生态网络中。因此本研究中在丰度变化上表现出正相关的微生物之间可认为存在互利共生或者协同作用。在T30 d时平均相对丰度大于0.05%的i-OTUs共计49个，而在MEnet中属于i-OTUs的共有37个，表明大多数的i-OTUs的相对丰度变化不仅和芘有关，还和它们之间可能存在的互利共生或协同作用有关。

寻找“关键种”是微生物生态领域的关键科学问题，基于微生物生态网络的拓扑属性(主要指点度中心度：在1个网络中与1个点直接相连的其他点的个数)寻找核心微生物是目前比较认可的方法^{[9, 14–15]}。基于点度中心度，本研究样品中的核心微生物包括Planctomyces, Erythrobacter等(表 3)。此外，代表互作关系的边是研究生态网络的另一个重要内容^[33]。为了评价2个不同属微生物之间的互作强度，我们引入NC值进行评价。NC值受以下2个方面因素的影响：每个属G所含OTUs的数目N；2个不同属微生物之间存在互作的倾向性P。N和P的值越高，NC值越高。一方面N值越高意味着属G在T30 d时所占比例越高，则属G在芘降解过程中越趋向于重要；另一方面P值越高两个属之间越倾向于互作。因此在生态网络中NC值越高的边越趋向于重要。值得注意的是，在g-MEnet中NC值较高的边(NC>5) 多连接到属Erythrobacter。综合点度中心度我们认为Erythrobacter为本研究中最重要的微生物类群，在我们的前期研究中发现，Erythrobacteraceae的微生物是重要的多环芳烃降解者^[12]。