中国北鳅谱系地理学分析

杨银盆 杨慧兰 王银肖 谭慧敏 陈咏霞

(河北大学生命科学与绿色发展研究院, 保定 071002)

生物多样性及物种的分布格局都与冰期及期间形成的避难所密切相关[1]。对当今物种分布影响较大的是始于3—2 Ma前, 结束于0.02—0.01 Ma前的第四纪冰期[2]。第四纪冰期的最大特点是冰期和间冰期往复震荡, 气候出现冷暖周期性变化, 强烈地影响了种群的历史动态和群落结构[3]。冰期来临时, 栖息于北部的温热生物向赤道方向迁移, 或沿着冰缘带迁移, 或躲藏于避难所; 此时, 海平面下降, 大陆架裸露, 为不同区系物种交流提供了通道[4]。间冰期时, 随着气候变得温暖和湿润, 适应于冷水的生物则向北迁移, 并扩大其分布范围[5]。第四纪冰期的这种往复变动引起种群分布范围及种群数量的扩张和收缩, 种群遗传结构受其历史动态变动而发生改变。对现生生物种群遗传结构及其地理分布的研究, 了解第四纪冰期对现生生物种群结构的影响以及生物避难所的确定是当前生物学领域研究的热点问题, 可以为生物多样性保护提供依据。

北鳅(Lefua costata)隶属于鲤形目(Cypriniformes)、条鳅科(Nemacheilinae), 为古北区一种冷水性鱼类, 在中国、韩国、日本、蒙古及俄罗斯东部等均有分布[6], 其中, 我国主要分布于淮河以北的东部地区。北鳅属现知物种有7种, 我国仅北鳅1种。

北鳅为小型底栖鱼类, 迁徙能力差, 对环境依存度高, 多见于溪流的上游, 喜低温, 受人类活动影响小, 其分布格局可以很好地反映环境与地质演变过程, 是探讨冰期气候变迁历史的理想物种。已有研究探讨了北鳅属的起源、分化及属内物种的亲缘关系。基于线粒体全基因组序列, 北鳅属大约在26 Ma前与高原鳅属等条鳅科鱼类发生分化[7]。Mihara等[1]认为北鳅属鱼类的祖先起源于亚欧大陆北部, 向南经北海道扩散至朝鲜半岛, 基于线粒体D-loop区序列分析认为北鳅属的分歧时间在(1.91±0.64) Ma前; 其中, 北鳅和短体北鳅(Lefua nikkonis)的分歧时间在(1.50±0.58) Ma前, 且二者的亲缘关系最近。而同样基于D-loop区序列分析, Shedko等[8]则认为北鳅属的分歧时间在11—12 Ma前, 北鳅的分歧时间在3 Ma前。目前, 对于北鳅遗传结构及生物地理学格局研究比较匮乏。基于线粒体D-loop区序列, 本文分析我国北鳅的遗传多样性, 并对冰期避难所及北鳅间冰期可能的扩散路径进行相关推测。

1 材料与方法

1.1 样本采集

样本于2016—2019年采自9个水系18样点(图1),共211尾, 所分析的样本涵盖了我国北鳅的主要分布范围, 按地理分布划分为7个组: 淮黄河(黄河水系和淮河水系, HH种群)2个样点: 山东省西营镇锦绣川水库、河南省汝州市蟒川镇汝河, 28尾; 海滦河(滦河水系和海河水系, HL种群)7个样点: 河北省闪电河、河北省府河、河北省拒马河、天津张庄水库、河南省狮豹头水库、河南省小南海水库、河南省塔岗水库, 共89尾; 辽河水系(LH种群)3个样点: 内蒙古赤峰市西辽河、辽宁省彰武县柳河、辽宁省朝阳市大凌河, 共29尾; 辽东半岛水系(LD种群)2个样点: 辽宁省瓦房店复州河、辽宁省普兰店沙河, 共33尾; 大洋河水系(DY种群)1个样点: 辽宁省圆台大沙镇大洋河, 共19尾; 鸭绿江水系(YL种群)2个样点: 辽宁省凤城和吉林省临江, 共5尾; 松花江水系(SH种群)1个样点: 黑龙江省桦南县八虎力河, 共8尾, 后文中不同种群均简称表示。标本保存于河北大学博物馆, 储存于-20℃冰箱保存。

