泥火山是由泥浆形成的，只因它的喷发形式和喷发后形成的堆积物在外形上与火山有些相似，因此称其泥火山^[1]。新疆乌苏具有亚洲最大的泥火山群，其特殊的环境是研究微生物资源的最佳场所^[2]。土壤微生物是土壤的重要组成部分，土壤微生物群落的变化在一定程度上可以反映土壤质量的变化趋势^[3]，而细菌是土壤中最丰富、分布最广泛的微生物类群^[4]，广泛参与土壤营养物质循环过程^[5-6]。传统的土壤微生物培养方法存在着耗时、工作量大的缺点，并且自然界微生物多是未培养的，从而导致无法全面解释微生物群落结构。近年来，随着高通量测序技术的不断发展^[7-9]，人们对土壤中微生物的研究日渐深入^[10]，它被广泛应用于环境微生物菌群结构与演替规律研究，其中，Illumina公司开发的MiSeq方法避免了通量低、准确率低、操作复杂等问题^[11]，与454焦磷酸测序相比，通量更高，操作简便、成本更低，更适合分析环境中群落结构多样性^[12]。本研究通过Illumina MiSeq测序技术分析新疆乌苏泥火山微生物的群落结构，这将为深入研究泥火山生态系统奠定基础，为合理利用和开发泥火山微生物资源提供指导。

乌苏泥火山位于新疆天山北麓乌苏市南约 40 km白杨沟镇东南约2 500 m处，属于天山北坡山前丘陵地带，四周群山环抱，属大陆干旱气候，空气干燥，降雨稀少，夏季炎热，光照充足，年均降水158 mm，平均气温7.3 ℃，其环境条件较为特殊。本课题以泥火山口为参照，选取具有代表的4种不同生境(命名为A、B、C、D)为研究对象，各个生境基本情况见表 1。

以4种不同生境(A、B、C、D)为研究对象，取样时间为2014年4、7、11月，以正方形5点取样法在样地内设置样方，草本采用1 m×1 m的样方，5次重复；灌木采用5 m×5 m的样方，3次重复。每个点取样量大体一致，去除土样里的植物根系和砾石，然后将重复的土样混合均匀，做好标记装入灭菌封口聚乙烯袋，低温保存带回实验室。将土样分成两份，一份在室内风干，将风干土分别过1 mm和0.15 mm筛子，用于理化因子的研究；一份放入-20 ℃，用于后续的分子方面的研究。最终，混合后的土壤样品共12个，根据采样季节不同，分别标记为A4、A7、A11，B4、B7、B11，C4、C7、C11，D4、D7、D11。

引物，上海生工合成；PCR和酶反应体系中DNA磁珠法回收纯化试剂盒，志昂生物科技有限公司。PCR仪，安捷伦科技有限公司；电泳仪和电泳凝胶成像系统，北京君意东方电泳设备有限公司。

1.4.1 土壤理化性质的测定：土壤含水量采用烘干测定法；pH和电导率采用pH电导率仪；根据水土比2.5:1和5:1进行测定；有机质采用重铬酸钾容量法-外加热法测定；土壤全磷采用HClO₄-H₂SO₄法测定；土壤全氮采用半微量凯氏法测定；全钾含量测定采用NaOH熔融-火焰光度计法^[13]。

1.4.2 土壤微生物基因组DNA的提取：DNA的提取采用美国(MOBIO PowerSoil DNA Isolation Kit)强力土壤DNA提取试剂盒。

1.4.3 土壤细菌16S rRNA基因PCR扩增：采用细菌16S rRNA基因V3-V4区的通用引物(336F：5′-GTACTCCTACGGGAGGCAGCA-3′；806R：5′-GTGGACTACHVGGGTWTCTAAT-3′)扩增泥火山土壤细菌16S rRNA基因片段。50 μL PCR反应体系：DNA模板1 μL，正反向引物(10 μmol/L)各 2 μL，dNTPs (2.5 μmol/L) 4 μL，10×Pyrobest Buffer 5 μL，Pyrobest DNA Polymerase (2.5 U/μL) 0.3 μL，ddH₂O 35.7 μL。PCR程序为：95 ℃ 5 min；95 ℃ 30 s，56 ℃ 30 s，72 ℃ 40 s，共25个循环；72 ℃ 10 min。将PCR产物用1% (质量体积比)的琼脂糖凝胶进行检测，参照电泳初步定量结果，使用Qubit荧光定量系统对PCR产物进行检测定量，之后按照每个样本的测序量要求，根据精确定量结果进行相应比例的混合。随后进行MiSeq文库构建，Illumina上机测序。

