氨氧化古菌(Ammonia-oxidizing archaea，AOA)是自然环境中一类具有将氨氧化成亚硝酸盐功能的古菌，在氮素的生物地球化学循环中起着重要的作用^[1]。2005年，Könneke等从西雅图热带海洋水族馆的海水中成功分离培养到第一株具有氨氧化功能的古菌Nitrosopumilus martimus SCM1，开辟了微生物氨氧化作用研究的新领域，自此氨氧化古菌作为新的氮循环微生物成为研究热点^[2]。2008年，科学家们通过系统发生学和酶学分析，将包括氨氧化古菌在内的部分泉古菌类群重新划分成一个新门——奇古菌门(Thaumarchaeota)^[3-4]。amoA基因是编码氨氧化过程中的关键酶——氨单加氧酶(Ammonia monooxygenase，AMO) α亚基的基因，对AOA的系统发育分析与16S rRNA基因有着高度的一致性^[5]，与氨氧化细菌(Ammonia- oxidizing bacteria，AOB)基因组中amoA基因存在多拷贝不同^[6]，目前已发现的AOA基因组中amoA基因均为单拷贝。目前，amoA基因作为分子标记物广泛地用于海洋、河口、沉积物、土壤等生境中氨氧化微生物的生态学研究^[7-9]。相较于保守的16S rRNA基因，amoA基因序列差异程度更大，作为分子标记分析氨氧化古菌多样性分辨率更高，可以更好地区分16S rRNA基因相似性较高的不同种。此外，由于AOA amoA基因为单拷贝，并不会出现类似16S rRNA基因由于多拷贝和基因组内部异质性导致的对多样性分析的偏差^[10]。但amoA基因作为分子标记分析氨氧化古菌多样性的过程中存在两大难点：(1)目前没有基于amoA基因序列构建的参考数据库(Taxonomic database)；(2)进行操作分类单元(Operational taxonomic unit，OTU)聚类分析时，由于AOA amoA基因与16S rRNA基因在进化速率上存在差异，AOA amoA基因在种水平上的相似性阈值并不清楚，amoA基因并不能类似于16S rRNA基因，将97%的序列相似度作为划分种水平OTU (Species-level OTU)的阈值^[11]，即对amoA基因测序结果聚类分析过程中没有明确地划分种水平OTU的阈值。

本文以南海氨氧化古菌的多样性分析为例，以amoA功能基因作为分子标记，基于高通量测序技术，利用MOTHUR软件^[12]进行数据分析，建立基于amoA基因高通量测序分析环境中氨氧化古菌多样性的研究方法。

选取南海S28站位(北纬18°00′00″，东经113°30′00″)，通过船载CTD (Conductivity- Temperature-Depth sensors)上配备的尼斯金采样瓶(Niskin bottle)采集不同深度水层水体(75、100、150、200、500、800、1 000 m)，量取2 L海水样品经0.22 μm聚碳酸酯滤膜(Millipore，USA)进行过滤，收集海洋水体中的微生物。将过滤后的滤膜装入冻存管中，置于液氮中冷冻保存。

PowerWater^® DNA Isolation Kit，Qiagen公司；GoTaq^® Colorless Master Mix，Promega公司；TaKaRa MiniBEST Agarose Gel DNA Extraction Kit，宝生物工程(大连)有限公司；PCR引物由Invitrogen公司合成。PCR仪，ABI公司；NanoDrop LITE核酸定量仪，Thermo公司；Mupid^®-One电泳仪，Advance公司；凝胶成像系统，Bio-Rad公司。

按照强力水样DNA提取试剂盒操作说明书提取不同深度水层样本滤膜上的微生物宏基因组DNA，使用特异性引物Arch-amoAF (5′-STAATGG TCTGGCTTAGACG-3′)和Arch-amoAR (5′-GCGG CCATCCATCTGTATGT-3′)扩增amoA基因，其中上游引物前端通过加8碱基的Barcode标记进行修饰，使用带有不同Barcode标记的上游引物扩增不同样品，PCR反应体系及反应条件参照Francis等^[13]。利用凝胶成像系统的定量工具初步测定PCR产物的相对丰度，根据各样本目标测序深度对各样本PCR产物进行相应比例混合，利用胶回收试剂盒对混合后的PCR产物进行切胶纯化回收。最后，将纯化后的PCR产物−80 ℃保存，委托苏州金唯智生物科技有限公司进行高通量测序。

