多基因解析长江与珠江鳤的遗传结构

陈蔚涛，段辛斌，高 雷，李新辉，杨计平，汪登强

1. 中国水产科学研究院珠江水产研究所/农业农村部珠江中下游渔业资源环境科学观测实验站，广东 广州 510380

2. 中国水产科学研究院长江水产研究所，湖北 武汉 430223

鳤 (Ochetobius elongatus) 隶属于鲤形目、鲤科、鳤属，曾是广泛分布于我国南方长江、珠江等诸多水系的重要经济鱼类[1-3]。受到涉水工程、江河环境污染、栖息地破坏、过度捕捞等因素的影响，其资源量萎缩极其严重，在很多河流江段已多年未被监测到[4-5]。2016年修订的《中国脊椎动物红色名录》已将鳤的保护等级上升至极度濒危[6]；近期颁布的《长江保护法》还特别针对鳤提出了重要指示，因此亟需给予重视和保护。然而，由于鳤样本较难获得，目前其相关研究鲜见报道，零星研究仅集中在大阶元的系统发育进化关系[7-8]、小样本的遗传变异分析[5]和局域分布鳤种群的遗传多样性[9]方面，这极大制约了有效保护策略的制定。

珠江和长江是我国的两大重要水系，在鱼类生物多样性和种质资源孕育方面起着重要作用[3,10]。南岭山脉作为珠江和长江两大重要水系的分水岭，是阻隔两个水系鱼类迁移扩散的天然屏障。已有大量遗传进化研究表明，南岭山脉导致珠江和长江两个水系共同分布的鱼类发生了遗传分化甚至形成高度分化的遗传谱系[11-13]。在此背景下，长江鳤种群与珠江鳤种群是否形成显著的遗传分化，遗传分化的程度如何，若存在遗传分化是在什么时间节点开始的，驱动分化的环境因素是什么，这些关键问题均尚未找到答案。

利用遗传进化方法揭示物种的遗传结构对于物种的管理保护具有重要意义[14]。本文使用珠江和长江两个水系鳤群体的2个线粒体基因 [细胞色素b (Cytb) 和酰胺腺嘌呤二核苷酸 (NADH) 脱氧酶亚单位2 (ND2)] 和2个核基因 [肌球蛋白重链6(myh6) 和重组激活基因2 (RAG2)] 的序列，利用系统发育分析、单倍型网状图、分化时间估算等方法分析了两个水系鳤种群的遗传结构，估算了分化时间，探讨了导致其分化的可能环境因素，以期为鳤的管理和保护提供科学依据。

1 材料与方法

1.1 样品采集

长江的7尾鳤样本 (体长21.0～24.2 cm，体质量74.5～112.0 g) 由中国水产科学研究院长江水产研究所于2020年在长江湖北公安江段常规渔业调查采集获得 (图1和表1)。剪取少量胸鳍妥善保存于无水乙醇中，用以基因组DNA的提取。处理后的鳤样本带回中国水产科学研究院长江水产研究所，用甲醛固定制作成标本用以科研和科普。另外，在珠江广东肇庆江段采集鳡 (Elopichthys bambusa)1尾，作为系统发育分析的外类群。

表1 采样信息和GenBank序列号Table 1 Sampling information and GenBank No.

图1 采样示意图Fig. 1 Map of sampling sites

1.2 基因组DNA提取、扩增与测序

采用高盐法提取基因组DNA，并利用质量分数为1.2%琼脂糖凝胶电泳检测DNA质量。本文分别选取了2个线粒体基因 (Cytb和ND2) 和2个核基因 (myh6和RAG2)。PCR反应体积为30 μL，其中 Taq PCR Master Mix (上海生工) 15 μL，包括Taq DNA 聚合酶、dNTP、PCR Buffer、PCR stabilizers、Gel loading和核酸染液，模板 DNA 1 μL，正反向引物各 1 μL (浓度为 10 μmol·L−1)，加灭菌超纯水至总体积30 μL。鳤和鳡4个基因扩增与测序的引物信息和PCR反应体系参考杨计平等[9]。PCR产物经质量分数为1.2%琼脂糖凝胶电泳检测，对于扩增效果良好的产物交由测序公司双向测序。

