马景雪张培玉王宗兴郑明刚高 萍曲凌云王 波郑风荣*

(1.青岛大学 环境科学与工程学院,山东 青岛 266071;2.自然资源部 第一海洋研究所,山东 青岛 266061;3.自然资源部 海洋生态环境科学与工程实验室,山东 青岛 266061)

中国拥有丰富的海洋渔业资源[1],海水养殖在海洋渔业中占有重要的地位。近年来,我国水产养殖业发展迅猛,其中贝类养殖已经发展成为我国海水养殖的支柱产业[2],据渔业年鉴统计,2019年全国贝类的海水养殖面积为120.42×104万m2,占全国海水养殖总面积60.45%;2019年全国贝类海水养殖产量为1 438.97万t,占全国海水养殖总产量的69.67%,其中山东省位居全国沿海省份贝类养殖产量首位[3]。

养殖规模的增加,加速了向沿海生态系统的养分排放[4-7],使近海养殖水域的环境恶化,导致病害频发,给贝类养殖产业造成了巨大的经济损失。养殖自身污染问题成为养殖海域环境恶化的主要原因之一[8],加剧了海水富营养化程度,破坏了生物多样性,在不同浓度有机质的影响下,微生物群落的生理机能也随之发生改变[9-10],而每天过滤大量海水的滤食性动物——双壳贝类则更容易受到水中细菌群落变化的影响[11],从而诱发贝类病变。

位于青岛崂山湾海域的青山村扇贝养殖区规模庞大,其中栉孔扇贝的筏式养殖占据主要比例。该养殖区位于河口附近,河流的输入会导致近海岸呈现显著的盐度梯度,p H、温度和营养盐也发生明显变化,从而形成动态复杂的水环境系统[12]。前人对养殖区浮游植物群落结构多样性颇有研究[13-20],而对近海扇贝养殖区浮游细菌群落结构研究鲜少。2019年本研究选取养殖活动频繁的5月至10月,对该海域环境理化因子进行监测,利用基于16S rRNA基因的Illumina Miseq高通量测序方法对其细菌群落多样性及其在空间和时间尺度上的变化进行分析,并探讨了环境因子与菌群结构的相互作用,以期为该水域的水产养殖及近海海域生态系统的健康发展提供支持。

1 材料与方法

1.1 样品采集时间与地点

本研究于2019年5至10月在青岛市崂山湾青山村近海扇贝养殖海域采样,每月采集1次。采样区(图1)分为养殖密区(120°41′40″E,36°09′21″N)、养殖疏区(120°42′04″E,36°09′28″N)、对照区(120°42′24″E,36°09′35″N)三个区,分别命名为M 区、S区和D 区,在每个区域均等地选3个站位采样(间距＞500 m)。

图1 样品采集地理位置示意图Fig.1 Geographic location of sample collection

1.2 样品采集与处理

用手持式柱状采水器在指定点采取1 L表层水,在低温(4 ℃)条件下运回实验室;采用0.22μm 的滤膜负压过滤300 m L水样,将处理好的滤膜用液氮速冻,放-80 ℃冰箱保存。每份水样过滤重复3次。

1.2.1 环境因子的测量

本研究采用WZB-170型便携式浊度计测量水样浊度,JENCO 9030M 便携式溶氧(S)测定仪测量溶氧(Dissolved Oxygen,DO)质量浓度和温度(θ),556MPS 便携式多参数水质测量仪测量盐度,TURNER Designs Triligy实验室荧光仪测量叶绿素a(Chl-a)质量浓度,其他参数在实验室中均用标准方法测量。

1.2.2 流式细胞术测定细菌丰度

各站位取4 m L海水至冻存管中,随后加入444μL戊二醛固定,在低温避光的条件下带回实验室,液氮处理后,-80 ℃保存。采用细胞流式仪进样前,取180μL 样品和20μLSYBRTMGreenⅠnucleic acid gelstain(invitrogen,USA)染料混匀,染色完成后,上机测样。使用低流速(14μL/min),激发电压为SSC 450、FL1 500、FL2 510和FL3 661。

1.2.3 Illumina Miseq高通量测序

