基于高通量测序和形态学鉴定的马来西亚马塘红树林区浮游动物多样性分析

2022-01-10张金勇高养春AmirShahRuddinMdSahAileenTANShauHwaiNorlailaBintiMohdZanuri李宏俊

大连海洋大学学报 2021年6期

张金勇，高养春,Amir Shah Ruddin Md Sah, Aileen TAN Shau-Hwai, Norlaila Binti Mohd Zanuri，李宏俊*

(1.国家海洋环境监测中心海洋生态室，辽宁大连 116023； 2.School of Biological Sciences，Universiti Sains Malaysia，Penang 11800，Malaysia; 3. Centre For Marine & Coastal Studies，Universiti Sains Malaysia，Penang 11800，Malaysia)

浮游动物自主活动能力较弱，生命周期短，种类组成繁杂，数量大且分布广[1]。海洋浮游动物作为生态系统中重要的次级生产者，其群落结构的动态变化可以通过上行效应影响鱼、虾和贝类等海洋渔业资源的结构和总量，也可通过下行效应制约着浮游植物的初级生产力变化[2-3]。因此，浮游动物是海洋生态系统物质循环、能量流动的关键环节，承担着枢纽的作用，对整个海洋生态系统产生重要的调控作用[4]。

浮游动物种类的准确鉴定是研究浮游生态学的基础，而浮游动物种类多而复杂，传统依托形态特征的鉴定方法对鉴定人员的专业技能要求较高，并且不同鉴定结果可比性不佳；此外，采样方式、物种自身形态特征(遗传变异性、不同生活史阶段、性别限制和隐存种等)及样品保存都会给物种形态学鉴定工作带来巨大挑战[5]。因此，为满足日益增多的监测和评估需求，亟待建立一种速度快、精确度高、可自动化及标准化的物种鉴定检测方法。近年来，高通量测序DNA宏条形码技术的出现使得评估环境样品的生物多样性及复杂的群落结构成为可能[6]。在浮游动物DNA条形码研究中应用较多的分子标志基因包括核糖体28S rDNA、18S rDNA、内转录间隔区(ITS)和线粒体细胞色素氧化酶亚基Ⅰ(COⅠ)等[7]。选择分子标志基因时应区分其适用范围，统筹考虑研究目的和待测生物类群。18S rDNA存在可变区和保守区且序列进化速率较慢，通用引物易于覆盖较广物种类别的问题，适用于利用宏条形码研究海洋浮游生物的多样性[8-11]。

马塘红树林保护区(Matang Mangrove Reserve)位于马来西亚半岛西海岸的霹雳州(Negeri Perak)，面积超过40 466 hm2[12]。马来西亚政府早在1902年就对这片红树林进行了有效保护，距今已有超过百年的管理历史经验。红树林保护区作为陆地和海洋之间重要的过渡地带，为多种生物提供了觅食、繁殖与栖息的场所，是近海生物多样性的摇篮，是一个独特而复杂的生态系统[13]。浮游动物是红树林生态系统食物链和生产力的关键环节，在红树林生态系统结构、功能和生物资源等要素循环中起着重要作用。浮游动物群落结构特征能及时和准确地反映红树林生态系统的健康状况及水域生态环境的优劣[14]。本研究中,在马塘红树林保护区采集了10个站位的浮游动物样品，分别采取DNA宏条形码方法和传统上基于形态特征鉴定方法分析保护区内浮游动物群落结构特征，并将两种方法进行对比，以期为评价基于高通量测序的DNA宏条形码鉴定方法在海洋生物多样性业务化监测中的推广与应用提供参考。

1 材料与方法

1.1 调查站位的设置

2018年8月在马来西亚马塘红树林保护区布设10个采样站位(SO1～SO10)(图1)，各采样站位的具体经纬度及采样水深如表1所示。

图1 马来西亚马塘红树林保护区采样站位Fig.1 Location of sampling stations in Matang Mangrove Reserve, Malaysia

