莱州湾近岸海域浮游纤毛虫群落组成和季节变化*

李雯璐 徐大鹏 季道德①

(1. 烟台大学海洋学院 山东烟台 264005; 2. 厦门大学近海海洋环境与科学国家重点实验室 海洋与地球学院 海洋微型生

物与地球圈层研究所 福建厦门 361104)

浮游纤毛虫是一类广泛分布于海洋中上层、营浮游生活的纤毛虫, 粒级大小一般在5～200 μm 之间。它们是pico-级(0.2～2 μm)和nano-级(2～20 μm)浮游生物的主要消费者(Suzukiet al, 2007), 同时也是后生动物和鱼类幼虫的食物来源(Stoeckeret al, 1987;Gómez, 2007), 在连接微食物网和传统食物链的物质循环和能量流动中发挥着重要作用(Pierceet al, 1992;Weisseet al, 2022)。

海洋浮游纤毛虫主要包括砂壳纤毛虫和无壳寡毛类, 也包括一些前口类、盾纤类等类群。依照Lynn(2008)的分类系统, 砂壳类纤毛虫隶属于旋毛纲(Spirotrichea)、环毛亚纲(Choreotrichia)、砂壳目(Tintinnida), 其重要特征是具有管状或瓶状的外壳,因其近岸丰度高、采样易获得的特点, 目前世界各海区已经鉴定约1 000 种(Dolan, 2010)。壳的口径和长度、壳的形状、黏着颗粒的有无是砂壳类重要的分类学依据(Capriuloet al, 1982; Dolanet al, 2002)。无壳寡毛类纤毛虫隶属于旋毛纲(Spirotrichea)、寡毛亚纲(Oligotrichia)、寡毛目(Oligotrichida)。目前已报道无壳寡毛类约200 余种(Liuet al, 2017), 然而据Foissner等人(2008)和Agatha (2011)估计, 无壳寡毛类在世界范围内约有830～1 280 种, 其多样性可能有83%～89%以上是未知的。

传统的形态学由于分辨率低、费时费力等弊端,难以避免地对形态上不可区分的物种(即隐性多样性)产生遗漏。高通量测序(high-throughput sequencing,HTS)因其通量高、耗时短、相对成本低等特点, 对微型生物生态学研究产生了革命性影响。对环境DNA 的高通量测序可以揭示纤毛虫的多样性、分布模式和长时间尺度上的动态变化, 从而获得更为全面的群落信息(Slapetaet al, 2006; De Schepperet al,2019)。然而研究显示基于HTS 得出的纤毛虫多样性大大超过已知形态物种数量(Bachyet al, 2013), 甚至高达6 倍之多(De Vargaset al, 2017), 因此基于形态和分子的分布模式是否一致尚未可知。近年来,越来越多的研究者尝试将多种方法相结合对纤毛虫进行生态调查: Santoferrara 等人(2016)对Lugol’s 固定样品镜检的同时结合HTS, 对大西洋砂壳类纤毛虫的群落组成和分布模式进行研究。Zhao 等(2017)利用HTS 比较了环境DNA 和RNA 两种方法在纤毛虫群落垂直分布上的差异, 显示出RNA 获得的多样性高于DNA。在调查南海浮游纤毛虫多样性和垂直分布时, Sun 等(2019)将HTS 与定量蛋白银染色(QPS)相结合, 结果显示两种方法在揭示较低分类等级的群落组成时存在差异。另外, 有研究者将QPS 的形态学方法与DNA 及其cDNA 的HTS 三者相结合, 调查沉积物中纤毛虫的多样性(黄平平等, 2017)。综合来看, 利用分子测序与传统形态学方法所揭示的群落组成及其变化趋势较为一致, 尤其是常见/优势物种, 而高通量测序因其测序深度高可发现大量稀有物种序列。因此将两种方法相结合能够在一定程度上弥补两者的弊端, 从而获得更为全面的纤毛虫群落信息。