图1 北鳅样本采集地点Fig.1 Map of the sample collection sites of Lefua costata

1.2 DNA提取、PCR扩增和测序

样本取自鱼的胸鳍、背鳍或肌肉组织, 剪碎后使用基因组提取试剂盒DP304(北京天根生化科技有限公司)进行提取。D-loop区扩增的引物序列为DL1 (3′-ACCCCTGGCTCCGAAAGC-5′)和DH2 (3′-ATCTTAGCATCTTCAGTG-5′)[9]。PCR反应条件为: 95℃预变性3min; 95℃变性30s, 退火30s温度为51℃, 72℃延伸40s, 循环30次; 最后72℃延伸10min,4℃停止反应。PCR产物用1.00%—2.00%琼脂糖凝胶电泳对其条带进行检测, 送往通用生物系统(安徽)有限公司完成序列测定。

1.3 数据处理和分析

利用BioEdit 7.0.9[10]软件进行DNA序列比对和拼接; 利用MEGA 5.05[11]软件对序列变异情况、碱基组成进行分析, 并计算转换/颠换比率等相关系数, 基于Kimura-2参数模型计算各地理群体间的遗传距离。利用DNAsp 5.1[12]软件进行遗传多样性相关分析。利用Arlequin 3.5[13]软件计算群体间分化系数(FST), 以及中性检验及错配分析等种群历史动态分析。

以细尾高原鳅(Tripiophysa stenura)、董氏须鳅(Barbatula toni)、短体北鳅(L.nikkonis)、斑北鳅(L.echigonia)和L.sp.作为外来群构建北鳅的系统发育树, 下载序列GenBank登录号见表1。利用jModel-Test 0.1.1[14]软件筛选构建系统发育树的最佳进化模型, 为HKY+G模型。利用PhyloSuite[16]中的MrBayes程序[17]构建单倍型贝叶斯(BI)系统发育树, 4条马尔科夫链同时运行1×106代, 通过后验概率评估每个节点, 每100代对系统树进行抽样。利用IQtree[15]构建单倍型最大似然(ML)系统发育树, 对每个节点的置信度对评估通过有1000个重复的bootstrap进行。利用PopART[18]软件采用TCS法构建单倍型网络关系图。利用BEAST 2.0[21]软件推算分歧时间, 采用花鳅属(Cobitis)鱼类化石(约23 Ma前)[19]和高原鳅属(Tripiophysa)鱼类化石(约13 Ma前)[20]进行时间校准, 下载序列GenBank登录号见表1。种群扩张时间的评估依据线粒体D-loop区每百万年1.00%的平均进化速率计算[8], 公式为τ=2ut[22],u为每条序列每个世代的变异速率, 根据u=2μk计算,μ为平均进化速率,k为序列长度。北鳅属鱼类为一次产卵类型,因此按1世代/年计算[23]。

表1 GenBank下载序列Tab.1 List of the sequences from GenBank

2 结果

2.1 序列特征与遗传多样性

9条水系211条北鳅D-loop区序列平均长度为925 bp: 保守位点823个, 约占序列总长度的88.97%;变异位点102个, 约占序列总长度的11.02%; 简约信息位点61个, 约占序列总长度的6.70%。序列转换颠换比(Ti/Tv)为1.88。碱基含量T为32.50%, C为18.61%, A为35.90%, G为12.99%, 表现为G偏移。A+T含量(68.40%)明显高于C+G含量(31.60%)。

