中国北方沿海竿虾虎鱼(Luciogobius guttatus)隐存多样性与群体历史动态*

邓颖达 薛东秀 刘进贤

(1. 中国科学院海洋研究所 海洋生态与环境科学重点实验室 青岛 266071; 2. 中国科学院大学 北京 100049;3. 青岛海洋科学与技术国家实验室海洋生态与环境科学功能实验室 青岛 266071)

由于地理隔离造成的基因交流限制和环境变化带来的选择压力, 在近缘物种或种内谱系间, 分子水平的显著分化可能先于其形态学差异而发生。然而,这些从形态学无法或难以区分的隐存单元可能预示着新物种的形成。线粒体DNA是海洋鱼类系统发育研究中常用的一种分子标记。由于其突变速率较高,且在脊椎动物中单亲遗传、极少发生同源重组, 故线粒体 DNA多态性位点丰富且连续, 对历史环境变化响应迅速, 在近缘鱼类物种间具有较高的解析度。竿虾虎鱼属(Luciogobius)是虾虎鱼科(Gobiidae)、拟虾虎鱼亚科(Gobionellinae)中的一个特化属, 已知分布于中、日、韩沿海(伍汉霖等, 2008)。早在1932年, 我国学者Chen Johnson(1932)首先报道了竿虾虎鱼属在我国沿海的分布记录, 并从形态学角度对其分类地位作出初步探讨。竿虾虎鱼属物种成体的第一背鳍完全退化消失, 无鳞, 侧线系统退化, 且脊椎数目增加以致躯干延长(伍汉霖等, 2008; Yamada et al, 2009),使得本属从形态上容易与其它属相互区分。但在本属内部的谱系甚至物种之间, 由于栖息地适应性的趋同进化, 从形态上难以有效进行分类。线粒体是细胞能量代谢的主要场所, 对于从寒温带到亚热带广泛分布的竿虾虎鱼属而言, 更新世东亚沿海气候和地理环境的变化可能在其线粒体DNA中留下了深刻的印迹; 尤其在形态分类不一致的情况下, 线粒体DNA标记有助于揭示竿虾虎鱼属物种的隐存多样性,从而促进明确一致的系统分类, 并反演群体的近期历史动态。

竿虾虎鱼(Luciogobius guttatus)是竿虾虎鱼属的常见种之一。竿虾虎鱼一龄鱼即性成熟, 怀卵量370—1542粒, 冬季产卵; 二龄鱼体长 4—6cm(伍汉霖等, 2008), 是一种底内生活型(infaunal)的虾虎鱼;加之其成体常栖息于暖水潮间带砾石之下, 捕食底栖无脊椎动物, 生殖季节时雄性个体在石下筑巢护卵, 故竿虾虎鱼亦是一种生态功能相对独立且迁移能力较弱的潮间带捕食者, 可能会通过下行控制效应和生物扰动效应影响其栖息地底栖生物群落的物种多样性。日本学者Mukai等(2004)从分子层面(线粒体16S rRNA)报道了竿虾虎鱼(L. guttatus)在日本沿海的谱系组成和地理分布, 发现竿虾虎鱼种内存在三个谱系。此后, 这一结果被Yamada等(2009)用核基因外显子和线粒体DNA标记(cytb)进一步证实。但这三个谱系在本属的线粒体基因树上并不总为单系群,且其中一个分布较广并且最先形成的谱系(根据Yamada等(2009)的命名, 为A谱系)与另外两个谱系(B和C)的遗传距离(P-distance)在cytb上都超过了2%的 DNA条形码种级经验阈值(Avise et al, 1999;Hebert et al, 2003)。然而, 我国分布的竿虾虎鱼隐存多样性的研究仍未见报道。为查明中国沿海竿虾虎鱼的谱系组成, 并深入比较谱系间历史动态差异, 我们根据《中国动物志·硬骨鱼纲·鲈形目·虾虎鱼亚目》的记载, 采集了竿虾虎鱼在我国北方沿海的群体样本,选用两个线粒体 DNA片段(cytb 和控制区 mcr), 和一个核基因外显子(钙离子通道蛋白脑亚型, 亦即：兰尼碱受体 3, ryr3)为分子标记, 并结合已发表的日本群体的基因序列, 对我国北方沿海竿虾虎鱼的隐存多样性与群体历史动态进行了初步研究。研究结果将加深对我国海岸带鱼类多样性的认识, 并为我国沿海竿虾虎鱼的保护提供分子依据。

