多棘海盘车(Asterias amurensis)是掠食性海洋底栖无脊椎动物, 属于棘皮动物门(Echinodermata)、海星纲(Asteroidea)、钳棘目(Forcipulatida)、海盘车科(Asteriidae)、海盘车属(Asterias)。多棘海盘车是北太平洋的常见种, 垂直分布于潮间带至200 m水深海域, 在中国渤海、黄海、美国阿拉斯加和阿留申群岛、北冰洋、白令海峡、鄂霍茨克海、日本海、鞑靼海峡均有分布(徐思嘉等, 2018)。同时, 在1986年多棘海盘车首次被记录为入侵种出现在澳大利亚塔斯马尼亚(Tasmania)岛, 并确定建立种群(Buttermore et al, 1994; Ward et al, 1995)。

多棘海盘车有“海底蝗虫”之称, 在包括我国和日本在内的北太平洋海域, 以及澳大利亚塔斯马尼亚海域等多个海域发生暴发事件, 造成了重大的经济损失和生态灾害。1954年日本东京湾仅因多棘海盘车摄食养殖贝类导致的经济损失高达4亿多日元(Kim, 1968); 1974年日本鸟取县发生多棘海盘车暴发性聚集, 导致增殖放养的东风螺几乎被捕食(孙光, 1983); 七年之后, 多棘海盘车再次在日本聚集, 对当地放流的魁蚶(Scapharca broughtonii)造成重大损失(孙光, 1983)。多棘海盘车在澳大利亚海域的首个记录发生在1986年10月的塔斯马尼亚岛的德温特河河口, 自此便持续发生因其聚集危害当地贝类养殖的事件(Buttermore et al, 1994; Ward et al, 1995), 被认为是本地海洋生物的严重害虫。它与塔斯马尼亚极度濒危鱼类粗体澳洲躄鱼(Brachionichthys hirsutus)急剧减少有关(Stevens, 2022)。自2000年起, 随着我国北方养殖业的发展, 海星暴发问题日趋严重, 严重危害当地养殖业。2006年7月多棘海盘车暴发, 密度高达300个/m² (周书珩等, 2008)。2007年3月胶州湾的16万亩菲律宾蛤仔(Ruditapes philippinarum)中60%遭到多棘海盘车侵害(周书珩等, 2008)。2021~2022年连续两年在山东青岛胶州湾海域出现多棘海盘车聚集的现象, 相比2021年2月的多棘海盘车聚集事件, 2022年的多棘海盘车暴发呈现新的特点, 一方面2022年分别在2月和7月共暴发了两次大规模聚集, 另一方面2022年7月的聚集灾害系多棘海盘车与软体动物经氏壳蛞蝓(Philine kinglipini)同时暴发。对于多棘海盘车大规模暴发, 目前尚无有效的检测方法和有效的治理措施, 使用最多的治理办法仍是在暴发之后进行捕捞清理。

尽管多棘海盘车分布广泛、且暴发性增殖对生态环境和水产养殖能够产生严重的负面影响, 但针对多棘海盘车遗传多样性的分子生物学研究相对较少。食品营养和药物方面的研究发现多棘海盘车含有多糖、脂类、甾醇、皂苷等多种活性成分, 对食品药品的功能性开发具有重要的研究价值(李敏晶等, 2017)。同时, 还有大量的研究集中于多棘海盘车的生物学特性和对菲律宾蛤仔、魁蚶等贝类的捕食特性(张秀梅等, 2014; 周珊珊等, 2014; 李淑芸等, 2014; 张天文等, 2015)。此外, 研究报道称海星等棘皮动物会吸收海水中的碳, 以无机盐形式形成外骨骼, 死亡后体内大部分含碳物质留在海底, 从而减少海洋中的碳进入大气, 在海洋碳循环中起着重要作用(Lebrato et al, 2010)。

迄今对于多棘海盘车的遗传多样性分析主要集中在日本海样本。采自日本岩手县的多棘海盘车样本的线粒体基因组最早得到解析(Matsubara et al, 2005)。最近, 采集自日本东北部青森县、宫城县女川和岩手县的16个多棘海盘车样本和南部冈山县濑户内市牛窓的10个多棘海盘车样本的线粒体基因组得到比较分析, 结果显示东北部的样本和南部的样本分别聚类为两个支系(Inoue et al, 2020)。针对我国多棘海盘车遗传多样性的研究较少，对于我国胶州湾海域, 已有的研究中仅有三个多棘海盘车幼虫的cox1序列报道(王敏晓等, 2011)。