在211个样本中检测到55个单倍型, 地理种群间共享单倍型有5个, 单倍型多样性(h=0.9304)和核苷酸多样性(π=0.0087)均较高(表2)。7个种群中LD种群单倍型多样性和核苷酸多样性均最高(h=0.8920,π=0.0074), SH种群(h=0.2500,π=0.0006)最低。7个种群的群体间遗传分化指数(FST)为-0.0026—0.6870(表3), 除LH种群与SH种群之间FST不显著外, 其余种群间均有一定程度的分化,FST均为显著(P≤0.01), 其中LD种群与其他地理种群之间FST均较高, 为0.4853—0.6679。7个种群间遗传距离(K2P)为0.0019—0.0158, 其中LD种群与其他地理种群间的遗传距离均大于其他地理种群间的遗传距离, 为0.0139—0.0158, 其次是YL种群与其他地理种群间遗传距离较大, 为0.0079—0.0139,DY种群与其他地理种群间遗传距离最小, 为0.0030—0.0158; 种内遗传距离最大的是LD种群 (0.0074),其次是YL种群(0.0053), SH种群最小(0.0006)。7个种群内的遗传差异(66.61%)大于种群间(33.39%);根据LD种群与其他6个种群间的分化程度, 将其与其他6个种群划分为2个组分析显示, 组间遗传差异(30.12%)大于组内遗传差异(29.49%), 而种群内遗传差异最大, 为70.51%(表4)。

表2 北鳅遗传多样性、中性检验、错配分析与扩张时间Tab.2 Genetic diversity, neutrality tests, mismatch distribution and population expansion time of Lefua costata

表3 北鳅群体间的遗传分化系数FST(对角线下方)、群体间的遗传距离(对角线上方)及群体内遗传距离(对角线上)Tab.3 Pairwise FST (below diagonal), genetic distances among populations (above diagonal) and genetic distances within populations (on diagonal) of Lefua costata

2.2 系统发育与单倍型网络图分析

以细尾高原鳅和董氏须鳅为外来群, 构建的ML和BI单倍型系统发育树, 其拓扑结构基本一致(图2)。我国北鳅9条水系的55个单倍型以高置信度(ML=96, BI=1)聚为一支, 并与短体北鳅形成姐妹群关系。9条水系的单倍型未能各自聚类, 而是嵌套分布, 形成3大分支(分支A、B和C)。分支A主要以LD种群的15个单倍型为主, 嵌套分布有HH种群单倍型Hap1和LH种群单倍型Hap3, 位于北鳅系统树的根部; 其余单倍型以略高的置信度(ML=91,BI=0.75)聚为一支, 形成分支B和分支C两个姐妹群,其中分支C主要以HH种群和HL种群为主, 嵌套分布于有LD种群单倍型Hap41和YL种群单倍型Hap5; 分支B则包含所有种群的单倍型。

表4 群体内和群体间变异的分子方差分析(AMOVA)Tab.4 Analysis of molecular variance (AMOVA) among populations and within the population of Lefua costata

单倍型网络图(图3)显示55个单倍型整体呈复杂的网状结构, LD种群Hap39、Hap36位于网状中心, 表明7个种群以LD种群为中心经历过复杂的基因交流, 结合系统发育树的结果, 将其分成3大分支。分支A以LD种群为主, 呈星状或辐射状, 单倍型Hap31位于中心位置, 表明此单倍型的祖先经历过种群扩张; 分支B包含所有种群的单倍型和种群间共享单倍型Hap3、Hap6、Hap8和Hap9, 形成复杂的网状结构, 共享单倍型均处于网状结构中心,表明7个种群经历过复杂的基因交流; 分支C以HH种群和HL种群为主, 呈星状或辐射, 单倍型Hap11位于中心位置。基于化石校准点推算出分支A与分支(B+C)的分歧时间约为2.6241 Ma (95%HPD:0.9264—4.6394 Ma)前, 分支C和分支B的分歧时间约为1.5808 Ma (95%HPD: 0.4686—3.0695 Ma)前(图4)。

图2 北鳅单倍型ML系统发育树(节点上方为ML支持率, 下方为BI后验概率)Fig.2 ML phylogenetic tree constructed by Lefua costata haplotype (bootstrap support for ML above the branch and posterior probability for Bayesian inferences below branch)

图3 北鳅单倍型网络进化图Fig.3 Haplotype evolution network diagram of Lefua costata

2.3 种群动态分析

