长江上游中华纹胸鮡遗传多样性及遗传结构研究

张力文，田辉伍，陈大庆，刘绍平，段辛斌，汪登强

(1.上海海洋大学水产科学国家级实验教学示范中心，上海 201306；2.中国水产科学研究院长江水产研究所，农业农村部长江中上游渔业资源环境科学观测实验站，武汉 430223)

纹胸鮡属(Glyptothorax)起源于喜马拉雅山及横断山脉中段，随着喜马拉雅山运动，沿着江河向各处扩散和迁移[1]，是鮡科鱼类中种类最多、分布最广的一个属。根据《中国动物志》和《四川鱼类志》的记录，在我国境内分布的纹胸鮡属共有22种[2]，主要分布在伊洛瓦底江、怒江、澜沧江、雅鲁藏布江、元江和长江及其附属水系。中华纹胸鮡(Glyptothoraxsinensis)和福建纹胸鮡(G.fukiensis)是该属分布在低海拔的鱼类，主要分布在长江及以南河流。谢仲桂等[3]和王汨等[4]从形态特征角度对这两种鱼进行对比后，认为二者并没有明显的性状分化，应为同种异名，因此，本研究将这两种鱼类统称为中华纹胸鮡。

中华纹胸鮡为小型底栖鱼类，喜栖息于流水环境中，产粘性卵，受精卵粘附在石头上发育孵化。中华纹胸鮡分布较为广泛，栖息在多种类型生态环境，既可生活在长江上游高山峡谷河流中，也可在长江中下游平原湖泊及河流中，形成了不同的生态群[5]，是研究生物地理学和环境适应性的良好对象。长江上游渔业资源调查显示，中华纹胸鮡是上游干支流主要渔获物之一[6，11，12]。目前长江上游建设或规划101座水电站[7]，导致中华纹胸鮡等鱼类栖息地缩小和破碎化，阻碍其基因交流，引起资源量下降[8]。研究中华纹胸鮡群体遗传多样性和遗传结构，对评估水电开发对鱼类的影响和制定资源保护措施具有重要意义。

目前有关中华纹胸鮡的研究主要集中在物种分类[9]、资源分布[10]和系统发育[11，20]方面，未见群体遗传学研究报道。本研究采用Cytb和COI基因对长江上游中华纹胸鮡群体进行遗传多样性分析，旨在阐明其遗传背景、不同地理群体的遗传结构，为资源保护提供科学依据。

1 材料与方法

1.1 样本采集

样本采集于2017年至2018年，共9个采样点341尾样本。采样点包括长江上游干流5个[巴南(BN)、江津(JJ)、合江(HJ)、泸州(LZ)、南溪(NX)]、岷江1个[犍为(QW)]、金沙江中游1个[攀枝花(PZH)]、金沙江下游1个[巧家(QJ)]及金沙江支流黑水河1个[宁南(NN)]。采样点和样本量信息见图1和表1。样本经鉴定[5]后剪取少量鳍条，保存于无水乙醇中备用。

图1 中华纹胸鮡采样位点Fig.1 Sampling sites of G.sinensis

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

样本经双蒸水冲洗后浸泡2 h以上除去乙醇，采用盐析法提取基因组DNA[12]。用线粒体DNACytb和COI基因序列进行分析，其中Cytb基因扩增引物可参考文献[13]，COI基因扩增引物参考文献[14]。两个基因采用相同的PCR体系和反应程序，PCR扩增体系为5 0 μL，含2×Master Mix酶(武汉天一辉远生物科技有限公司)25 μL，10 μmol/L引物各2 μL，DNA模板2 μL，灭菌去离子水补至50 μL。PCR扩增反应程序如下：94 ℃预变性3 min；94 ℃变性30 s，55 ℃退火30 s，72 ℃延伸1 min，35个循环；72 ℃再延伸8 min。扩增产物用1%琼脂糖凝胶电泳检测，将扩增效果良好的扩增产物送至武汉天一辉远生物科技有限公司进行双向测序，测序引物与PCR引物相同。

1.3 数据处理和分析