用DNeasy Powerwater Pro Kit试剂盒(QIAGEN,USA)提取DNA,用1%的琼脂糖凝胶检测提取的DNA 质量。分别用细菌通用引物338F-806R[21]和弧菌特异性引物V_169F-V_680R[22-23]两对引物对DNA 进行PCR 扩增(表1)。将同一样本的PCR 产物混合后使用2%琼脂糖凝胶回收PCR产物,利用Axy Prep DNA Gel Extraction Kit(Axygen Biosciences,Union City,CA,USA)回收产物纯化,2%琼脂糖凝胶电泳检测,并用QuantusTMFluorometer(Promega,USA)对回收产物进行检测定量。使用NEXTFLEX Rapid DNA-Seq Kit建库。

表1 引物序列Table 1 Primer sequences in this study

1.2.4 数据处理

利用均值法处理环境因子的数据及细菌细胞的丰度;使用Fastp软件[24]和FLASH 软件[25]分别对原始测序序列进行质控和拼接,利用UPARSE 软件根据97%的相似度对序列进行OTU 聚类并剔除嵌合体[26-27]。利用RDP classifier[28]对每条序列进行物种分类注释。利用R 语言(version 3.3.1)进行物种组成分析、环境因子与群落结构的相关性分析及冗余分析。利用SPSS 22.0 的单因素方差分析(one-way Analysis of Variance,one-way ANOVA)进行组间差异显著性检验。

2 结果

2.1 环境因子理化指标

养殖环境是水产养殖中一个不可忽视的部分,养殖用水的理化性质是养殖活动的重要保证,表2展示了崂山湾青山村近海养殖区在2019年5月至10月的环境理化指标变化,养殖水质符合《海水水质标准》(GB 3097-1997)[29]第二类水质标准要求。

表2 2019年崂山湾青山村扇贝养殖区理化性质指标Table 2 Physical and chemical properties of scallop culture area in Qingshan Village,Laoshan District from May to October in 2019

续表

2.2 流式细胞仪检测细菌丰度结果

从海水中细菌单细胞的集聚分布情况(图2)所示,侧向散射光(Side Scattering,SSC)与细胞的内容物含量有关,FL表示颗粒被染上的荧光的数量。通过SSC对FL1和FL3对FL1双参数图对细菌细胞进行区分和计数。R1和R2分别为在不同的荧光强度下根据聚集情况划分的细菌细胞门。由流式细胞仪计数结果(图2)可知,每毫升海水的细胞数的变化呈现先上升后下降的趋势。养殖密区及疏区均在9月达到最高点,对照区的最高点出现在8月,其最低点均出现在5月(图3)。

图2 细胞流式图Fig.2 Cell flow diagram

图3 2019年5月至10月养殖区海域细菌细胞数量Fig.3 The number of bacterial cells in the aquaculture area from May to October in 2019

2.3 水体菌群多样性及其对环境因子响应

2.3.1 微生物Alpha多样性分析

通过聚类操作,根据97%的相似度水平,对所有序列进行OTU 划分,青山村扇贝养殖区及对照区细菌的香农-威纳指数(Shannon-Weiner Index)呈现波动的变化趋势,养殖疏区呈现先升高后下降的趋势。弧菌的香农-威纳指数均呈现M 型变化。2019年9月,细菌多样性升高而弧菌多样性降低,且细菌细胞数也达到峰值(图4a和4b)。

图4 2019年5月至10月水环境细菌和弧菌多样Fig.4 The diversities of bacteria and Vibrio in water environment from May to October in 2019

2.3.2 水环境群落组成

本研究共检测出52 个细菌门,其中存在11 个优势门(丰度占比大于1%),分别为变形菌门(Proteobacteria)、蓝细菌门(Cyanobacteria)、拟杆菌门(Bacteroidetes)、放线菌门(Actinobacteria)、厚壁菌门(Firmicutes)、Patescibacteria、疣微菌门(Verrucomicrobia)、Marinimicrobia_SAR406-clade、绿弯菌门(Chloroflexi)、Epsilonbacteraeota、酸杆菌门(Acidobacteria),其中变形菌门(Proteobacteria)为优势菌门(图5)。

