崆峒岛海域(Kongtong Dao, KTD), 位于山东省烟台市芝罘区北部, 属于北黄海；东港海域(Donggang, DG)位于山东省日照市东港区南部, 属于南黄海。黄海是一个半封闭的陆架浅海, 以山东半岛的成山头至朝鲜的长山串之间的连线为界划分成南、北两部分(Liu et al, 1998)。其中, 北黄海位于山东半岛、辽东半岛和朝鲜半岛之间, 西有经渤海海峡流入的黄河物质, 南北有南沿岸流、鲁北沿岸流、西朝鲜沿岸流等沿岸流系汇入(杨青等, 2012)；南黄海位于山东半岛成山头与朝鲜半岛长山串连线以南、长江口北岸启东嘴与济州岛连线以北, 除了最南部的长江汇入以外, 西与山东半岛南部、江苏北部相邻, 有灌河、射阳河、淮河等河流入海(王俊杰等, 2018)。由于两处海域特殊的地理位置和环境影响, 使得其生态环境较为复杂, 是研究浮游生物种类及群落结构的热点海区。

真核浮游生物主要包括浮游植物和浮游动物。其中, 浮游植物是海洋生态系统中的初级生产者, 是水域生产力的基础(左涛等, 2008; 张婷等, 2014)。浮游动物作为海洋生态系统的次级生产者, 是经济水产动物, 特别是上、中层鱼类和一些幼鱼的饵料基础(吴瑞等, 2016)。二者都是海洋生态系统能量流动和物质循环的重要环节, 也是滤食性贝类饵料的重要组成部分(刘瑞玉, 1992)。一般浮游生物可分为六类:巨型浮游生物, 大于1cm, 如海蜇；大型浮游生物, 5—10mm, 如大型桡足类、磷虾类；中型浮游生物, 1—5mm, 如小型水母、桡足类；小型浮游生物, 50μm—1mm, 如硅藻、蓝藻；微型浮游生物, 5—50μm, 如甲藻、金藻；超微型浮游生物, 小于5μm, 如细菌(王先锋等, 2017)。其中, 微型也是较大浮游生物种类生长周期中某一阶段(高亚辉等, 1989)。

KTD、DG两海域分别位于黄海南、北两端, 都是山东省重要贝类养殖区, 养殖有大量的滤食性海洋双壳贝类。双壳经济贝类食物营养级较广, 对符合其口径大小的真核浮游生物皆有摄食, 其中包括了大部分微型和微微型的浮游植物及部分微型的浮游动物(Davenport et al, 2000; Lehane et al, 2006)。由于微微型、微型浮游生物体型微小, 传统的显微镜检测所需的样本量又较大, 准确判断滤食性贝类的养殖区饵料成分较为困难。本研究为了更好的对养殖区浮游生物物种类和多样性进行分析, 以18S rDNA V4区为目标基因, 采用高通量测序技术, 对山东半岛南北两侧KTD、DG海域的真核浮游生物群落结构逐月进行了调查, 同时, 测定每个月份的环境因子, 并对环境因子和浮游动植物丰度进行关联分析。研究结果对了解双壳经济贝类养殖区饵料组成及其在时空的变化, 以及引起饵料组成变化的原因, 对海岸带食物网、生态基础管理和海洋经济贝类养殖生产等方面提供数据支撑。

1 材料与方法 1.1 样品采集及环境因子调查

在KTD和DG海域近岸养殖区设置站位, 采样区域见图 1。自2017年11月开始至2018年7月, 逐月采集水样, 水样经0.2μm滤膜过滤, 将滤膜转移至1.5mL无菌离心管中, 置于–20℃条件下冷冻保存运输。

环境因子调查与生物样品取样同期进行, 包括氨氮(NH₄⁺-N)、亚硝酸盐氮(NO₂^--N)、硝酸盐氮(NO₃^--N)、硅酸盐(SiO₃^2--Si)、活性磷酸盐(PO₄^3--P)、叶绿素a(chl a)、溶氧(DO)、水温(SST)、盐度(S)、透明度(TRA)等, 各类营养盐通过各站点采集1L海水水样, 装入至聚乙烯瓶中, 并加入分析纯级氯仿进行固定；叶绿素a采用分光光度法测定, 通过各站点采集500mL海水, 并加入5mL碳酸镁, 经过0.22μm孔径的微孔滤膜抽滤获得；水温用表层温度计测定, 透明度使用透明度盘测定, 样品的采集、保存及分析均按中华人民共和国国家标准《海洋调查规范》(GB12763-2007)进行。

