基于线粒体控制区序列的黄海日本鳀(Engraulis japonicus)的群体遗传结构*

邹聪聪, 王丽娟, 吴志昊, 尤锋

1. 中国科学院实验海洋生物学重点实验室, 海洋大科学中心(中国科学院海洋研究所), 山东 青岛 266071;

2. 青岛海洋科学与技术试点国家实验室实验海洋生物学与生物技术实验室, 山东 青岛 266237;

3. 中国科学院大学, 北京 100049

遗传多样性是影响种群适应能力的关键因素,尤其是在多变和极端环境下(Ketchum et al, 2016),而群体结构则反映了种群在进化过程中因地理或生态变化等所产生的环境适应性分化(高天翔 等,2020)。它们是研究种群动态变化的基础, 也是生物分子系统地理学研究的重要内容, 能为可持续渔业的开发、物种管理和保护提供重要建议(Montes et al,2017)。因此, 为了解海洋生物遗传多样性的现状、系统地理结构和物种进化过程, 许多海洋生物开展了生物分子系统地理学研究。借助分子标记, 研究鱼类群体遗传变异和多样性可以有效地了解目前该鱼种的资源状况和追溯群体的历史变化(Menezes et al, 2012; Zhang et al, 2014)。线粒体DNA具有母系遗传、进化速度快等特点, 已被广泛应用于群体遗传、历史动态分析研究(Guo et al, 2014; Munian et al,2015)。其中, 线粒体控制区是一段非编码基因序列,由高变位点组成, 因其受自然选择压力小, 突变率高, 被广泛应用于群体遗传变异分析(Han et al,2008; Myoung et al, 2014)。

日本鳀(Engraulis japonicus)隶属于鲱形目(Clupeiformes)、鳀科(Engraulidae), 是一种小型远洋型的群居鱼种, 广泛分布于西北太平洋的沿岸(Ohata et al, 2018)。在中国其主要分布在渤海、黄海和东海, 是黄海和东海单种鱼类资源生物量最大的鱼种之一(朱琨 等, 2020)。日本鳀属暖水性中上层鱼类, 主要以浮游生物为食, 其在幼虫阶段被称为“海蜒”, 是海洋食肉动物的主要猎食对象, 如海洋哺乳动物和海鸟等(Checkley et al, 2017)。因此, 日本鳀在海洋食物链中起着承上启下的作用(胡东方,2009), 是黄海生态系统动力学研究的资源关键种(唐启升 等, 2000)。日本鳀具有很高的营养价值, 被称之为“鱼中珍品”, 主要用于食用和动物饲料等。20世纪90年代之前, 其捕获方式主要是由沿岸小型渔具和近海底拖网作业兼捕所得; 而在90年代之后, 由于大规模地开发利用, 由兼捕对象转为了主捕对象, 特别是在中国、日本和韩国(薛利建 等,2005)。因此, 对日本鳀的分子系统地理学研究迫在眉睫。目前, 有关我国沿海鳀鱼群体变异的研究已有些报道, 但大多集中于较大地理尺度的海域间的群体分析, 且仅限于遗传多样性和群体结构研究,对于其历史动态变化分析甚少, 如Liu等(2006)发现中国沿海包括黄海日本鳀表现出很高的单倍型多样性; Kim等(2004)利用线粒体12S核糖体序列研究了黄海南部、南海东部和济州岛西岸日本鳀的群体遗传结构, 结果显示各群体间的遗传距离并不显著。基于线粒体DNA的COI序列和cyt b序列对黄海和东海日本鳀的群体遗传变异及遗传结构进行的分析发现, 这两个群体无显著的群体遗传结构(Yu et al, 2005; Zheng et al, 2015), 后者推测其可能与浮游幼虫的扩散行为和洋流有关。而有关各海域内特别是黄海的日本鳀群体遗传结构和多样性研究较少。因此, 本文对黄海4个日本鳀地理群体进行遗传变异分析, 探究其遗传多样性、群体结构和历史动态变化。

1 材料与方法

1.1 样品采集

