DNA复合条形码技术(metabarcoding)是一种结合了DNA条形码技术和高通量测序技术的快速鉴定混合样本中物种的方法^[1-3]。DNA条形码技术通过一段短的DNA序列来鉴定物种，由Hebert等人于2003年提出^[4]。传统的Sanger测序只能对单一的DNA分子进行测序，而高通量测序技术可以分别对混合DNA中的每一条DNA分子进行测序。两种技术的结合运用可以获得混合样本中大量的含有分类学信息的DNA片段的序列，然后将这些序列与合适的数据库比对即可确定其代表的类群。随着高通量测序技术的持续发展与分析软件的不断完善，在分析生物多样性时，DNA复合条形码技术可以作为代替传统方法的一个选择^[2-3]。

DNA复合条形码技术已普遍使用于微生物群落的研究，并且正越来越多的应用到动植物多样性研究中^{[1, 5]}。在土壤动物方面：Douglas等人构建了7个节肢动物群落，通过对这些节肢动物群落的DNA复合条形码研究证明了该方法可以精确估计群落内的系统发育多样性和群落间的差异性^[2]；Bienert等人对2个不同地点的土壤样品进行了研究，发现这一方法能很好地区分这2个地点的蚯蚓群落组成，并且与使用传统方法得到的结果一致^[6]；Yang等人对基于该技术分析土壤小型动物多样性的3种数据处理流程进行了分析比较，结果表明3种流程的处理结果基本一致^[7]；魏攀利用该技术对小五台山国家级自然保护区1m²草地中无脊椎动物多样性进行了研究，共调查得到125种无脊椎动物^[8]；Yang等人使用该方法对采自土壤和落叶层的样本与采自地上的样本进行了比较，结果表明，在研究生物多样性格局时，两类样本可以得到相似的生态学信息^[9]。此外，还有许多关于土壤线虫的研究^[10-12]。

目前, 科学界已经认识到土壤是地球上生物多样性最丰富的生境之一^[13]。虽然土壤生物多样性的研究已有很长时间，但是，关于土壤生物多样性本身的很多基础科学问题仍然不清楚。从某种意义上说，对土壤生物多样性的研究还处于初期阶段^[14]。

我国对土壤动物的研究主要集中在东北、亚热带地区、热带雨林地区以及新疆地区，对于中西部地区的土壤动物研究还很少^[15-16]。太白山为秦岭主峰，位于动物古北界与东洋界的交汇和过渡地带，是华北、华中和横断山脉3个植物区系的交汇处，自然植被带完整，生态系统复杂多样，生物资源丰富，在海拔3400m以上存有完整的第四纪古冰川遗迹，是一个保存完好的自然生态系统，是我国生物多样性保护的关键区域^[17]。

本文选取18S核糖体小亚基RNA(18S rRNA)基因和细胞色素C氧化酶亚基I(COI)基因的部分序列作为分子标记，应用DNA复合条形码技术对秦岭太白山地区5种不同生境(高山草甸，针叶林，落叶小叶林，落叶阔叶林，农田)的中小型土壤动物多样进行研究，并讨论该方法目前应用于中小型土壤动物研究的优缺，为进一步探讨秦岭中小型土壤动物多样性格局以及秦岭生物多样性保护提供科学依据。

1 研究方法 1.1 样本的采集和分离

于2014年6月上旬在太白山自然保护区及周围选取高山草甸、针叶林、落叶小叶林、落叶阔叶林和农田5种生境进行取样。每种生境选取10m×10m的3个样方，在每个样方内随机采取3份土样，采样深度为0—15cm。将土样带回实验室，用Tullgren漏斗进行分离(40W白炽灯分离48 h)。将收集到的土壤动物保存在100%乙醇中，用于总DNA的提取。样地详细信息见表 1。

1.2 总DNA提取

将一个样方分离到的所有土壤动物混合放入一个1.5mL离心管，使用血液/细胞/组织基因组DNA提取试剂盒(天根生化科技有限公司)进行总DNA提取，得到15份总DNA。将每种生境的3份总DNA进行等浓度混合，得到用于后续PCR的5份总DNA。

1.3 PCR及测序

PCR反应体系为25μL，包括2×Taq PCR StarMix with Loading Dye 12.5μL，灭菌超纯水9.5μL，正向引物(10μmol/L)1.0μL，反向引物(10μmol/L)1.0μL，总DNA 1.0μL。

18S rRNA基因部分序列(大约350bp)的扩增使用正向引物#3：5′-GYGGTGCATGGCCGTTSKTRGTT-3′和反向引物#5_RC：5′-GTGTGYACAAAGGBCAGGGAC-3′^[18]。为了同时对多个样本进行测序，在每个样本的正向引物5′端加上了7bp的barcode序列，用于后续分析中区分序列来自哪个样本。PCR的循环程序为：95℃10min；10个循环：95℃10s，57℃30s(每个循环降低0.5℃)，72℃60s；30个循环：95℃10s，55℃30s，75℃60s；最后在72℃下延长5min。