表1 采样站位经纬度及水深Tab.1 Latitude, longitude and water depth at sampling stations

1.2 方法

1.2.1 样品采集使用浅水Ⅱ型浮游生物网(网目孔径0.160 mm，网口内径31.6 cm，网长140 cm)，参照《海洋监测规范》(GB 17378.7—2007)采集浮游动物样品。将每个站位采集的样品进行充分混匀后均分成两份，分别进行传统形态学鉴定和分子技术鉴定。形态学样品用体积分数5%的福尔马林溶液固定保存；分子技术鉴定样品利用同孔径筛绢过滤掉海水后，将截留浮游动物的筛绢放入新的采样瓶，加入体积分数100%乙醇固定样品，用于DNA提取[15]。

1.2.2 样品种类鉴定形态学鉴定在体视显微镜下对浓缩后的浮游动物样品进行分类鉴定和计数。分子技术鉴定分3次从用乙醇保存的每个站位的样品(总体积100 mL)中吸取浮游动物样品，每次吸取10 mL，将3次吸取的样品均匀混合后用网目20 μm的筛绢过滤，去除样本中的残留浮游藻类及乙醇等[15]。

1.2.3 基因组DNA的提取使用DNeasy Blood and Tissue Kit (Qiagen, Hilden, 德国)提取样品中的总DNA，DNA提取完成后用NanoDrop 2000(Thermo Fisher Scientific, USA)检测其浓度，用10 g/L琼脂糖凝胶电泳检测DNA的完整性，检测合格的DNA置于-20 ℃冰箱中保存备用。

1.2.4 PCR扩增和Illumina测序以稀释后的基因组DNA为模板，使用浮游动物通用引物Uni18SF(5′AGGGCAAKYCTGGTGCCAGC 3′)和Uni18SR(5′GRCGGTATCTRATCGYCTT 3′)对18S rDNA 基因的V4 区进行扩增[6]。将引物添加8碱基的Barcode 标签用以区分不同样品，在上海美吉生物医药科技有限公司进行合成[16]。PCR反应体系(共20 μL)：5×FastPfu Buffer 4 μL，2.5 mmol/L dNTPs 2 μL，上、下游引物(5 μmol/L)各0.8 μL，FastPfu Polymerase 0.4 μL，DNA模板10 ng，用ddH2O补足至20 μL。使用ABI GeneAmp®9700型PCR系统进行扩增。PCR反应程序：95 ℃下预变性3 min；95 ℃下变性30 s，55 ℃下退火45 s，72 ℃下延伸30 s，共进行32个PCR循环；最后在72 ℃下再延伸10 min。PCR产物经纯化、浓度检测合格后，使用NEB Next UltraTM DNA Library Prep Kit for Illumina(New England Biolabs)进行文库构建，构建好的文库经过Qubit 定量和文库检测，使用MiSeq PE300 测序仪(Illumina Inc., San Diego, CA, USA)进行测序。

1.3 数据处理

测序得到的原始数据使用Flash 1.2.11进行数据拼接，参照Qiime(Quantitative insights into microbial ecology 1.9.1)[17]质量控制流程将拼接后的序列去重过滤、序列分类注释、Beta多样性距离计算。采用Uparse 7.0[18]对有效数据进行OTUs聚类(阈值设定相似度97%以上)，利用Usearch 8.1进行OTUs统计处理，将期望误差阈值(EER)设为1.0。使用RDP Classifier 2.11方法与NCBI Nr数据库对OTUs代表序列进行序列分类注释，在此基础上再对OTUs进行丰富度、α多样性计算等，同时对物种注释在各分类水平上进行群落结构的组成分析。浮游动物的形态学鉴定需要依据门、目和属进行专业技术人员鉴定。应用Mothur 1.30.2构建样品稀释性曲线，使用R语言Vegan包完成热点图(Heatmap图)分析，热点图是以颜色梯度来表征二维矩阵或表格中的数据大小，并呈现群落物种组成及物种的丰富度信息。将OTU相对丰富度数据进行标准化处理后，基于Bray-Curtis距离矩阵，进行CLUSTER聚类分析。