本研究重点分析了采集于2022年7月在胶州湾多棘海盘车暴发时的12个样本, 以及3个采集自连云港海域的样本, 系统构建了其线粒体基因组、核糖体基因簇以及多种通用分子标记, 并开展了比较分析, 尝试开发了针对多棘海盘车的高分辨率分子标记。本研究不仅首次测序组装了我国胶州湾海域致灾物种多棘海盘车的多种分子序列, 开发的高分辨率分子标记也将为实时监控胶州湾以及其他海域多棘海盘车的多样性特征和时空变化提供技术支持。

1 材料与方法 1.1 样品采集与处理

2022年7月6日搭乘青岛胶州湾渔民的渔船采用680目底拖网在胶州湾海域(图 1a)底拖获得海星样本, 挑选12个健康个体(表 1)暂养于自然海水, 1小时内送达实验室。样品到达实验室后, 在灭菌海水中使用毛刷迅速清理, 使用灭菌手术剪剪取每个个体的一条腕足, 收集外壳于−80 ℃冰箱迅速冷冻。2022年8月5日搭乘科学3号科考船在黄海连云港海域底拖网采集了3个海星样本(表 1), 现场海水冲洗, 使用灭菌手术剪剪取每个个体的一条腕足, 收集外壳于−80 ℃冰箱迅速冷冻。此外, 为了进行系统的比较研究, 本文下载了已发表的25个日本多棘海盘车样本二代基因组测序的数据(Inoue et al, 2020) (表 1)。

1.2 DNA文库制备与测序

样本基因组DNA采用传统的CTAB (Cetyltrimethylammonium Bromide)方法进行提取, 所提取的DNA经Covaris S220超声破碎仪(美国)进行片段化, 获得读长为350 bp的DNA序列; 经末端修复、加A尾、加测序接头、纯化、PCR及产物富集。文库质量分别采用NGS3K/Caliper和real-time PCR (Qubit 3.0 fluorometer, Invitrogen, 美国)进行文库大小分布评估和定量分析。质量合格的文库采用深圳华大DNBSEQ-T7进行测序。

1.3 线粒体基因组组装与注释