COI基因部分序列(大约313bp)的扩增使用正向引物mlCOlintF：5′-GGWACWGGWTGAACWGT WTAYCCYCC-3′^[19]和反向引物jgHCO2198：5′-TANACYTCNGGRTGNCCRAARAAYCA-3′^[20](用N替换了原序列中的I)。同样在每个样本的正向引物5′端加上了7bp的barcode序列。PCR的循环程序为：16个循环：95℃10s，62℃30s(每个循环降低1℃)，72℃60s；25个循环：95℃10s，46℃30s，72℃30s；最后在72℃下延长5min。

每份总DNA扩增3次，然后将3份PCR产物等浓度混合。PCR产物的测序使用Illumina MiSeq平台，采用250PE测序策略，由上海派森诺生物科技有限公司完成。

1.4 数据处理 1.4.1 原始数据处理

步长为1的5bp窗口从序列起始位置开始滑动，要求窗口中碱基平均质量≥Q20，低于该值则截断序列。然后使用FLASH(http://ccb.jhu.edu/software/FLASH/)对通过质量过滤的序列进行拼接。使用Qiime^[21]对拼接后的序列进一步过滤，舍弃含有N碱基、长度< 150bp、barcode有错配和引物错配碱基数> 1的序列。使用mothur^[22]去除嵌合体序列，得到用于后续分析的序列。

1.4.2 OTU聚类及注释

使用Qiime将序列按相似度0.97聚类为OTU(Operational Taxonomic Unit)^[23]，选取每个OTU最长的序列作为代表序列。将18S rRNA基因和COI基因的OTU代表序列分别与Silva^[24]数据库(115版)和BOLD^[25]数据库(截止到Release 5.00-v1)进行比对，获得每个OTU的分类学信息。使用两个标记的OTU列表生成物种列表，然后合并，删去不属于土壤动物的类群。使用Excel绘制目水平上的群落组成百分比堆积图。使用R软件VennDiagram程序包绘制Venn图。

1.4.3 α多样性分析

使用R软件vegan程序包对各样本的Chao1丰富度指数，ACE丰富度指数，Shannon-Wiener多样性指数和Simpson多样性指数进行计算。

1.4.4 群落相似性分析

使用R软件vegan程序包对各样本之间的Bray-Curtis距离进行计算，然后使用ape程序包基于Bray-Curtis距离对样本进行主坐标分析(Principal Coordinates Analysis; PCoA)。

1.4.5 聚类分析

根据物种列表，使用R软件pheatmap程序包将科水平上的分类信息进行聚类分析并绘制Heatmap图。

2 结果 2.1 原始数据处理与OTU聚类

经最初的质量过滤与拼接后得到611959序列，其中18S rRNA基因309167条，COI基因302792条。经Qiime与mothur进一步过滤后得到用于后续分析的序列，共437313条，其中18S rRNA基因211436条，COI基因225877条。合格序列聚类后，18S rRNA基因共得到1653个OTU，COI基因共得到2138个OTU。

2.2 种类组成

本文的土壤动物分类参照《中国土壤动物检索图鉴》^[26]。调查共鉴定到土壤动物3门9纲28目199科，其中18S rRNA基因鉴定到2门8纲23目101科，COI基因鉴定到3门7纲16目128科，两者共同鉴定到2门6纲12目31科。高山草甸包含21目136科，针叶林包含20目141科，落叶小叶林包含21目122科，落叶阔叶林包含21目138科，农田包含23目132科。各样本所包含的类群及各类群在各样本中所含序列条数见附表。

对各样本做Venn图(图 1)可得，各样本在科水平共有的类群数为84，在高山草甸、针叶林、落叶小叶林、落叶阔叶林、农田各自总类群数中所占比例分别为61.76%、59.57%、68.85%、60.87%、63.64%。高山草甸独有的类群数为9，占其总类群数6.62%；针叶林独有的类群数为13，占其总类群数9.22%；落叶小叶林独有的类群数为7，占其总类群数5.74%；落叶阔叶林独有的类群数为10，占其总类群数7.25%；农田独有的类群数为7，占其总类群数5.30%。

如图 2所示，在目水平上，真螨目、弹尾目和蜘蛛目的丰度在各个样本中都较高，其中蜘蛛目在落叶小叶林中所占比例达64.62%。鳞翅目在高山草甸、针叶林和农田中的丰度较高，而在落叶小叶林和阔叶林中丰度较低。半翅目在高山草甸、针叶林和落叶小叶林中丰度较高，而在落叶阔叶林和农田中丰度较低。双尾目在落叶阔叶林中所占比例达61.04%，而在其他生境中含量较低。条马陆目在针叶林中含量较高，鞘翅目在农田中含量较高。其余类群在各生境中含量均较低。

2.3 α多样性分析

各生境的土壤动物群落Chao1丰富度指数为：针叶林>落叶小叶林>落叶阔叶林>高山草甸>农田；ACE丰富度指数为：针叶林>落叶阔叶林>高山草甸>落叶小叶林>农田；Shannon-Wiener多样性指数为：针叶林>高山草甸>农田>落叶阔叶林>落叶小叶林；Simpson多样性指数为：针叶林>高山草甸>农田>落叶阔叶林>落叶小叶林。各生境多样性指数见表 2。

2.4 群落相似性分析

Bray-Curtis距离可以作为衡量群落间相似性的指标，Bray-Curtis距离越小，表示两个样本的群落组成越相似。表 3为各样本之间的Bray-Curtis距离。根据各样本之间的Bray-Curtis距离，可以做PCoA分析。样本的相似性越高，它们在PCoA图中的距离就越近。图 3为根据各样本间的Bray-Curtis距离，选取贡献率最大的主坐标组合所做的PCoA图。图 3显示土壤动物群落组成在高山草甸、针叶林、农田3种生境中的相似性较高，落叶小叶林和落叶阔叶林的土壤动物群落组成与这3种生境的相似性较低，落叶小叶林与落叶阔叶林之间的土壤动物群落组成相似性也较低。

2.5 聚类分析

Heatmap图可以直观地将数据值的大小以颜色的深浅表示出来，同时可根据需要将数据进行物种或样本间的相似性聚类。图 4为使用欧式距离和最长聚类法所作的Heatmap图，展现了科水平各类群在各样本中的丰度分布，同时也在科水平对各样本进行了聚类。可以看到，高山草甸与针叶林首先聚类，然后与农田聚类，然后与落叶阔叶林聚类，最后与落叶小叶林聚类。

3 结论与讨论

从类群列表来看，两个分子标记鉴定到的类群没有完全重叠，两者共同鉴定到的类群只有31科。而18S rRNA基因单独鉴定到70科，COI基因单独鉴定到97科，远大于它们共同鉴定到的类群数。即使是在目水平上，两者鉴定到的类群也有很大的区别。这是由于不同分子标记的通用引物在不同类群中的扩增能力不同与数据库的覆盖度不同造成的^{[18-19, 27]}。因此单独使用一种分子标记来鉴定样本中含有的类群会造成对类群数的低估，同时使用多种分子标记可以更全面的鉴定样本中所含类群^{[10, 28]}。考虑到数据库的可用性，对于后生动物来说，较好的分子标记有COI、18S rRNA和28S rRNA基因^{[18, 25, 29-30]}。

尽管采样地点的海拔跨度高达2500m以上，生境类型发生了剧烈的变化，但从各样本群落组成来看，在目水平，真螨目、弹尾目和蜘蛛目在各生境中所占的比例都较大，在科水平，5种生境的共有类群达84种，在各生境类群数中所占比例大多超过60%(在针叶林中为59.57%)，这表明太白山中小型土壤动物中大多数类群为广布类群。另外，各生境均含有类群数低于10%的特有类群，表明海拔与生境的变化会对太白山中小型土壤动物的群落组成造成一定的影响。这与以往的研究结果一致^[31-33]。

在α多样性分析中，虽然Chao1指数与ACE指数同为丰富度指数，但它们的排序结果却有细微的差别，落叶小叶林中的土壤动物群落Chao1指数排在落叶阔叶林和高山草甸之前，而ACE指数却排在了落叶阔叶林和高山草甸之后。这是因为ACE指数的计算公式中加入了一个定义“优势”类群的阈值，这个阈值在vegan程序包中为类群包含10条序列，而Chao1指数并没有这个规则^[34]。因此可以推断落叶小叶林土壤动物群落中包含较多的低丰度类群(所含序列数< 10)。在DNA复合条形码技术中，Shannon-Wiener和Simpson多样性指数是基于一个类群所含序列数计算的，而不是基于一个类群的个体数。Carew等人^[28]指出，在一个样本中，一个物种的序列数和个体数是正相关的，但序列的频率并不能完美的反映物种的频率。Piñol等人^[35]的研究也表明该技术的定量潜力有限，只有定性的研究是可靠的。因此，Shannon-Wiener和Simpson多样性指数仅能作为一个参考，不能通过这两个指数很好的判断土壤动物群落的多样性。定量问题是目前该技术用于土壤动物研究的局限性。

从各样本之间的Bray-Curtis距离与PCoA图来看，针叶林与高山草甸的土壤动物群落组成最为相似，针叶林与农田次之。针叶林、高山草甸与农田三者的土壤动物群落组成相似性较高，而落叶小叶林和落叶阔叶林的土壤动物群落组成与这3种生境的相似性都较低，落叶小叶林和落叶阔叶林的土壤动物群落组成相似性也较低。聚类分析得到了与PCoA分析相符的结果。这种相似性的形成可能与各生境中的气温和土壤的理化条件相关^[32]，有待后续的实验。

根据以往的土壤动物研究可知，在一次调查中会得到成千上万的土壤动物，且土壤动物包含的类群很多，鉴定这些土壤动物需要多人长时间作业。对于其它的环境样本，也存在同样的问题。而通过DNA复合条形码技术，则可以快速简便的完成大量样本的分类鉴定，与传统方法相比，省时省力，这也是该技术的优势所在。随着高通量测序技术的不断发展和分析方法的不断完善，DNA复合条形码技术的高成本效益和全面性将使更容易的大尺度监测生物多样性趋势成为可能^[36-37]。