鳗弧菌和迟钝爱德华氏菌对短蛸肝脏菌群结构的影响

朱兴华，王金龙，袁廷柱，冯艳微，李 赞，徐晓辉，王卫军，孙国华，杨建敏

( 1.上海海洋大学，水产科学国家级实验教学示范中心，上海 201306； 2.山东省海洋生态资源重点修复实验室，山东省海洋资源与环境研究院，山东 烟台 264006； 3.上海海洋大学，上海水产养殖工程技术研究中心，上海 201306； 4.长岛综合试验区海洋经济促进中心，山东 烟台 265800； 5.鲁东大学 农学院，山东 烟台 264025 )

水产养殖中细菌性疾病普遍存在，而且危害极大。弧菌病和爱德华氏菌病是养殖动物常患的细菌性疾病。鳗弧菌(Vibrioanguillarum)是一种能够感染贝类和鱼类等水产动物的条件性致病菌，其通过损伤皮肤继而感染水产动物，发病快、传染性强、死亡率高[1]；迟钝爱德华氏菌(Edwardiantarda)是对水产动物最具危害的病原菌之一，其含有大量致病因子，可经血液传播到各组织器官，引起感染，严重时可导致个体死亡。

短蛸(Octopusocellatus)主要栖居于温带偏北海域，是我国北方沿海渔业重要的经济物种之一[2]，具有较高的经济价值和药用价值，市场前景广阔[3]。目前，短蛸的人工繁育和养殖处于初步阶段，尚未形成一定的规模化[4-5]。在养殖过程中易感染弧菌属、爱德华氏菌属、气单胞菌属(Aeromonas)和假单胞菌属(Pseudomonas)等细菌[6-7]，给短蛸养殖业带来了巨大挑战。研究鳗弧菌、爱德华氏菌等细菌对短蛸的致病机理，对养殖过程中病原菌的诊断与防治具有重要意义。

肝脏(肝胰腺)作为软体动物中重要的代谢和免疫相关组织，在营养摄取、激素合成和碳氮代谢中起主要作用，同时参与病原体清除、抗原加工以及病原菌感染宿主所引起的代谢变化[8]。肝脏在头足类中发育良好，占胸腔的绝大部分[9]。研究发现，当短蛸感染了鳗弧菌和副溶血弧菌(V.parahaemolyticus)后，2种弧菌均可破坏短蛸肝脏的组织结构和细胞结构，导致肝功能损害，加速短蛸死亡[10]。微生物及其代谢物质在维持宿主生理生命活动稳态方面具有重要作用，它形成了一个动态平衡系统，可抵抗病原菌的侵染，有利于机体正常生命活动的进行[11]。当病原菌入侵宿主时，会抑制或损害宿主共生微生物群落中的某些成员，对微生物菌群结构产生一定的影响[12]，导致生物体的自我调节能力发生紊乱，产生病变。病原菌对宿主微生物菌群影响方面的研究主要集中在肠道组织[13-16]方面。而无脊椎动物的肝脏中也存在着非致病性的共生菌，对宿主的抗病性和免疫功能有一定的影响[8]。目前病原菌对肝脏菌群的影响在文蛤(Meretrixmeretrix)中有过相关研究[8]，而在短蛸中尚未见报道。

笔者利用16S rRNA基因扩增子技术检测短蛸肝脏的菌群结构，分析鳗弧菌和迟钝爱德华氏菌感染后肝脏菌群结构的变化，以期为今后短蛸养殖过程中细菌性疾病的准确诊断与防控提供基础资料。

1 材料与方法

1.1 试验材料

短蛸于2019年11月采自山东烟台东方海洋基地附近海域，共66只，体质量(60±5) g，体长(15±2) cm，暂养于70 cm×60 cm×50 cm恒温水箱中，水温21 ℃，每日更换新鲜海水，空腹7 d待稳定后进行试验。

1.2 病原菌分离与培养