通过COPE软件(Connecting Overlapped Pair-End，V1.2.3.3) ，对于Illumina测序的测序结果进行分析，去掉测序质量不好的序列，并对序列长度进行筛选，删掉长度小于200 bp的序列，根据Barcode将序列确定为本样品的序列，并去除Barcode和引物序列，最终得到有效的序列文件。在0.97相似度下利用Qiime (V1.8.0) 软件将其聚类为用于物种分类的OTU (Operational taxonomic units)，统计各个样品每个OTU中的丰度信息，OTU的丰度初步说明了样品的物种丰富程度。

利用QIIME软件计算样品包括Chao1^[14]指数、Shannon^[15]指数和Simpson指数^[16]的Alpha多样性值，其中Chao1指数表示群落丰富度，其值越高表示群落物种的丰富度越高；而Shannon指数和Simpson指数可以反映样品的多样性程度，其值越高表示群落物种的多样性越高；经过UniFrac算法利用系统进化的信息来比较样品间的物种群落差异，进行Beta多样性(Beta diversity)分析；通过SPSS 17.0^[17]软件进行细菌多样性与环境因子的相关性分析。

由表 2可知，不同生境土壤理化性质存在显著性差异(P＜0.05) ，pH、含水量、电导率、有机质、全磷、全氮和全钾的含量波动范围分别为7.41-10.20、2.72%-33.56%、0.12-7.97 ds/m、9.88-163.91 g/kg、0.33-1.06 g/kg、0.57-15.12 g/kg和3.24-6.27 g/kg。其中除了pH、电导率，在4种生境中，生境C的含水量、有机质、全磷、全氮和全钾含量最高；生境D的全氮、全磷、全钾含量相对较低；从时空动态变化方面看，pH、含水量、有机质(除了生境C)的变化趋势均是4月份最高；全磷的变化趋势均是11月份含量最高，7月份最低；全氮的变化趋势(除了生境A)是11月份含量最高；全钾的变化趋势均是4月份最高。

通过对泥火山土壤细菌V3-V4区测序，12个样品原始序列条数为302 439，过滤掉低质量的序列后，有效序列的总数为287 467。根据Barcode标签进行样品序列拆分，在对初始序列进行去冗余处理后，获得16S rRNA基因Unique Reads，并在97%相似度下将其聚类为用于物种分类的OTU (Operational taxonomic units)，统计得到各个样品在不同OTU中的丰度信息，12个样品共产生29 005个OTU，平均测序读长在420-480 bp之间。每个样品的有效序列数量如表 3所示。

稀释曲线(rarefaction)反映了样品的取样深度，可以用来评价测序量是否足以覆盖所有类群。从稀释曲线(图 1)中可知，12个样品稀释曲线有一些基本趋于平缓，例如样品C11，曲线已趋于平缓，即再增大数据量对OTU的发现也没有影响，可以说样本C11的OTU的覆盖度已达到饱和；而有一些仍未进入饱和状态，例如样品A4、C4，说明随着测序量的增加可能会发现新的细菌种系。

图 1　相似度为97%条件下各土壤样品的稀释曲线 Figure 1　Rarefaction curves of each soil sample at cutoff level of 3%

图选项

2.3.1 Alpha多样性分析：Alpha多样性指数是对单个样品中物种多样性的分析。由表 4可知，从生境类型来看，生境D的平均Chao1指数3 666.25和平均Shannon指数9.65最高；生境A的平均Simpson指数0.99最高；最低的是生境C，Simpson值0.86。从季节变化方面看，OTU数、丰度指数Chao1、多样性Shannon和Simpson指数均是7月份最高，依次是2 694、3 240.40、88.94、0.98。

2.3.2 Beta多样性分析：Beta多样性(Beta diversity)分析是用来比较一对样品在物种多样性方面存在的差异大小。UniFrac是通过系统进化的信息来比较样品间的物种群落差异，其计算结果可以作为一种衡量Beta多样性的指数，它考虑了物种间的进化距离，该指数越大表示样品间的差异越大。由图 2可知，基于UniFrac的加权和非加权的主坐标分析发现其第一主成分、第二主成分和第三主成分的贡献率分别为39%、21%和12%，20%、14%和12%，生境A、生境B、生境C、生境D之间明显有分离的现象，但又存在一定的交叉。这种现象说明，有些分类操作单元(OTU)在不同生境之间均存在，但含量不同。其中明显可以看出生境A的4月、7月、11月份的3个样品聚集在一起，说明其组成较为相似，也就是说生境A的微生物群落在季节动态变化方面并无太大差异。

