三个地理种群的董氏须鳅遗传多样性及种群历史动态

陈 浩 杨银盆 张 慧 陈咏霞

(河北大学生命科学学院, 保定 072000)

董氏须鳅(Barbatula toni), 隶属于鲤形目(Cypriniformes)、条鳅科(Nemacheilidae), 喜栖息于水质清澈、具砂石底质的山涧溪流中, 是一类小型喜冷水的底栖鱼类, 分布于俄罗斯、蒙古及我国的东北地区[1,2]。须鳅鱼类肉质鲜美、营养丰富, 经济价值高, 在我国北方尤其是东北地区深受欢迎, 也是我国北方河流的主要渔获物之一, 具有重要养殖潜力的淡水鱼类[3,4]。据调查, 由于过度捕捞以及栖息环境的破坏等影响, 须鳅鱼类自然资源近年来不断下降。

线粒体DNA是分析种群遗传结构、系统地理格局与历史的常用分子标记[5]。线粒体Cytb和D-loop是最常采用的分子标记, 其进化速率较快, 适用于种间及种内分析, 在高原鳅属(Triplophysa)、花鳅属(Cobitis)、荷马条鳅属(Homatula)等鱼类的多样性研究中均有应用[6—9]。关于须鳅鱼类的研究, 国内外主要集中在形态学分类研究[1,2,10], 仅有少数学者利用DNA分子标记对其不同群体进行研究, 例如Jasminca等[11]利用线粒体COⅠ基因作为遗传标记, 对欧洲须鳅属鱼类Barbatula barbatula的上游群体以及下游洞穴内群体进行分析, 发现尽管2个群体有着明显的形态差异, 但属同一物种, 并且洞穴群体具有较低的遗传多样性。目前对于亚洲地区须鳅鱼类的群体遗传结构及遗传多样性的研究尚未见报道。

本研究将采用线粒体Cytb和D-loop基因序列,分析分布于我国北方地区的董氏须鳅的地理格局;以鳅科鱼类Cytb的进化率, 估算种群分歧时间及种群扩张时间, 探讨种群演化历史及扩散方式, 为我国鳅超科鱼类多样性研究及资源保护提供科学依据。

1 材料与方法

1.1 研究材料

样本于2016—2017年分别采自内蒙古赤峰市阴河(42°34′N, 118°71′E)和克什克腾旗木希嘎河(43°73′N, 117°66′E), 黑龙江省哈尔滨市松花江(45°81′N, 126°47′E)和穆棱市乌苏里江支流穆棱河(44°92′N, 130°53′E), 共计4个采样点(图 1), 共119尾。根据分布水系, 董氏须鳅分为3个地理种群： 西辽河群体、松花江群体和穆棱河群体。采集工具包括定置地笼和拖网。在样本采集后, 依据形态特征, 参照Prokofiev[1]、Kottelat[2]和Cao等[10]对物种进行鉴定。新鲜样本用95%乙醇固定, 置于-20℃冰箱中保存备用, 保存于河北大学博物馆。

1.2 DNA提取及PCR扩增

取每尾样本左侧胸鳍剪碎, 使用组织基因组提取试剂盒DP304 (北京天根生化科技有限公司)进行提取。Cytb基因扩增和测序引物为L14724(3′-GACTTGAAAAACCACCGTTG-5′)和H15915(3′-CTCCGATCTCCGGATTACAAGAC-5′)[12]; D-loop基因为DL1(3′-ACCCCT GGCTCCGAAAGC-5′)和DH2(3′-ATCTTAGCATCTTCAGTG-5′)[13]。PCR反应条件为： 95℃预变性3min; 95℃变性30s,退火30s温度为56℃, 72℃延伸40s, 循环30次; 最后72℃延伸10min, 4℃停止反应; D-loop基因退火温度为51℃, 其他条件相同。PCR产物用1%琼脂糖凝胶电泳检测, 条带清晰明亮产物送往通用生物系统(安徽)有限公司完成序列测定。

1.3 数据处理和分析