近岸环境与人类社会的关系最为密切, 且同时受到人为因素和自然因素的双重影响而具有一定的特殊性。一方面, 三山岛近岸海域受到多重人类生产活动的影响, 海水富营养化严重、水质下降; 另外, 受到风暴潮等自然因素的影响, 水体扰动频繁使得浑浊度较高。这些因素可能对该生境中微食物网成员(如纤毛虫)的丰度、物种多样性和群落结构产生影响。本研究首次利用传统形态学和分子测序技术相结合的手段,对莱州湾三山岛近岸海域的浮游纤毛虫进行周年生态调查, 旨在了解该海域的浮游纤毛虫的多样性、群落组成和时空变化, 以及对环境因子的响应。

1 材料与方法

1.1 调查站位、样品采集和预处理

本研究于2020 年9 月至2021 年8 月在莱州湾三山岛近岸海域设置3 个采样点, A、B、C 站点分别位于码头、离岸300 m 处和养殖厂排水口(图1), 使用渔船进行每月两次的浮游样品采集。由于海况原因,缺失3 个站点的1、2、7、8 月上旬和5 月下旬的样品及数据, 以及离岸的B 站位11、1 和4 月下旬样品及数据。因此, 秋季(9～11月)共采集样品17 份, 冬季(12～2月)样品11 份, 春季(3～5月)样品14 份, 夏季(6～8月)样品12 份, 共计54 个样品。

图1 莱州湾三山岛近岸采样点示意图Fig.1 Locations of the sampling sites in Sanshan Island, Laizhou Bay

使用便携式水质分析仪(YSI, USA)测量原位水温、盐度、溶解氧和pH 值。使用卡扣式采水器采集表层水样(深度20 cm), 每一个站位采集两个水样各1 L, 混合为一个2 L 的样品。用孔径为200 μm 的筛绢进行预过滤, 再使用聚碳酸酯滤膜 (Millipore,Ireland) (直径45 mm、孔径0.22 μm)对1 500 mL 水样进行抽滤收集, 设置三组重复(每张膜500 mL)。将抽滤得到的含有环境DNA 的滤膜保存于冻存管中, 立即置于－80 °C 冰箱中保存。取100 mL 表层原位海水,经滤膜孔径为0.45 μm 的针头过滤器过滤后收集在聚乙烯瓶中, 并置于－20 °C 冰箱中保存, 用于营养盐分析。氨态氮(NH4＋-N ), 硝态氮(NO3--N ), 亚硝态氮(NO2--N ), 磷酸盐(PO34-)和硅酸盐含量(SiO23-)使用营养盐自动分析仪(Seal, Germany)进行分析。水样的总碳(TC)和总无机碳(TIC)含量使用Vario 微型立方元素分析仪(Elementar, Germany)进行测定。使用400 目浮游生物网(孔径37 μm)在A、B 两个站点水平拖网采集镜检样品并用广口采样瓶收集, 获得35 个样品用于纤毛虫的形态学分类鉴定。

1.2 高通量测序方法

环境总DNA 使用Fast DNA SPIN Kit for Soil DNA快速提取试剂盒(Mpbio, USA) 进行提取。使用NanoDrop One (Thermo-Fisher, USA)对DNA 的浓度和纯度进行测定后, 于－80 °C 保存。

利用真核生物通用引物528F (5′-GCGGTAATT CCAGCTCCAA-3′)和 706R (5′-AATCCRAGAATTT CACCTCT-3′)(Cheunget al, 2010) 对18S rRNA 基因的V4 区进行PCR 扩增, 于北京诺禾致源生物信息科技有限公司利用Illumina NovaSeq 测序平台对文库进行双末端PE250 测序。

1.3 数据处理及分析

