DNA复合条形码技术(DNA metabarcoding)是在DNA条形码基础上结合高通量测序技术和分子分类学技术兴起的一种用于生物多样性检测的方法。该技术既可以在短时间内对大量样本的物种组成进行识别,也可以对大批量的多物种混合环境样本进行快速的DNA条形码鉴定,避免了单物种DNA条形码鉴定的低效率^[1-3]。随着高通量测序技术的进步、条形码参考数据库的完善和分析软件的开发,DNA复合条形码技术已经成为生物多样性检测最常用的方法。

复合条形码技术从诞生以来广泛应用于多方面的研究,如Yang等采用18S rDNA作为标记基因分析了土壤小型动物的多样性^[4]；Hajibabaei等对格兰德河水生生物多样性进行了研究^[5]；Nolte等比较了奥地利一个湖中不同季节原生动物多样性^[6]；Porazinska等比较了哥斯达黎加不同热带雨林微生境中线虫类的生物多样性^[7]。此外,DNA复合条形码在微生物和真菌多样性^[8-9]、动物排泄物以及食性结构^[10]和古细菌^[11]的研究上都有很广泛的应用。

自然界内江河、湖泊、水库众多,水生生物种类繁多,其中水生动物是最主要的类群,包括浮游动物和底栖动物。浮游动物包括原生动物、轮虫、枝角类和桡足类等,底栖动物包括扁形动物门、环节动物门、节肢动物门等^[12]。由于广泛的水域面积和复杂的水生环境,关于水生动物多样性的研究,虽然有很多发现,但仍然有一些科学问题亟待解决。秦岭是横亘于我国中部东西走向的巨大山脉,是我国的南北分界线^[13]。其内部山谷众多,大小溪流贯穿其中,水生生态环境复杂,但目前还缺乏对水生动物区系和群落的研究。因此,本文选取COⅠ和18S rRNA基因作为条形码标记,用DNA复合条形码技术对秦岭水生动物区系组成和群落结构组成进行了初步研究,同时用α和β多样性以及聚类分析来比较不同生境和取样环境下,水生动物多样性组成的差别,为秦岭生物多样性的保护和研究提供数据支持。

1 材料与方法 1.1 样品采集与准备

从秦岭中选择金龙峡、石砭峪、五台山、子午峪和沣峪口5个淡水水域,分上下游取样,在上下游分别设置3个采样点,每个采样点用50mL冻存管取30mL水样,同时用毛笔和镊子在采样点的石头上或底层采集肉眼可见底栖动物,放入50mL冻存管中。采样地点和环境如表 1所示。

将采得的水样带回实验室用真空抽滤机浓缩,一个样地3个取样点的水样浓缩为一个浮游动物水样。然后将浓缩水样和底栖动物用100%乙醇保存,放入-20℃的冰箱,用于提取总DNA。

1.2 总DNA的提取

将浓缩水样和底栖动物,混合放入一个经高压灭菌的1.5mL离心管中,用DNeasy Tissue and blood kit(Qiagen Inc)试剂盒提取总DNA,得到10个DNA样品,编号如表 1所示。

1.3 PCR扩增以及测序

根据Hajibabaei等^[5]和Nolte等^[6]类似研究采用的引物,实验扩增COⅠ序列的正向引物为LepF1：5′-ATTCAACCAATCATAAAGATATTGG-3′,反向引物为EPT-long-univR： 5′-AARAAAATYATAAYAAANGCGTGN ANNGT-3′,扩增长度约为130bp；扩增18S rRNA序列的正向引物为fw：5′-ATTAGGGTTCGATTCCGGAGAGG-3′,反向引物为rv：5′-CTGGAATTACCGCGGSTGCTG-3′,扩增长度约为180-200bp。为了能够同时对多个样品进行测序,以及在后续分析中区分每个样品序列,在每个样品的5′端都加了一个7bp的标签序列。

实验中PCR体系为25μL,包括12.5μL 的2×Taq PCR StarMix with Loading Dye,1.0μL的DNA模版,上游引物(10μmol/L)和下游引物(10μmol/L)各1μL,9.5μL的双蒸水。COⅠ基因的扩增循环程序是：95℃ 5min预热,35个循环：95℃ 40s,43.5℃ 1min,72℃ 30s,72℃ 5min延伸,4℃保温；18S rRNA基因的扩增循环程序是：98℃ 2min预热,35个循环：98℃ 30s,65℃ 45s,72℃ 90s,72℃ 7min延伸,4℃保温。