采用GetOrganelle软件(Jin et al, 2020)从头组装获得线粒体基因组序列, 其中使用的组装软件为SPAdes (3.10.1) (Bankevich et al, 2012), 组装过程中使用的seed参考序列为NCBI登录序列A. amurensis AB183559.1 (Matsubara et al, 2005)。组装后的线粒体基因组序列通过BWA v0.7.17软件的MEM运算法对所构建的线粒体基因组序列进行质检(Li, 2011), 并通过IGV v2.8.12进行可视化比对分析(Robinson et al, 2011)。本研究采用在线软件MITOS (Bernt et al, 2013)、MFannot (https://megasun.bch.umontreal.ca/RNAweasel/)和ORF finder (Rombel et al, 2002)对线粒体基因组进行注释。线粒体圈图采用在线软件OGDraw (Greiner et al, 2019)进行绘制。

1.4 核糖体基因簇组装与注释

核糖体基因簇的组装、质检方法同1.3, 组装过程中使用的seed参考序列为NCBI登录序列蓝指海星(Linckia laevigata) LC505033.1 (Hiruta et al, 2020)。核糖体基因簇的注释采用Mega 7将组装序列与LC505033.1进行序列比对。

1.5 系统发育分析

针对分析的序列采用在线mafft软件(Madeira et al, 2022)进行比对, 然后在Mega 7 (Valach et al, 2014)中进行手动裁剪, 使用最大似然法进行系统发育树构建, bootstrap值=1 000。

1.6 基因组序列比较分析

采用Dotter软件(Sonnhammer et al, 1995)对序列进行两两比对分析, 分析的对象是以样本编号为CNS01590的线粒体基因组为参照, 分别与日本个体On1、U1和胶州湾个体CNS01589进行比较分析, 所有参数为默认参数。

1.7 分化时间估算

以已发表线粒体基因组的所有海星纲物种以及棘皮动物其他纲物种(外类群)为对象, 采用Phylosuite v1.2.2 (Zhang et al, 2020)对其线粒体基因组的13个蛋白编码基因的核酸序列进行提取、mafft比对、trimAL裁剪、Concatenate串联序列, 然后使用IQ-tree进行系统发育树构建并获得树形。分化时间估算采用PAML v4.8a的MCMCTree软件包(Yang, 2007)。首先, “Branch lengths”、“gradient (g)”和“Hessian (H)”的估算采用最大似然估算(maximum likelihood estimates, MLE)和GTR + G substitution model (model=7), 同时采用independent rates clock model (clock=1)。分化时间的估算设置三个估算时间点(calibration points), 包括棘皮动物门(Echinodermata)分化时间[501~542 Million years ago (Ma)], 游移亚门(Eleutherozoa)分化时间(482~521 Ma)和海星亚门(Asterozoa)分化时间(> 479 Ma)。系统发育树的每个节点通过FigTree v1.4.3采用95%的HPD(highest posterior density interval)进行可视化展示。

1.8 基因组位点多样性分析

以编号为CNS01590的多棘海盘车个体的线粒体基因组序列为参考序列, 采用BWA v0.7.17分别将25个日本多棘海盘车样本的DNBSEQ-T7测序结果clean reads与其进行比对。然后使用SAMtools对比对结果进行分析获得每个样本的单核苷酸多样性位点(single-nucleotide variants, SNVs)信息(homozygous support > 85%) (Koboldt et al, 2012)。所有样本的SNV位点通过内部Python脚本进行分析, 以CNS01590个体的线粒体基因组序列为参考, 以400 bp为滑移窗口尺寸进行连续计算。所有株系的每400 bp窗口的SNV密度和灵敏度通过CIRCOS进行可视化展示(Krzywinski et al, 2009)。

2 结果 2.1 物种鉴定及生物地理分布 2.1.1 形态描述

于2022年7月在胶州湾多棘海盘车暴发时采集的12个多棘海盘车样本个体大小各异, 形态特征基本相同(图 1b)。样本颜色与日本海域采集到的多棘海盘车样本基本相同, 口面呈橙色, 反口面呈紫红色, 棘刺呈白色。本次胶州湾采集的多棘海盘车样本与以前在我国渤海、黄海海域采集到的样本的特征基本一致, 盘扁平, 具有5个腕, 基部宽, 腕扁平。背板结合紧密, 背板上有背棘, 具有小刺和不明显沟槽; 背棘基部具有交叉叉棘。腹侧板上有少量皮鳃和直形叉棘; 侧步带棘交替排列, 末端尖细无小刺。根据形态分析, 判定12个胶州湾海星样本均为多棘海盘车。

2.1.2 系统发育分析特征

通过对测序结果进行线粒体组装, 获得所有12个海星样本的cox1全长序列, 与已发表的序列A. amurensis MZ435266.1、MZ435267.1进行在线Blastn同源比对, 相似性达到99.10%和99.36% (Coverage 100%)。根据Inoue等(2020)的研究, NCBI SRA已公开了25个日本海域多棘海盘车样本的Illumina测序数据, 本研究采用同样的线粒体组装和注释方法获得所有样本的线粒体基因组和cox1序列。同时, 下载了NCBI现已发表的全球多棘海盘车cox1序列, 基于此51个样本的cox1序列构建系统发育树, 结果如图 1c所示, 51个样本在发育树中聚类为两个主要的分支, 包括Clade Ⅰ和Clade Ⅱ。Clade Ⅰ的样本全部采集自日本海域, 包含日本东北部的青森县、宫城县女川和岩手县的16个样本。Clade Ⅱ包含了日本海乌苏里湾、日本南部冈山县濑户内海、我国青岛胶州湾(2006年和2022年)、我国连云港、澳大利亚塔斯马尼亚岛的样本(图 1a)。Clade Ⅱ内部聚类形成两个分支, 包括Clade ⅡA和Clade ⅡB, 这两个分支的样本显示没有明显的地域偏好性(图 1c)。

基于形态特征和系统发育特征, 可以确定2022年夏季胶州湾海星暴发样本为多棘海盘车。比较分析发现多棘海盘车具有较高的种内遗传多样性。

2.2 多棘海盘车的遗传进化分析

为了解析中国胶州湾多棘海盘车与全球其他海域多级盘车之间的遗传进化关系, 本研究对其线粒体基因组和核糖体基因簇进行了组装、注释和比较分析。

2.2.1 基于线粒体基因组的多棘海盘车种内多样性分析

本研究测序组装了采集自胶州湾的12个和连云港的3个多棘海盘车样本的线粒体基因组, 共15个样本的线粒体基因组均完成环化。比较分析发现这些线粒体结构基本一致, 无大片段结构差异。以CNS01590样本的线粒体图谱为例(图 2a, 表 2), 均包含13个蛋白编码基因、2个rRNA和22个tRNA, 长度为16 419~16 421 bp。将CNS01590的线粒体基因组分别与图 1a中三个分支中的代表样本On1、CNS01589和U1的线粒体基因组进行序列比对分析(图 2b~2d), 发现两两比较的样本之间序列相似度较高, 没有结构差异和重复序列等。

为进一步探究多棘海盘车线粒体基因组的种内变异, 本研究以2022年于胶州湾采集的样本CNS01590的线粒体基因组序列为参考, 计算了包括其他11个2022年于胶州湾采集的样本、3个2022年于连云港采集的样本和25个于日本海域采集的样本的39个多棘海盘车样本(表 1)的单核苷酸变异位点的密度(图 3a), 结果显示主要存在两个高变异位点区, 第8 000~8 600 bp和第10 900~11 300 bp, 分别位于nad5~nad6和rns与rnl之间的tRNA基因和非编码区。基于此, 本研究构建了一个长度为3 700 bp的分子标记。该分子标记可以将所有不同的个体充分区分开来, 具有很高的分辨率(图 3b~3d)。基于本研究构建的高分辨率分子标记的系统发育树(图 3b)显示, 除了U1与U2之间无法区分, 其他样本之间得以区分, 该分辨率与线粒体基因组全长序列的分辨率(图 3c)相当。对于常用的cox1分子标记, 尽管全长序列也不能区分多组样本(图 3d), 譬如A3、A5、A9与A10之间, 而本研究建立的高分辨率分子标记分辨能力远高于cox1全长序列的分辨率。因此, 该高分辨率分子标记可用于多棘海盘车的种内遗传多样性分析。

2.2.2 多棘海盘车与海星纲物种线粒体基因组的种间多样性

为了探究多棘海盘车与其他海星纲物种间的进化特征, 本研究比较分析了海星纲物种间的线粒体结构组成和共线性(图 4), 结果发现海星纲物种的基因组成极为保守, 均包含13个蛋白编码基因、2个rRNA和22个tRNA, 并且基因的组成顺序和在基因组中的方向基本一致, 仅Anthernca aspera MT476596.1的trnQ基因位于基因组的反义链(按照cox1基因位于正义链), 与其他物种相反。

基于线粒体基因组的13个蛋白编码基因, 进行了包括多个多棘海盘车样本在内的海星纲物种分化时间估算(图 5)。海星纲物种形成大约在322.1 Ma, 海盘车科(Asteriidae)物种约在161.3 Ma分化形成, 海盘车属(Asterias)约在36.6 Ma分化, 多棘海盘车物种约在8.0 Ma分化形成。多棘海盘车种内的变异可能发生在1.7~4.2 Ma。

2.2.3 基于核糖体基因簇的多棘海盘车种内多样性分析

为进一步揭示海星纲物种种间与多棘海盘车种内的多样性特征, 本研究构建了40个多棘海盘车样本和NCBI SRA数据库中已提交测序数据的海星纲其他21个物种的核基因组基因核糖体基因簇。经统计, 海星纲物种的核糖体基因簇组成如表 3, 海星纲物种的18S rDNA、ITS和28S rDNA的长度范围分别为1 815~1 832 bp、966~1 658 bp和4 039~ 5 062 bp。ITS长度变化主要在于ITS1和ITS2, 5.8S rDNA的序列长度和序列组成较保守(165~168 bp)。对于多棘海盘车而言, 该物种种内的组成相当保守, 仅CNS01591、A4和On1三个样本与其他37个样本的ITS2存在1~2 bp的长度差异。

根据基于18S rDNA、ITS、28S rDNA以及18S-ITS-28S rDNA的系统发育树(图 6), 18S rDNA、ITS、28S rDNA以及18S-ITS-28S rDNA均能够对海星纲物种种间的系统发育关系有效区分, 对于多棘海盘车种内多样性的分辨极低, 换言之, 多棘海盘车核糖体基因簇进化保守, 变异率极低。

3 讨论 3.1 基于不同分子标记的多棘海盘车种内多样性比较分析

本研究通过构建线粒体基因组(包括单基因通用分子标记cox1)和核糖体基因簇, 分别获得了中国胶州湾的多棘海盘车样本的细胞器和核相关基因的系统发育特征。根据18S rDNA、ITS、28S rDNA以及18S-ITS-28S rDNA的系统发育树, 多棘海盘车的核糖体基因簇序列种内极为保守, 仅仅CNS01591、A4和On1三个样本与其他37个样本的ITS2存在1~2 bp的长度差异, 18S rDNA、ITS1、5.8S rDNA、28S rDNA序列组成完全一致。Petrov等(2016)在比较日本海彼得大帝湾三个多棘海盘车样本的ITS序列发现种内的变异位点仅在ITS2序列, 我们的结果与其一致。相比于核糖体基因簇, 多棘海盘车的cox1序列具有较高的种内多样性。对于cox1的全长序列, Clade Ⅰ与Clade Ⅱ之间的序列相似度在96.9%~98.7%内; 即使在每个分支内部仍显示了不同个体间的明显差异, Clade Ⅰ和Clade Ⅱ内部序列相似度分别在97.8%~ 100%和98.4%~100%范围内。由于相比cox1全长序列, 线粒体基因组随着序列长度的增加, 具有更多的变异位点, 因此种内的分辨率更高, 基本上实现对本研究中的多棘海盘车样本的种内分辨(除了U1与U2之间的区分)。因此, 对于多棘海盘车而言, 线粒体基因组相比核糖体基因簇具有明显的进化变异, 线粒体基因组适合于多棘海盘车的进化研究, 核糖体基因簇更适用于物种的种间鉴定。

此外, 本研究比较分析了40个多棘海盘车个体的线粒体基因组多样性, 并在此基础上构建了一个约3 700 bp的高分辨率分子标记, 该分子标记的分辨率水平能够达到线粒体基因组全长序列的分辨率, 可以后续应用于基于三代扩增子技术的宏条形码分析, 从而能够实现对环境样本eDNA中多棘海盘车遗传多样性的快速检测。

3.2 多棘海盘车的遗传进化特征

基于线粒体基因组的比较分析, 我们发现多棘海盘车种内线粒体基因组结构完全一致, 甚至海星纲的物种间几乎没有结构的变异, 除了Anthernca aspera MT476596.1 (Quek et al, 2021)的trnQ基因所在正负链的位置与其他物种相反。

然而, 多棘海盘车线粒体基因组种内的变异十分显著。在本研究分析的40个样本中, 仅有两个样本的线粒体基因组序列组成完全一致。我国胶州湾和连云港的多棘海盘车样本与日本牛窓海域的样本聚为一枝, 说明我国胶州湾、连云港以及日本牛窓海域间多棘海盘车存在基因交流, 这可能与多棘海盘车具有较强的适应能力有关。根据海星纲物种的分化时间估算, 多棘海盘车物种的形成大约在8.0 Ma, 胶州湾样本的种内变异大约发生在4.2 Ma。根据线粒体基因组的系统发育关系, 澳大利亚样本(MZ435267)与胶州湾样本亲缘关系更近。20世纪80年代初多棘海盘车可能通过船舶的压舱水从北太平洋跨越热带而引入澳大利亚塔斯马尼亚东南部, 之后向北扩散, 1995年发现维多利亚州菲利普湾港已建立种群, 现已在菲利普港外又发现了四个种群(Inverloch、San Remo、Tidal River和Gippsland Lakes; 均在维多利亚州境内), 表明该物种目前正在扩大范围(Richardson et al, 2015)。据报道, 多棘海盘车已经成为入侵南北极的潜在物种(Byrne et al, 2016)。适应是生命进化多样化的核心特征和主要解释, 因其复杂性和超长的时间尺度而难以开展研究(Battlay et al, 2023)。生物入侵是快速进化的典型模型(Reznick et al, 2001), 遗传多样性常用于指示进化适应的潜力(Lee, 2002)。促使入侵物种适应新环境的遗传变化来自新突变还是遗传变异是进化生物学中的一个关键问题, 这个问题的解答与理解入侵性的快速进化同样重要(Hermisson et al, 2005)。Battlay等(2023)追踪一种侵略性入侵的杂草Ambrosia artemisiifolia时发现, 在富含平行适应特征的基因组区域中识别出具有多个特征的单倍型块, 这些单倍体块赋予其对当地气候的平行适应的能力, 与快速适应性状相关, 并在空间和时间上表现出显著的频率变化。多棘海盘车显著的变异能力为研究进化适应提供了潜在的实验材料。

4 结论

本研究比较研究了全球多个海域多棘海盘车样本的单基因通用分子标记cox1、线粒体基因组和核糖体基因簇的系统发育特征, 不仅揭示了中国海域多棘海盘车的遗传多样性, 同时也认识了全球多棘海盘车的遗传特点。多棘海盘车种内具有极高的遗传多样性, 变异水平可能远超现有实验样本所体现的水平, 这可能与其强大的适应能力有关。后续仍需要密切跟踪多棘海盘车的遗传多样性特征, 下一步工作的重点是建立可用于多棘海盘车环境样本eDNA宏条形码分析的分析方法和数据平台。

致谢 感谢中国科学院海洋研究所海洋大数据中心的支持; 感谢在2022胶州湾样本采样过程中青岛市胶州市东营码头渔民胡大山在采样工作中给予的支持和便利。