图 2　Beta多样性的主坐标分析(PCoA)图 Figure 2　Principal coordinate analysis (PCoA) scores plot of beta diversity

图选项

2.3.3 细菌多样性与理化因子的相关性分析：从表 5可看出，Chao1指数、Shannon指数均与含水量呈极显著负相关关系，与全磷呈显著负相关关系；Simpson指数与含水量、有机质、全氮均呈极显著负相关关系，与全磷呈显著负相关关系。总体来看，土壤细菌多样性指数与理化因子之间呈负相关关系，且均与含水量、全磷相关性较高。

2.4.1 细菌门分类水平的比较：根据图 3，在门分类水平上共得到38种细菌类群，在这些样本中含量最高的门是变形菌门(Proteobacteria)，在4类生境中相对丰度分别为22.27%、45.58%、72.44%、32.74%，其次是放线菌门(Actinobacteria)相对丰度分别为44.52%、19.16%、41.12%、22.21%，拟杆菌门(Bacteroidetes)相对丰度分别为11.54%、7.98%、3.59%、7.73%，而其他菌类比例相对较低，例如芽单胞菌门(Gemmatimonadetes)、厚壁菌门(Firmicutes)、酸杆菌门(Acidobacteria)、绿弯菌门(Chloroflexi)、疣微菌门(Verrucomicrobia)等。这表明在泥火山土壤细菌中，变形菌门、放线菌门、拟杆菌门为三大优势菌群，而其他细菌类群为非优势菌群。

图 3　细菌门分类水平的比较 Figure 3　Comparison of bacteria groups at phylum level

图选项

2.4.2 细菌属分类水平的比较：根据图 4，在属分类水平上共得到72种细菌类群，其中含量最高的是未分类细菌，在4类生境中相对丰度分别为81.38%、74.34%、81.55%、64.12%；其次是Other、Halomonas (盐单胞菌属)、Kaistobacter、Rhodoplanes (红游动菌属)、Pseudomonas (假单胞菌属)、Streptococcus (链球菌属)、Rubrobacter (红色杆菌属)、Staphylococcus (葡萄球菌属)等，这些菌属的丰度比例相对较低。

图 4　细菌属分类水平的比较 Figure 4　Comparison of bacteria groups at genus level

图选项

2.5.1 OTU水平的PCA图：PCA (Principal component analysis) 分析，即主成分分析，是一种分析和简化数据集的技术。通过分析不同样品OTU (97%相似性)组成可以反映样品的差异和距离，如果两个样品距离越近，则表示这两个样品的组成越相似。从图 5可以看出，第一成分的差异贡献值是23.21%，第二成分的差异贡献值是15.35%，其中生境A的3个样品距离最近，说明生境A的OTU群落结构组成在季节动态方面是相似的，这与PCoA分析的结果一致；而生境D的这一组较为分散，说明在季节变化方面OTU组成差异较大，其中生境C和生境D距离较近，说明两者在一部分OTU上有相似的组成成分。

图 5　基于OTU丰度的PCA分析 Figure 5　Principal component analysis based on OTU level

图选项

2.5.2 丰度较高的前20个OTU heatmap图： Heatmap图通过颜色变化来反映二维矩阵或表格中的数据信息，将高丰度和低丰度的物种分块聚集，通过颜色梯度及相似程度来反映多个样品在各分类水平上群落组成的相似性和差异性。图 6是在属分类水平上12个样品的群落组成的Heatmap图，每小格代表所在样品中某种属的相对丰度。颜色越红代表相对丰度越高。由图 6可知，在属分类水平上，Shigella flexneri (X96963) (0.006%-64.470%)、OTU116 (0.005%-14.010%)、Halomonas (AY268080) (0.02%-10.90%)、OTU130 (0.022%-10.850%)、Pseudomonas (AB021404) (0.011%-5.170%)相对丰度较高。从群落组成上来看，生境A的7月份和11月份组成相似；生境C 4月份和生境D 11月份组成相似；生境B的4月和7月组成相似；生境B11月份和生境D的7月份组成相似，说明不同样本之间有交叉，且有相似之处。此外，可以明显看到heatmap图中出现大面积的白色方块，这是因为这些菌的含量非常低，导致没有数值。

图 6　基于丰度比较高的前20个OTU的Heatmap图 Figure 6　Heatmap based on the top 20 OTU abundance of bacteria genera