从下机数据中拆分出各样本数据, 截去barcode和引物序列后使用FLASH 软件(V 1.2.11, http://ccb.jhu.edu/software/FLASH/)(Magočet al, 2011)对样本的测序结果进行拼接, 得到低质量原始测序数据。随后使用fastp 软件对得到的低质量原始测序数据进行质控, 得到高质量数据。利用QIIME 和Mothur 软件包(Schlosset al, 2009; Caporasoet al, 2010)去除嵌合体, 利用UCLUST 以97%的相似性水平划分可操作分类单元(operational taxonomic units, OTUs), 去除只有一条代表序列的OTU, 利用原生生物核糖体参考数据库(Protist Ribosomal Reference, PR2)数据库进行物种分类注释, 得到最终的OTU 表格。抽取纤毛门数据用于纤毛虫群落的分析, 将各样品的序列相对丰度数据按照季节划分整合。

为了保证样品间α 多样性比较的可靠性和准确性, 以所有样品中最小序列数为标准, 对样品数据进行标准化, 然后计算α 多样性指数(Chao 1、Richness、Shannon 和 Simpson 指数)。通过单因素方差分析(ANOVA)对浮游原生生物α 多样性的时空差异进行统计检验。

基于Bray-Curtis 距离矩阵进行浮游原生生物群落的β 多样性的分析。使用PRIMER 6.0 软件对浮游原生生物群落进行非计量多维尺度分析(non-metric multidimensional scaling, NMDS); 利 用 R 中 的ggplot2 包进行主坐标分析(principal co-ordinates analysis, PCoA), 并利用ANOSIM 分析检验三个站点不同季节的浮游原生生物群落结构的差异情况; 使用SPSS 软件对浮游原生生物主要类群和环境因子进行Spearman 相关分析; 使用R 中的vegan 包进行冗余分析(redundancy analysis, RDA), 寻找影响浮游原生生物群落结构的主要环境因子。

1.4 形态学方法

将拖网样品充分混匀后利用解剖镜进行镜检观察, 依据文献对浮游纤毛虫进行鉴定和计数(Kofoidet al, 1929; 宋微波等, 2009; 张武昌等, 2012), 并拍照记录生物学特征及出现频次, 以获得浮游纤毛虫多样性信息。将镜检鉴定后的纤毛虫用灭菌水清洗数次后转移至灭菌后的1.5 mL 离心管中, 用于提取DNA。引物为: EukA (5′-AACCTGGTTGATCCTGCCA GT-3′)、EukB (5′-TGATCCTTCTGCAGGTTCACCT AC-3′)、Euk82F (5′-GAAACTGCGAATGGCTC-3′)和U1492R (5′-GGTTACCTTGTTACGACTT-3′)(Medlinet al, 1988; López-Garcíaet al, 2001; Dotyet al,2005)。PCR 扩增程序依照(Wanget al, 2019)进行。使用TIANgel Midi Purification Kit (天根生化科技有限公司, 上海) 试剂盒纯化PCR 产物。将纯化后的扩增产物插入质粒载体 pTZ57R/T (InsTAcloneTM PCR Cloning Kit, Fermentas, Canada), 并将重组质粒转移至感受态细胞Trans1-T1 (全式金生物技术股份有限公司, 中国北京)中。摇菌培养1 h 后接种至琼脂糖培养基上过夜培养。随机挑选6 个重组克隆, 利用引物M13-20 和M13-26 对菌样进行PCR 验证。将满足条件的菌液送至青岛生工测序公司进行测序。将测序后的DNA 序列与数据库中的序列进行比对, 相似度99%以上则认为是同一物种。

2 结果

2.1 环境因子

12 个环境因子(包括温度、盐度、溶解氧、pH、总有机碳、总无机碳、总碳、磷酸盐、氨态氮、亚硝态氮、硝态氮、硅酸盐含量)时空变化如图2 所示。三个站点间的所有环境因子均无显著差异, 有6 个环境因子表现出明显的季节模式: 水温(ANOVA,P＜0.05)、总碳(ANOVA,P＜0.01) 、总无机碳(ANOVA,P＜0.01)和总有机碳含量(ANOVA,P＜0.01)一致呈现为夏秋季＞冬春季的趋势; 盐度(ANOVA,P＜0.01)和溶解氧(ANOVA,P＜0.05)呈现冬春季＞夏秋季的趋势。三个站点间盐度和有机碳含量的季节变化折线图变化趋势相似。

图2 不同站点与季节的环境因子Fig.2 Environmental factors of different sites and seasons.

A 站点水温年变化为1.7～27.9℃, 平均15.4 °C,其中最高温出现在2020 年9 月上旬; 水温最低出现在2021 年1 月下旬; B 站点水温年变化为4～27 °C, 平均14 °C, 其中最高温出现在2021 年8 月下旬; 水温最低出现在2020 年12 月和2021 年2 月下旬; C 站点水温年变化为10.9～30.4 °C, 平均22.2 °C, 其中最高温出现在2021 年7 月下旬, 水温最低出现在2021 年2 月下旬。三个站点的溶解氧在夏秋两季含量较为一致且最低值出现在夏季, 而冬季A 站点溶解氧含量最高(27.14 mg/L), 其次是B 站点(21.52 mg/L)。A、B两个站点的pH 值在秋冬两季的变化幅度较大。A、B站点的磷酸盐含量在春夏两季较高, 氨态氮和亚硝态氮含量在春夏季交替时升高, 硝态氮则相反, 在春夏季交替时含量降低, 硅酸盐含量在夏秋两季最高,而在春季达到最低值。C 站点的营养盐含量季节变化较大, 氨态氮、硝态氮和亚硝态氮含量均呈现出夏季最低, 秋冬季高的特点, 最低值出现在7 月上旬和8月下旬, 最高值出现在9 月上旬和12 月下旬。磷酸盐含量呈现夏季最低, 春秋季最高, 最低值出现在6月下旬, 最高值出现在4 月上旬。

在三个站点间, A、B 站点环境因子的全年变化拟合程度较高, 但pH 在秋冬季变化趋势相反, 磷酸盐浓度在春夏季呈现明显不同的变化。总的来看, C 站点海水理化因子变化较为平缓而营养盐浓度变化浮动较大, 与A、B 站点变化拟合程度较差。

2.2 群落组成与相对丰度

高通量测序结果表明, 旋毛纲(Spirotrichea)是最优势类群, 占纤毛虫序列的64.10%, 其次是寡膜纲(Oligohymenophorea, 2.17%) (图 3a)。将旋毛纲(Spirotrichea)在目(order)水平上进一步细分如图3b 所示, 无壳寡毛类(急游虫目和环毛目共占序列总丰度的40.14%)和砂壳类(砂壳目, 10.67%)是浮游纤毛虫的主要优势类群, 并且显示出明显的季节模式。其中急游虫目(Strombidiida)序列的相对丰度为 27.27%,环毛目(Choreotrichida)为12.88%。急游虫目在全年的变化较为稳定并且具有明显的季节变化模式和站点差异。季节变化表现为秋冬＞春夏, 且在春夏季样品中呈现A＞B＞C, 在秋冬季样品B＞A＞C。环毛目在站点间无明显差异, 而显示出秋＞夏＞冬＞春的季节模式,春季相对丰度明显降低。砂壳目在站点间差异显示出A(4.34%)＞B(4.10%)＞C(2.23%), 且A 和B 站点砂壳目在夏秋两季的相对丰度相近, 而在冬春两季的相对丰度显示出相反的变化模式。

寡毛类和砂壳类在科(family)级水平上的群落组成如图 3c 所示, 寡毛类中相对丰度较高的有Strombidiidae_M (16.62%)、Strombidiidae (14.63%)、Strombidiidae_R (9.77%) 和Strobilidiidae_G (8.00%);砂壳类中相对丰度较高的是Tintinnida_X (7.59%)、Favellidae (2.53%) 和TIN_07 (2.32%)。

在属阶元上对寡毛类和砂壳类进行进一步划分如图 3 d, 寡毛类的相对丰度主要由急游虫属(St ro m b i d i u m, 1 9.0 1%)和急游虫科未定类群(Strombidiidae_X, 12.24%)组成, 相对丰度显示出先下降后上升的趋势, 并在冬季达到最大值。砂壳类主要由砂壳目未定类群(Tintinnida_XX, 3.84%)、拟铃虫属(Tintinnopsis, 3.27%)和网纹虫属(Favella, 1.32%)组成, 其相对丰度在春夏季达到最大值后下降, 冬季的相对丰度最低。

图3 浮游纤毛虫在纲(a)、目(b)阶元上的相对丰度, 以及主要优势类群在科(c)、属(d)阶元上的相对丰度Fig.3 Relative sequence abundance of planktonic ciliates at class (a), order (b), family (c), and genus (d) level

站点间的物种组成存在明显差异, 值得注意的是无论处于哪个分类阶元, C 站点的夏季样品在物种组成和相对丰度上均与其他分组样品有明显的不同,主要表现在总的物种多样性低, 以及除寡毛类和砂壳类外的纤毛虫相对丰度与其他样品有明显差别。在纲水平上, 旋毛纲相对丰度明显降低, 而叶咽纲(Phyllopharyngea)、叶口纲(Litostomatea)、寡膜纲(Oligohymenophore)具有明显高的相对丰度。在旋毛纲中, 管口目(Cyrtophoria)、前口目(Haptoria)、缘毛目(Peritrichia)的相对丰度明显高于其他站点, 而CONThreeP_X、Tintinnida 和Strombidiida 明显低于其他站点。A、B 站点物种组成的季节变化较为相似(图3b), 不同之处主要在于A 站点的寡毛类纤毛虫相对丰度在春季明显大于B站点, 随着季节推移二者差距减小, 直至秋季相对丰度相当。冬季则呈现出B 站点的急游虫目纤毛虫相对丰度大于A 站点。而砂壳类纤毛虫的相对丰度变化趋势相反, 在春季B 站点大于A 站点, 随后差距缩小, 直至冬季A 站点砂壳类纤毛虫的相对丰度大于B 站点。

2.3 α 多样性

计算54 个样品的α 多样性指数, 并进行站点和季节间的ANOVA 方差分析检验, 结果如表 1 所示,浮游纤毛虫群落的丰富度和多样性指数均未显示出站点间和季节间差异(P＞0.05)。Chao 1 指数显示出秋季＞春季＞冬季＞夏季的季节分布, 而除Chao 1 指数外的α 多样性指数均显示出夏季＞秋季＞春季＞冬季, 其中春秋两季的α 多样性指数均最为接近。

为了进一步明确α 多样性指数和环境因子之间的相关性, 我们进行了Spearman 相关性分析, 结果如表2 所示。Chao 1 指数与总碳和总无机碳含量显著正相关(R＞0,P＜0.05)。Richness 指数与所有环境因子均不呈现相关性。Shannon 指数和Simpson 指数均与盐度呈现显著相关关系(R＜0,P＜0.05)。

表 1 各季节和各站点的浮游纤毛虫α-多样性指数Tab.1 α-diversity of planktonic ciliates in different samples

表2 浮游纤毛虫α 多样性指数与环境因子的Spearman相关性分析Tab.2 Spearman’s correlation coefficients between α diversity indices and environmental factors

2.4 β 多样性

基于Bray-Curtis 距离矩阵对分组后的浮游纤毛虫群落进行非度量多维尺度分析(non-metric multidimensional scaling, NMDS)。如图4 所示, 每个季节的浮游纤毛虫样品聚合性良好且和其他季节样品有明显分离, 说明纤毛虫群落具有季节性。秋冬两季样品有部分交叉并且和春夏季样品显著分开。而不同站点的样品混合散乱排布且重叠明显, 纤毛虫群落未显示出明显的水平分布差异。

图4 不同季节(a)和不同站点(b)的浮游纤毛虫群落结构的非度量多维尺度图Fig.4 Nonmetric multidimensional scaling plots showing the community structure of planktonic ciliates in season (a) and in region (b)

对不同组间浮游纤毛虫群落结构差异性进行ANOSIM 检验可知(表3), 不同站点的纤毛虫群落结构没有显著性差异(ANOSIM,P＞0.05)。而不同季节的纤毛虫群落结构差异显著(ANOSIM,P＜0.01)。

表3 站点和季节变化对浮游纤毛虫群落β-多样性的ANOSIM 分析Tab.3 ANOSIM results on changes of plankton ciliates β-diversity in season and in region

?

2.5 与环境因子间的相关关系

为了解影响浮游纤毛虫群落多样性差异的环境因子, 对12 个环境因子以及54 个样品进行了冗余分析(RDA), 结果如图5 所示。总有机碳和溶解氧是影响三山岛近岸海域浮游纤毛虫群落结构的主要环境因子。总有机碳含量与B 站点的部分样品呈现正相关性; 与C 站点的夏秋季样品的群落结构呈现正相关性。溶解氧含量与A 站点样品和B 站点绝大多数样品的群落结构呈现正相关, 且总有机碳与溶解氧含量之间显示出负相关关系。

图5 浮游纤毛虫群落与环境因子的冗余分析Fig.5 The redundancy analyses (RDA) plot showing correlations between ciliate communities and environmental factors

为了解纤毛虫主要类群与环境因子的相关性,对属(genus)阶元下的无壳寡毛类和砂壳纤毛虫与环境因子进行了Spearman 相关性分析(表4), 与环境因子无相关关系的类群未呈现在表格中。结果显示无壳寡毛类和砂壳类纤毛虫与氨态氮和磷酸盐均无相关性。可以看到在无壳寡毛类中相对序列丰度最高的Choreotrichida_XX 与盐度(R＜0,P＜0.05)和碳含量(R＜0,P＜0.01)呈现显著负相关关系;Leegaardiella与硅酸盐呈显著负相关(R＜0,P＜0.01)。砂壳类中相对序列丰度最高的Tintinnopsis与温度(R＜0,P＜0.01)呈显著负相关, 而与溶解氧(R＞0,P＜0.01)、pH (R＞0,P＜0.05) 和硅酸盐(R＞0,P＜0.01)呈显著正相关关系。

2.6 形态学分类鉴定

对A、B 两个站点的35 个拖网样品进行活体观察和单细胞18S rRNA 基因测序, 累计分离鉴定17 种纤毛虫(图6; 表5), 分别隶属于7 个属: 拟铃虫属(Tintinnopsis)、类铃虫属(Codonellopsis)、网纹虫属(Favella) 、 真 铃 虫 属(Eutintinnus) 、 前 口 虫 属(Frontonia) 、 栉 毛 虫 属(Didinium) 、 旋 游 虫 属(Spirostrombidium), 有9 种已测序定种。其中砂壳类有13 种, 11 种属于黏着壳, 仅有两种透明壳物种。如表5 所示, 砂壳类不仅物种多样性较高, 在检获频率上也远大于其他类群。由于少数物种检获个数和频次少, 未能获得清晰的显微照片或单细胞测序结果, 因此仅能够鉴定到属级阶元或未能呈现。

表5 三山岛近岸海域浮游纤毛虫种名录Tab.5 List of planktonic ciliates species in Sanshan Island

图6 浮游纤毛虫的活体显微照片Fig. 6 Microscopic photograph of planktonic ciliates

与高通量测序结果相比, 形态学结果无论是α 多样性还是β 多样性均较低, 而砂壳类的季节变化趋势在两种方法的结果中较为一致。形态学结果显示拟铃虫属在春夏秋三季均大量出现, 出现频率约为纤毛虫总检获频次的69%, 其中主要在秋季大量检获。这与分子测序结果相似: 拟铃虫属春夏秋三季的序列相对丰度占全年总丰度的80%, 春季的相对丰度最高。形态和分子结果均显示透明壳种类仅在夏秋季高水温时出现。另外, 站点间差异在两种方法的结果中也显示出一致性, 二者均显示A 站点寡毛类和砂壳类多样性大于B 站点。

3 讨论

3.1 浮游纤毛虫群落组成的时空变动

迄今为止在莱州湾进行的浮游纤毛虫生态研究很少, 仅有陈雪等人(2014, 2015)针对砂壳类纤毛虫丰度的时空分布调查。本研究发现无壳寡毛类和砂壳类纤毛虫是浮游纤毛虫的优势类群, 均呈现明显的季节变化。其中无壳寡毛类在秋冬季的序列相对丰度最高, 主导了冬季纤毛虫的多样性和相对丰度, 因此对季节模式具有较高的贡献度。无壳寡毛类纤毛虫在冬季的生存优势可能与它们的营养模式有关。研究发现有些纤毛虫不仅可以进行吞噬营养, 还可以利用内共生藻类或滞留叶绿体进行光合营养(Bernardet al,1994; Pettigrossoet al, 2016), 据Stoecker 等人(2017)估计, 约30%～40%的寡毛类能够进行混合营养。光合作用可以为混合营养纤毛虫提供大量的有机碳,因此相比于专性异养的纤毛虫, 混合营养纤毛虫能够通过营养模式的转变, 提高摄食效率、增加碳源,并降低对氧气、营养盐和食物浓度的依赖程度, 从而更好的适应多变的环境(Estebanet al, 2010)。本研究发现冬季氨态氮含量显著提高, 这种结果在一定程度上与Schoener (2013)的推测相符: 混合营养纤毛虫能够为藻类提供氨态氮, 而非与藻类竞争无机氮。因此推测大多数寡毛类纤毛虫能够降低冬季营养盐变动和食物来源减少带来的不利影响, 从而在冬季样品中占据优势地位, 并能够通过固碳和固氮效应参与浮游生态系统的物质循环。

在其他温带海湾, 如渤海湾(于莹等, 2018, 2022)、Hiroshima 湾(Kamiyamaet al, 1996)、胶州湾(Chenet al, 2009; Jianget al, 2011; 于莹等, 2011)、桑沟湾(于莹等, 2013)等的周年调查与本研究中利用HTS 得到的结果呈现较高的一致性: 砂壳类纤毛虫的丰度和物种组成呈现明显的季节变化, 且在春夏季达到峰值。另外, 本研究通过镜检观察发现大多为黏着壳类物种, 少数为透明壳物种, 这可能与近岸海水浑浊度高有关(Dolanet al, 2013)。

在C 站点的夏季样品中, 不仅物种多样性低而且相对丰度与其他两站点显著不同。其中旋毛纲(Spirotrichea) 相对丰度全年最低, 而寡膜纲(Oligohymenophorea)、叶口纲(Litostomatea)和叶咽纲(Phyllopharyngea)的相对丰度显著高于其他两站点。对环境因子的全年监测显示C 站点的氮磷元素浓度在夏季出现骤降, 推测氮磷等营养盐含量能够影响浮游纤毛虫的物种组成和优势种的相对丰度。

无论是对浮游纤毛虫种群的α 多样性分析还是群落的β 多样性分析均未显示出站点间差异, 原因可能是站点设置距离太近。然而在群落组成上, 砂壳目在A 站点的相对丰度最高。浮游纤毛虫的空间分布同样会受到环境因子的驱动(Yanget al,2020), 本研究中溶解氧含量可能是A 站点浮游纤毛虫群落的重要控制因素, 前人研究也显示出溶解氧对浮游纤毛虫总生物量的影响(Zhanget al,2015)。Capriulo 等(1982)的研究表明在近岸和上升流区域, 水体扰动能够使一部分沉降的砂壳类克服重力影响而分散在表层水体中。另外, 不同水团也会对近岸浮游纤毛虫群落生态特性的显著影响(Guet al, 2021)。因此我们推测由于A 站点位于码头, 船只靠岸和海浪冲击导致的水体扰动会明显影响近岸环境中砂壳类纤毛虫的水平分布模式。

除环境因子对浮游纤毛虫群落结构产生直接影响外, 生物间相互作用(捕食、寄生、共生等) (Dhibet al, 2013; Elloumiet al, 2014; Fuet al, 2020)或与环境因子的叠加效应能够共同作用于纤毛虫群落(Huanget al, 2021), 从而产生一定的时空分布模式。

3.2 传统形态学与分子生态学方法相结合

在以往对于莱州湾纤毛虫的研究中, 大多利用传统形态学的方法对特定类群进行的丰度和生物量调查。传统的镜检方法较为费时费力, 且镜检样品的体积受到时间和精力的限制(约50 mL), 往往体型较大物种(＞40 μm)容易被检获并鉴定到种水平, 并且鉴定过程需要观察者具备一定的分类学基础和经验。与传统形态学研究方法相比, 基于基因测序的生态调查能够揭示更高的生物多样性(隐性种), 目前对18S V4 区进行高通量测序已经广泛应用于真核微生物的生态调查中(Majanevaet al, 2018; Gran-Stadniczeñkoet al, 2019; Sunet al, 2019)。

本研究使用高通量测序和显微镜检方法所揭示的季节变动和群落组成基本一致, 且砂壳类纤毛虫群落的β 多样性并未受到两种方法样本量不同的影响, 这与之前利用两种方法的调查研究相似(Santoferraraet al, 2016)。另外, 高通量测序结果显示全年寡毛类(40.14%)与砂壳类(10.67%)相对丰度的比例大致为4∶1, 这与胶州湾的调查结果十分接近(寡毛类占总纤毛虫丰度的59.91%, 砂壳类占21.81%)(Jianget al, 2011), 可能暗示了温带海湾浮游纤毛虫群落具有相似的群落组成模式。

对两种方法得到的结果进行比较, 我们发现砂壳类序列的相对丰度较低(10.67%), 而镜检中的出现频率远高于寡毛类, 导致这种差异可能的原因主要有: (1) 拖网采样对较小个体的遗漏以及皮膜脆弱个体的损坏; (2) 分类学经验不足, 对体型小于40 μm个体的忽略; (3) 高通量测序仅能得到相对丰度数据,由于类群间rRNA 基因拷贝数的差异(寡毛类的rRNA基因拷贝数可能较高), 从而放大了某些物种的相对丰度(Medingeret al., 2010)。

另外, 我们对物种的分类注释结果中显示出大量未知种序列, 如寡毛类中相对丰度最高的Choreotrichida_XX 和砂壳类中的Tintinnida_XX, 由于参考序列数据库中基于形态鉴定的寡毛类和砂壳类的序列较少, 数据库不完善, 因此可能会导致OTU 分配出现错误, 降低了分类的广度和质量。另外,高通量测序不能区分活体和胞囊信息, 因此会大大高估环境中纤毛虫多样性。同时结合高通量测序和显微镜检的研究显示, 砂壳类纤毛虫群落的β 多样性并未受到两种方法样本量不同的影响(Santoferraraet al,2016), 但由于OTU 信息很可能无法完全对应形态学鉴定的物种(Santoferraraet al, 2020), 该结果仅适用于对纤毛虫群落整体的研究, 而不针对某一类群。

目前我们仍迫切需要环境中的OTU 数据来补充纤毛虫生物多样性信息(Warrenet al, 2017), 但由于分子测序方法具有一定的局限性, 例如仅能获得物种相对丰度数据、难以确定获得的多样性是真实存在的还是被人为放大, 同时测序结果也会受到测序过程和数据处理方法的影响, 因此应该与形态学方法结合使用、相互补充。未来应利用形态学方法补充完善数据库中参考序列信息, 连接基因种和形态种, 以实现更精确的分类识别。另外由于该海域相关研究的欠缺, 未来也需要进行长时间尺度和更大范围的重复采样, 以获得较为全面的纤毛虫多样性和分布信息, 并通过一定程度的标准化来连接形态学和高通量测序结果。

4 结论

本研究首次利用传统形态学和分子生态学方法相结合的手段, 对莱州湾三山岛近岸海域浮游纤毛虫进行周年生态调查。高通量测序结果显示无壳的寡毛类和有壳的砂壳类是浮游纤毛虫的主要优势类群,且具有显著的季节变化模式。寡毛类占总纤毛虫序列的40.14%, 在冬季多样性最高, 其中主要优势类群是急游虫(Strombidium)。砂壳类相对丰度为10.67%,并在春夏季出现峰值。A 站点(码头)黏着壳砂壳类的相对丰度明显多于其他站点, 推测水体扰动会影响近岸环境中砂壳类纤毛虫的物种组成和水平分布模式。然而站点间水平分布差异不显著, 原因可能是站点设置距离过近导致的。C 站点(养殖厂排水口)夏季样品中寡毛类和砂壳类纤毛虫多样性和相对丰度均低于其他站点, 可能的原因是养殖厂排放的低溶氧尾水会对浮游纤毛虫的群落结构产生影响。在全年12 个环境因子中, 总有机碳和溶解氧含量是驱动浮游纤毛虫群落结构变化的主要环境因子。通过传统形态方法共检获17 种纤毛虫, 分别隶属于7 个属, 多数为砂壳类, 其中拟铃虫属(Tintinnopsis)检获次数最高。本研究主要基于相对丰度数据, 因此某些类群的相对丰度可能会受到其他分类单元绝对丰度变化的影响。目前高通量测序和显微镜检两种方法均有不足,需要结合使用以获得更为全面的环境生态信息。