1 材料与方法

1.1 样品采集

于2013—2015年共在我国北方沿海10个位点采集到 125条竿虾虎鱼个体(图 1)。所有个体均是退潮时在潮间带手工捕捉采集, 活体带回实验室, 无水乙醇保存备用。其中, 2013年采样点为： 青岛市汇泉湾(HQ,4月, 12个个体)、石老人海水浴场(SLR, 8月, 8个个体)、崂山仰口(LS, 6月, 12个个体); 2014年10月采样点为：荣成俚岛(RC, 11个个体)、大连黑石礁(DL, 21个个体);2015年4月青岛采样点： 汇泉湾(12个个体)、石老人海水浴场(12个个体)、唐岛湾(td, 5个个体)、金沙滩(js, 7个个体)、栈桥(zq, 8个个体)、第三海水浴场(sy, 9个个体)、麦岛(md, 8个个体)。

图1 竿虾虎鱼谱系的已知分布格局Fig.1 Distribution of lineages of L. guttatus

1.2 DNA提取

剪取鳍条组织提取总 DNA, 提取过程按照苯酚-氯仿-异戊醇法(Sambrook et al, 1989)进行。所提取的总DNA溶于TE溶液, -20°C保存。

1.3 PCR扩增

使用已发表的引物及自主设计的引物(表 1)进行目标DNA片段的PCR扩增, 具体PCR反应体系为：2×PCR-mix (东盛)12.5μL, H2O 9.5μL, 正反向引物各1μL(浓度 0.1mmol/L), 模板总 DNA 1μL, 共 25μL;PCR反应程序为： 94°C预变性10min, 94°C变性1min,退火1min (温度见表1), 72°C延伸1min; 循环35次;然后在72°C终末延伸10min。

1.4 序列测定

采自2013—2014年的样本至少在一个基因标记上被测序(表 2), 各片段皆为双向测序。2015年采自汇泉湾和石老人海水浴场的样本仅测定了控制区mcr片段序列。测序由上海英潍捷基公司完成。2015年采自唐岛湾、金沙滩、栈桥、第三海水浴场以及麦岛的样本仅通过一对谱系特异性引物(fDloopThrF和LgCA, 引物LgCA根据mcr中的一处长20bp的谱系C特异性插入序列自主设计)用PCR-琼脂糖凝胶电泳方法检查其是否属于谱系C。

表1 本研究使用引物及其来源、用途、退火温度和序列Tab.1 Source, application, annealing temperature, and sequence of primers used in this study

表2 测序群体在每个标记上的测序个体数Tab.2 The number of individuals sequenced at each marker per population

1.5 序列比对

使用DNAStar软件包中的SeqMan程序查看测序峰图, 除去引物区, 舍弃双峰区段并对稀有变异位点逐个检查。在MegAlign程序中用ClustalV算法进行初步比对, 结果保存为paup文件, 后者再经DNASP软件转为fasta文件。在MEGA5(Tamumraet al, 2011)中对 fasta格式的序列进行再次比对后, 根据日本学者已发表的同源序列, 对标记cytb和ryr3截至相同长度用于后续分析; 对线粒体控制区mcr保留测序全长。考虑到不同标记的突变速率、碱基替代模型和选择压力可能存在差异, 我们没有把不同标记拼接成联合片段。本研究所提交序列的 GenBank登录号为KT598738-KT598761(ryr3), KT598808-KT598860(cytb),KT598970-KT599022(mcr); 本研究引用序列登录号为AB503858-AB504047(cytb), AB504202-AB504232(ryr3), KF486456-KF486460(ryr3)。

1.6 系统发育分析