采用Coverage指数(C)估算样本文库的覆盖率，采用Shannon指数(H)估算样本中物种多样性水平，采用丰富度指数(SChao1)评价群落中物种的丰富程度，采用Pielou均匀度指数(E)度量群落中相对物种丰富度，采用Ace指数(SAce)估计群落中OTU数目的指数，该指数在生态学中同样作为度量物种丰富度的指标，计算公式[19-22]为

C=1-n1/N，

(1)

(2)

(3)

E=H/Hmax,

(4)

(5)

式中，

其中：n1为只含有一条序列的OTU数目；N为抽样中出现的总序列数；Sobs为实际观测到的OTU数目；ni为第i个OTU所含的序列数；SChao1为估计的OTU数；n2为只含有两条序列的OTU数目；H为Shannon指数；Hmax为在物种丰富度相同的情况下，能够达到的最大Shannon指数(即当群落中所有物种丰度完全一致时)；Srare为含有“abund”条序列或者少于“abund”的OTU数目；Sabund为多于“abund”条序列的OTU数目。

2 结果与分析

2.1 高通量测序浮游动物稀释性曲线

18S V4区引物测序结果显示，对样本的原始测序数据经过拼接和质量过滤后，总共获得571 865条序列，每个站点样品的有效序列片段数量范围为46 335～71 624条，平均为57 186.5条序列。经去除重复序列、OTU聚类(97%的阈值)等步骤，所有站位样品总共获得249个OTUs。根据OTUs数目所做的样品稀释曲线(图2)可以看出，所有站位样品OTUs数目都趋于饱和，表明各样品测序数据量足够覆盖大部分的种类。

图2 各采样站位浮游动物样品OTUs数目稀释性曲线Fig.2 Rarefaction curves of the number of OTUs in zooplankton samples at each sampling station

2.2 基于高通量测序的浮游动物多样性指数

由表2可看出，各样品文库的覆盖率较高，表明样本序列检出率高，测序结果真实可靠。基于DNA宏条形码分子鉴定方法计算得到的各站位样品α多样性指数显示，SO4、SO1、SO2的Shannon指数值处于前3位，SO3和SO5的Shannon指数值处于后2位；SO4、SO1、SO2的Ace指数值也处于前3位，SO7和SO10的Ace指数值处于后2位；SO2、SO4和SO1的Chao1指数值处于前3位，SO7和SO9的Chao1指数值处于后2位；SO9、SO4和SO10的Pielou均匀度指数值处于前3位，SO3的Pielou均匀度指数值最低(0.27)(表2)。

表2 基于宏条形码鉴定方法的浮游动物多样性指数统计Tab.2 Zooplankton diversity index statistics based metabarcoding identification methods

2.3 基于高通量测序的浮游动物组成分析

热点图可以在门、目和属水平上呈现样本之间的关系。根据门和目分类水平上热点图样本聚类树可以看出，10组样本大致可划分为3类(图3(a), 图4(a)),第Ⅰ类包括SO1、SO2、SO3和SO6，其样本间丰度组成的相似性较高，第Ⅱ类是样本间较为相似的SO4、SO7、SO8、SO9和SO10，而第Ⅲ类SO5与其他样本差别最大，单独为一类。在属分类水平上，根据热点图样本聚类树依然可将10组样本划分为3类，但站位SO4从门和目分类水平的第Ⅱ类站位划归到第Ⅰ类站位中，即第Ⅰ类包括SO1、SO2、SO3、SO4和SO6(图5(a))。

图3 各站位浮游动物样品门水平上组成分析Fig.3 Analysis of the composition of zooplankton samples at each sampling station at phylum level

图4 各站位浮游动物样品目水平上组成分析Fig.4 Analysis of the composition of zooplankton samples at each sampling station at order level