1.2 基因组DNA的提取

将0.22μm滤膜剪碎置于1.5mL离心管中, 采用CTAB法提取滤膜上所有物种的宏基因组, 利用琼脂糖凝胶电泳和紫外分光光度计检测DNA的纯度和浓度, 检测合格后置于–20℃冰箱保存备用。

1.3 18S V4区PCR扩增

将检测合格的基因组DNA交由诺禾致源生物信息科技有限公司, 以稀释后的基因组DNA为模板, 使用带Barcode的特异引物、New England Biolabs公司的Phusion^® High-Fidelity PCR Master Mix with GC Buffer和高效高保真酶进行18S V4区DNA的PCR扩增。V4引物序列为528F (5′-GCGGTAATTCCAG CTCCA-3′)和706R (5′-AATCCRAGAATTTCACC TCT-3′)。PCR产物使用2%浓度的琼脂糖凝胶进行电泳检测合格后, 根据PCR产物浓度进行等量混样, 使用1×TAE浓度2%的琼脂糖胶电泳纯化PCR产物, 剪切回收目标条带。

1.4 文库构建和上机测序

基于IonS5TMXL测序平台, 利用单端测序(Single-End)的方法, 使用Thermofisher公司的Ion Plus Fragment Library Kit 48 RXNS建库试剂盒进行文库的构建, 构建好的文库经过Qubit定量和文库检测合格后, 使用Thermofisher的Ion S5TMXL进行上机测序。

1.5 数据分析 1.5.1 测序数据分析

测序得到的原始数据(Raw data)使用Cutadapt软件过滤和按barcode拆分样本后得到有效数据(Clean data)。为了研究样品的物种组成多样性, 对所有样品的有效数据(Clean reads)进行聚类, 以97%的一致性(Identity)将序列聚类成为OTUs (Operational Taxonomic Units)。然后对OTUs的代表序列进行物种注释, 并在各分类水平上进行群落结构的统计, 得到对应的物种信息和基于物种的丰度分布情况。同时, 对OTUs进行丰度、Alpha多样性计算等分析, 以得到样品内物种丰富度和均匀度信息、不同样品间的共有和特有OTUs信息等。另一方面, 通过主坐标分析(PCoA, Principal Co-ordinates Analysis)展示各样品之间的相似与差异。

1.5.2 理化指标及环境因子相关性分析

环境因子相关数据采用SPSS 17.0软件进行单因素方差分析(One-way ANOVA), 差异显著则进行Duncan’s多重比较检验, 显著水平设定为0.05。统计数据均以平均值±标准差(Mean±SD)的形式表示。对环境因子与物种及物种丰富度之间进行(Canonical Correspondence Analysis, CCA)典范对应分析(Sheik et al, 2012)。用CCA典范对应分析可以检测环境因子、样品、种群三者之间的关系或者两两之间的关系, 可得到影响样品分布的重要环境驱动因子。

2 结果与分析 2.1 各样品有效序列数及物种OTUs数目统计

对测序得到的有效序列数进行OTUs聚类, 每个OTUs代表一类物种。本研究对聚类得到的OTUs数目进行统计, 并根据OTUs的数目做样品稀释曲线发现, 曲线显示所有样品OTUs数目都趋于饱和, 表明样品测序数据量较为合理, 再多的测序数据量只会产生少量的物种(图 2)。对聚类得到的OTUs数目进行统计发现:所有月份样品中种水平总共455个OTUs。其中KTD海域的OTUs数目总体高于DG海域, KTD海域浮游植物OTUs数目在2018年的1月、2月、3月、4月较高, 最低是2018年的6月；浮游动物OTUs数目最高是在2018年的4月、5月, 最低是在2018年的6月。DG海域浮游植物OTUs数目最高在2017年的11月、12月份, 最低是在2018年的4月；浮游动物OTUs数目最高是在2017年的12月, 最低是在2018年的4月(表 1)。

2.2 KTD、DG海域各月份样品群落多样性指数(Alpha多样性)统计分析

多样性指数(Alpha Diversity)可用于分析样品内的生物群落多样性, 通过此分析我们可以对各样品内生物群落的丰富度和多样性进行检验(Li et al, 2013)。本研究应用QIMME软件(Caporaso et al, 2010)计算得到各海水样品中的生物群落多样性指数见表 2。由表 2可见, 各样品文库的覆盖率(Coverage)指数较高, 表明样本中序列没有被测出的概率较低, 测序结果能够反映样品中生物的真实情况。ACE指数与Chao1指数是用来估计种群丰度的的指数, 它们的值越大, 表明样本的生物种类越多, 丰富度越高。Shannon多样性指数是用来估算样品中生物多样性的指数, 其值越高, 说明群落多样性越高。本研究对各样品Alpha多样性指数统计显示: KTD海域2017年11月、2018年1月、2月、3月物种多样性较高；2018年6月、7月物种多样性最低。DG海域2017年12月物种多样性较高, 2018年2月、4月、5月物种多样性偏低。