图5 细菌门水平群落结构组成Fig.5 Compostion of community structure at bacterial phylum level

在纲水平,共检测出115个纲,其中存在13个优势纲(丰度占比大于1%),相对丰度较高的是γ-变形菌纲(Gammaproteobacteria)、α-变形菌纲(Alphaproteobacteria)、拟杆菌纲(Bacteroidia)和放线菌纲(Actinobacteria)(图6)。

图6 细菌纲水平群落结构组成Fig.6 Composition of community structure at bacterial class level

弧菌属在细菌属水平丰度中占据重要的比例(图7),养殖密区在2019年6月和7月出现较高的丰度比例,对照区在2019年7月的丰度比例最高。

图7 养殖区弧菌丰度比例分布Fig.7 Diagram of abundance proportion of Vibrio in culture area from May to October in 2019

对各区不同月份的样品及密区、疏区和对照区的样品进行优势门的组间差异显著性检验分析,结果表明,各区样品在不同月份Bacteroidetes、Verrucomicrobia、Marinimicrobia_SAR406_clade、Chloroflexi、Acidobacteria具有差异(P＜0.05)或显著差异(P＜0.01),而密区、疏区和对照区的样品在门水平上微生物组成没有显著差异(图8)。

图8 组间差异显著分析Fig.8 Analysis of significant differences between groups from May to October in 2019

2.3.3 水环境群落组成与环境因子的相关性

RDA 分析揭示了不同区域样品中微生物组成与环境因子的相关性,采用bioenv函数检验了每个环境因子的显著性,在养殖密区、养殖疏区和对照区的RDA 分析中,第一坐标轴和第二坐标轴共贡献解释量分别为59.25%、47.1%和45.27%。结果表明,浊度(P=0.027)对养殖密区OTU 水平的物种组成和群落结构具有显著的影响;硝酸盐(P=0.017)、浊度(P=0.033)对养殖疏区OTU 水平的物种组成和群落结构具有显著的影响,磷酸盐(P=0.004)、DO(P=0.001)、p H(P=0.004)对其具有显著的影响;DO(P=0.018)对对照区OTU 水平的物种组成和群落结构具有显著的影响(图9)。

图9 细菌OTU 水平环境因子RDA 分析Fig.9 RDA analysis on the correlation between environmental factors and bacterial OTU level

通过Spearman 相关性计算,利用求平均的方式对环境因子和物种层级进行聚类,选取门水平上丰度为前11位的物种进行分析,相关性Heatmap图和环境因子聚类树分析结果表明:

1)在门水平上,养殖密区的细菌群落物种与环境因子均具有不同程度的相关性,其中,拟杆菌门与盐度(Spearman相关系数r=-0.776,P＜0.001)呈现极显著负相关关系;绿弯菌门与浊度(Spearman相关系数r=0.744,P＜0.001)呈现极显著正相关关系。养殖疏区的放线菌门与浊度呈现极显著正相关关系(Spearman相关系数r=0.808,P＜0.001);Marinimicrobia_SAR406_clade与盐度呈现极显著正相关关系(Spearman相关系数r=0.856,P＜0.001),厚壁菌门分别与DIN(Spearman相关系数r=0.760,P＜0.001)和硝酸盐(Spearman相关系数r=0.751,P＜0.001)呈现极显著正相关关系(图10和图11)。

图10 养殖密区细菌门水平环境因子相关性Heatmap图Fig.10 Heatmap diagram of correlation between bacterial phylum level and environmental factors in intensive culture area

图11 养殖疏区细菌门水平环境因子相关性Heatmap图Fig.11 Heatmap diagram of correlation between bacterial phylum level and environmental factors in sparse culture area

2)在门水平上,对照区的物种与环境因子具有不同程度的相关性,与养殖区不同,其中,酸杆菌门和厚壁菌门与叶绿素分别呈现显著负相关关系(Spearman相关系数r=-0.612,P＜0.01)和极显著正相关关系(Spearman相关系数r=0.736,P=0.001),Patescibacteria与温度呈现显著正相关关系(Spearman相关系数r=0.650,P＜0.01),绿弯菌门和拟杆菌门分别与浊度具有显著正相关(Spearman相关系数r=0.642,P＜0.01)和极显著负相关关系(Spearman相关系数r=-0.722,P＜0.001),放线菌门与溶解氧呈现显著正相关关系(Spearman相关系数r=0.607,P＜0.01)(图12)。