首先基于最大似然法(ML method), 对 2013—2014年测序所获得的所有竿虾虎鱼的cytb片段和ryr3片段, 在MEGA5(Tamumraet al, 2011)中筛选最佳碱基替代模型(BIC信息标准值最低)并构建最大似然树, 自举1000次。另外, 分别基于cytb和ryr3两个基因片段, 应用Beast程序(Drummondet al, 2012)中的贝叶斯方法(同时考虑遗传距离、替代模型和碱基顺序)并加入外群物种(Yamadaet al, 2009), 构建了贝叶斯系统发育树。突变速率设为正态分布, 其均值在cytb为2.7%每碱基每百万年(Yamadaet al, 2009),ryr3的突变均速按总体P遗传距离在两个标记间的比例设定; 其标准差的设定保证 95%可能的突变速率落入 0.1%—10%每碱基每百万年的取值区间内; 各标记的突变模型与相应最大似然树所用模型相同;马尔科夫链长200万代, 采用Speciation： Yule Process假设, 每 1000代建立一棵系统发育树; 舍弃前 10%后将所有分枝平均后验概率最高者选出, 在 FigTree中显示最后所得树形图(贝叶斯系统发育树)。

1.7 群体分化与历史动态分析

针对单一谱系 A 或 C, 在群体分化方面, 应用DNASP软件(Libradoet al, 2009)比较各分子标记的多态性(核苷酸多态度、多态位点数目、单倍型数目);通过Arlequin3.5软件(Excoffieret al, 2010)计算cytb和mcr基因片段所揭示的群体间遗传分化指数FST及其显著性概率P值。另外, 为比较中日大群体(所有中国沿海单倍型和所有日本沿海单倍型)间的遗传分化, 亦应用 Arlequin3.5软件构建竿虾虎鱼谱系间的单倍型最小进化网络图。最后, 应用 Beast软件包中的Tracer程序, 选择与本文系统发育分析同样的突变速率分布和碱基替代模型, 但采用 Coalescent：Bayesian Skyline假设与linear stepwise溯祖模型, 对每个谱系进行了有效种群大小的历史动态反演, 其结果以Bayesian skyline plots(BSP)方式呈现。

此外, 为探讨竿虾虎鱼谱系 A中群体内的遗传变异是否在时间上具有稳定性, 我们对所获得的汇泉湾和石老人海水浴场的2013年群体样本和2015年群体样本的 mcr基因片段序列进行比对分析, 应用MEGA5软件构建ML进化树。

2 结果与分析

2.1 DNA分子标记多态性比较

在本研究中, 所获得的cytb的608bp片段中共有103个多态性位点, 在52个测序个体中检测到22个单倍型, 核苷酸多样性 π值为 0.069; 所获得的 mcr的917bp片段中共有88个多态性位点与27个indel位点, 在53个测序个体(不包括2015年采集的汇泉湾和石老人群体样品)中检测到40个单倍型, 核苷酸多样性π值为0.036(考虑indels： 0.047); 所获得的核基因ryr3的505bp片段中共有12个多态性位点, 在24个测序个体中区分出 10个单倍型, 核苷酸多样性 π值为0.0063。

2.2 竿虾虎鱼隐存多样性及其地理格局

如图1所示, 中国北方沿海分布有偏南、偏北两个竿虾虎鱼的谱系, 分别对应竿虾虎鱼日本沿海群体的谱系A和C(Yamada et al, 2009)。其中, 偏南谱系(谱系A)集中分布在青岛市石老人海域及其以南(除一个个体混居于大连群体中, 两个个体混居于荣成群体外); 而偏北谱系(谱系C)只分布在崂山沿海及其以北。

2.3 谱系间分子系统发育比较

结合已发表的 L. guttatus的同源基因序列(Yamada et al, 2009), 我们发现从线粒体DNA片段和核基因外显子片段上分开的竿虾虎鱼谱系能够相互对应, 但分支顺序不同(图2)。且当加入外群时, 基于ryr3构建的贝叶斯树揭示竿虾虎鱼谱系C最先分出,谱系B其次, 谱系A处于树的顶端; 而基于线粒体基因cytb构建的贝叶斯树则表明谱系B、C先聚成一枝,谱系A单成一枝并先与其他竿虾虎鱼属物种相聚(图3)。