取养殖水样通过2216E固体培养基进行细菌培养，挑取单克隆扩大培养后进行16S rRNA基因测序和菌种鉴定，得到鳗弧菌和迟钝爱德华氏菌2种菌株。将2种菌株分别置于2216E液体培养基中培养24 h，细菌沉淀离心后用磷酸缓冲盐溶液冲洗3次，最后稀释成1.01×1010cfu/mL的鳗弧菌悬液和1.47×109cfu/mL的迟钝爱德华氏菌悬液，备用。

1.3 注菌及取材

试验分为A、B、C组，每组15只短蛸。A组每只腹腔注射鳗弧菌悬液200 μL，B组每只腹腔注射等剂量的迟钝爱德华氏菌悬液，C组每只腹腔注射等剂量磷酸缓冲盐溶液作为对照。每日观察各组短蛸的行为和生理变化，记录死亡情况，找到半致死的时间点。

另设鳗弧菌组和迟钝爱德华氏菌组，每组8只短蛸，分别在腹腔注射鳗弧菌悬液和迟钝爱德华氏菌悬液，达到半致死时间点时，各组随机挑选存活短蛸3只，取肝脏组织放入-80 ℃冰箱中保存。对照组无死亡，随机取3只短蛸的肝脏组织保存备用。

1.4 16S rRNA基因V3～V4高变区扩增

利用CTAB法提取肝脏组织DNA。选择16S rRNA基因的V3～V4区域作为扩增和测序的目的片段，引物序列为：

341F：CCCTACACGACGCTCTTCCGATCTG；

805R：GACTGGAGTTCCTTGGCACCCGAG AATTCCA。

合成带有标签序列的特异性引物。PCR扩增共两轮，第一轮反应体系为20 μL：2×Taq master Mix 15 μL，上、下游引物各1 μL，模板DNA取10～20 ng，加无菌水定容至20 μL。PCR反应条件：94 ℃预变性3 min；94 ℃变性10 s，55 ℃退火30 s，72 ℃延伸30 s，25个循环；72 ℃延伸5 min。第二轮的反应体系为30 μL：2×Taq master Mix 15 μL，上、下游引物各1 μL，上一轮的PCR产物20 ng，加无菌水定容至30 μL。PCR反应条件：94 ℃预变性3 min；94 ℃变性20 s，55 ℃退火20 s，72 ℃延伸30 s，5个循环；72 ℃延伸5 min。PCR产物用2.0%琼脂糖凝胶进行电泳，选择长度在400～450 bp的条带，使用GeneJET胶回收试剂盒进行回收纯化。

1.5 文库构建及Illumina MiSeq测序

使用NEB DNA建库试剂盒进行DNA文库构建，构建文库完成经Qubit定量和检验合格后，采用Illumina MiSeq测序平台进行上机测序。

1.6 数据处理及分析

将测序得到的原始数据剔除非特异性扩增序列及嵌合体之后获得有效序列；按照有效序列之间的距离进行聚类，根据序列间的相似性划分不同的运算分类单元(OTU)，通常情况对97%相似水平下的运算分类单元进行生物信息统计分析；利用R的Venn Diagram功能制作韦恩图展现样品的运算分类单元数目的相似性和组间重叠情况。使用Mothur软件计算每组样品的Alpha多样性指数；利用R软件制作稀疏性曲线图。利用RDP Classifier贝叶斯算法对运算分类单元进行物种分类，统计每个样品在不同分类层级水平上的群落组成，利用R进行作图；利用PICRUSt软件预测短蛸肝脏菌群功能。

2 结 果

2.1 序列分析

测序共产生217 287条质控后序列，其中对照组、鳗弧菌组和迟钝爱德华氏菌组的质控后序列分别为78 155、59 159和79 973条，序列长度为421.80～428.53 bp，除去嵌合体以及非特异性扩增序列后，对照组、鳗弧菌组和迟钝爱德华氏菌组的剩余序列数目分别为74 082、58 759和72 464条(表1)。

表1 样品测序数据

2.2 运算分类单元分布韦恩图