2.3 真核浮游生物组成分析

用RDP Classifier方法(Wang et al, 2007)与Sliva数据库(Quast et al, 2013)对各样品OTUs代表序列进行物种注释, 并分别在各个分类水平:界(kingdom)、门(phylum)、纲(class)、目(order)、科(family)、属(genus)、种(species)统计各样本的群落组成。其中, 在门分类水平上, KTD海域9个月共检测出36个门类。DG海域共检测出34个门类；在种分类水平上, KTD海域四个季度共检出物种424种, DG海域四个季度共检出物种365种。选取每个样品相对丰度在门水平排名前10位的类群, 生成各门相对丰度柱形累加图, 选取相对丰度在属水平上排名前30位的类群, 生成各属相对丰度柱形累加图(图 3, 图 4)。

由图 3可见, 各个月份样品之间门水平上差异明显。其中, 有待进一步明确其分类地位的门类占16%—61%。绿藻门(Chlorophyta)、硅藻门(Diatomea)是两养殖海域浮游植物中整体丰度最高的门类, 在KTD海域, 绿藻门各月丰度在3.0%—21.3%之间, 其中7月(K1807)最高, 达到了21.3%；硅藻门各月丰度在2.0%—16.6%之间, 其中2月(K1802)最高, 达到了16.6%。DG海域, 绿藻门各月丰度在2.0%—12.3%之间, 其中11月(D1711)最高, 达到了12.3%；硅藻门各月丰度在2.0%—47.0%之间, 其中1月(D1801)最高, 达到了47.0%。浮游动物在两地海水群落中占了很大比重, 其中占地位的主要是节肢动物门类的物种, 其每月丰度分别在6.0%—38.9%和7.6%—48.6%之间。在KTD海域的2017年11月、12月和2018年的1月、2月、3月、4月、7月以及DG海域的2017年11月、12月和2018年的1月、3月、4月、5月、7月的样品中占比较高。

由图 4可见, 除却有待明确分类地位的一些真核生物外, 丰度在属水平上排名前30位的类群中, 浮游动物主要有桡足类的哲水蚤亚目(Calanoida)、纤毛虫类的海游虫(Pelagostrobilidium)及Parastrombidinopsis、环节动物门的叶须虫(Phyllodocida)、双壳类的贻贝目(Mytiloida)和牡蛎目(Ostreoida)、海鞘类的复鳃目(Stolidobranchia)等。浮游植物中丰度排名前30的属较多, 主要有甲藻门的夜光藻属(Noctiluca sp.)、异甲藻属(Heterocapsa sp.)、双鞭甲藻类(Haplozoon)、沟鞭甲藻类(Gyrodiniellum)、营寄生的鞭毛虫类(Haplozoon sp.)；绿藻门的葡萄藻属(Bathycoccus sp.)、微球藻属(Ostreococcus sp.)、微胞藻(Micromonas sp.)、塔胞藻属(Pyramimonas sp.)；硅藻门的海链藻属(Thalassiosira sp.)、角刺藻属(Corethron sp.)、楔形藻属(Licmophora sp.)、细柱藻属(Leptocylindrus sp.)、舟形藻属(Navicula sp.)以及隐藻门的全沟藻属(Teleaulax sp.)等。

2.4 KTD、DG海域各月份样品浮游动、植物群落比较分析

利用PCoA方法探明了各样品之间的关系(图 5), PCoA点聚集程度表明两地样品中, 不同月份之间存在相似和差异, 其中KTD海域样品2017年11月(K1711)与12月(K1712)相似程度较高；KTD海域2018年3月(K1803)样品与DG海域2018年2、3、4月(D1802、D1803、D1804)样品相似程度较高；KTD海域2018年1月(K1801)样品与DG海域2017年12月(D1712)样品相似程度较高, 其余月份之间差异较大。

2.5 环境因子数据及相关性分析 2.5.1 环境因子数据分析