PCR产物送往上海派森诺生物科技有限公司进行高通量测序,采用Illumina MiSeq平台,进行双末端测序。最后将得到的测序数据返回,进行数据处理和生物多样性分析。

2 数据处理和分析 2.1 原始序列处理：质量过滤和嵌合体的去除

首先,将返回的双末端Fataq序列用fastqc软件对序列文件进行质量检测；然后用trimmomatic 0.33软件进行质量过滤,要求步长为1的5bp窗口从第1个碱基开始移动,窗口中碱基的平均质量≥Q20(即碱基准确率为99%),从第1个低于Q20处截断序列,最终要求序列长度≥150bp,且不容许有模糊碱基；然后利用Flash软件对质量过滤的序列进行拼接,要求两个序列数据的最大重叠片段≥10bp,且不容许有碱基错配。

由于PCR扩增可能会产生嵌合体,测序过程中会产生点突变等测序错误,因此还要对拼接后的序列进一步的过滤和去除嵌合体。实验用Qiime^[14]软件进行序列过滤,用mothur^[15]软件去除嵌合体序列。

2.2 OTU分类和注释

实验用Qiime软件中的uclust^[16]对所得到的优质序列按照0.97的相似度进行OTU聚类,选取每个OTU中最长序列为代表序列；用Qiime软件中的blast,将得到的代表序列与数据库比对,获得OTU分类学信息。其中18S rRNA基因与Silva^[17]数据库比对,COⅠ基因与BOLD数据库比对。

然后,根据注释得到的精简OTU列表,用excel和mothur软件做目水平上的群落组成图和样品间venn图。

2.3 群落多样性分析

根据物种丰度,用mothur软件中的summary.single命令,求出每个样品的Chao,ACE,Shannon和Simpson指数。其中Chao和ACE指数是群落丰富度(Community richness)指数,二者的值越大,说明群落丰富度越高；Shannon和Simpson指数是群落多样性(Community diversity)指数,Simpson指数越大,说明群落多样性越低,Shannon指数越大,说明群落多样性越高。

用Qiime软件,根据各样品的物种进化和丰度信息,进行Unifrac分析^[18-19],得到样品间差异距离矩阵,然后进行PCoA分析。

2.4 聚类分析

根据精简后的OTU列表,用R软件中的pheatmap程序包进行属水平上的聚类分析,并绘制出heatmap^[20-21]图。

3 结果 3.1 测序结果

原始序列经过滤和嵌合体的去除,两个标记基因共得到1120240条有效序列,971319条优质序列,其中18S rRNA的有效序列有767441条,优质序列有640700条,而COⅠ的有效序列有352799条,优质序列有330619条,它们的优质序列所占比例分别为：83.49%和93.71%。

3.2 区系组成和群落结构

将优质序列与数据库比对,共得到水生动物10门48纲89目,其中18S rRNA鉴定了9门42纲52目,COⅠ鉴定了5门11纲36目。其中18S rRNA鉴定的各样品组成情况是：J1包含31目55科,J2包含23目41科,D1包含28目49科,D2包含24目44科,W1包含27目45科,W2包含25目41科,Z1包含23目47科,Z2包含20目46科,F1包含20目35科,F2包含19目30科。

根据OTU列表做venn图,如图 1和图 2所示,COⅠ鉴定得到：J1与J2的共有物种数是218,特有物种数分别为212和152；D1与D2的共有物种数是236,特有物种数分别为95和116；W1与W2的共有物种数是88,特有物种数分别为214和170；Z1与Z2的共有物种数是122,特有物种数分别为153和208；F1与F2的共有物种数是265,特有物种数分别为80和123；所有取样点下游与上游的共有物种数是628,特有物种数分别是272和155；D、F、J、W和Z的总物种数分别是：683、733、800、560和605；18S rRNA鉴定得到：J1与J2的共有物种数是154,特有物种数分别为65和47；D1与D2的共有物种数是158,特有物种数分别为74和28；W1与W2的共有物种数是131,特有物种数分别为56和67；Z1与Z2的共有物种数是160,特有物种数分别为49和68；F1与F2的共有物种数是139,特有物种数分别为43和53；所有取样点下游与上游的共有物种数是305,特有物种数分别是11和10；J、D、W、Z和F的共有物种数是140,特有物种数分别为3、2、9、2和0,D、F、J、W和Z的总物种数分别是：418、374、420、375和437。