图5 各站位浮游动物样品属水平上组成分析Fig.5 Analysis of the composition of zooplankton samples at each sampling station at genus level

群落柱形图可以在门、目和属水平上直观展示各样本的浮游动物组成。在门分类水平上(图3(b))，第Ⅰ类站位中刺胞动物门Cnidaria、棘皮动物门Echinodermata和节肢动物门Arthropoda的平均reads数较多，处于前3位，相对丰度分别占总丰度的10.64%、7.50%和6.69%；第Ⅱ类站位(SO4、SO7、SO8、SO9和SO10)中平均reads数占比最高的是节肢动物门，达58.74%，除节肢动物门外，Ⅱ类站位reads数占比较高的门类有刺胞动物门(6.11%)和环节动物门(5.78%)；第Ⅲ类SO5中半索动物门Hemichordata、节肢动物门和刺胞动物门的平均reads数所占比例处于前3位，相对丰度分别占64.44%、27.18%和0.91%。

在目分类阶元上(图4(b))，第Ⅰ类站位中心形海胆目Spatangoida、哲水蚤目Calanoida和硬水母目Trachymedusae的平均reads数占比处于前3位，相对丰度分别占7.48%、6.20%和3.66%；第Ⅱ类所有站位中哲水蚤目平均reads数占比最高，相对丰度占比达到51.21%；第Ⅲ类SO5中哲水蚤目平均reads数占比亦是最高，相对丰度达20.33%。

在属分类阶元上(图5(b))，3大类站位浮游动物的相对丰度也存在较大差异，第Ⅰ类站位中石笔海胆属Archaeopneustes和球水母属Sphaeronectes的平均reads数较高，相对丰度分别占6.35%和3.05%；第Ⅱ类站位中矮隆哲水蚤属Bestiolina的平均reads数占比最高，相对丰度占比达到47.40%；而第Ⅲ类SO5中腺头虫属Glandiceps的平均reads数占比最高，相对丰度占比达到64.44%。

2.4 基于高通量测序和形态学特征鉴定得到的浮游动物组成比较

借助光学显微镜，在各站位样品中共检测出的浮游动物包括哲水蚤目、剑水蚤目Cyclopoidea、猛水蚤目Harpacticoida、十足目Decapoda、甲壳类幼虫Nauplii及鱼卵等。基于形态特征鉴定的各站位浮游动物丰度平均值为129.8 ind./m3(不包括鱼卵)，范围为12～307 ind./m3，站位SO5平均丰度最高(307 ind./m3)，站位SO6平均丰度最低(12 ind./m3)(表3)。

基于形态学特征鉴定的结果显示，SO10、SO2和SO5的Shannon指数处于前3位，SO6的Shannon指数值最低(0.28)；SO10、SO9和SO6的Pielou均匀度指数值处于前3位，SO4的值Pielou均匀度指数最低(0.60)(表4)。

表3 基于形态学特征鉴定的马塘红树林保护区浮游动物丰度分布Tab.3 Distribution of abundance of zooplankton identified by morphological characteristics in Matang Mangrove Reserve

表4 基于形态学鉴定方法的浮游动物多样性指数统计Tab.4 Zooplankton diversity index statistics based morphological identification methods

将采用DNA宏条形码方法获得的OTUs序列与NCBI Nr数据库进行物种比对，成功注释浮游动物25门88目138属；基于形态学鉴定，共鉴定出浮游动物2门6目13属。本研究发现，利用形态学方法鉴定出来的大部分类别(门100%、目83%、属54%)同时也能用宏条形码方法鉴定出来；而利用宏条形码分子方法鉴定出的大部分类别(门92%、目94%、属95%)却未能利用形态学方法鉴定出来。此外，宏条形码分子方法鉴定结果显示，矮隆哲水蚤属的reads数在Ⅱ类站位SO7、SO8、SO9和SO10中占比都是各自采样点最高的，相对丰度占比分别达到40.73%、64.95%、36.98%和46.92%；而形态学方法鉴定结果显示，丰度最高的属是大眼剑水蚤属Corycaeus，达到24.07%(表3)。