调查海域中, 监测到KTD海域9个月的水温最高是2018年7月的21.2℃, 最低是同年2月的4.5℃；DG海域8个月最高水温是2018年7月的28.4℃, 最低水温是同年1月的5℃。其间监测KTD海域叶绿素a浓度冬春季节(2018年1—5月)显著增加, 最高是2018年2月份的3.4μg/L；DG海域叶绿素a浓度高水平出现在2017年的12月和2018年的1月及4月, 最高值是2017年12月份的5.4μg/L, 也较符合冬春季节含量较高的特征。五项营养盐中氨氮(NH₄⁺-N)、活性磷酸盐(PO₄^3--P)含量较为稳定, 两海域各月份之间差异不显著；亚硝酸盐氮(NO₂^--N)、硝酸盐氮(NO₃^--N)、硅酸盐(SiO₃^2--Si)、变化波动较大, 其中, KTD海域硝酸盐氮(NO₃^--N)、亚硝酸盐氮(NO₂^--N)、硅酸盐(SiO₃^2--Si)含量整体呈现出一个秋冬季到春季降低, 夏季升高的态势。DG海域除2018年4月份硝酸盐氮(NO₃^--N)、亚硝酸盐氮(NO₂^--N)略有增高外, 整体情况与KTD海域相似, 详见表 3。

2.5.2 环境因子与浮游真核生物群落数据的相关性分析

通过envfit函数筛选出具有显著影响的环境因子与浮游生物群落数据做CCA典范对应分析(图 6)。CCA分析图中, 箭头表示环境因子, 箭头连线的长度代表环境因子与群落分布和种类分布间相关程度的大小, 箭头越长, 说明相关性越大, 反之越小。箭头连线和排序轴的夹角代表某个环境因子与排序轴的相关性大小, 夹角越小, 相关性越高, 反之越低。环境因子之间的夹角为锐角时表示两个环境因子之间呈正相关关系, 钝角表示呈负相关关系。由此图可以发现各月份的水温、溶氧(DO)、pH、硅酸盐、硝酸盐氮, 对各浮游生物门类群落结构的改变有较大影响。水温与pH、硝酸盐氮和硅酸盐之间呈正相关, 与其他环境因子呈负相关。CCA排序图上各门物种的位置反映了各物种对于不同环境资源的依赖程度。由此图可以看出, 海水样品中丰度较高的硅藻门和绿藻门浮游植物与溶氧、硝酸盐、硅酸盐之间呈现了较大的相关性。浮游动物中的纤毛虫类、鞭毛虫及后生动物门类的种则与水温、溶氧、pH等相关性较大。

3 讨论 3.1 真核浮游生物种类鉴别方法的选择

高通量测序技术的广泛应用, 极大地推动了浮游生物多样性的研究, 尤其是对显微镜难以分辨的微型及微微型浮游生物(Howard et al, 2006; Amaral-Zettler et al, 2010; 王丽平等, 2015)。于杰(2014)采用高通量测序技术, 以18S rDNA V9区为目标基因检测了渤海山海关海域微型浮游生物多样性；宋伦等(2016)以18S rDNA V4区为目标基因, 结合高通量测序技术, 对渤海长兴岛海域微型和微微型浮游植物多样性进行了检测, 研究结果表明V4引物在微型、微微型浮游植物群落鉴定方面更为高效。王先锋等(2017)通过高通量测序技术, 扩增了18S rDNA V4区DNA, 研究了平岛海域长牡蛎的真核食物组成。本研究以18S rDNA的V4区为目标基因, 通过高通量测序技术, 检测了海水样品中的真核浮游生物的组成, 相对于传统镜检的方法, 该技术具有简单快速、测序通量高、错误率低和成本低等特点, 节省了大量的人力物力, 能够很好地对两海域近岸养殖区浮游生物群落多样性进行检测。

3.2 KTD和DG两海域环境因子分析

本研究发现, KTD海域叶绿素a浓度冬春季节(2018年1—5月)显著增加, 最高是2018年2月份；DG海域叶绿素a浓度高水平出现在2017年的12月和2018年的1月及4月, 最高值是2017年12月份, 也呈现出冬春季节含量较高的特征。认为叶绿素a和溶氧的增高是受到浮游植物丰度升高和水产经济作物养殖量减少的影响。冬季前后, 正是滤食性经济贝类陆续大量上市的季节, 对水体中浮游植物的摄食量大大减少, 且冬春季节温度降低, 贝类呼吸代谢及摄食率也大大降低(董波等, 2000)。