将DNA序列用BioEdit 7.0.9软件进行人工校正比对并拼接。利用MEGA 5.05[14]软件分析序列中碱基组成、变异情况、计算转换/颠换比率, 基于Kimura-2参数模型计算群体间及群体内遗传距离。使用DNAsp 5.1软件[15]对所有样本的单倍型组成进行分析, 计算群体间核苷酸歧义度Dxy, 并分析核苷酸多样性(π)及单倍型多样性(h)。使用Arlequin 3.5软件[16]对数据进行中性检验及错配分析,种群分化指数(FST)以及τ值计算。使用MrBayes 3.1.1程序[17]和IQtree程序[18]分别构建贝叶斯(BI)和最大似然(ML)单倍型系统发育树, 以条鳅科高原鳅属的武威高原鳅(Triplophysa wuweiensis, KT2243 65)和怒江高原鳅(Triplophysa nujiangensis, KT213 598)作为外类群, 并以jModelTest 0.1.1软件[19]选择最适合的进化模型, 最佳模型为GTR+G。BI分析运用MrBayes 3.1.1程序, 4条马尔科夫链同时运行1×106代, 每100代对系统树进行抽样, 最终构建一致树, 通过后验概率评估每个节点。ML分析采用IQtree软件, 并通过具有1000个重复的bootstrap分析评估每个节点的可信度。使用Network软件[20],采用中接法(Median-joining)构建单倍型网络关系图。根据线粒体Cytb每百万年0.68%—0.84%的进化速率[21], 利用公式T=Dxy/2μ估算各类群间分歧时间[22], 其中μ为平均进化速率。利用公式τ=2ut来估算种群扩张时间[23],u为每条序列每个世代的变异速率, 根据u=2μk计算,μ为平均进化速率,k为序列长度。关于须鳅属鱼类性腺发育的研究显示, 须鳅一年仅在3月底至4月中旬产卵一次[24], 因此按1世代/年来计算。

2 结果

2.1 核酸序列与单倍型分析

图 1 采样点示意图Fig. 1 Map of the sample collection sites

线粒体Cytb序列总长度为1140 bp： 保守位点1074个, 约占序列总长度的94.2%, 变异位点66个,约占序列总长度的5.8%; 碱基组成T为28.9%, C为28.5%, A为24.6%, G为18.1%, 119个个体中共确认单倍型44个, 单倍型多样性(h)介于0.705—0.918, 核苷酸多样性(π)介于0.00138—0.00411。D-loop序列总长度896 bp： 保守位点867个, 约占序列总长度的96.8%, 变异位点29个, 约占序列总长度的3.2%; 碱基组成T为33.3%, C为18.7%, A为34.1%, G为13.8%, 119个个体中共确认单倍型20个, 单倍型多样性(h)介于0.353—0.844, 核苷酸多样性(π)介于0.00041—0.00324。相比较下, D-loop序列相对更为保守, 具有较少的变异位点, 但A+T含量明显高于Cytb序列。

序列联合后总长度为2036 bp, 共确认单倍型57个(表 1)。整体单倍型多样性(h)为0.960, 核苷酸多样性(π)为0.00496。其中, 西辽河群体具有最高的单倍型多样性(h=0.926,n=22)和核苷酸多样(π=0.00458,n=22), 穆棱河群体次之(h=0.923,n=77;π=0.00124,n=77), 松花江群体最低(h=0.805,n=20;π=0.00095,n=20)。

2.2 遗传分化

基于联合基因的AMOVA分析结果显示, 群体间变异占总变异的79.45%, 显著高于群体内变异的20.55% (表 2)。3个地理群体间, 遗传分化系数最大值为松花江和穆棱河群体间, 遗传变异固定指数(FST)为0.876 (表 3); 西辽河和松花江群体间以及西辽河和穆棱河群体间的遗传分化系数FST分别为0.716和0.721, 各群体间遗传分化系数差异均显著(P＜0.01)。3个地理群体间, 穆棱河和松花江群体间遗传距离最大, 为1.0% (表 3), 西辽河和穆棱河群体间最小, 为0.8%; 各群体内遗传距离0.2%—0.5%,其中西辽河群体内遗传距离最大, 为0.5%。

2.3 分子系统发育关系分析

基于联合基因, 以BI法和ML法分别构建的单倍型系统发育树拓扑结构基本一致(图 2)： 董氏须鳅以高置信度(BI=1.00, ML=100)聚为一支, 可分为3个谱系： 西辽河群体、穆棱河群体和松花江群体。在这3个谱系中, 松花江群体位于系统发育树的最外部, 其中单倍型Hap22与群体其他单倍型并未聚为一支, 而是形成并系群; 西辽河群体和穆棱河群以较高置信度成姐妹群(BI=0.76, ML=99)。以Cytb基因为参考, 根据Dxy估算分支分歧时间： 西辽河群体和穆棱河群体的分歧时间约为0.872—0.669 Ma前, 松花江群体与其他2群体的分歧时间约为1.082—0.876 Ma前(表 4), 群体分歧时间处于更新世中期。

表 1 董氏须鳅线粒体基因遗传多样性、中性检验及错配分析Tab. 1 Genetic diversity, neutrality tests and mismatch distribution based on mtDNA of B. toni