从GenBank数据库(https://www.ncbi.nlm.nih.gov/genbank/)中下载目前已通过分离或富集培养获得菌株的amoA基因序列作为种子序列，同时下载FunGene数据库^[14] (http://fungene.cme.msu.edu/)中收集的环境样品amoA基因序列，将获得的所有序列进行整合，用MOTHUR软件(https://mothur.org/wiki/MiSeq_SOP)对序列进行筛选：(1)以种子序列作为参比序列进行序列对齐；(2)对环境样品序列进行去冗余，提取非冗余序列；(3)对序列进行长度控制，删除长度小于555个碱基的序列；(4)对剩余序列进行嵌合体检测，去除含有嵌合体片段的序列。利用筛选处理后的序列生成基因序列文件(Template file)，基于amoA基因的序列号从NCBI数据库中搜集各序列对应的种属信息，生成种属信息文件(Taxonomic file)。基因序列文件和对应的种属信息文件构成功能基因的参考数据库^[15]。

收集NCBI数据库(https://www.ncbi.nlm.nih.gov/)中amoA基因以及16S rRNA基因序列信息均完整的氨氧化古菌菌株，利用MEGA 7.0分别计算菌株间amoA基因以及16S rRNA基因的相似度。以16S rRNA基因相似度作为X轴，功能基因amoA基因相似度作为Y轴，用GraphPad Prism7.00进行散点图绘制，通过回归分析计算菌株间2个基因相似度的相关关系。基于回归方程得出97%的16S rRNA基因相似度所对应的amoA基因相似度，从而得到amoA基因在种水平上的相似性阈值(Species-level cut-off value)。

利用MOTHUR (V.1.39.5)软件对原始序列进行分析处理，包括序列质量控制和相关生物信息学分析。对初始序列进行质控过滤，通过长度筛选以及嵌合体检测，去掉过短序列和嵌合体序列(Chimera sequence)，序列对齐后，去除Barcode和引物序列，得到有效序列用于进一步分析。通过与参考数据库进行比对，对有效序列进行物种注释，绘制种群结构组分图(属水平)，分析南海不同深度水层中氨氧化古菌的种属组成；以2.5中确定的相似性阈值进行OTU聚类分析，将相似度处于阈值范围内的序列归为同一个OTU；去除稀有OTU (序列数 < 总序列数的0.1%)，绘制稀释曲线(Rarefaction curve)进行取样深度验证，评价测序量是否足以覆盖样品中的所有类群；选取主要OTU (序列数 > 总序列数的1%)的代表序列构建系统发育树，分析样品中优势物种进化关系及分类地位，同时基于主要OTU的丰度信息绘制热图(Heatmap)，反映不同样本间物种分布的相似性与差异性；基于非稀有OTU (序列数 > 总序列数的0.1%)在每个样品中丰度分布状况(格式为Shared的文件)，计算样品中的Chao1指数^[16]、Shannon指数^[17]等多样性指数，从而分析样品中的α多样性(Alpha-diversity)；基于主坐标分析(Principal coordinates analysis，PCoA)计算样品间的多样性差异，进行β多样性(Beta-diversity)分析。

将得到的AOA amoA基因高通量测序结果提交至NCBI SRA (Sequence read archive)数据库获取登录号。

构建的序列种子库中，共收集了迄今通过分离或富集培养获得的34株氨氧化古菌的amoA基因序列信息，基于种子序列构建的系统发育树见图 1。可以看出，目前已知的氨氧化古菌可分为4大簇，包括海洋簇(Group I.1a)、土壤簇(Group I.1b)、嗜酸簇(Group I.1a-associated)、嗜热簇(THAOA)。从FunGene数据库中共下载到47 792条环境样品序列。经筛选处理并整合后，得到AOA amoA基因序列文件(Template file)。表 1简要总结了参考数据库中序列的质量信息，基因序列文件中共收集了26 091条amoA基因序列，所有序列中均不含模糊碱基，但含有3−5个碱基长度的同聚物，序列长度控制在555−635个碱基的范围。基因序列文件与种属信息文件均由2个纵列组成，基因序列文件第1列为序列号，第2列为碱基序列；种属信息文件第1列与基因序列文件的第1列一致，第2列为对应的分类信息，其中不同分类级别用分号分隔。基因序列文件和种属信息文件构成amoA基因参考数据库。

图 1　已知的氨氧化古菌amoA基因系统发育分析 Figure 1　Phylogenetic analysis of recognized species of archaeal amoA genes 注：基于邻接法，节点上的数值为自展值(Bootstrap)，低于50%的数值(1 000次重复)没有显示，序列选自NCBI数据库，括号内为含有amoA基因的序列登录号. Note: Based on neighbor-joining method, numbers at the node are the bootstrap values (%), bootstrap values (1 000 replicates) less than 50% are not shown. Sequences are obtained from NCBI database, accession numbers of sequences of amoA gene are shown in parentheses.

图选项

通过菌株间的两两比对，得到成对菌株间16S rRNA基因相似度和amoA基因相似度(Pairwise distances)，从回归分析结果(图 2)可以看出，这两者之间有着很好的相关性，线性拟合R²达到0.84，而多项式拟合R²达到0.91。表明16S rRNA基因与amoA基因在进化过程中有很高的相关性，但16S rRNA基因比amoA基因在进化速率上表现得更保守。选择R²值更高的多项式拟合曲线作为16S rRNA基因与amoA基因相关关系的回归方程可以得到，16S rRNA基因97%的序列相似度对应的amoA功能基因的相似度为89.4%。因此，基于功能基因amoA序列对氨氧化古菌的多样性分析，选取以89%的序列相似度作为在种水平上划分OTU的阈值。

图 2　amoA基因相似度与16S rRNA基因序列相似度相关性分析 Figure 2　Correlation of amoA gene similarity versus 16S rRNA gene similarity 注：黑色直线为线性拟合，红色曲线为二项式拟合. Note: Black line represents linear fitting, red curve represents binominal fitting.

图选项

AOA amoA基因高通量测序结果在NCBI SRA数据库获取的登录号为SRP131845。

运行MOTHUR软件中的classify.seqs语句，基于构建的参考数据库对质量控制后的有效序列进行分类学注释。利用运行结束后生成的分类结果统计文件，计算氨氧化古菌属水平下的相对丰度分布情况，图 3为基于OriginPro 2017绘制的不同深度水层AOA在属水平上的相对丰度分布图。从图 3可以看出，随着深度变化，南海水体中的AOA菌群分布有明显的变化趋势：南海75 m和100 m的浅水层中的AOA主要属于Nitrosopelagicus属；500、800和1 000 m的深水层样品序列匹配到的参考序列主要为未培养(Uncultured)的环境样品序列，未明确分类地位菌属主要分布在深层水体；而150 m和200 m中层水体的氨氧化古菌种属分布表现出一个明显的过渡段。此外，测序产生的所有序列均来自于奇古菌门的amoA基因，说明本研究选择的扩增引物对AOA的amoA基因具有很强的特异性。而分析结果中仅有很小一部分序列(3.7%)被划分到了Unclassified，即仅有3.7%的序列无法从参考数据库中找到匹配序列，说明构建的参考数据库收集的序列信息较为完整，库容量足够支撑环境中氨氧化古菌的多样性分析。

图 3　南海不同深度水层氨氧化古菌的属水平分布 Figure 3　Relative abundance of the AOA at genus level in different depth in the South China Sea

图选项

利用MOTHUR平台上的dist.seqs命令语句计算有效序列间的相似度，生成距离矩阵，将相似度≥89%的序列归为一个OTU。挑取主要OTU的代表序列，得到了8个主要OTU的代表序列，利用MEGA 7.0进行系统进化树的构建(图 4)，分析主要OTU间的系统进化关系；基于OTU在每个样品的丰度分布及进化关系进行OTU热图(Heatmap)绘制(图 5)，通过颜色梯度变化及相似程度来反映不同深度水层中AOA种属组成的相似性和差异性。系统发育分析结果表明，所有样品序列均属于奇古菌中的Group I.1a部分，其中丰度最高的OTU01与从海洋中通过富集培养得到的Candidatus Nitrosopelagicus brevis CN25^[18]有很近的亲缘关系，序列经两两比对(BLAST)后相似度为99%。结合图 4和图 5来看，样品序列被划分为两大簇，其中包含有OTU01、OTU04、OTU07的一簇主要分布在表层及中层(75−200 m)，在底层水体中丰度极低；而另一簇(OTU02、OTU03、OTU05、OTU06、OTU08)主要分布在中层和底层(150−1 000 m)，在表层水体中几乎没有分布。氨氧化古菌在南海海洋水体中表现出很高的垂直分布差异。

图 4　基于主要OTU构建的系统进化树 Figure 4　Phylogenetic tree of dominant OTUs in the South China Sea 注：基于邻接法，Bootstrap值选自1 000重复，黑点标志为参考序列，参考序列选自NCBI数据库，圆括号里为参考序列的登录号，方括号内为各主要OTU序列数占总序列数的比例. Note: Based on neighbor-joining method, numbers on branches are the supporting percentage by 1 000 replicates. Reference sequences obtained from NCBI database are marked by black circle, accession number of reference sequences are shown in parenthesis, the percentage of each dominant OTU in total sequences are shown in bracket.

图选项

图 5　基于主要OTU绘制的热图 Figure 5　Heat map of dominant OTUs in the South China Sea 注：热图左侧为样品名称和样本聚类树，上面和下面分别为物种聚类树和主要OTU名称. Note: The left side of the heat map are the sample name and the sample clustering tree, the upper and the bottom side are the species clustering tree and dominant OTU name, respectively.

图选项

通过随机抽样的方法提取数据，以抽取的序列数作为横轴，以这些序列所包含的OTU数目作为纵轴，绘制稀释曲线(Rarefaction curve)。稀释曲线反映了样品的取样深度，可以用来评价测序量是否足以覆盖所有类群。南海不同深度水层样品的稀释曲线如图 6所示，可以看出稀释曲线已趋于平缓，表明当前的测序量已经能够很好地反映出样品中AOA的群落结构和多样性，测序数据量合理，再增大数据量对发现新OTU的影响很小。

图 6　操作分类单元稀释曲线 Figure 6　Rarefaction curves for OTUs

图选项

利用Mothur软件对南海不同深度水层样品进行多样性指数(Alpha-diversity)统计分析(表 2)，测序结果覆盖率均在99.55%以上，表明数据有效可靠。从整体上来看，150 m和200 m的样品中OTU数目最多，Shannon指数最高，Simpson指数最低；而75 m和100 m的样品中OTU数目最少，Shannon指数最低，Simpson指数最高。这表明南海氨氧化古菌在表层水体中的物种丰富度、种群多样性以及物种组成均匀度最低，而在中层水体中表现最高，深层水体次之。

基于OTU聚类结果，利用MOTHUR软件中的thetayc算法计算样品间的差异，运行pcoa语句进行主坐标分析(Principal coordinates analysis)，结果见图 7。由图 7可知，主坐标分析其第1主成分和第2主成分的贡献率分别为68.64%和24.85%，其中，75 m和100 m距离较近，而代表 500、800和1 000 m样品的点很好地聚集在一起，代表样品的点按深度聚集成3簇。由此可以得出，南海氨氧化古菌在浅层(75、100 m)、中层(150、200 m)、深层(500、800、1 000 m)的海洋水体中分布差异明显，氨氧化古菌在南海水体的分布呈现垂直分层现象。

图 7　南海不同深度水层样品中氨氧化古菌种群结构主坐标分析 Figure 7　Principal coordinates analysis for AOA community in different depth in the South China Sea

图选项

基于MOTHUR软件平台，利用构建的amoA基因参考数据库以及将89%的相似度作为在种水平上划分OTU的阈值，可以有效分析环境样品中的氨氧化古菌多样性。对南海不同深度水层水体样品中amoA基因高通量测序结果的分析表明，本文建立的分析方法可以有效揭示南海水体中氨氧化古菌的垂直分布差异。

本研究构建了含有26 091条序列的参考数据库，其中大部分amoA基因序列来源于未培养微生物，仅是通过分子生物学的手段从环境样品中获得，其对应的种属信息和分类地位都不明确。因此，已鉴定物种过少以及分类信息不完整仍是目前针对环境样品中氨氧化古菌多样性进行分析的制约因素。2012年Biller等分析了NCBI中来自不同生境的8 296条AOA amoA基因序列，按照85%的相似度在种水平上进行物种划分时自然环境中AOA的物种数在138个左右，当按照90%的相似度进行物种划分时AOA的物种数为315个^[19]。而由于氨氧化古菌生长缓慢，代时长，目前通过分离或富集培养得到并鉴定的AOA仅有二十几个物种，这就意味着自然界中将近85%−95%的AOA还未被人类所认识，因此对氨氧化古菌新物种的分离与鉴定仍然是目前微生物氨氧化研究的重点。

参照Purkhold等^[20]确定的氨氧化细菌(AOB) amoA基因在种水平上相似性阈值的方法，Pester等^[21]基于当时16S rRNA基因与amoA基因序列信息均完整的10株氨氧化古菌菌株，提出将amoA基因序列相似度大于85%的AOA定义为同一个物种。本研究中增加了近年新鉴定的氨氧化古菌菌株，利用两两比对后的561组数据重新计算了AOA amoA基因在种水平上的相似性阈值，将该值从85%提升至89%。而随着分离鉴定的氨氧化古菌菌株的增加，数据量更加丰富，阈值的选取也将更加真实准确。基于目前已发现的氨氧化古菌的全基因组，并未发现AOA基因组中16S rRNA基因与amoA基因存在多拷贝的现象，因此本文对阈值的选取过程均是基于基因组中单拷贝的16S rRNA基因与单拷贝的amoA基因比对得出。但有研究发现，古菌基因组中16S rRNA基因存在多拷贝现象，而且多拷贝的基因序列并不完全一致，即16S rRNA基因在基因组内部存在异质性。16S rRNA基因的多拷贝及异质性对基于16S rRNA的多样性分析会造成分析结果的偏差，会引起一定程度的高估^[10]。但对于功能基因相似性阈值的选取，由于绝大多数原核微生物基因组中多拷贝的16S rRNA基因序列通常是相同或是相似度非常高的(序列差异小于1%)^[22]，而对AOA amoA基因的相似性阈值的确定是基于16S rRNA基因相似度与amoA基因相似度在统计学上表现出的相关关系得到的，因此16S rRNA基因的多拷贝及异质性对相似性阈值的选取影响较小。

目前已发展了多种技术手段用于环境微生物组研究，包括基因芯片、高通量测序、宏基因组测序等。高通量测序技术以其低成本、高通量、流程自动化的优点备受青睐，功能基因的高通量测序是继16S rDNA/18S rDNA/ITS等特定区域片段PCR产物测序后的又一项高靶向性测序技术。选择对功能基因区域覆盖度和特异性最佳的扩增引物，通过PCR扩增技术获得功能基因扩增产物，利用功能基因的碱基差异来反映环境中功能微生物物种之间的差异。功能基因具有的高特异性与高分辨率，使其作为分子标记不仅能更加准确地靶定环境微生物群中具有特定功能的类群微生物，也能更加清晰地反映菌群中特定类群微生物的物种组成及多样性；不仅可以区分生态功能不同但亲缘相近的微生物种群，也可以鉴定功能相同但亲缘关系较远的微生物种群。此外，随着测序技术的迅猛发展，测序成本的不断降低，尤其第三代测序技术的推广，让更长的功能基因用于功能微生物多样性分析成为了现实。因此，采用功能基因的高通量测序技术，利用本文中提出的分析方法，可以针对自然环境中的特定功能微生物进行有效的多样性分析。