五项营养盐中氨氮(NH₄⁺-N)、活性磷酸盐(PO₄^3--P)含量较为稳定, 两海域各月份之间差异不显著；亚硝酸盐氮(NO₂^--N)、硝酸盐氮(NO₃^--N)、硅酸盐(SiO₃^2--Si)、变化波动较大, 其中, KTD海域硝酸盐氮(NO₃^--N)、亚硝酸盐氮(NO₂^--N)、硅酸盐(SiO₃^2--Si)含量整体呈现出一个秋冬季到春季降低, 夏季升高的态势。DG海域除2018年4月份硝酸盐氮(NO₃^--N)、亚硝酸盐氮(NO₂^--N)略有增高外, 整体情况与KTD海域相似。认为是受到气候因素和地理位置影响:两处海域都是近岸海域, 有沿岸多条径流汇入；秋冬季节降水减少, 河流径流量也随之减少；夏季则降雨增多, 雨水及河流径流的补充, 使得两处海域营养盐含量得以补充。总体来说, 近岸海域尤其是距海岸带不远的海域, 受地理位置因素和气候因素干扰较大。

3.3 浮游真核生物的群落数据与环境因子相关性分析

综合分析浮游真核生物的群落数据与环境因子的相关数据发现, 真核浮游生物各门类受溶氧(DO)、pH、水温、硅酸盐和硝酸盐氮的影响较大, 且不同样品之间的差异较为显著。水温与pH及硝酸盐氮之间呈正相关, 与海水中溶氧含量呈负相关；硅酸盐和硝酸盐之间呈正相关, 与水温、pH及溶氧呈负相关。认为在春季至夏季水温升高的过程中, 降雪降雨的增加及河流流入带来水体中硝酸盐氮和硅酸盐含量的增高, 促进了浮游植物丰富度的增加。同理, 硝酸盐氮与硅酸盐含量的增高使浮游植物丰富度增加的同时, 也促进了水体中溶氧的增高(图 6)。而叶绿素a的浓度可以作为浮游植物生物量的指标, 可以反映海水中浮游植物现存量的变化(张学雷等, 2004), 通过观测叶绿素a浓度观测值的变化规律也可以得出相同结论。

3.4 影响浮游动植物优势种的因素分析

本研究结果显示, 两地养殖海域中的浮游植物类群主要优势种以硅藻、绿藻和甲藻为主, 浮游动物主要是纤毛类、鞭毛类和一些后生动物门类种。从两地海域浮游动植物的优势种分析来看, 两处海域饵料资源较为丰富, 生产力水平较高, 适合发展水产养殖业。

从浮游生物群落多样性指数的分析结果来看(表 2), 秋冬季节物种多样性较春夏两季要高, 且检测出来的总OTUs个数也保持一个较高的数目, 这与谢斌等(2017)研究的结论—海州湾浮游植物种类和丰度季节变化的表现为秋季 > 夏季 > 春季的趋势相似。谢斌等(2017)及高春梅等(2015)认为, 无机氮在春、夏、秋3季均由沉积物向水体进行迁移, 在很大程度上增加了水体中各类营养盐的含量, 直接或间接促进了各类浮游生物丰度的变化。对KTD和DG两海域而言, 两地水产养殖业较为发达, 主要养殖种类有海湾扇贝(Argopecten irradias)、长牡蛎(Crassostrea gigas)和贻贝(Mytilus edulis), 其摄食饵料范围包括大部分微型、微微型浮游植物以及一部分大小与贝类摄食口径合适的纤毛虫、鞭毛虫类。冬季前后, 正是滤食性贝类大量上市的季节, 贝类的大量产出, 减少了对作为饵料的浮游植物和一些浮游动物的摄食, 而每年的4—6月, 又是贝类性成熟繁殖阶段, 春到夏季, 大量稚贝又回到海区, 对海水中浮游植物的摄食量又大大增加, 这也符合本调查中叶绿素a浓度的变化规律, 根据此次针对浮游生物种类及丰度的调查, 结合我国北方沿海浮游动物摄食率的研究以及贝类摄食机制的特征(Wang et al, 1998; 董婧等, 2003; 张学雷等, 2004), 综合认为, 除环境因子对浮游生物种群的影响, 贝类摄食是对养殖区浮游生物种类及丰度造成改变的主要原因。

4 结论

(1) KTD和DG海域浮游生物种类较丰富, 两养殖海域浮游植物中的优势种群主要是硅藻门、绿藻门、甲藻门的物种；浮游动物在两地海域生物群落中占优势地位的主要是节肢动物门、纤毛虫门类的物种。

(2) 两海域的真核生物群落组成受水温、溶氧、硝酸盐氮和硅酸盐的影响较大, 贝类摄食是影响两海域浮游生物群落多样性的重要因素。