使用DNASTAR软件包中的SeqMan对测序结果进行拼接，序列比对由CLUSTL X完成。使用MEGA v6.06[15]将序列进行翻译，以检查是否存在测序错误，并计算碱基组成和转换/颠换比。群体遗传多样性参数，如变异位点、单倍型数量、单倍型(Hd)和核苷酸(Pi)多样性使用软件Dnasp version 5.10[16]分析。应用Arlequin version 3.1.1[17]进行分子方差分析(AMOVA)估算群体间遗传变异及遗传分化指数(FST)，分析群体遗传结构。群体间基因流用公式计算。为了检验群体间地理距离与遗传分化是否符合距离隔离(Isolation by distance，IBD)[18]模式，用SPSS 22.0计算距离和遗传分化指数的相关性。

单倍型网络结构用NETWORK 5.0[19]完成。采用MEGA 6.0构建基于Kimura-two-Parameter(K2P)单倍型邻接系统发育树(NJ树)，boostrap设为1000，检验分支支持率。为了确定本研究所采集样本中是否存在隐存种或亚种，采用目前常被用作物种鉴定的DNA条形码的COI基因进行Automated barcode gap detection(ABGD)分析，由在线软件服务器ABGD web(https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html)完成，以P=0.01作为种数参考指数[20]。采用Structure 2.3.4 软件分析群体的遗传结构，设置运行参数K=2-9，每个K值10次重复，将运算结果提交到在线工具Structure Harvester(http://taylor0.biology.ucla.edu/struct_harvest/)[21]计算△K，以△K的峰值，判断群体的最佳K值。

为评估群体历史动态，首先用Dnasp完成中性检验(Tajima′s D,Fu′s Fs)和碱基错配分析。其次用BEAST version 2.4.3[22]进行Bayesian skyline plot(BSP)分析，推断中华纹胸鮡的群体动态。以上三种方式仅使用Cytb序列，BSP分析采用鱼类线粒体Cytb的平均突变率(0.65%每百万年)[23]，保证ESS值大于200，结果的图形化用软件Tracer version 1.5绘制。

2 结果与分析

2.1 遗传多样性

序列对比后得到341条长度为1 057 bp的Cytb序列和109条长度为607 bp的COI序列。序列中的转换数明显高于颠换数，Ts/Tv比分别为6.328和14.331。Cytb序列和COI序列A+T含量分别为58.02%和54.3%，高于C+G含量(41.95%和45.7%)，与硬骨鱼类一致[24]。

Cytb序列共检测到80个变异位点，其中简约信息位点52个，单一变异位点28个。341条序列定义了65个单倍型，其中单倍型HB_2频率最高，为35.78%，其次是HB_19为12.32%，HB_9为8.21%。其中单倍型HB_2和HB_9为长江干流和岷江6个群体共享。HB_19为QJ和NN 2个群体共享。总样本Hd和Pi分别为0.845和0.007(表1)。各群体中，BN群体Hd最高为0.867，PZH群体Pi最高为0.005，NN群体最低，分别为0.251和0.002。

COI序列共检测到31个变异位点，其中13个简约信息位点，18个单一变异位点。109条序列共定义了23个单倍型，其中单倍型HC_1和HC_4频率最高，分别为38.53%和25.69%，为LZ、NX和QW群体共享。单倍型HC_7频率为9.17%，仅出现在NN群体。总样本Hd和Pi分别为0.774和0.005(表1)。两种标记得到的遗传多样性指数相近，金沙江的3个群体遗传多样性低于长江上游群体。

表1 中华纹胸鮡线粒体DNA群体遗传多样性和中性检验Tab.1 Genetic diversity parameters and neutral test of G.sinensis populations based on mtDNA

2.1.2 遗传结构

分子变异方差分析(AMOVA)显示，无论是Cytb序列还是COI序列，遗传变异主要来自群体间(Cytb，59.62%；COI，55.24%)，群体间出现显著遗传分化(Cytb，FST=0.60，P<0.01；COI，FST=0.55，P<0.01)(表2)。

基于Cytb序列的两两群体间FST分析结果见表3，长江干流5个群体和岷江的犍为群体间FST均小于0.05，这6个群体与QJ、NN和PZH群体之间的FST都大于0.7。QJ和NN群体间FST也小于0.05，它们与PZH群体间的FST为0.459和0.644。基因流分析显示长江干流和岷江群体间Nm值(或绝对值)均大于8，它们与QJ、NN和PZH群体间对群体间FST(表3)与地理距离(表4)进行Pearson相关性分析表明，这两个变量间呈线性正相关关系(r=0.696，P<0.001)。

表3 基于Cyt b的中华纹胸鮡群体间两两FST(对角线下方)和基因流 Nm(对角线上方)Tab.3 Pairwise FST(below diagonal)and Nm values(above diagonal)among populations of G.sinensis based on Cyt b sequences

表4 中华纹胸鮡群体间的地理距离Tab.4 Pairwise geographical distances among G.sinensis populations km

利用Median Joining方法构建Cytb单倍型网络结构图(图2)显示，中华纹胸鮡群体形成三个明显分支：Clade I、Clade II和Clade III，分支间突变步骤为12-22个，分支内相邻单倍型间突变步骤最多为5个(图2 a)。Clade I个体主要来自于BN、JJ、HJ、LZ、NX、QW等群体；Clade II主要来自NN和QJ等群体；Clade III主要来自PZH等群体(图2 a)。COI单倍型网络结构图形成两个分支：Clade A、Clade B，Clade A主要来自LZ、NX、QW等群体，与Cytb的Clade I对应；Clade B主要来自宁南等群体，与Cytb的Clade II对应(图2 b)。

以老挝纹胸鮡(G.laosensis)和扎那纹胸鮡(G.zanaensis)(登录号：HM636516.1,HM636515.1，HQ898002.1,HQ593565.1,DQ514349.1,KM610 755.1)为外类群，构建NJ系统发育树见图2。系统发育树具有相似的拓扑结构，其中Cytb系统发育树呈现三个具有高支持率的分支，COI树有两个高支持率分支，分别与单倍型网络结构图的分支对应。基于遗传距离的ABGD分析结果表明，P=0.01时，样本被划分为一组。

对全部样本进行Structure聚类分析，根据△K结果显示K=2时最佳。当K=2时，长江上游6个群体和金沙江3个群体被分为2个不同群体(图3a)；当K=3或4时(图3 b、c)，样本被划分为3个聚类群，长江上游、QJ和NN、PZH群体各为一个聚类群，与网络结构图(图2)和系统发育树(图2)相符合。

图3 基于Cyt b序列的中华纹胸鮡群体Structure聚类分析Fig.3 Structure clustering conducted based on Cyt b within populations of G.sinensis

图2 基于Cyt b序列(a)和COI序列(b)构建的中华纹胸鮡单倍型网络结构图(右)和邻接系统发育树(左)Fig.2 Haplotype network(right)and NJ phylogenetic tree(left)of G.sinensis based on(a)Cyt b and(b)COI sequences单倍型网络结构图中一个圆圈代表一个单倍型，圆圈大小代表样本频率.

2.1.3 种群历史

对Cytb序列的9个群体以及Clade I、II(图4)分别进行Tajima′s D[25]和Fu′s Fs[26]中性检验，结果显示Clade I(-59.18，-2.18)均为显著负值，长江上游群体也均为负值(表1)。碱基错配分析显示，Clade I、Clade II呈多峰结构(图4)。BSP分析结果表明，Clade I群体在距今0.1～0.04 Ma(百万年)经历了扩张事件。

图4 基于Cyt b序列的中华纹胸鮡错配分布分析和BSP分析Fig.4 Graphs of mismatch distribution and extended BSP for G.sinensis based on Cyt b

3 讨论

3.1 遗传多样性

目前已开展群体遗传学研究的纹胸鮡属鱼类有怒江扎那纹胸鮡(Hd=0.851,Pi=0.014)[27]和澜沧江老挝纹胸鮡(Hd=0.299,Pi=0.000 32)[24]，与这两种鱼比较，中华纹胸鮡总体遗传多样性(表1)与扎那纹胸鮡相似，高于老挝纹胸鮡。但从不同群体看，各群体遗传多样性差异较大，其中金沙江的3个群体(PZH、NN、QJ)遗传多样性明显低于长江上游及岷江的6个群体(表1)。金沙江中华纹胸鮡群体低水平的遗传多样性可能是由于奠基者效应造成。受第四纪冰期影响，冰期期间中华纹胸鮡可能仅存于长江中下游，随着冰期结束，金沙江各江段气候随海拔从低到高逐渐回暖，中华纹胸鮡从长江中下游逐渐往上迁移而形成目前的格局。长江和金沙江的泉水鱼(Semilabeoprochilus)[28]和裸体异鳔鳅鮀(Xenophysogobionudicorpa)[29]也有类似的情况。

3.2 群体遗传结构

AMOVA分析、network网络结构、系统发育树以及Structure聚类分析的结果都表明，长江上游及金沙江中华纹胸鮡群体存在显著遗传分化，至少形成3个遗传谱系，分别为金沙江中游(PZH)、金沙江下游(QJ、NN)、长江上游(BN、JJ、HJ、LZ、NX、QW)群体。Pearson相关性分析显示中华纹胸鮡群体间FST与地理距离呈正相关关系(r>0.6)，说明中华纹胸鮡群体的遗传结构是距离隔离造成，符合IBD模式。距离和生境的差异是中华纹胸鮡产生群体隔离的可能原因。中华纹胸鮡为底栖鱼类，产粘性卵，活动范围窄，本研究采样点最大距离达540 km(表4)，攀枝花到重庆巴南江段海拔落差达920 m左右，金沙江段内地形极为复杂，众多高山深谷相间并列，流域内气候时空变化较大，垂直差异十分显著。这些因素阻碍了中华纹胸鮡迁移，基因流分析结果也表明了三个谱系间几乎没有基因交流。

历史研究表明中华纹胸鮡存在多个生态类群[5]，本研究结果显示中华纹胸鮡存在多个遗传谱系，为了揭示中华纹胸鮡是否存在隐存种或亚种，本研究用两种方法进行检验，即遗传距离和ABGD法。首先，根据李博[30]基于Cytb序列计算三种纹胸鮡属鱼类[中华纹胸鮡、扎那纹胸鮡、穴形纹胸鮡(G.cavia)]之间遗传距离，种间遗传距离在0.067～0.098之间；王利华等[31]基于Cytb序列计算的鲌属6种鱼类的种间遗传距离在0.1～7.5。本研究中中华纹胸鮡9个群体间的遗传距离在0.003～0.020，单倍型之间遗传距离最大为0.027(HB_21-HB_61)，远小于鲌属和纹胸鮡属种间的遗传距离。其次，通常用DNA条形码ABGD作种类鉴定时，一般以P=0.01显示组数作为种类数量的依据[20]。本研究COI基因的ABGD分析结果表明，所分析的样本在P=0.01下只有一个组即为同一个种。综合以上两种方法结果表明，长江上游及金沙江中华纹胸鮡群体中并无隐存种或亚种存在。

本研究的采样江段中建有两座电站，即向家坝水电站和溪洛渡水电站，位于宜宾和巧家之间，AMOVA分析显示水电站上游(BN、JJ、HJ、LZ、NX、QW)、下游(QJ、NN、PZH)群体有明显的遗传结构差异(表2)，但这种差异是种群历史造成，而不是大坝阻隔的影响。Sturcture分析显示Clade I(主要是大坝下游群体)包含有大坝上游的少量个体，而Clade II和III完全没有大坝下游个体(图3)，目前大坝运行距采样时间仅5年，难以对群体差异产生如此大的影响。并且，攀枝花江段与巧家江段之间目前未受大坝阻隔，但两群体间也发生了显著遗传分化，因此，大坝建设是否会对群体基因交流造成影响，有待今后长期监测。

3.3 种群历史

BSP分析表明，中华纹胸鮡长江上游群体(Clade I)在晚更新世(0.13～0.01 Ma)发生了种群扩张事件，此时为末次冰期后期，长江上游气候回暖，有利于中华纹胸鮡群体增长。长江上游一些特有鱼类，如红唇薄鳅(Leptobotiarubrilabris)、长薄鳅(L.elongata)[32]、裸体异鳔鳅鮀、异鳔鳅鮀(X.boulengeri)[29]等也在这一时期检测到种群扩张事件。金沙江下游群体(Clade II)在0～0.0125 Ma也检测到扩张现象(图4)，扩张时间比Clade I的要晚，说明金沙江出现适宜于中华纹胸鮡气候的时间更晚，与其较高海拔相符。

本研究结果显示，长江上游及金沙江中华纹胸鮡群体发生了显著的遗传分化，将金沙江中游、金沙江下游、长江上游群体分别视为一个进化显著单元(ESU)[33]，应对三个单元分别加以保护。巧家、宁南、攀枝花江段中华纹胸鮡群体遗传多样性较低，对环境变化的适应能力较低，在金沙江水电工程开发过程中需要特别加以关注。