中性检验分析显示, Tajima’sD值仅LH种群呈显著负值(P＜0.05, 表2); Fu’sFS值仅LD种群和DY种群呈显著负值, 表明LH种群、LD种群和DY种群可能经历过种群扩张事件。错配分析中,Hri检验结果均为不显著(P＞0.05),SSD检验结果仅LH种群显著偏离了种群扩张模型, 表示7个种群均经历过种群扩张事件。错配分布图显示(图5), DY种群呈单峰分布曲线; HH种群呈波浪多峰曲线; HL种群、LH种群、LD种群、YL种群和SH种群呈锯齿状下降曲线, 表明DY种群经历过种群扩张事件,其他地理种群未出现种群扩张现象。参考北鳅属鱼类D-loop区序列1.00%/百万年的平均进化速率,推算出北鳅不同种群扩散时间为2.3784—0.0477 Ma前(表2)。

图4 基于松散分子钟构建的北鳅分歧时间Fig.4 Divergence times among Lefua costata derived from the Bayesian relaxed-molecular

3 讨论

3.1 北鳅种群遗传多样性和种群分化

现生生物种群的遗传结构和分布区系与第四纪冰期和间冰期的往复更迭造成的冷暖气候变化密切关联, 冰期和间冰期引起的海进和海退为不同区系物种迁徙以及基因交流提供通道。通常在冰期降温时, 温热性鱼类往往会缩小其分布范围, 而冷水性鱼类则会扩大其分布范围, 而间冰期气温回暖时, 温热性鱼类的分布范围会扩大, 冷水性鱼类的分布范围则会缩小[24]。种群遗传结构分析能反映生物对自然选择及环境变迁的适应能力和演化潜力[25], 种群遗传多样性越高, 适应环境能力越强,承受的生存压力越小。群体遗传多样性衡量指标则有赖于核苷酸多样性指数和单倍型多样性指数[26]。种群经历过小种群事件, 如奠基者效应或瓶颈效应,种群遗传结构通常表现为低的遗传多样性[27]; 而种群经历过长期积累遗传分化, 如种群处于物种起源中心或冰期避难所, 种群遗传结构则表现为高的遗传多样性[7,28]。在本研究中, 核苷酸多样性和单倍型多样性最高的是辽东半岛种群(LD), 距离辽东半岛越远的种群, 其种群的单倍型多样性与核苷酸多样性也普遍越低, 如距离辽东半岛较近的大洋河种群(DY)和鸭绿江种群(YL)遗传多样性较高, 仅次于辽东半岛种群(LD), 距离辽东半岛较远的松花江种群(SH)遗传多样性则最低。种群分歧时间推算显示我国北鳅7个地理种群中, 辽东半岛种群(LD)与其他地理种群分化时间最早, 这表明我国北鳅种群起源地之一为辽东半岛或辽东半岛种群(LD)受第四纪冰期影响较小, 该种群经历过长期遗传分化的积累, 其他地理种群则是以辽东半岛种群为中心在第四纪冰期和间冰期时向南或向北扩散类群。

从种群间的遗传分化程度来看, 辽东半岛种群(LD)与其他地理种群的遗传分化程度最为显著, 且距离辽东半岛近的种群与其他地理种群的分化程度高, 而距离辽东半岛远的种群与其他地理种群的分化程度低, 这表明我国北鳅以辽东半岛为中心在末次冰期向南扩散的种群和间冰期向北扩散的种群积累的遗传变异最少。北鳅为底栖小型鱼类, 迁移能力弱, 且为冷水性鱼类, 多生活在山涧溪流中,对环境依存度高, 在末次冰期后不同地理群体之间无基因交流, 表现为群体内遗传变异大于群体间遗传变异, 符合距离隔离模型(Isolation by distance,IBD)理论[29], 即远距离分开的不同个体群由于有限距离的基因迁移而产生的遗传分化[30]。

3.2 北鳅种群历史动态

北鳅属鱼类的祖先种最早起源于亚欧大陆北部, 始于渐新世约26 Ma前[7], 经东亚北部穿过对马海峡和鞑靼海峡迁移至日本相关岛屿, 向南经北海道扩散至朝鲜半岛[1]。在我国, 与朝鲜半岛相邻的吉林省延边图们江有较北鳅物种原始的新记录种——史氏北鳅(Lefua pleskei) (待发表), 及我国北鳅种群中与朝鲜半岛距离最近的辽东半岛种群(LD)和鸭绿江种群(YL)具有较高的遗传多样性, 而越往北和往南北鳅种群遗传多样性降低, 表明我国北鳅属鱼类的起源地是由朝鲜半岛迁徙而来。