2.4 中日沿海群体分化与有效种群扩张

图2 基于cytb基因片段(上)和ryr3基因片段(下)构建的竿虾虎鱼ML进化树Fig.2 ML trees of L. guttatus based on cytb (upper) and ryr3(nether)

对单一谱系从群体间遗传分化指数 FST值判断,在中日群体间, 竿虾虎鱼谱系A在线粒体cytb标记上的遗传分化较大(FST= 0.12, P＜0.05), 但谱系C在中日大群体间的分化更具显著性(FST= 0.31, P＜0.01)。单倍型网络图(图 4)更直观地揭示了中日大群体间的遗传分化。基于线粒体cytb我国所有测序群体除大连-荣成外, 其余两两群体间皆存在显著的遗传分化(0.17＜FST＜0.23, P＜0.05); 使用线粒体控制区标记 mcr亦检测到较大的群体间遗传分化(-0.02 ＜ FST＜ 0.96),且汇泉-石老人(FST=0.08)，大连-荣成(FST=-0.02)，大连-崂山间(FST=0.00)分化不显著(P ＞ 0.05)，其余两两群体间的遗传分化均达到显著水平(P＜0.03)。另一方面, 基于cytb的BSP结果表明, 谱系A、C都呈现出3—5万年以来的有效种群扩张, 且谱系C的变化幅度更大; 用核基因片段ryr3并未检测到近期的种群扩张, 但谱系C的总体瓶颈效应更为明显(图5)。

2.5 竿虾虎鱼群体遗传组成时间上的稳定性

基于线粒体 mcr基因片断对分别于 2013年和2015年采自汇泉湾与石老人海水浴场的四个群体进行群体遗传学分析, 结果显示同一地理群体不同采样年份的样本间 FST值均较低(HQ2013-HQ2015, 群体间 FST= -0.01, P ＞ 0.05; SLR2013-SLR2015, 群体间FST= 0.05, P ＞ 0.05)。如图6所示, 采自汇泉湾和石老人海水浴场的个体并未在所构建的 ML树上分出与年份对应的枝系, 两个群体在同一年份的样本都是分散的, 且以遗传距离衡量其分枝长度在 1.5%以内。此外, 在混合这两个群体时, ML树仍呈现出以地理群体为主导的分枝结构, 汇泉湾群体更多地处于树的基部, 而石老人群体更向顶端分枝集中(图 6,左)。这一结果表明, 竿虾虎鱼具有较强的本地化栖息习性, 群体内世代间的遗传分化较小, 整体遗传结构主要受空间因素影响。

图3 基于ryr3(a)和cytb(b)的竿虾虎鱼属贝叶斯进化树Fig.3 Bayesian trees of Luciogobius based on ryr3 (a) and cytb (b)

图4 竿虾虎鱼各谱系中国大群体(黑色)与日本大群体(白色)的cytb单倍型网络图Fig.4 Haplotype network of L. guttatus lineage A and C based on cytb

3 讨论

在我国和日本沿海, 竿虾虎鱼谱系A、C都存在同域分布的情形, 而我国竿虾虎鱼谱系A、C的分布比重改变于崂山西南侧的石老人与东侧的仰口之间。崂山是我国海岸线上海拔最高的山峰, 第四纪冰川在崂山(尤其东侧海岸)留有多处遗迹(Kusky et al,2011)。且崂山附近海底地貌复杂、坡度较大, 以致夏季表层水温在青岛沿海区域性急剧下降(王彬, 2010),且东(北)风下的表层海流在崂山附近形成涡旋(王彬,2010)。而根据近海遥感资料, 悬浮物浓度、叶绿素a浓度、硅氮比等生态指标都能在青岛近海沿南北方向形成主要梯度(窦勇, 2012)。这些历史、地理和生态背景可能是阻碍竿虾虎鱼谱系A、C在崂山附近海岸相互渗透的重要原因。而根据 GenBank上已有的竿虾虎鱼线粒体DNA能够推断, 竿虾虎鱼谱系A在北海道至本州岛南部沿海以及对马海峡西岸都有分布(Yamada et al, 2009; Jeon et al, 2012), 本研究也进一步佐证了这一海洋型谱系(Hashimoto et al, 2014)在我国北方沿海相对于偏淡水型谱系 C的广布性(没有来自韩国沿海谱系C的分子信息)。谱系B已知分布于西南日本沿海及琉球群岛, 在本研究中未能检测到。