图12 对照区细菌门水平环境因子相关性Heatmap图Fig.12 Heatmap diagram of correlation between bacterial phylum level and environmental factors in control area

3 讨论

3.1 海水水样细菌群落多样性

通过Illumina Miseq高通量测序的方法对养殖区海水水样进行门、纲水平的群落结构分析,结果表明,在水样样品中共检测到52个门、115个纲,香农-威纳指数为4.762～5.539,Sobs指数为676～2 140。2019年8月和9月的细菌多样性指数较高,香浓指数总体上随温度的变化而变化,但2019年7月的养殖密区和对照区的多样性指数相对低于趋势线,这可能是由在该月份取样时出现连续的雷雨天气所致。根据图5可知,不同区域不同月份的样品中的优势群落均为变形菌门(Proteobacteria),占总群落50%以上,主要包括α-Proteobacteria和γ-Proteobacteria。其他的优势菌门排列依次是蓝细菌门(Cyanobacteria)、拟杆菌门(Bacteroidetes)、放线菌门(Actinobacteria)、厚壁菌门(Firmicutes)、Patescibacteria、疣微菌门(Verrucomicrobia)、Marinimicrobia、绿弯菌门(Chloroflexi)、Epsilonbacteraeota、酸杆菌门(Acidobacteria)。检测结果与近年来的研究结果相似,杨铭玉等[30]对东海舟山渔场近岸细菌群落结构分析得出,该养殖区主要的细菌群落类群为变形菌(62.99%)、拟杆菌群(10.23%)、蓝藻细菌(9.69%)、放线菌门(8.94%)。辽东湾沿海的仿刺参、对虾混养养殖水体中主要是变形菌门(51.5%),其次是放线菌门(23.6%)、拟杆菌门(11.3%)、厚壁菌门(9.1%)[31]。王飞飞[32]利用高通量测序、分子生物学方法等手段对黄海北部海域贝类养殖区进行细菌群落分析,也证明养殖区的海水中变形菌门占绝对优势。

在纲水平中,γ-变形菌纲在群落中占据了主导地位,且养殖区丰度占比高于对照区,其中包含许多重要的病原菌,如沙门氏菌属(Salmonella)、耶尔辛氏菌属(Yersinia)、弧菌属(Vibrio),这些菌群会引起不同程度的生物感染,严重者会导致大面积流行病的爆发。变形菌门中σ-变形菌纲的比例相对较小,具有硫氧化、固碳、C1利用和异养能力[33],在2019年8、9和10月份活动较为频繁。SAR406在2019年9和10月份占比较高,该菌门参与厌氧氨氧化与磷循环的过程[34-35],且能够影响无机碳的固定[36]。

在属水平,弧菌属是我们重点关注的对象,研究发现2019年6月的弧菌属比例在养殖区处在最高点,弧菌多样性也处于较高值,说明单一弧菌种的数量较低,且细菌细胞数量也较低,因此患病的几率大幅降低。2019年7月弧菌的相对比例较高,细菌总数下降,弧菌多样性降低,因此单一弧菌种数量增加。2019年7月、8月和9月处于高温时期,适宜弧菌的快速繁殖,同时会降低扇贝的活性。因此该时期易爆发弧菌病。2019年9月份虽然弧菌属占比不高,但是细胞总数处于峰值,且细菌多样性增大而弧菌多样性减少。因此2019年9月相较于其他高温月份(2019年7月和8月)更具有爆发弧菌病的条件。邓欢等[37]也曾报道过在2019年9月份大连爆发了扇贝弧菌病,造成了扇贝群体突发性死亡。

对各区不同月份的样品及密、疏、对照区的样品进行优势门的组间差异显著性检验分析,发现养殖活动对海水中菌群优势门的影响不大,组间没有明显差异,表明从养殖区到海区,浮游细菌群落优势门季节变化大于空间尺度差异。Lü 等[38]研究表明黄河三角洲演替潮滩的细菌群落结构存在季节性差异。而MA等[39]认为黄海沿岸由于河口的存在,从河流到深海,浮游细菌群落的空间变化大于季节尺度。