根据OTU列表,可以用样品在各分类水平(门、纲、目、科、属)上的物种组成比例情况,来反应样品在不同分类水平上的群落结构。COⅠ从目水平上分析结果如图 3所示：D1、D2、F1、F2、J1、J2和W1中双翅目的含量相对较高,所占比例分别为45.2%、46.3%、46.9%、27.1%、83.1%、53.9%和48.9%；W2中毛翅目的含量较高,所占比例为56.6%；Z1中基眼目的含量较高,所占比例为29.3%；Z2中鞘翅目的含量相对较高,所占比例为45.3%。18S rRNA从门水平上分析结果如图 4所示：D1、D2、J1、J2、W1、Z1和Z2中节肢动物门的含量相对较高,所占比例分别为59.2%、84.5%、92.7%、65.8%、91.6%、53.5%和66.9%；F1和F2中软体动物门的含量相对较高,所占比例分别为55.2和58.6%；W2中扁形动物门的含量相对较高,所占比例为54.5%。可见,COⅠ调查表明：样本中含量相对较高的类群是双翅目、毛翅目、基眼目、鞘翅目、蜉蝣目等；而18S rRNA调查表明：样本中主要有节肢动物门、软体动物门和扁形动物门三个大门。

3.3 α多样性分析

COⅠ基因的Chao和ACE指数表明J1、J2和F2的群落丰富度相对较高,Z1、W2和Z2的则相对较低；Simpson和Shannon指数表明J1、Z2和F2的群落多样性相对较高,而J2、D1和W2的则相对较低。18S rRNA基因的Chao和ACE指数表明J1、D1和J2的群落丰富度相对较高,而F1、F2和W1的则相对较低,Simpson和Shannon指数表明J1、Z1和W2的群落多样性相对较高,而F2、F1、和D2的则相对较低。整体来看,J、F、D的群落丰富度和群落多样性相对较高,Z、W的则相对较低。

3.4 β多样性分析

用样品间物种进化和丰度信息,进行Unifrac分析,得到样品间差异距离矩阵,然后用差异距离矩阵进行PCoA分析,得到PCoA图,如图 5所示,图中两点之间的距离越近说明两者之间的生物群落差异性越小。18S rRNA基因的PCoA图显示：Z1、D2和J1,J2和D1之间的生物群落结构差异性较小；COⅠ基因的PCoA图显示：D1、Z1、D2和J2,Z2和W2之间的生物群落差异性较小。

3.5 聚类分析

18S rRNA基因的聚类情况显示Z2、J2和D2,F2和F1,W1和W2,Z1和D1在属水平上群落相似性相对较高(图 6)；COⅠ基因的聚类情况显示J2、D2、D1和Z2,J1和W1,F1和F2,Z1和W2在属水平上的群落相似性相对较高(图 7)。

4 分析与讨论 4.1 环境因素与结果讨论

由于水环境的复杂性和多变性,使得水生生物多样性受多方面影响,其可以分为自然因素和人为因素两个方面^[12]。自然因素主要包含上下游、水体深浅以及流动性等,而人为因素主要是指人类活动对生物多样性的影响。本文主要从以上这些方面对秦岭水域的水生动物多样性做了初步调查。

对于水流来说,水的流动性会让上游的营养物质向下游流动,使得下游的营养高于上游,从而导致下游的物种丰富度高于上游^[22]。从venn图可以看到,各样地下游的物种丰度普遍高于上游；群落结构图显示,各样地下游的群落结构比上游的复杂,并且下游群落结构中优势物种的含量也高于上游；从α多样性指数可以看到,虽然个别样品上游的群落丰富度和群落多样性较下游的高,但绝大多数样品下游的群落丰富度和群落多样性还是较上游的高。可见,本次调查也表明下游的物种丰富度要高于上游。

实验选择了秦岭中5个淡水水域：沣峪口、五台山、子午峪、金龙峡和石砭峪,其中沣峪口,金龙峡和石砭峪在景区内,人为干扰较大,而子午峪和五台山相对来说人为干扰较小,基本上处于自然状态。从实验结果可以看出各样地之间的不同：如各样地物种丰度的大小顺序是：W＜Z＜D＜F＜J,表明受人为影响较大的J、F和D 3个取样地的物种丰度大于受人为影响较小的W、Z两个取样地；群落结构图显示,D1、D2、F1、F2、J1和J2的主要类群是双翅目,W1和W2的主要类群是双翅目和毛翅目,Z1和Z2的主要类群是基眼目和鞘翅目；α多样性指数表明,J、F和D的群落丰富度和群落多样性都相对较高,而Z和W的则相对较低。这些都表明一个样地的生物多样性会受到人为活动的影响。“放牧效应”^[23-24]在一般生态系统中都是适合的,一定程度的外来干扰会对一个地方的生物多样性有明显的提高,如本实验中有人为干扰的3个取样地的物种丰度、群落丰富度和群落多样性都比自然状态下的两个取样地要高。但是当外来干扰超过一定程度之后,就会使该地的生物多样性急剧下降,对原来的生态系统造成严重的破坏。