图5 核苷酸错配分布图Fig.5 The nucleotide of mismatch distribution

北鳅分化时间约在3 Ma前[9], 处于晚第三纪上新世晚期, 我国北鳅鱼类主要分布于朝鲜半岛邻近地区, 如鸭绿江和辽东半岛地区等地。在第四纪冰期时, 由于辽东半岛东侧地势由东北向西南渐低的长白山余脉千山山脉减轻了冰期气候对生物的影响, 为生物提供了避难场所, 因此分布于辽东半岛的北鳅鱼类受冰期影响小, 能够长期积累遗传分化,具有较高的单倍型多样性和核苷酸多样性; 其他地理种群随冰缘带借助黄渤海平原[31]的河口向南及周围扩散, 南迁至海河、黄河等, 西迁至大凌河、滦河等。文献报道, 渤海是第四纪才下陷的, 在渤海下陷之前, 古黄河、古辽河及两河中间水系曾为同一水系[32], 渤、黄海经历几次海进海退, 为北鳅鱼类在冰期和间冰期迁徙及基因交流提供了通道。本文研究与地质历史事件相吻合, 以辽东半岛种群(LD)为主的分支A分化时间(4.6394—0.9264 Ma前)发生在第一次冰期期间[33], 为我国北鳅较原始的类群, 分支A所含有的淮黄河种群(HH)单倍型Hap4、辽河种群(LH)单倍型Hap27及海滦河种群(HL)共享单倍型Hap19, 表明辽东半岛种群(LD)往南往西扩散至黄河、滦河和辽河干流等。

冰期避难所存在多个[34], 北京、河北等地在第四纪冰川期经历过四次持续性且分布范围相似的冰期[32], 燕山山脉东侧由滦河堆积形成的雾灵山为末次冰期的避难所之一[33,35]。本文研究支持雾灵山为冰期避难所之一。以海滦河种群(HL)为主的分支C分化时间(3.0695—0.4686 Ma前)约在第一次间冰期和第二次冰期, 结合海滦河种群(HL)与辽东半岛种群(LD)的共享单倍型Hap19(位于分支A), 与大洋河种群(DY)、淮黄河种群(HH)的共享单倍型Hap3和Hap9(位于分支B), 及海滦河种群(HL)较高的核苷酸多样性和较低的单倍型多样性, 表明海滦河种群(HL)是在第一次冰期或间冰期由辽东半岛种群(LD)扩散而来, 在以后的几次冰期和间冰期中经历了瓶颈效应后种群迅速扩张; 同时与大洋河种群(DY)、淮黄河种群(HH))的基因交流, 表明海滦河种群曾发生过不同地理种群的二次接触。

末次冰期后, 北鳅往北扩大其分布范围才迁徙至西辽河、东辽河以及松花江水系。因为在早更新世和中更新世期间, 辽河古水系(东辽河和西辽河)及松花江古水系流向松辽平原北部的松辽古大湖[36,37], 晚更新世时, 由于松辽分水岭的上升和形成, 松辽古大湖消失[38—41], 末次冰期黄、渤海海平面下降, 大陆架裸露, 辽河干流袭夺了东辽河、西辽河, 形成辽河水系[31,42]; 另外松辽分水岭的形成使得原本作为西辽河上源的嫩江与第一松花江合并后在三江平原与黑龙江干流合并[31], 最终形成现今松花江水系格局。在末次冰期后, 随着辽河水系的形成和松花江的改道为北鳅往北迁徙提供了通道, 最终形成现今的辽河种群(LH)和松花江种群(SH)。淮河水系的北鳅起源于黄河水系, 是在黄河夺淮[43]事件后, 淮河水系才有了北鳅的分布, 由于隔绝时间短, 黄河水系和淮河水系北鳅仍存在基因交流。