表 2 基于线粒体Cyt b和D-loop基因联合计算群体内和群体间变异的分子方差分析Tab. 2 Analysis of molecular variance (AMOVA) among population and within population based on combined Cyt b and D-loop gene sequences

表 3 基于线粒体Cyt b和D-loop基因联合计算群体间的遗传分化系数FST (对角线下方)、群体间的遗传距离(对角线上方)及群体内遗传距离(对角线上)Tab. 3 Population pairwise FST (below diagonal), genetic distances among populations (above diagonal) and genetic distances within populations (on diagonal) based on combined Cyt b and D-loop gene sequences

图 2 基于线粒体Cyt b和D-loop基因联合构建的董氏须鳅贝叶斯(BI)单倍型系统发育树(节点上方为BI后验概率, 下方为ML支持率)Fig. 2 Haplotype bayesian tree (BI) based on combined Cyt b and D-loop gene sequences of B. toni (posterior probability for Bayesian inferences above branch and bootstrap support for ML below the branch)

2.4 单倍型网络分析

单倍型进化网络图与单倍型系统发育树结果一致(图 3), 形成3个明显的地理聚群。其中, 地理分布较近的松花江群体与穆棱河群体之间具有较大的遗传距离。在单倍型网络图中, 来自松花江的Hap16和来自穆棱河的Hap27位于各自星状网络结构类群的中央位置, 可能为各自地理群体的原始倍型。而西辽河的单倍型网络结构中部分缺失中央连接单倍型, 且并未表现出明显的中央单倍型。

2.5 种群动态分析

西辽河流域群体Tajima’sD和Fu’sFs值均呈负数且不显著(表 1), 表明群体未偏离原来假设, 符合中性进化假说。错配分析SSD和Hri均表明未偏离原来假设, 即符合群体扩张假说, 但其核苷酸错配分布图呈现多峰(图 4A)。穆棱河和松花江群体Tajima’sD和Fu’sFs值均呈负数且已经达到显著水平, 错配分析SSD和Hri均不显著, 表明这2个群体都可能经历过群体扩张, 且核苷酸错配分布图呈单峰模式(图 4B、4C)。

以鳅科鱼类Cytb基因进化率为参考, 根据核苷酸错配分布的τ值, 估算出松花江群体和穆棱河群体的体扩张时间(表 4)： 松花江群体的扩张时间为0.071—0.057 Ma前, 而穆棱河群体的扩张时间为0.058—0.047 Ma前。

3 讨论

3.1 遗传多样性及群体结构

王雪等[6]基于Cytb基因对赤水河两种荷马条鳅属鱼类的研究表明荷马条鳅属鱼类遗传多样性较为丰富, 单倍型多样性(h)介于0.821—0.936, 核苷酸多样性介于0.00493—0.01105; 杨成等[7]基于Cytb基因对黄河上游拟硬刺高原鳅的研究显示各地理群体单倍型多样性(h)介于0.650—0.890, 核苷酸多样性介于0.00100—0.00250。我们基于Cytb基因对董氏须鳅的3个群体研究显示, 董氏须鳅具有较高的单倍型多样性(h=0.705—0.918)和较低的核苷酸多样性(π=0.00138—0.00411), 这与其他鳅超科鱼类相同。其中, 松花江群体的单倍型多样性和核苷酸多样性在3个群体中均表现出最低值(h=0.705,π=0.00138), 总体上遗传多样性处于较低水平; 穆棱河群体的单倍型多样性处于中等水平, 但核苷酸多样性处于低水平(h=0.858,π=0.00159); 西辽河群体的单倍型多样性和核苷酸多样性均明显高于另外2个群体(h=0.918,π=0.00411)。已有研究表明, 庞大的种群数量及快速增多的生活习性是维持自然群体单倍型多样性高的基础[25]。我们的研究结果显示, 穆棱河群体及西辽河群体与松花江群体相比,存在较高的单倍型多样性, 推断可能是因为穆棱河群体及西辽河群体数量较大, 目前正处于群体成长阶段。而松花江群体具较低的单倍型多样性可能是因为松花江群体经受了环境胁迫, 我们采样调查了解到, 由于污染以及捕捞等原因, 须鳅的种群数量较过去有大幅度的减少, 很多河段及支流已不见须鳅的分布。因此致使松花江群体有效群体数量降低, 造成单倍型多样性降低。同时, 我们认为不能够排除由于样本量不充分而导致结果出现偏差。