对于所有采样站位，基于宏条形码方法鉴定的平均Pielou均匀度指数(0.45±0.13)低于基于形态学鉴定的平均Pielou均匀度指数(0.86±0.10)，但在站位SO4中发现基于分子技术鉴定的Pielou均匀度指数略高于基于形态学鉴定的结果；同时，基于宏条形码鉴定方法的多样性指数均高于形态学鉴定方法得出的多样性指数(表2,表4)。

3 讨论

3.1 DNA宏条形码法监测浮游动物群落的可行性

DNA宏条形码技术是结合DNA条形码和高通量测序的方法，这一技术的应用对浮游生物多样性的研究提供了强大的技术支撑[23]。但利用DNA宏条形码进行物种鉴定选择合适的标志基因是基本保障。由于标志基因对不同物种的通用性和分辨率存在差异，高灵敏度的通用引物仍然是DNA条形码技术发展的限制之一。目前，大部分浮游动物研究仍基于COⅠ、18S rDNA和28S rDNA基因。为探究浮游动物DNA宏条形码标志基因的扩增差异，Zhan等[24]比较研究了COⅠ、16S和18S等标志基因后发现，传统COⅠ基因难以提供稳定的PCR产物，建议选择18S和16S基因进行研究。然而，高旭等[25]的研究结果发现，虽然18S V9对浮游动物具有更高的物种覆盖度，16S具有更好的特异性，但COⅠ的物种覆盖度、物种识别敏感性和物种特异性都适中，更加适合浮游动物DNA 宏条形码监测。

本研究中选取的Uni18S引物已广泛应用于海洋及淡水生物多样性研究与评价[15]。虽然此引物并不能将所有物种鉴定到种的水平，但对大部分种类而言鉴定到科及属水平具有较高可信度[6]。利用此引物，本研究中鉴定出的浮游动物类别在门、目及属分类水平上(25门88目138属)明显优于形态学鉴定出的类别(2门6目13属)，且宏条形码分子鉴定的多样性指数(Shannon多样性指数)明显高于基于形态学特征鉴定出的结果，平均Pielou均匀度指数(0.45±0.13)低于形态学鉴定的平均Pielou均匀度指数(0.86±0.1)。但在浮游动物多样性研究与评价中，对于同一海域不同采样站位应采取相同的浮游动物鉴定方法，多样性指数比较亦是如此。同时，根据课题组前期工作，高养春等[15]估算了两种方法鉴定每份浮游动物样本所需的时间成本和费用成本(表5)，结果发现，与形态学鉴定方法相比，宏条形码分子鉴定具有鉴定效率高、费用低的优点。这些研究结果表明，DNA宏条形码分子鉴定方法在浮游生物多样性研究与评价中具有较高的应用潜力。

表5 基于两种方法鉴定每份浮游动物样本所需费用和时间成本Tab.5 Economical and time costs of analyzing one zooplankton sample based on both methods

3.2 DNA宏条形码方法监测浮游动物群落面临的问题

尽管DNA宏条形码技术的优势(快速、操作简便、人为干扰小等)已经被大家公认[10，26]，但该技术本身也还存在一些问题需要在未来浮游动物鉴定应用中加以解决。这些问题主要集中在3个方面，即如何选择合适的引物与引物偏好性问题、绝对定量问题和参考数据库不完整性问题。首先，在浮游动物研究中，有很多标志基因可以选择，但是没有哪一个标志基因能够同时扩增出所有物种并保证每一个物种的扩增效率是一致的。因此，针对不同的生物组成如何选择合适的标志基因是DNA宏条形码技术能否用于环境评价的限制之一[27]。其次，在PCR扩增及高通量测序过程中均会产生假阳性及假阴性的错误，OTUs的序列数可以近似反映生物量的多少，由于各种因素的干扰，很可能导致生物量和OTUs序列数之间并不是线性相关的，例如DNA的提取方法、PCR扩增效率、DNA混合池甚至是生物信息学分析都会干扰定量计算[28]。最后，利用DNA宏条形码技术进行物种鉴定最终要回归到测序序列和参考数据库的比对上。因此，参考数据库的完整性和质量直接决定了DNA宏条形码技术获得数据的可靠性。虽然近年来参考数据库中物种信息增长很快，GenBank、Barcode of Life (BOLD)等公共数据库收录了大量数据，但相对陆生生态系统，对水生生态系统收录的类群较少。数据库的不完善导致很多宏条形码数据没有办法进行有效注释，处于未知状态[29]。