样品间运算分类单元数目组成相似性及其重叠情况见图1。3组共有的运算分类单元为141个，对照组特有的运算分类单元为1834个，鳗弧菌组特有的运算分类单元为870个，迟钝爱德华氏菌组特有的运算分类单元为2207个。对照组与迟钝爱德华氏菌组共有的运算分类单元为373个，对照组与鳗弧菌组共有的运算分类单元为166个，迟钝爱德华氏菌组和鳗弧菌组共有的运算分类单元为200个。

图1 运算分类单元样品分布韦恩图Fig.1 Venn diagram of OTUs distribution among samples不同样品用不同颜色表示，图中数字代表特异或共有的运算分类单元数.Different samples are shown in different colors. The numbers in the diagram represent specific or common OTUs.

2.3 多样性指数分析

Alpha指数可以反映微生物群落的多样性和丰度，其中Chao1指数预估群落分布丰度，香农指数与辛普森指数预估样品中微生物的多样性。对照组的Chao1指数为21 460.22，迟钝爱德华氏菌组为9196.04，鳗弧菌组为5550.49，表明对照组群落丰度最大，鳗弧菌组群落丰度最小；香农指数与辛普森指数，对照组为4.78、0.03，迟钝爱德华氏菌组为4.09、0.09，鳗弧菌组为0.45、0.91，表明对照组的群落多样性最高，其次是迟钝爱德华氏菌组，鳗弧菌组最低(表2)。3组样品的文库覆盖率达到0.98，以香农指数构建的稀释性曲线趋于平缓(图2)，说明各组样品测序数据数量合理，可以较为全面地体现样本中微生物菌群的组成情况。

表2 Alpha多样性指数统计

图2 各组样品稀释性曲线Fig.2 Dilution curve of samples in each group

2.4 菌群结构组成分析

在门水平上3组样本的菌群结构组成见图3。对照组的核心菌门(相对丰度>1%)是变形菌门(46.2%)、拟杆菌门(13.32%)、软壁菌门(9.03%)、厚壁菌门(7.9%)、酸杆菌门(5.76%)、浮霉菌门(3%)和放线菌门(2.62%)；迟钝爱德华氏菌组的核心菌门为变形菌门(43.33%)、软壁菌门(31.33%)、厚壁菌门(16.5%)和放线菌门(2.59%)；鳗弧菌组的核心菌门是软壁菌门(95.95%)和变形菌门(2.2%)。对照组的变形菌门、拟杆菌门、酸杆菌门等相对丰度高于其他2组；而迟钝爱德华氏菌组的变形菌门相对丰度高于鳗弧菌组，厚壁菌门相对丰度高于其他2组；鳗弧菌组只有软壁菌门相对丰度远高于其他2组(图4)。

图3 门水平上样本的群落结构分布Fig.3 The distribution diagram of the community structure of all samples at the phylum level.

图4 门水平上样本相对丰度对比Fig.4 Comparison of relative abundance of samples at phylum level

属水平上3组样本的菌群结构组成见图5。对照组的核心菌属由支原体属(Mycoplasma,9.01%)、短波单胞菌属(Brevundimonas,5.74%)、Gp4(5.16%)、假单胞菌属(4.26%)、硫氧化湖沉积杆菌属(Limnobacter，2.29%)、甲基杆菌属(Methylobacterium,1.73%)和狭义梭菌属(Clostridiumsensustricto,1.07%)构成；迟钝爱德华氏菌组中丰度较高的属依次为支原体属(31.28%)、乳球菌属(Lactococcus, 9.49%)、肠杆菌属(Enterobacter, 6.07%)、短波单胞菌属(3.71%)、假单胞菌属(Pseudomonas, 3.44%)、克雷伯氏菌属(Klebsiella, 2.06%)；鳗弧菌组核心菌属为支原体属(95.94%)。迟钝爱德华氏菌组与鳗弧菌组均未发现Gp4属。对照组未分类占37.4%，迟钝爱德华氏菌组未分类占12.7%，鳗弧菌组未分类占1.15%。对照组唯有未分类菌属相对丰度明显高于其他2组；迟钝爱德华氏菌组的乳球菌属以及肠杆菌属相对丰度高于其他2组；而鳗弧菌组的支原体属相对丰度远远超过其他2组(图6)。