图5 基于cytb与ryr3基因片段构建的竿虾虎鱼A、C谱系的BSP图Fig.5 Bayesian skyline plots of L. guttatus lineage A and C based on cytb and ryr3

竿虾虎鱼属内物种分类的混乱与各谱系形态的趋同进化密不可分。本研究进一步证明, 对隐存多样性高的海洋鱼类物种(如： 竿虾虎鱼), 使用高多态性的 DNA标记(如： 线粒体 DNA)可以得到比形态分类更准确的系统发育关系。同时我们也看到, 从线粒体基因与核基因上分出的竿虾虎鱼谱系虽然能够对应起来, 但分支的顺序是不一致的。核基因ryr3编码一种内质网上的钙离子通道蛋白, 主要在神经系统中表达(Nakashima et al, 1997); 在高纬度沿海更大的年温差下, 谱系C可能具备更灵活的钙离子信号的调控以促进对季节间温度变化的适应; 而谱系间所受选择压力的差异可能是竿虾虎鱼在核基因与线粒体DNA上呈现出系统发育模式差异的一个潜在因素。另外, 本研究只选用了一个核基因片段, 所得结果可能受到片段特异性选择压力的影响。为更加完整地阐明竿虾虎鱼谱系间的遗传分化, 后续工作中应选用更多的核基因序列片段进行深入研究。

查明群体遗传组成的时空变化对于全面理解竿虾虎鱼谱系内群体间遗传分化至关重要。虽然竿虾虎鱼成体不适合长距离游泳, 但其幼体浮游期一般为20—35天, 且随降温而增加, 这在其 1—2年的生活史中, 与多数洄游型鱼类相比所占的比例并不小(Kitano et al, 2003); 加之其对海水、淡水环境均能迅速适应, 故竿虾虎鱼的幼体扩散能力可能要高于预期。因此, 竿虾虎鱼幼体可能在海流的作用下发生较远的空间迁移, 同一地理群体不同时间的补充群体的遗传组成也可能存在差异。本研究结果表明, 虽然采样时间并不在同一年的同一季节, 竿虾虎鱼谱系内部群体的遗传组成在时间上仍具有稳定性。此外, 竿虾虎鱼各谱系内部的我国各群体间虽然存在一定的遗传分化, 但部分群体间可能由于幼体扩散导致高水平的基因交流, 其遗传分化未能达到显著水平。

图6 基于mcr基因片段的汇泉-石老人两群体(左)、石老人(右上)与汇泉(右下)单群体年际ML树Fig.6 ML trees of L. guttatus sampled in different years from HQ-SLR (left), SLR (upper right) and HQ (lower right) based on mcr

本研究基于线粒体 DNA的分析结果表明, 竿虾虎鱼各谱系都经历了与末次冰期相关的有效种群扩张(开始于约3—5万年前), 以及开始于日本海南部海域东岸的群体分化(开始时间更早, 约10—30万年前)：其中谱系C的中日间群体分化早于谱系A, 且更可能受到近30万年来全球海平面涨落的影响(Lambeck et al, 2014)。本研究中的核基因片段多态性较低, 且序列较少, 未能反映出末次冰期后的有效种群恢复, 但基于核基因ryr3获得的BSP图进一步支持了谱系C在末次冰期前(约 12—30万年前)有效种群数量趋于减小, 而谱系A并未体现出末次冰期前的瓶颈效应。我们推断, 偏北分布的谱系C在末次冰期初期可能经受了更强的选择压力, 而谱系 A很有可能在间冰期时随海侵扩散到青岛沿海。综观各标记的分析结果,竿虾虎鱼可能是一个物种复合体(L. guttatus agg.)：其中谱系B、C起源于第四纪初期并在冰期发生纬向隔离(谱系B可能是L. ryukyuensis的同种异名); 谱系A起源于中新世末期并广泛分布于中间纬度, 虽可能与谱系 C同域分布, 但应该加以分子鉴别并单独保护。

