中国海洋湖沼学会主办。
文章信息
- 王静, 段泽林, 何子岩, 陈楠生. 2023.
- WANG Jing, DUAN Ze-Lin, HE Zi-Yan, CHEN Nan-Sheng. 2023.
- 2022年夏季胶州湾致灾物种多棘海盘车的分子解析
- MOLECULAR ANALYSIS OF ASTERIAS AMURENSIS FROM A STARFISH OUTBREAK IN JIAOZHOU BAY IN SUMMER 2022
- 海洋与湖沼, 54(6): 1656-1671
- Oceanologia et Limnologia Sinica, 54(6): 1656-1671.
- http://dx.doi.org/10.11693/hyhz20230300045
文章历史
-
收稿日期:2023-03-02
收修改稿日期:2023-07-14
2. 崂山实验室 海洋生态与环境科学功能实验室 山东青岛 266237;
3. 中国科学院海洋大科学研究中心 山东青岛 266071;
4. 中国科学院大学 北京 100049
2. Laboratory for Marine Ecology and Environmental Science, Laoshan Laboratory, Qingdao 266237, China;
3. Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China;
4. University of Chinese Academy of Sciences, Beijing 100049, China
多棘海盘车(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个/m2 (周书珩等, 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)。
编号 | 海域 | 经纬度 | 采样时间/(年.月.日) | 参考文献 |
CNS01589 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01590 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01591 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01908 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01909 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01910 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01911 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01912 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01913 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01914 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01915 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01916 | 中国山东, 胶州湾 | 36°6.30′ N, 120°12.017′ E | 2022.7.7 | 本文 |
CNS01988 | 中国江苏, 连云港 | 35°0.26′ N, 119°46.04′ E | 2022.8.5 | 本文 |
CNS01990 | 中国江苏, 连云港 | 35°0.26′ N, 119°46.04′ E | 2022.8.5 | 本文 |
CNS01991 | 中国江苏, 连云港 | 35°0.26′ N, 119°46.04′ E | 2022.8.5 | 本文 |
A1 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A2 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A3 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A4 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A5 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A6 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A7 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A8 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A9 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
A10 | 日本, 青森县 | 40°89′ N, 140°86′ E | 2017 | Inoue et al, 2020 |
On1 | 日本, 宫城县女川 | 38°44′ N, 141°46′ E | 2018 | Inoue et al, 2020 |
On2 | 日本, 宫城县女川 | 38°44′ N, 141°46′ E | 2018 | Inoue et al, 2020 |
On3 | 日本, 宫城县女川 | 38°44′ N, 141°46′ E | 2018 | Inoue et al, 2020 |
On4 | 日本, 宫城县女川 | 38°44′ N, 141°46′ E | 2018 | Inoue et al, 2020 |
On5 | 日本, 宫城县女川 | 38°44′ N, 141°46′ E | 2018 | Inoue et al, 2020 |
U1 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U2 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U3 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U4 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U5 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U6 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U7 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U8 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U9 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
U10 | 日本, 冈山县牛窓 | 34°60′ N, 134°14′ E | 2017 | Inoue et al, 2020 |
样本基因组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), 发现两两比较的样本之间序列相似度较高, 没有结构差异和重复序列等。
基因名称 | 类型 | +/−链 | 起始碱基 | 终止碱基 | 长度/bp | 密码子 | ||
起始密码子 | 终止密码子 | 反义密码子 | ||||||
cox1 | CDS | + | 1 | 1 554 | 1 554 | ATG | TAA | |
trnR | tRNA | + | 1 553 | 1 622 | 70 | TCG | ||
nad4L | CDS | + | 1 623 | 1 919 | 297 | ATT | TAA | |
cox2 | CDS | + | 1 920 | 2 609 | 690 | ATG | TAA | |
trnK | tRNA | + | 2 611 | 2 685 | 75 | CTT | ||
atp8 | CDS | + | 2 689 | 2 856 | 168 | ATG | TAA | |
atp6 | CDS | + | 2 847 | 3 533 | 687 | ATG | TAA | |
cox3 | CDS | + | 3 537 | 4 319 | 783 | ATG | TAA | |
trnS | tRNA | − | 4 318 | 4 388 | 71 | TGA | ||
nad3 | CDS | + | 4 414 | 4 764 | 351 | ATG | TAA | |
nad4 | CDS | + | 4 782 | 6 164 | 1 383 | ATG | TAA | |
trnH | tRNA | + | 6 315 | 6 383 | 69 | GTG | ||
trnS | tRNA | + | 6 385 | 6 452 | 68 | GCT | ||
nad5 | CDS | + | 6 453 | 8 369 | 1 917 | ATG | TAA | |
nad6 | CDS | − | 8 402 | 8 890 | 489 | ATG | TAG | |
cob | CDS | + | 8 906 | 10 099 | 1 194 | ATG | TAA | |
trnF | tRNA | + | 10 044 | 10 115 | 72 | GAA | ||
rns | rRNA | + | 10 115 | 11 009 | 895 | |||
trnE | tRNA | + | 11 009 | 11 075 | 67 | TTC | ||
trnT | tRNA | + | 11 078 | 11 147 | 70 | TGT | ||
rnl | rRNA | − | 11 640 | 13 246 | 1 607 | |||
nad2 | CDS | − | 13 270 | 14 331 | 1 062 | GTG | TAA | |
trnI | tRNA | − | 14 332 | 14 403 | 72 | GAT | ||
nad1 | CDS | − | 14 404 | 15 381 | 978 | GTG | TAA | |
trnL | tRNA | − | 15 382 | 15 454 | 73 | TAA | ||
trnG | tRNA | − | 15 487 | 15 555 | 69 | TCC | ||
trnY | tRNA | − | 15 560 | 15 630 | 71 | GTA | ||
trnD | tRNA | + | 15 634 | 15 704 | 71 | GTC | ||
trnM | tRNA | − | 15 707 | 15 779 | 73 | CAT | ||
trnV | tRNA | + | 15 811 | 15 882 | 72 | TAC | ||
trnC | tRNA | − | 15 882 | 15 950 | 69 | GCA | ||
trnW | tRNA | − | 15 952 | 16 019 | 68 | TCA | ||
trnA | tRNA | + | 16 030 | 16 101 | 72 | TGC | ||
trnL | tRNA | − | 16 103 | 16 174 | 72 | TAG | ||
trnN | tRNA | − | 16 175 | 16 246 | 72 | GTT | ||
trnQ | tRNA | + | 16 250 | 16 321 | 72 | TTG | ||
trnP | tRNA | − | 16 323 | 16 391 | 69 | TGG |
为进一步探究多棘海盘车线粒体基因组的种内变异, 本研究以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 | ITS1 | 5.8S rDNA | ITS2 | ITS | 28S rDNA | 18S-ITS-28S | Accession No. in NCBI |
Acanthaster planci DRR251129 | 1 816 | 398 | 166 | 655 | 1 219 | 4 089 | 7 124 | |
Anthenea aspera SRR12622300 | 1 818 | 386 | 166 | 591 | 1 143 | 4 230 | 7 191 | |
Archaster typicus SRR12622308 | 1 816 | 367 | 167 | 432 | 966 | 4 039 | 6 821 | |
Asterias amurensis A1 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A2 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A3 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A4 | 1 815 | 635 | 167 | 658 | 1 460 | 4 059 | 7 334 | |
Asterias amurensis A5 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A6 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A7 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A8 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A9 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis A10 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis CNS01589 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630824 |
Asterias amurensis CNS01590 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630825 |
Asterias amurensis CNS01591 | 1 815 | 635 | 167 | 657 | 1 459 | 4 059 | 7 333 | OR630826 |
Asterias amurensis CNS01908 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630827 |
Asterias amurensis CNS01909 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630828 |
Asterias amurensis CNS01910 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630829 |
Asterias amurensis CNS01911 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630830 |
Asterias amurensis CNS01912 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630831 |
Asterias amurensis CNS01913 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630832 |
Asterias amurensis CNS01914 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630833 |
Asterias amurensis CNS01915 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630834 |
Asterias amurensis CNS01916 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630835 |
Asterias amurensis CNS01988 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630836 |
Asterias amurensis CNS01990 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7332 | OR630837 |
Asterias amurensis CNS01991 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | OR630838 |
Asterias amurensis On1 | 1 815 | 635 | 167 | 657 | 1 459 | 4 059 | 7 333 | |
Asterias amurensis On2 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis On3 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis On4 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis On5 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U1 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U2 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U3 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U4 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U5 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U6 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U7 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U8 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U9 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias amurensis U10 | 1 815 | 635 | 167 | 656 | 1 458 | 4 059 | 7 332 | |
Asterias rubens ERR3312498 | 1 815 | 605 | 167 | 653 | 1 425 | 4 063 | 7 303 | |
Coscinasterias acutispina SRR19025778 | 1 815 | 575 | 167 | 611 | 1 353 | 4 072 | 7 240 | |
Cryptasterina hystera SRR12430566 | 1 817 | 391 | 166 | 575 | 1 132 | 4 119 | 7 068 | |
Culcita novaeguineae SRR12622305 | 1 817 | 394 | 166 | 677 | 1 237 | 4 094 | 7 148 | |
Euretaster insignis SRR12622299 | 1 823 | 842 | 168 | 577 | 1 587 | 5 062 | 8 472 | |
Goniodiscaster scaber SRR12622302 | 1 817 | 425 | 166 | 541 | 1 132 | 4 268 | 7 217 | |
Iconaster longimanus SRR12622306 | 1 827 | 426 | 167 | 682 | 1 275 | 4 177 | 7 279 | |
Linckia laevigata LC505033.1 | 1 816 | 461 | 167 | 603 | 1 231 | 4 263 | 7 310 | |
Luidia maculata SRR12622298 | 1 832 | 508 | 166 | 657 | 1 331 | 4 602 | 7 765 | |
Marthasterias glacialis ERR6054755 | 1 815 | 685 | 167 | 676 | 1 528 | 4 071 | 7 414 | |
Nepanthia sp SRR12622301 | 1 818 | 408 | 166 | 563 | 1 137 | 4 079 | 7 034 | |
Ophidiaster granifer SRR12622297 | 1 828 | 686 | 166 | 806 | 1 658 | 4 560 | 8 046 | |
Patiria miniata SRR343071 | 1 816 | 467 | 165 | 588 | 1 220 | 4 041 | 7 077 | |
Patiria pectinifera DRR157371 | 1 816 | 448 | 165 | 575 | 1 188 | 4 047 | 7 051 | |
Patiriella regularis ERR995408 | 1 816 | 476 | 166 | 630 | 1 272 | 4 125 | 7 213 | |
Pentaceraster mammillatus SRR12622307 | 1 817 | 397 | 166 | 689 | 1 252 | 4 125 | 7 194 | |
Plazaster borealis SRR16914471 | 1 815 | 561 | 167 | 649 | 1 377 | 4 060 | 7 252 | |
Protoreaster nodosus SRR12622304 | 1 817 | 404 | 166 | 670 | 1 240 | 4 131 | 7 188 |
根据基于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胶州湾样本采样过程中青岛市胶州市东营码头渔民胡大山在采样工作中给予的支持和便利。
王敏晓, 程方平, 李超伦, 等, 2011. 基于线粒体cox1片段序列的胶州湾浮游动物DNA条形码分析[J]. 海洋与湖沼, 42(5): 702-710. |
孙光, 1983. 海星类的生态及其在贝类增养殖中的危害和防除方法的研究[J]. 国外水产, 3(4): 16-20. |
李敏晶, 孙蔷薇, 王旭莹, 等, 2017. 多棘海盘车多肽的抑菌活性研究[J]. 应用化工, 46(10): 2073-2076. DOI:10.3969/j.issn.1671-3206.2017.10.051 |
李淑芸, 张秀梅, 聂猛, 等, 2014. 不同温度下多棘海盘车对菲律宾蛤仔的摄食选择性研究[J]. 水产学报, 38(7): 981-991. |
张天文, 刘广斌, 刘恩孚, 等, 2015. 多棘海盘车对魁蚶摄食量、选择性及昼夜摄食差异的初步研究[J]. 中国海洋大学学报自然科学版, 18(12): 24-29. |
张秀梅, 李淑芸, 刘佳, 等, 2014. 青岛近海多棘海盘车繁殖生物学的初步研究[J]. 中国海洋大学学报自然科学版, 44(11): 16-24. |
周书珩, 王印庚, 2008. 近海水域海星泛滥引起的反思[J]. 水产科学, 27(10): 555-556. DOI:10.3969/j.issn.1003-1111.2008.10.015 |
周珊珊, 张秀梅, 蔡星媛, 等, 2014. 温度对魁蚶稚贝潜沙能力及对多棘海盘车捕食魁蚶稚贝能力的影响[J]. 水产学报, 38(9): 1439-1446. |
徐思嘉, 肖宁, 曾晓起, 2018. 中国海域海盘车科(棘皮动物门, 海星纲)种类记述[J]. 海洋科学, 42(10): 53-63. |
BANKEVICH A, NURK S, ANTIPOV D, et al, 2012. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing[J]. Journal of Computational Biology, 19(5): 455-477. DOI:10.1089/cmb.2012.0021 |
BATTLAY P, WILSON J, BIEKER V C, et al, 2023. Large haploblocks underlie rapid adaptation in the invasive weed Ambrosia artemisiifolia[J]. Nature Communications, 14(1): 1717. DOI:10.1038/s41467-023-37303-4 |
BERNT M, DONATH A, JÜHLING F, et al, 2013. MITOS: Improved de novo metazoan mitochondrial genome annotation[J]. Molecular Phylogenetics and Evolution, 69(2): 313-319. DOI:10.1016/j.ympev.2012.08.023 |
BUTTERMORE R E, TURNER E, MORRICE M G, 1994. The introduced northern Pacific seastar Asterias amurensis in Tasmania[J]. Memoirs of the Queensland Museum, 36: 21-25. |
BYRNE M, GALL M, WOLFE K, et al, 2016. From pole to pole: the potential for the Arctic seastar Asterias amurensis to invade a warming Southern Ocean[J]. Global Change Biology, 22(12): 3874-3887. DOI:10.1111/gcb.13304 |
GREINER S, LEHWARK P, BOCK R, 2019. Organellar Genome DRAW (OGDRAW) version 1.3. 1: expanded toolkit for the graphical visualization of organellar genomes[J]. Nucleic Acids Research, 47(W1): W59-W64. DOI:10.1093/nar/gkz238 |
HERMISSON J, PENNINGS P S, 2005. Soft sweeps: molecular population genetics of adaptation from standing genetic variation[J]. Genetics, 169(4): 2335-2352. DOI:10.1534/genetics.104.036947 |
HIRUTA S F, ARAI M, CHAVANICH S, et al, 2020. Complete mitochondrial genome of a sea star, Linckia laevigata (Echinodermata, Asteroidea, Valvatida, Ophidiasteridae)[J]. Mitochondrial DNA Part B, 5(2): 1342-1343. DOI:10.1080/23802359.2020.1735271 |
INOUE J, HISATA K, YASUDA N, et al, 2020. An investigation into the genetic history of Japanese populations of three starfish, Acanthaster planci, Linckia laevigata, and Asterias amurensis, based on complete mitochondrial DNA sequences[J]. G3 Genes Genomes Genetics, 10(7): 2519-2528. |
JIN J J, YU W B, YANG J B, et al, 2020. GetOrganelle: a fast and versatile toolkit for accurate de novo assembly of organelle genomes[J]. Genome Biology, 21(1): 241. |
KIM Y S, 1968. Histological observations of the annual change in the gonad of the starfish, Asterias amurensis Lüken[J]. Bulletin of the Faculty of Fisheries Hokkaido University, 19(2): 97-108. |
KOBOLDT D C, ZHANG Q Y, LARSON D E, et al, 2012. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing[J]. Genome Research, 22(3): 568-576. |
KRZYWINSKI M, SCHEIN J, BIROL İ, et al, 2009. Circos: an information aesthetic for comparative genomics[J]. Genome Research, 19(9): 1639-1645. |
LEBRATO M, IGLESIAS-RODRíGUEZ D, FEELY R A, et al, 2010. Global contribution of echinoderms to the marine carbon cycle: CaCO3 budget and benthic compartments[J]. Ecological Monographs, 80(3): 441-467. |
LEE C E, 2002. Evolutionary genetics of invasive species[J]. Trends in Ecology & Evolution, 17(8): 386-391. |
LI H, 2011. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data[J]. Bioinformatics, 27(21): 2987-2993. |
MADEIRA F, PEARCE M, TIVEY A R N, et al, 2022. Search and sequence analysis tools services from EMBL-EBI in 2022[J]. Nucleic Acids Research, 50(W1): W276-W279. |
MATSUBARA M, KOMATSU M, ARAKI T, et al, 2005. The phylogenetic status of Paxillosida (Asteroidea) based on complete mitochondrial DNA sequences[J]. Molecular Phylogenetics and Evolution, 36(3): 598-605. |
PETROV N B, VLADYCHENSKAYA I P, DROZDOV A L, et al, 2016. Molecular genetic markers of intra- and interspecific divergence within starfish and sea urchins (Echinodermata)[J]. Biochemistry, 81(9): 972-980. |
QUEK Z B R, CHANG J J M, IP Y C A, et al, 2021. Mitogenomes reveal alternative initiation codons and lineage-specific gene order conservation in Echinoderms[J]. Molecular Biology and Evolution, 38(3): 981-985. |
REZNICK D N, GHALAMBOR C K, 2001. The population ecology of contemporary adaptations: what empirical studies reveal about the conditions that promote adaptive evolution[J]. Genetica, 112(1): 183-198. |
RICHARDSON M F, SHERMAN C D H, 2015. De novo assembly and characterization of the invasive Northern Pacific seastar transcriptome[J]. PLoS One, 10(11): e0142003. |
ROBINSON J T, THORVALDSDóTTIR H, WINCKLER W, et al, 2011. Integrative genomics viewer[J]. Nature Biotechnology, 29(1): 24-26. |
ROMBEL I T, SYKES K F, RAYNER S, et al, 2002. ORF-FINDER: a vector for high-throughput gene identification[J]. Gene, 282(1/2): 33-41. |
SONNHAMMER E L L, DURBIN R, 1995. A dot-matrix program with dynamic threshold control suited for genomic DNA and protein sequence analysis[J]. Gene, 167(1/2): GC1-GC10. |
STEVENS C, 2022. Asterias amurensis (northern Pacific seastar). CABI Compendium [DB]. doi: 10.1079/cabicompendium.92632.
|
VALACH M, BURGER G, GRAY M W, et al, 2014. Widespread occurrence of organelle genome-encoded 5S rRNAs including permuted molecules[J]. Nucleic Acids Research, 42(22): 13764-13777. |
WARD R D, ANDREW J, 1995. Population genetics of the northern Pacific seastar Asterias amurensis (Echinodermata: Asteriidae): allozyme differentiation among Japanese, Russian, and recently introduced Tasmanian populations[J]. Marine Biology, 124(1): 99-109. |
YANG Z H, 2007. PAML 4: phylogenetic analysis by maximum likelihood[J]. Molecular Biology and Evolution, 24(8): 1586-1591. |
ZHANG D, GAO F L, JAKOVLIĆ I, et al, 2020. PhyloSuite: An integrated and scalable desktop platform for streamlined molecular sequence data management and evolutionary phylogenetics studies[J]. Molecular Ecology Resources, 20(1): 348-355. |