图选项

本研究利用Illumina MiSeq测序技术通过对16S rRNA基因的V3-V4高变区的测序分析，成功检测了泥火山土壤细菌群落结构的多样性，获得了大量、全面且深入的菌群信息。12个样品共获得有效序列的总数为287 467，平均测序读长在420-480 bp之间。在97%的相似度水平下共产生有效OTU个数为29 005，通过与RDP数据库^[18]进行比对，在门水平上得到优势菌群均为Proteobacteria、Actinobacteria、Bacteroidetes，其中Proteobacteria的丰度最高，这与Roesch等^[19]对土壤微生物和Zhang等^[20]对污泥微生物群落的研究相一致；在属水平上，含量最高的是未分类细菌，这说明泥火山样品中蕴含着丰富的微生物资源，很多物种是我们现在还未知的，而这些未知的物种为我们探究物种的多样性提供了很好的资源。此外，还有一些功能菌属被检测出，例如，Sphingomonas (鞘氨醇单胞菌属)是一类丰富的新型微生物资源，可用于芳香化合物的生物降解，该属菌株凭借自身的高代谢能力与多功能的生理特性，在环境保护及工业生产方面具有巨大的应用潜力，但是由于对鞘氨醇单胞菌的认识较晚，该菌的生态价值及经济价值很少被关注，对它的研究也停留在初级阶段^[21]。这一类菌属的发现可以为我们提供一个新思路，在后续研究中可以通过筛选新的培养基及分离程序，进一步研究可培养微生物及其功能，从而进一步发现有利用价值的资源。

MiSeq测序结果通过丰度指数(Chao1) 和多样性指数(Shannon、Simpson)反映了细菌的群落结构多样性，其中生境D的Shannon指数和Chao1指数最高，生境B和生境C相对较低；在土壤理化因子测定结果中，生境C的含水量、有机质、全磷、全氮和全钾含量最高，生境D的土壤养分(全氮、全磷、全钾)含量相对较低，并且通过细菌群落多样性与理化因子的相关性分析，结果表明泥火山细菌群落多样性随着土壤养分的增加而基本呈下降趋势，说明了物种多样性指数与理化因子之间呈负相关关系，这与陈芳^[22]的研究结果一致；群落结构多样性与含水量、全磷的相关性较高，说明泥火山微生物活动对磷循环有着很大的影响，并且磷在生物遗传中也扮演着重要的角色，是影响微生物多样性的一个重要因素，这也曾在高慧琴等^[23]的研究中提及过。

OTU水平的PCA分析表明，除了生境A的群落组成在时空演替方面没有差异，生境B和生境C的群落结构均受到了季节变化的影响，尤其是生境D，其3个样品较为分散，说明微生物群落结构的差异较大，群落组成明显不同，而生境C和生境D距离较近，表明两者在一部分OTU上有相似的组成成分，说明这两个生境之间群落组成有相似之处。而且，在Heatmap图中，也更进一步验证了PCA的结果，并发现Shigella flexneri、OTU116、Halomonas、OTU116、Pseudomonas的相对丰度较高。

传统测序技术(Sanger测序)由于其读长长、准确率高，能根据测得的全长16S rRNA基因序列注释到种，但Sanger测序通量较低，大规模测序成本高，操作复杂。Illumina MiSeq高通量测序技术，它能够产生大量Reads，一般每个样品可以保证测定40 000-60 000个序列数，覆盖深度非常大，对物种多样性的分析十分有利，并且在一般规模实验中，它的成本比率更易掌控^[24]。但是高通量测序读长短，不可能将16S rRNA基因的9个可变区全部测序，所以往往只选择1-3个可变区作为测序区域。自从Woese等将16S rRNA基因作为原核生物系统发育和分子生态学研究的经典标准以来^[25-26]，16S rRNA基因便成了细菌群落结构分析常用的系统进化标记分子。然而，由于发现该基因不同区域存在不同程度的异质性，不同可变区的选择对菌群具有偏向性，因而利用不同区域的序列进行高通量测序会造成不同程度的多样性高估，而相关研究^[27]指出，基因异质性趋向于集中在特殊的区域，V1和V6区域的高估程度最高，V4和V5区域则显示了最低的高估程度。因此，在本实验中，结合相关方面的文献选择了泥火山细菌16S rRNA基因的V3-V4可变区作为研究对象，从而降低了泥火山细菌多样性分析的高估程度，使得测序结果更加全面和准确。