王 彬, 2010. 地形对西南黄海环流的影响. 青岛： 中国科学院海洋研究所硕士学位论文, 14—15

伍汉霖, 钟俊生. 2008. 中国动物志： 硬骨鱼纲、鲈形目(五)、虾虎鱼亚目. 北京： 科学出版社, 484—486

窦 勇, 2012. 基于RS、GIS调查资料的青岛市海岸带生态系统健康评价. 青岛： 中国海洋大学博士学位论文, 133—142

Avise J C, Walker D, 1999. Species realities and numbers in sexual vertebrates： Perspectives from an asexually transmitted genome. Proc Natl Acad Sci USA, 96(3)：992—995

Chen J T F, 1932. Note sur un nouveau poisson chinois appartenant au genre Luciogobius. Bull Mus Natl Hist Nat Ser 2, 4： 648—650

Cheng Y Z, Xu T J, Jin X X et al, 2012. Universal primers for amplification of the complete mitochondrial control region in marine fish species. Mol Biol, 46(5)： 727—730

Drummond A J, Suchard M A, Xie D et al, 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol, 29(8)： 1969—1973

Excoffier L, Lischer H E L, 2010. Arlequin suite ver 3.5： a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Res, 10(3)： 564—567

Hashimoto S, Koizumi I, Takai K et al, 2014. Different habitat salinity between genetically divergent groups of a worm-like goby Luciogobius guttatus： an indication of cryptic species.Environ Biol Fish, 97(10)： 1169—1177

Hebert P D N, Cywinska A, Ball S L et al, 2003. Biological identifications through DNA barcodes. Proc Roy Soc B Biol Sci, 270(1512)： 313—321

Jeon H B, Choi S H, Suk H Y, 2012. Exploring the utility of partial cytochrome c oxidase subunit 1 for DNA barcoding of gobies. Anim Syst Evol Divers, 28(4)： 269—278

Kitano T, Hatakeyama R, Akiyama N et al, 2003. Relationship between annual reproductive cycle of female flat-head goby,Luciogobius guttatus and the seasonal change in water temperature in estuaries along the northern Suruga bay.Aquacult Sci, 51(1)： 41—48

Kusky T, Guo L, Xiang S B et al, 2011. A critical examination of evidence for a quaternary glaciation in Mt. Laoshan, eastern China. J Asian Earth Sci, 40(1)： 403—416

Lambeck K, Rouby H, Purcell A et al, 2014. Sea level and global ice volumes from the Last Glacial Maximum to the Holocene. Proc Natl Acad Sci USA, 111(43)： 15296—15303

Librado P, Rozas J, 2009. DnaSP v5： a software for comprehensive analysis of DNA polymorphism data.Bioinformatics, 25(11)： 1451—1452

Mukai T, Nishida M, 2004. Intraspecific mitochondrial DNA phylogeny of a Japanese brackish water goby, Luciogobius guttatus. Jpn J Ichthyol, 51(2)： 157—161

Nakashima Y, Nishimura S, Maeda A et al. 1997. Molecular cloning and characterization of a human brain ryanodine receptor. FEBS Lett, 417(1)： 157—162

Sambrook J, Fritsch E F, Maniatis T, 1989. Molecular Cloning： A Laboratory Manual. 2nd ed. Plainview, NY： Cold Spring Harbor Laboratory Press, 3： E3—E4

Tamumra K, Peterson D, Peterson N et al, 2011. MEGA5：molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol, 28(10)： 2731—2739

Wang B, Li Y, Yuan D L, 2013. Effects of topography on the sub-tidal circulation in the southwestern Huanghai Sea(Yellow Sea) in summer. Acta Oceanol Sinica, 32(3)： 1—9

Yamada T, Sugiyama T, Tamaki N et al, 2009. Adaptive radiation of gobies in the interstitial habitats of gravel beaches accompanied by body elongation and excessive vertebral segmentation. BMC Evol Biol, 9： 145