图5 属水平上所有样本群落结构分布Fig.5 The distribution diagram of the community structure of all samples at the genus level

图6 属水平上样本相对丰度对比Fig.6 Comparison of relative abundance of samples at genus level

2.5 菌群分布组成差异

根据3组样本中丰度最高的前50个物种的分类信息绘制的聚类分析见图7。对照组和迟钝爱德华氏菌组的物种丰度较接近，鳗弧菌组差异较大。对照组与迟钝爱德华氏菌组聚成一支，鳗弧菌组单独聚成一支，说明对照组和迟钝爱德华氏菌组菌群分布较类似，距离较近。

图7 属水平物种丰度聚类分析热图Fig.7 The heat map of species abundance clustering at the genus level

2.6 菌群功能预测

PICRUSt分析表明，对照组预测的功能基因可注释到KEGG数据库中4条第一层级功能通路以及17条第二层级功能通路，其中包括：新陈代谢(相对丰度为49.86%)类通路中的氨基酸代谢(10.35%)、碳水化合物代谢(9.89%)、能量代谢(5.96%)、辅助因子和维生素代谢(4.22%)、脂质代谢(3.71%)、核苷酸代谢(3.42%)、外源生物降解和代谢(3.16%)、聚糖的生物合成和代谢(2.18%)、萜类化合物与聚酮类化合物的代谢(2.10%)、其他氨基酸代谢(1.90%)等；遗传信息处理(17.09%)类通路中的复制与修复(7.50%)，翻译(4.79%)，折叠、分选与降解 (2.43%)，转录(2.38%)；环境信息处理(12.47%)类通路中的膜转运(9.98%)、信号转导(2.31%)；细胞进程(4.43%)类通路中的细胞运动(3.49%)(图8)。

图8 基于KEGG的功能结构分布Fig.8 The histogram of function structure distribution based on KEGG

鳗弧菌组与迟钝爱德华氏菌组注释到了与对照组相同的第一层级与第二层级功能通路，但各功能丰度普遍低于对照组正常水平。其中鳗弧菌组细胞进程类通路的细胞运动功能丰度为172 676，与对照组对比下降92.62%；新陈代谢类通路的外源生物降解和代谢功能丰度为252 854，下降88.06%，氨基酸代谢功能丰度下降79.64%，环境信息处理类通路的信号转导功能丰度下降92.78%，遗传信息处理类通路中的转录功能丰度下降69.74%；迟钝爱德华氏菌组新陈代谢类通路的氨基酸代谢功能丰度下降25.36%，细胞进程类通路的细胞运动功能丰度下降26.09%(图9)。

图9 基于KEGG的功能组成和丰度对比Fig.9 The histogram of function composition and abundance comparison based on KEGG

3 讨 论

随着分子生物学技术的不断发展，以16S rRNA基因扩增子测序技术为代表的高通量测序逐渐取代了传统的分子生物学分析手段，广泛应用到各种组织菌群的研究[17-18]，为全面探索微生物菌群结构及群落差异指明了新方向。笔者首次采用16S rRNA基因扩增子测序法探究短蛸肝脏菌群多样性与功能，并对迟钝爱德华氏菌和鳗弧菌感染后肝脏菌群的结构变化进行了分析，以期为短蛸养殖过程中细菌性疾病的准确诊断与防控提供基础资料。

3.1 菌群丰度与多样性影响