AMOVA分析结果显示, 3个地理群体的变异主要来源是群体间变异(79.45%), 3个群体间的间遗传分化系数远大于0.25(FST=0.71564—0.87618), 表明3个群体间已呈现高度分化[26,27]。系统发育树和单倍型网格图也表明3个群体间呈高度分化, 较远的地理距离等地理屏障能够阻碍群体间的基因交流, 并且由于它们的小生境适应性和迁移能力差的原因, 逐渐加大群体间的变异, 形成具有遗传差异的不同种群。遗传距离分析显示, 3个地理群体间遗传距离为0.8%—1.0%, 群体内遗传距离为0.2%—0.5%。Avise[26]提出以线粒体基因作为遗传标记时,种内遗传距离往往小于1.0%, 极少高于2.0%。陈咏霞等[9]对北方地区花斑鳅3个亚种的研究表明, 河北地区花斑鳅与俄罗斯3个亚种间的遗传距离为1.7%—2.6%, 各亚种种内遗传距离为0.1%—0.6%。因此, 从遗传距离来看, 董氏须鳅3个地理种群间的遗传分化应属于种内分化。从形态特征来看, 董氏须鳅3个地理种群在须鳅属的鱼类鉴别特征上[10,28],如鳞片、鼻孔位置、下唇侧叶、吻长等特征上无明显差异, 可见董氏须鳅3个地理种群尽管分化程度较高, 但并没有达到种间分化的程度, 属于种内分化。值得注意的是, 在这3个地理群体中, 地理距离较近的松花江群体和穆棱河群体间的遗传分化系数及遗传距离最大, 表明这2个群体间的亲缘关系更远。系统发育树和单倍型网络分析同样也暗示着松花江和穆棱河2个群体相比西辽河流域群体,有着更大的群体间差异, 这可能与董氏须鳅的扩散和迁移路径、冰期所带来的影响以及栖息环境差异等有着密切关系, 使这2个群体间累积了更多的群体间遗传变异。

表 4 基于Cyt b序列的3个地理群体分歧时间及扩张时间Tab. 4 Divergent time and population expansion time of three populations based on Cyt b gene sequences

图 3 基于线粒体DNA Cyt b基因和D-loop基因联合分析的单倍型进化网络图Fig. 3 Haplotype network based on combined Cyt b and D-loop gene sequences of mtDNA

图 4 基于线粒体DNA Cyt b基因和D-loop基因联合计算不同地理群体的董氏须鳅核苷酸错配分布图Fig. 4 The mismatch distribution based on combined Cyt b and D-loop gene sequences for different populations of B. toniA. 西辽河群体(XLH); B. 松花江群体(SHJ); C. 穆棱河群体(MLH)A. population of West Liao River (XLH); B. population of Songhua River (SHJ); C. population of Muling River (MLH)

3.2 群体历史动态分析

3个群体的中性检验Tajima’sD和Fu’sFs均呈负值, 可用于推断种群扩张[29]： 其中西辽河群体中性检验结果不显著(P＞0.05), 说明未检测到群体扩张迹象, 可能未发生过群体扩张事件; 松花江和穆棱河群体中性检验结果显著(P＜0.01), 表明这2个群体可能经历过群体快速扩张事件。3个群体错配分析Hri和SSD检验均不显著, 均未偏离群体扩张模型。通过绘制歧点分布图表明(图 4), 西辽河群体错配分布图呈多峰模式, 暗示没有发生过群体扩张事件; 松花江和穆棱河群体核苷酸错配分布图明显呈单峰模式, 表明这2个群体近期都经历了快速的群体扩张。

根据Cytb基因作为分子钟推算种群分歧时间和扩散时间, 3个群体的分歧时间在1.082—0.669 Ma前, 受到第四纪更新世早期的几次冰期循环的影响[30]。李郊[31]对黑龙江地区林蛙的生物地理学研究表明, 我国的松花江中游以及乌苏里江流域在早更新世晚期时存在生物避难所, 这为松花江群体和穆棱河流域的董氏须鳅提供了存活下来的条件。同样, 我们认为西辽河流域也可能存在这样一个避难所。考虑到董氏须鳅迁移能力差, 其种群分化极有可能是随水位变化而进行迁移。因此, 我们认为董氏须鳅在冰期期间, 冰川融水使得河流水位上升形成临时的水系间的联系, 通过这种临时的水系间的连接, 迁移至不同避难所, 经过长时间逐渐累积差异, 最终形成具有较大遗传差异的不同种群。根据Cytb计算出的核苷酸错配分布τ值估算,松花江和穆棱河群体近期发生扩张时间在0.071—0.047 Ma前, 处于更新世晚期的末次冰期, 有研究表明该时期东北地区存在未被冰川完全封闭的水体[30], 而在此时期东北地区的董氏须鳅发生了快速群体扩张事件可能暗示了末次冰期内该地区发生的区域气候变化。