对一个生态系统来说,影响生物多样性的因素是多方面的,除了前面所述的上下游和人为活动的影响外,本实验还关注了水的深浅、流动性以及阳光的有无对生物多样性的影响。如取自深水的D2和F1样品与取自浅水的D1和F2样品,跟其他都取自浅水的样品相比,有相对较多的共有物种； PCoA图分析表明不同环境因素下样品的群落结构差异性较大,而相似环境因素下样品的群落结构差异性则相对较小,如J2和D2都是取自较阴暗的上游流动水,取样环境相似,因此它们之间的群落结构差异性相对较小。此外,Heatmap图也显示,环境相似性较高的样品首先聚集在一起,并且群落相似性相对较高,如Z1和D1都是取自较浅的下游静水,取样环境相似,因此它们首先聚类在一起,表明它们之间的群落相似性相对较高。

4.2 实验方法与结果讨论

实验结果显示,两个分子标记所得的实验结果并不完全一致。这主要是由于每个基因的通用引物在不同类群中的扩增能力和在数据库中的覆盖范围不同导致的。对基于高通量测序的复合条形码技术来说,由于样品的复杂性和多重性,在研究时如果选择多个分子标记基因可能会得到更全面的实验结果,如本论文中18S rRNA基因鉴定了9门42纲52目,COⅠ基因鉴定了5门11纲36目,而两个分子标记基因共同鉴定了10门48纲89目。Carew^[25]等人的研究结果也表明,使用两个或两个以上的分子标记可以更全面的鉴定到样本中所包含的物种。这使得一些依赖于类群数计算的各种生物多样性指数,将变得更加准确。但是在使用多个分子标记进行分析时,由于每个分子标记在进化速率和单独鉴定的类群上不同,可能使得构建的系统发育树之间有区别,从而影响了依赖于系统发育树计算的unifrac分析结果。

目前,应用DNA复合条形码技术对生物多样性的研究越来越多,但是关于复合条形码技术的一些本质问题仍然没有得到很好的解决。其中一个最主要的问题就是：复合条形码技术只能做定性分析,而不能很准确的进行定量分析。在研究中,一直用序列数的多少来表示物种个体数量的多少,虽然有些研究认为序列数与物种的个体数之间存在相关性,但也仅仅是它们的变化趋势相对一致,而并不是很完美的契合。因此,在分析的过程中,只能用序列数大致的估计物种的个体数,而不能完全用序列数来代替物种的个体数。Piñol^[26]等人的研究也表明复合条形码技术的定量研究能力有限,只有定性分析的结果才是可靠的。对于本研究来说,有些分析会不可避免的应用到物种个体数的多少,如群落结构图中各样本中类群的含量,以及α多样性中的Shannon指数等,所以对于这些结果,只能定性来看,如果想要得到更加精确的结果,还须进一步的验证。当然,在本实验中也避免了一些用物种个体数的分析,如没有使用Bray-Curis指数进行β多样性的分析,而是使用了依赖于物种数和进化信息的jaccard和unweighted_unifrac指数进行β多样性和后续的分析。但是,生态学研究中的群落组成和群落结构以及群落之间的关系,都是依赖于物种的个体数量,所以DNA复合条形码技术还不能很准确的应用到生态学研究中的每一个方面。

虽然复合条形码技术还不能更精确的进行生物多样性研究,但相比于传统耗时耗力以及单一的物种鉴定和生物多样性研究,复合条形码技术有明显的优势：它能在短时间对大量复合样本进行物种鉴定和生物多样性研究,既省时又省力。随着测序技术的进步和提高,条形码数据库的全面和完善,生物信息学软件的开发和使用,复合条形码技术将会成为生物多样性研究的主要方法,并且准确性也会大幅度提高,更重要的是它将能够进行更大尺度和更大规模的生物多样性研究。