1.3 数据处理与分析

测序所得序列使用Lasergene软件包 (DNASTAR Inc., Madison, WI) 来检测测序质量并生成一致序列。所得的一致序列使用MUSCLE[14]进行比对，然后剪掉首尾两端噪音序列至同样长度。在Gen-Bank核苷酸序列数据库下载了珠江52尾鳤样本已公布的4个基因序列，具体采集信息参考杨计平等[9](图1和表1)。这些公布序列将与新获得序列比对后一起用于后续分析。为了获得更多的变异信息，分别合并了2个线粒体基因和2个核基因用以后续分析。线粒体DNA单倍型和核基因等位基因均使用DnaSP 6.12[15]进行划分。利用DnaSP 6.12统计单倍型数、计算单倍型多样性和核苷酸多样性，并估算了两个水系鳤群体的遗传分化状况。使用Network 10.2 (http://www.fluxus-engineering.com)构建了Cytb+ND2的单倍型中介网状图[16]和myh6+RAG2的等位基因中介网状图[16]。

为了揭示珠江和长江鳤种群的系统发育关系，本文通过MrBayes 3.2.7[17]构建贝叶斯树，并以鳡作为外类群。合并了4个基因进行系统发育树的构建。贝叶斯树使用马尔科夫链蒙特卡洛方法(Markov chain Monte Carlo, MCMC) 分析运行108万代，每1 000代取样1次，并舍弃初始的25×105代。另外，本文使用同样的参数设置对Cytb+ND2和myh6+RAG2也分别进行了贝叶斯树的构建。不同基因组合的最优核苷酸替代模型分别利用MrModeltest 2.3[18]计算获得 (表2)。最后，基于Kimura双参数模型 (Kimura-2-parameter,K2P)[19]，利用MEGA 7.0[20]计算两个水系群体之间的遗传距离。

表2 系统发育树构建和分化时间估算的最优模型Table 2 Optimal model of phylogenetic tree construction and differentiation time estimation

本文联合4个基因，使用BEAST v.1.8.1[21]软件中的贝叶斯马尔可夫链蒙特卡洛方法估算两个水系鳤群体的分化时间，采用严格分子钟模型、Yule先验模型和GTR+I+G最优核苷酸替代模型。由于缺乏有效的化石记录，本文采用固定进化速率来估算分化时间。联合基因的进化速率由Cytb的进化速率估算获得。由于联合基因的K2P平均遗传距离 (0.29%) 与Cytb的K2P平均遗传距离 (0.43%) 的比值为0.674，因此利用Cytb的进化速率乘以比值便可获得联合基因的进化速率。1%～2%是鲤科鱼类Cytb基因普遍接受的一个进化速率范围[22-24]，因此联合基因的进化速率为0.674%～1.348%。分别利用0.674%和1.348%的进化速率进行运算，每次运算共运行1×109代，每1 000代进行1次抽样，舍弃初始5×108代。运行结束后，利用Tracer 1.6[25]评价每个参数的稳定性状况。最后使用TreeAnnotator v.1.8.1[21]总结出一致树，在Figtree v.1.3.1 (http://tree.bio.ed.ac.uk/software/figtree) 显示节点的分化时间和后验概率。

2 结果

2.1 序列信息

本文分析数据包括新获得的长江7尾鳤的全部4个基因片段以及杨计平等[9]52尾珠江样本的对应4个基因片段，共计59尾样本。序列比对后共获得895 bp的Cytb序列，包含了30个变异位点；935 bp的ND2序列，包含了26个变异位点；763 bp的myh6序列，包含了3个变异位点；857 bp的RAG2序列，包含了5个变异位点。在4个基因片段中均未发现缺失和插入。此外，成功获得了鳡的4个基因序列，GenBank序列号分别为MW657589、MW657597、MW657605和 MW657613。

2.2 遗传多样性、系统发育分析与遗传结构

基于线粒体联合基因和核基因联合基因估算的单倍型数目、单倍型多样性和核苷酸多样性见表3。基于线粒体基因估算的两个水系鳤的单倍型多样性为0.950，核苷酸多样性为0.494%。基于核基因计算的单倍型多样性和核苷酸多样性明显低于线粒体联合基因，这主要是与核基因较慢的进化速率有关。基于线粒体基因组合发现两个水系之间的遗传分化系数 ФST为0.863 (P<0.001)，基于核基因组合获得两个水系之间的遗传分化系数FST为0.518(P<0.001)。

表3 鳤种群的单倍型与遗传多样性指数Table 3 Haplotypes and genetic diversity indexes of O. elongatus populations

线粒体基因组合共界定了29个单倍型，其中珠江群体26个 (H1—H26)，长江群体3个 (H27—H29，图2-a)。对于核基因组合，共鉴定了13个等位基因型 (A1—A13)，其中珠江和长江群体分别界定了5和8个等位基因型 (图2-b)。基于线粒体基因组合的单倍型中介网状图表明，珠江和长江的鳤群体形成显著的遗传分支结构，并且两个分支之间相差22个突变位点 (图2-a)。基于核基因组合的等位基因中介网状图虽然未形成明显的遗传分支结构，但两个水系群体之间却各自形成了特有的等位基因型 (图 2-b)。

图2 单倍型/等位基因中介网状图注：圆的大小代表单倍型频率，颜色代表所属群体，空心圆代表未检测到的单倍型。Fig. 2 Haplotype/allele median joining networkNote: The size of the circle represents the haplotype frequency; the color represents the group, and the hollow circle represents the undetected haplotype.

基于4个基因组合和线粒体基因组合获得的系统发育树拓扑结构一致，均支持珠江和长江形成了显著的支系结构 (图3-a和图3-c)。基于核基因组合构建的贝叶斯树并未形成可信的支系结构 (图3-b)。基于4个基因组合、线粒体基因组合和核基因组合估算的两个水系群体间的遗传距离分别为0.87%、1.52%和0.14%。

图3 基于不同基因组合的贝叶斯系统发育树Fig. 3 Bayesian phylogenetic trees based on different gene combinations

2.3 分化时间估算

Beast运算的有效取样大小 (Effect sample size,ESS) 参数均大于200，表明运算是收敛的。基于不同进化速率估算的分化时间显示，珠江和长江的鳤群体形成了2个进化支系，并获得100%的后验概率，分化时间大概在 0.38～0.76 Ma (百万年前) (图4)。

图4 基于线粒体基因组合的分化时间结果Fig. 4 Estimation of divergence time based on combination of two mitochondrial genes

3 讨论

3.1 谱系分化

基于线粒体基因组合的贝叶斯系统发育树、单倍型网状图和基于4个联合基因的系统发育树均发现，两个水系的鳤群体存在高度分化的遗传谱系，且形成了严格的珠江-长江分化格局。此外，两个水系群体之间的遗传距离在线粒体基因层面上也呈现较高水平 (1.52%)。这与范凤娟和章群[5]的研究结果相类似，其获得了4尾长江样本和2尾珠江样本，利用线粒体Cytb基因发现两个水系样本各自形成了两个独立遗传谱系，且谱系间的遗传距离为1.5%。这些研究结果表明两个水系鳤群体经历了长期的地理隔离。作为珠江和长江水系的分水岭，南岭山脉可能是导致其分化的重要推动力。类似的例子也发生在珠江和长江共存的其他鲤科鱼类上，比如䱗 (Hemiculter leucisculus)[26]、马口鱼(Opsariichthys bidens)[27]。基于核基因组合获得的系统发育树和等位基因网状图发现，两个水系鳤群体虽未形成独立的分支结构，却形成了特有的等位基因型，说明两个水系群体之间不存在基因交流现象。然而，由于本文中涉及的样本量较少，后期需要使用更多样本量来验证这个观点。

草鱼 (Ctenopharyngodon idella)、青鱼 (Mylophayngodon piceus)、鲢 (Hypophthalmichthys molitrix) 和鳙 (H. nobilis) 是我国重要的经济鱼类，在珠江和长江均有广泛分布[1-3]。基于草鱼[28]、青鱼[29]、鲢[30]和鳙[30]的相关遗传分化研究发现，虽然珠江和长江野生群体形成了一定程度的群体分化，但均未形成独立的遗传谱系。本文发现两个水系的鳤群体形成了两个深度分化谱系，表明鳤的分化程度显著高于两个水系共同分布的这几种重要经济鱼类。鳤有限的分布范围和较小的种群数量可能是造成这种差异遗传模式的原因[1,3]。早期的地质学和鱼类区系研究认为，珠江和长江水系的部分河流相隔距离很近，甚至有少数河流直接连通[31-33]。四大家鱼分布广泛、种群数量庞大，促使两个水系群体之间可以通过一些相邻的河流、江段进行迁移扩散，降低了群间的遗传分化。另外，四大家鱼人工繁育种群数量庞大，不可避免的养殖逃逸和人工增殖放流也会导致两个水系四大家鱼的遗传种质趋于同质。相比之下，鳤主要分布在大型河流江段与湖泊，种群规模极其有限，加上人为干预极少，因此两个水系群体之间进行相互扩散的概率较小。

已有大量研究表明古地质历史事件是促进物种分化的重要驱动因子[26-27,34-35]。本文发现珠江和长江鳤群体的谱系分化大概发生在更新世中期。这一时期刚好与上新世后期开始的青藏高原快速隆升阶段第二幕的昆仑-黄河运动时间段吻合 (1.1～0.6 Ma)。研究表明，青藏高原在上新世后期进行快速抬升，共分为3个阶段[36-37]：第一阶段是青藏运动，时间介于3.6～1.7 Ma；第二阶段是昆仑-黄河运动，时间介于1.1～0.6 Ma；第三阶段是共和运动，时间大致在0.15 Ma。青藏高原的快速隆升重塑了东亚地区的景观特征，包括河流和山脉。早期研究也证明了青藏高原的快速抬升导致了一些鱼类的分化甚至成种[26,38-40]。因此，两个水系鳤的谱系分化很可能与青藏高原的快速隆升阶段相关。

3.2 保护建议

针对特定物种展开遗传结构解析能够有效辅助物种的管理和保护[26,41-42]。鳤的资源极其缺乏，近几年关于鳤的零星发现都能引起社会的广泛关注。由于鳤的野外资源稀少以及人工繁殖尚未成功，细化保护单元并展开就地保护是一个重要举措。因此，通过遗传进化分析解析鳤在不同水系的遗传结构对于细化其保护单元非常必要。针对濒危或稀有物种的保育，许多学者认为演化显著单位 (Evolutionary Significant Unit, ESU) 为适当的保护管理单元[43]。ESU是指一个种群在线粒体基因或者核基因水平已经形成独立的单系群 (Monophyly)[44]。本文发现珠江和长江鳤群体形成了独立的线粒体基因遗传谱系和特有的等位基因型，因此可以看作两个不同的ESU，在保护策略制定和实施上应该区别对待。然而，由于本研究样本量和分子标记有限，今后需要通过扩大研究范围和样本量，结合更多的标记 (如核基因组数据) 来全面评估鳤的种群遗传结构，以期为鳤的保育策略的制定提供更系统的支撑。