试验结果显示：迟钝爱德华氏菌组与对照组共有的运算分类单元为373个，鳗弧菌组与对照组共有的运算分类单元为166个；聚类分析热图中，对照组与迟钝爱德华氏菌组聚成一支，鳗弧菌组单独聚成一支。由此可知，迟钝爱德华氏菌组与对照组的菌群组成较为相近，而鳗弧菌组与其他两组差异较大，表明鳗弧菌相比爱德华氏菌，对肝脏菌群产生了更大的影响。对照组肝脏菌群的丰度和多样性最高，迟钝爱德华氏菌组次之，鳗弧菌组最少。这种变化可能是由于病菌的入侵改变了细菌群落的组成与种间相互作用，从而改变微生物介导的功能，影响了宿主的免疫应答，宿主的自我调节能力紊乱，发生病变。细菌感染后，宿主共生菌多样性降低与宿主健康受损的相关性在许多研究中均有发现[19-20]。

3.2 不同水平上优势菌群影响

门水平上细菌分类结果显示，变形菌门、厚壁菌门和拟杆菌门是短蛸肝脏的核心菌门。这与文蛤肝脏中的优势菌群组成相似，文蛤肝脏的优势菌群依次为放线菌门、拟杆菌门、厚壁菌门、变形菌门[8]。以上细菌同时也是白甲鱼(Onychostomasima)、凡纳滨对虾(Litopenaeusvannamei)、拟穴青蟹(Scyllaparamamosain)等多种水产动物肠道中的优势菌群[21-23]。本试验中，变形菌门是短蛸肝脏中最优势的细菌类群，其中包含了γ-变形菌纲的短波单胞菌属、假单胞菌属，α-变形菌纲的甲基杆菌属，β-变形菌纲的硫氧化湖沉积杆菌属等。此外，厚壁菌门的狭义梭菌属也是短蛸肝脏的优势菌属。短波单胞菌属和假单胞菌属常被报道为各种对虾肠道中共有的优势菌属[24]；狭义梭菌属被报道是鳢属(Channa)鱼类和团头鲂(Megalobramaamblycephala)等肠道的优势菌属[25-26]。说明短蛸肝脏菌群和以上水生动物肠道菌群存在一定的共性。

属水平上细菌分类结果显示，对照组的核心菌属中支原体属占9.01%，经过菌感染后，试验组支原体属丰度激增，在鳗弧菌组中占95.94%。在弧菌感染文蛤的试验中同样出现了支原体激增的状况[8]。前期研究发现，短蛸在注射鳗弧菌后，肝脏组织出现了部分细胞膜溶解，细胞界限模糊，脂肪颗粒及细胞质外溢，甚至细胞核和部分细胞器破裂溶解的现象[13]。笔者推测，原因可能是激增的支原体通过黏附作用与宿主细胞结合，获取细胞膜中的脂质和胆固醇，破坏细胞膜结构。

3.3 PICRUSt功能预测影响

PICRUSt预测结果显示，短蛸肝脏菌群主要的功能与新陈代谢类功能有关，包括氨基酸代谢、碳水化合物代谢、能量代谢、核苷酸代谢、脂质代谢等，表明肝脏菌群能够参与短蛸的代谢活动。狭义梭菌属是短蛸肝脏的优势菌属。研究表明，狭义梭菌属细菌可以利用发酵糖类和碳水化合物，提供机体所需能量，促进生长、增强免疫功能和抗癌功能[27-28]。本试验结果表明，肝脏菌群在短蛸的生理活动中发挥着十分重要的作用。注射病菌后，各代谢通路以及细胞其他基础生命活动功能通路丰度都有所降低，表明当菌群结构遭到破坏后，肝脏的生理机能也相应受损。

4 结 论

短蛸肝脏具有稳定的菌群结构，在短蛸的日常代谢活动中发挥一定的作用。当受到鳗弧菌与爱德华氏菌侵染后，肝脏菌群的多样性与丰度均发生变化。与对照组相比，爱德华氏菌组菌群结构变化较小，鳗弧菌组变化较大，表明鳗弧菌对短蛸肝脏的菌群结构破坏更大。研究结果可丰富和拓展短蛸肝脏菌群的基础研究，并为短蛸养殖过程中细菌性疾病的准确诊断与防控提供基础资料。