3.3 宏条形码方法与形态学方法监测马塘红树林保护区浮游动物群落结构差异性

近年来，已有国内研究人员对福建漳江口、雷州半岛及海南东寨港红树林自然保护区大型、小型底栖动物进行过相关研究，有关浮游动物的研究仅见黄勃课题组对东寨港红树林区的分析和研究[13-14]。而早在20世纪80年代初，国外就有关于红树林林区浮游动物的研究报道。目前，相关研究主要集中在盐度、水温等理化因子对浮游动物分布的影响等方面，而有关红树林区生态系统中浮游动物群落结构的研究却较少[30-31]。本研究中，采用宏条形码和形态学鉴定两种技术方法对马塘红树林区10个站位浮游动物样品进行浮游动物群落组成和结构分析，结果表明，宏条形码分子方法能得到更多的浮游动物种类和丰富度信息，基于形态学特征鉴定出的浮游动物在门、目和属分类水平上明显少于分子鉴定出的类别，且发现了大量基于形态学特征无法鉴定出的种类。由此结果可以推知，基于18S rDNA V4区扩增的宏条形码方法具有较高的物种分辨率和通用性，鉴定效率大大高于形态学鉴定方法。另外，两种方法鉴定到的丰度最高的属在各站位是不同的，宏条形码方法检测出丰度最高的属为矮隆哲水蚤属，且各站位矮隆哲水蚤属所占的比例沿采样点SO1至SO10有逐渐升高的趋势。基于形态学特征鉴定出丰度最高的属是大眼剑水蚤属，且在大部分采样点该属所占的比例均较高，少部分采样点(如SO2和SO8)矮隆哲水蚤属所占的比例较高。出现此结果，一方面由于形态学鉴定对人员的专业技能要求较高，不同鉴定人员的结果可比性较差；另一方面限于目前宏条形码分子方法存在的一些技术问题(引物偏好性问题、绝对定量问题和参考数据库不完整性问题)，导致部分站位宏条形码数据未能有效鉴别。黄进等[32]采用宏条形码分子方法结合显微镜计数的方法研究了鸭河口水库微型浮游生物群落组成，比较两种方法，宏条形码分子方法能检测微型后生动物Metazoa、双鞭毛虫Dinoflagellata、丝足虫Cercozoa等8个门，而形态学鉴定方法在浮游动物检测中仅观察到轮虫动物门Rotifera 1个门。施军琼等[33]对真核浮游生物18S rDNA V4区进行高通量测序，探讨了大宁河不同水华期真核浮游生物群落组成，结果表明，V4区引物在微微型、微型浮游植物群落鉴定方面更为高效，能有效弥补光学显微镜观察中丢失的一些种类信息。宏条形码分子方法能检测到更丰富的浮游生物群落信息，国外的研究中也普遍发现类似的结果[34-35]。

通过对比高通量测序的宏条形码方法和传统光学显微镜鉴定数据发现，高通量测序不但具有鉴定效率高、费用较低等优点，而且能得到更多的浮游动物种类信息。随着宏条形码技术体系的不断完善，数据库不断更新，使得DNA条形码技术的准确性更高，将为海洋浮游生物多样性快速评估提供强大的技术支撑。