试验材料来自2016年8—9月渔业资源调查船“北斗号” 2016003航次在黄海海区16个站位利用底层单拖网方式采集的样品。如图1所示, 将从16个站位的采集的131尾日本鳀划分成北黄海北部(NYN,n=43, 体重=19.3±6.1g)、北黄海南部(NYS,n=36, 体重=11.6±7.5g)、南黄海北部(SYN,n=30, 体重=7.8±4.2g)和南黄海中部(SYC,n=22, 体重=3.8±2.0g) 4个地理群体。采集每个个体的肌肉组织,并于-20℃保存, 用于总DNA的提取。实验流程均符合中国科学院海洋研究所实验动物福利伦理和动物实验安全审查制度要求。

图1 黄海日本鳀样本采集的地理分布图图中表示北黄海北部群体(NYN);表示北黄海南部群体(NYS);表示南黄海北部群体(SYN);表示南黄海中部群体(SYC)。该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2020)4629的标准地图制作Fig. 1 Map of sampling locations of Engraulis japonicus in the Yellow Sea., north of the Northern Yellow Sea population (NYN);, south of the Northern Yellow Sea population (NYS);, north of the Southern Yellow Sea population (SYN);, central part of the Southern Yellow Sea population (SYC)

1.2 试验方法

将采集样品根据形态特征, 参考鱼类检索(成庆泰 等, 1987)和Fishbase (http://www.fishbase.org)进行物种鉴定。使用DNA提取试剂盒(中国, 天根)提取日本鳀个体的总DNA, 用1%琼脂糖凝胶电泳对提取的DNA进行检测。使用Nanodrop 2000分光光度计(美国, Thermo)测定核酸浓度和纯度, 然后于-20℃下保存。根据NCBI日本鳀线粒体控制区序列(GenBank序号: NC_003097.1), 利用Primer 3软件(https://www.yeastgenome.org/primer3), 设计了一对引物(EJ-CR-F1: CAAAGGGAAGAGGCTTGAAC和EJ-CR-R1: TGCTTACGGAACTTTTCAGGGT)用于PCR扩增。PCR扩增的反应体系为50μL, 包括25μL 2× Es taq酶(中国, 康为世纪), 2mM正向引物/反向引物, 50～100ng DNA模板, 超纯水补充加至50μL。PCR程序为: 预变性94℃ 2 min, 变性94℃ 30 s, 退火55℃ 30 s, 延伸72℃ 1min, 其中预变性的循环次数为1个, 变性、退火和延伸的循环次数为35个。最后在72℃下进行终延伸5min。PCR产物在1.2%琼脂糖凝胶上进行检测, 利用胶回收试剂盒(美国,Omega)进行纯化。纯化后的PCR产物送至青岛擎科梓熙生物技术有限公司, 利用扩增引物进行测序,并设计中间引物双向测通。

1.3 数据处理

利用ClustalX 2.1 (Larkin et al, 2007)进行序列拼接和比对, 并在NCBI上进行同源检测。用Arlequin 3.5 (Excoffier et al, 2010)分析群体的核苷酸组成、单倍型数量、单倍型多样性(h)、多态性位点数、碱基颠换、转换、插入、核苷酸多样性(π)和遗传分化系数Fst值。将4个群体按地理位置进行种群划分: 一个基因池(所有群体), (NYN、NYS)和(SYN、SYC)两个基因池, (NYN、NYS、SYN)和(SYC)两个基因池, 通过分子方差分析(Analysis of Molecular Variances, AMOVA)评估群体间的遗传变异。中性检验推测群体是否发生种群扩张事件。利用R语言进行主成分判别分析(Discriminant Analysis of Principal Components, DAPC)。基于赤池信息准则(Akaike Information Criterion, AIC)准则,在jModeltest 2.1 (Darriba et al, 2012)中选择最佳拟合的核苷酸替代模型, 模型参数为HKY (Hasegawa-Kishino-Yano)+I+G, 其中I=0.835、G=0.409。通过MrBayes 3.2 (Ronquist et al, 2012)构建贝叶斯系统发育树。利用PopART 1.7 (Leigh et al, 2015)绘制单倍型网络图。根据公式:τ=2tμ, 式中τ是扩张参数,μ是控制区的突变率(单位: Ma-1),t是自扩张以来的时间(单位: Ma), 在错配分布中评估种群的扩张时间。由于日本鳀的突变速率未曾报道,根据鱼类的突变速率(Donaldson et al, 1999), 本研究选取了3%·Ma-1作为日本鳀的突变速率。同时利用Beast 2.5 (Bouckaert et al, 2014)绘制群体的贝叶斯天际线图, 以评估有效种群大小随时间的变化, 使用Beast 2.5的Yule模型评估谱系之间的分歧时间。

2 结果

2.1 遗传多样性与系统发育分析

对4个地理群体的131尾鱼的线粒体控制区序列与基因组序列进行比对, 获得了长度为1024bp的全部线粒体控制区片段。在日本鳀线粒体基因组中, 该片段位于15652～16675bp区域。131个序列中,碱基的组成分别为A 29.40%、T 30.78%、C 23.56%和G 16.26%, 其中, 单倍型数量为126个, 多态位点数为108个, 碱基转换数为85个, 碱基颠换数为38个和碱基插入数为4个(表1)。4个群体的单倍型多样性均很高, 其中NYS群体的最低(0.995±0.009),SYC群体的最高(1.000±0.014)。而核苷酸多样性均较低, 为0.010±0.005～0.011±0.006。NYN群体的多态位点数最多, 为71个; SYC群体的多态位点数最少, 为50个。在4个群体中, 仅NYS和SYC群体有碱基插入情况。根据单倍型网络图可知, 126个单倍型被明显分成了两个分支, 并且它们与地理位置无显著关系(图2)。在系统发育树中, 黄海日本鳀的单倍型同样分化成了2个谱系(图3)。

图3 黄海日本鳀126个单倍型的贝叶斯系统发育树Fig. 3 Bayesian inference phylogenetic tree of 126 haplotypes of Engraulis japonicus in the Yellow Sea

表1 日本鳀4个地理群体的遗传多样性参数Tab. 1 Genetic diversity parameters in four Engraulis japonicus populations

图2 黄海日本鳀的单倍型网络图分支上的短线表示单倍型间的突变Fig. 2 The haplotype network of Engraulis japonicus in the Yellow Sea. The short lines on the branches indicate mutation step between haplotypes

2.2 种群分化与结构

日本鳀4个地理群体两两间的Fst值均很小(0.010～0.018)(表2)。在AMOVA分析中, 4个群体按地理位置进行了3个类型的种群划分。结果显示, 3个类型的种群划分类型均无明显的群体结构, 且分子变异主要来源于群体内(表3)。DAPC分析结果与上述结果一致, 4个群体比较集中, 无明显的群体分离现象(图4)。

图4 日本鳀4个地理群体的主成分判别(DAPC)分析Fig. 4 DAPC analysis of four Engraulis japonicus populations in the Yellow Sea

表3 分层结构的分子方差分析Tab. 3 Analysis of Molecular Variance (AMOVA) in hierarchical structure

2.3 历史动态分析

如表4所示, 4个群体的Tajima’sD值均为负值,且p值均大于0.05。谱系1和谱系2的Tajima’sD值为负值, 但p值均小于0.05。群体和谱系的Fu’sFs值从–24.790到–13.580,p值均为0.000。Ramos-Onnsis and Roza’s检验结果显示, 无论是群体还是谱系,R2的p值均为0.000。以上结果表明, 4个群体和2个谱系可能发生过群体扩张。

表4 日本鳀4个地理群体的中性检验Tab. 4 Neutrality tests of four Engraulis japonicus populations in the Yellow Sea

群体扩张的错配分布结果显示, 仅谱系1的扩张后群体大小θ1(99999.000)值远大于初始群体大小θ0(0.228), 且群体和谱系的偏差平方和SSD值和Harpending’s Raggedness指数HRI均无显著差异(表5)。根据扩张参数τ值评估群体的扩张时间,结果表明, 谱系1和谱系2分别于0.135Ma和0.164Ma前发生了群体扩张, 这可以追溯到更新世时期。在空间扩张模型下, SYC、SYN和谱系1的M值远大于θ。群体和谱系的SSD值和HRI也无显著差异。如图5显示, 谱系1和谱系2在错配分布中均呈现为单峰。贝叶斯天际线图显示, 根据线粒体控制区序列, 黄海日本鳀的有效群体数量随着时间的推移并不稳定, 在0.150Ma前出现了明显的下降(图6)。

图5 谱系1 (a)和谱系2 (b)的错配分布图Fig. 5 Mismatch distributions of lineage 1 (a) and lineage 2 (b)

图6 黄海日本鳀的贝叶斯天际线图实线表示中值, 阴影区域表示95%的最高后验密度Fig. 6 Bayesian skyline plot of Engraulis japonicus populations in the Yellow Sea. The solid line represents median values,and the shading shows the 95% highest posterior density

表5 群体扩张和空间扩张模型下的错配分布Tab. 5 Mismatch distribution under demographic expansion and spatial expansion models

3 讨论

群体的遗传多样性和遗传结构是群体动力学研究的重要内容(Domingues et al, 2018), 也是渔业管理和物种保护的基础(Nykänen et al, 2020)。日本鳀广泛分布于中国沿海地区, 是黄海单种鱼类资源生物量最大的鱼种之一(朱琨 等, 2020), 也是其生态系统动力学研究的资源关键种(唐启升 等, 2000)。但目前针对黄海海域日本鳀的群体遗传研究较少。本研究利用线粒体控制区序列, 对黄海日本鳀进行了遗传多样性和群体结构分析。结果显示, 黄海日本鳀具有较高的遗传多样性, 4个地理群体的总单倍型多样性达到了0.999±0.001, 这与Yu等(2005)利用线粒体DNA的COI序列和cyt b序列检测的黄海日本鳀遗传多样性水平结果一致, 其单倍型多样性水平很高, 也达到0.934±0.061～0.981±0.031。这也与许多海洋鱼类如分布在西大西洋、印度洋和红海的丝鲨(Carcharhinus falciformis)(Clarke et al, 2015)、智利巴塔哥尼亚和阿根廷海岸的智利黍鲱(Sprattusfuegensis)(Canales-Aguirre et al, 2018)和马六甲海峡、中国南海和西里伯斯海的青甘金枪鱼(Thunnus tonggol)(Kasim et al, 2020)等一样具有很高的单倍型多样性。不仅如此, 在我国沿海海域的许多鱼种也表现出很高的单倍型多样性, 如小黄鱼(Larimichhys polyacis)(Zhang et al, 2017)、黑鲷(Acanthopagrus schlegeli)(赵爽 等, 2010)等。除此之外, 鲱形目的许多鱼种也具有很高的单倍型多样性,如分布在苏鲁西里伯斯海的沙丁鱼(Sadinella lemuru), 单倍型多样性达0.9981±0.0050～1.0000±0.0041 (Pedrosa-Gerasmio et al, 2015), 西北太平洋沿岸的鳓(Ilisha elongata)平均单倍型多样性达0.98 (Wang et al, 2016)。有研究表明, 群体内维持较高的单倍型多样性可能与庞大的群体数量, 环境的异质性和生活史有关(Nei, 1987)。尽管本研究中4个地理群体均表现出很高的单倍型多样性, 但反映群体中所有单倍型对之间平均差异的核苷酸多样性这一敏感指标(Ma et al, 2010; Sun et al, 2019)却都很低, 仅为0.010±0.005。根据Grant和Bowen分类表,这种高单倍型多样性和低核苷酸多样性可能与种群的扩张有关(Grant et al, 1998), 种群数量因环境变化急剧增加, 同时碱基的突变使得单倍型数量和单倍型多样性增加, 而核苷酸多样性没有获得足够的时间积累(梁述章 等, 2019), 故鳀推测黄海日本群体可能经历了群体扩张。然而, 虽然本研究结果显示, 黄海日本鳀的单倍型多样性水平很高, 说明目前资源仍较丰富, 但有研究表明, 黄海日本鳀的常规生物量正在下降(Jung et al, 2016)。因此, 合理科学的物种保护措施仍需实施。

AMOVA分析结果和Fst值显示, 黄海4个日本鳀群体间无明显的群体分化。在群体的动态变化中,许多因素都可能影响海洋种群的群体结构, 如海洋动物的生活史、洋流和历史变化等(Zhang et al, 2014;Cheng et al, 2018)。在海洋环境中, 许多物种都需经历较长的浮游阶段, 而较强的浮游能力可以促进群体间的基因交流, 从而使得原有的地理距离被克服,产生遗传上的同质性(Otwoma et al, 2018)。尤其是浮游鱼类, 在其早期生命阶段, 卵和幼体都会被动地被洋流所运送(Xing et al, 2020)。有研究指出, 日本鳀的卵和仔鱼会受洋流影响漂浮至较温暖的沿海和近海水域(Hayashi et al, 2016)。而且大约在产卵后2周内, 大部分产的卵和分布在日本南部黑潮周围新孵出的仔鱼会从黑潮及其下游区域朝东北方向输送至过渡区(Takahashi et al, 2005)。Zheng等(2015)发现, 从黄海和东海采集的日本鳀样品间无明显的群体结构, 推测可能是浮游幼体强大的扩散能力和洋流系统所致。黄海水域中分布着多种表层流, 如黄海沿岸流、苏北海流和韩国沿岸流等(Yu et al,2020)。而表层环流受风的影响呈现出季节性的变化特点, 这使得浮游幼体的被动迁移也出现季节性的变化, 因此, 浮游物种群体间的联系在一定程度上取决于产卵的时间(Sefc et al, 2020)。日本鳀属于洄游性鱼类, 会在产卵场(3—6月)、索饵场(7—10月)和越冬场(12—2月)之间有节律地做季节性洄游(赵宪勇, 2006; Pandey et al, 2017)。每年初春时分, 成体鳀鱼就会开始从越冬场迁移至沿海的产卵场(Kwon et al, 2012)。当海水温度适宜时, 就会促使日本鳀迁移(Zhou et al, 2015)。黄海海域日本鳀的产卵季节较长, 从4月持续至10月, 5月到7月为产卵的盛期(Yu et al, 2020)。而黄海中部近海域是黄海主要的产卵场之一, 也是黄海日本鳀产卵场的重要组成部分(李秀梅 等, 2016)。因此, 这些原因都可能使得黄海海域的日本鳀表现出遗传上的均一性, 本研究也证实了这一点。另外, Liu等(2006)、Yu等(2005)和Zheng等(2015)都对日本鳀展开了群体遗传研究, 研究海域主要集中于较大地理尺度的海域间群体分析。虽然他们也分析了黄海的样品, 但采样点及样品数量少, 并且Yu等(2005)和Zheng等(2015)使用的分子标记是COI和cyt b。本文使用的分子标记是序列变异更高的线粒体DNA控制区序列, 广泛应用于群体遗传和历史动态分析。本文研究结果显示, 黄海海域日本鳀无明显的群体结构。Yu等(2005)和Zheng等(2015)与Liu等(2006)也分别发现, 黄海与东海日本鳀、中国沿海与日本海域日本鳀同样无显著的群体结构。综上所述, 西北太平洋日本鳀可以作为一个渔业单位进行相关资源的管理和保护。

根据贝叶斯进化树的结果, 黄海日本鳀分化出了两个谱系, 虽然两个谱系的置信值比较低, 但是单倍型网络图中也产生了两个明显的谱系。中性检验分析结果显示, 黄海日本鳀2个谱系均具有显著的Tajima’sD、Fu’sFs和Ramos-Onnsis and Roza’sR2值, 说明其可能经历过突然的种群扩张、基因流动、方向选择或平衡选择等(Wang et al, 2016)。结合错配分布图呈现出的单峰图, 本研究推测谱系1和谱系2可能发生过群体扩张, 而扩张时间分别为0.135Ma和0.164Ma前, 即在更新世。并且两个谱系的分化时间为0.701Ma前, 也可以追溯到更新世(2.60Ma—0.012Ma前)。研究发现, 许多海洋生物在更新世都会经历群体瓶颈、群体扩张等变化, 如东海近岸的带鱼(Trichiurus japonicus)(吴仁协 等, 2019)、马来西亚沿海的长尾金枪鱼(Thunnus tonggol)(Kasim et al, 2020)和印度—西太平洋中部的蓝圆鲹(Decapterus maruadsi)(Jamaludin et al, 2020)等。Liu等(2006)也发现, 在更新世期间, 中国沿岸和日本沿岸的日本鳀和澳大利亚鳀(Engraulis australis)发生了群体分化。同样地, 我们的研究结果也发现黄海日本鳀的谱系分化发生在更新世。更新世期间的显著特点就是冰期和间冰期的循环交替, 导致海平面和气候发生巨大变化, 从而使得群体间发生分化(Pillans et al, 1998; Inoue et al, 2014)。因此, 更新世对于海洋动物包括日本鳀的进化、分布以及群体结构会造成很大影响(Liu et al, 2007; Shen et al, 2011)。除此之外, 本研究还发现在0.150Ma前, 日本鳀有效群体数量开始下降, 这可能与更新世期间海平面降低导致沿岸海洋物种栖息地的减少有关。