3.2 环境因子相关性分析

环境因子对微生物群落结构的影响一直是学者们关注的热点。不同的地理位置、环境因子与海域浮游细菌的群落结构组成之间会进行相互作用。在黄海沿岸,盐度是导致微生物群落结构向河流、河口两部分空间分化的主要因素[39]。而在黄河三角洲附近,温度和亚硝酸盐构成菌群结构季节变化的主要驱动力[38]。王彩霞等[9]对渤海三湾环境因子与菌群结构进行Heatmap相关性分析得出,温度、磷酸盐、硝酸盐和铵盐与门水平上物种分布有显著正相关关系(P＜0.01),而溶解氧、纬度和盐度与门水平上物种组成有明显的负相关关系(P＜0.05)。养殖活动是否会影响环境因子与近海细菌群落结构的相互作用,大家鲜少关注。通过对环境因子与菌群结构的Heatmap相关性分析,我们发现对照区与养殖区所呈现的相关性结果之间还是有明显的差异。在门水平上,养殖区与对照区的相同的菌群所受到的不同环境因子影响,不但表现在相关性的程度上,甚至还改变了其作用方向,例如疣微菌门,在养密区与硝酸盐呈现负相关关系,而在对照区与其呈现微弱的正相关关系。绿弯菌门在养殖密区与温度(Spearman相关系数r=0.636,P＜0.01)呈现显著正相关关系,养殖疏区与温度(Spearman相关系数r=0.678,P＜0.01)呈现显著正相关关系,而在对照区绿弯菌门与温度(Spearman相关系数r=0.555,P＜0.05)的相关性区域性降低。绿弯菌门在海洋生态系统中参与光合作用的过程,利用3-羟基丙酸途径来固定二氧化碳,养殖区中的绿弯菌门更易受温度升高而增加,加强光合作用,释放更多的溶解氧,而扇贝是高耗氧生物,集聚性养殖会消耗大量溶解氧,绿湾菌门在养殖区与对照区中的相关性差异结果也证实养殖活动会影响环境因子与菌落结构的相关程度,这与先前的研究结果一致[10]。研究近海养殖区海域微生物群落结构多样性和环境因子相关性为保护近海养殖生态环境可持续发展以及扇贝健康养殖提供了更好的理论基础。

4 结论

本文利用细胞流式仪和高通量测序技术,对青岛崂山湾近海扇贝养殖水样进行多层次的分析,集中分析扇贝高发病月份的细菌群落变化情况,并综合考虑环境因子的变化,对样品菌群结构与环境因子相关性进行了分析,得到的结论如下:

1)青岛崂山湾近海扇贝养殖水样中的优势菌门为变形菌门,其中,γ-变形菌纲在群落菌门中占据主导地位,养殖区γ-变形菌纲丰度占比高于对照区,由于γ-变形菌纲中含有多种致病菌,丰度占比增长的同时也增加了养殖区扇贝感染流行病的潜在风险。9月处于高温时期,细菌总数增大,养殖区的细胞多样性上升且弧菌多样性指数降低,使得扇贝在9月爆发弧菌病的可能性骤增。

2)在该近海区,养殖活动对海水细菌群落结构的影响主要以季节变化为主,在门水平上,浊度对养殖密区的菌群结构具有显著影响,DO 对对照区的菌群结构具有显著影响。

3)养殖活动不仅影响了环境因子与菌群结构相关性的程度,甚至还会改变其相互作用的方向。

研究结果表明,青岛崂山湾近海扇贝养殖活动影响了正常的近海菌落的结构分布,集聚式养殖不同程度上增加了某些病原菌的丰度,增加了该区近海生物疾病爆发的风险,因此在近海养殖活动中要做到合理布局、科学喂养、适当降低养殖密度,这些措施对降低养殖活动对菌落结构多样性的影响具有重要意义。本文的研究结果可为实现近海科学养殖提供理论依据。