海南岛宽额鳢(Channa gachua)群体遗传变异与生物地理过程

陈佳琪,李 潮,张雯君,李 炜,高天扬,赵 俊

广州市亚热带生物多样性与环境生物监测重点实验室,广东省水产健康安全养殖重点实验室,广东省水产优质环保养殖工程技术研究中心,华南师范大学生命科学学院,广州 510631

海南岛是我国纬度最低、面积最大的热带岛屿,该岛北以琼州海峡与广东划界,西临北部湾与越南遥望。由于火山作用和海平面上升,海南岛最初在大约2—2.5百万年前从大陆中分离出来[1- 2],而地质证据表明,在更新世冰期期间,海南岛与大陆有数次的连接与分离[3- 4]。海南岛的地形复杂,以五指山和鹦哥岭两座山脉为隆起核心,从中部和南部地区陡峭上升,并向北延伸到一个宽阔的平原,水系从中部山区向周围延伸组成辐射状水系。南渡江为最大河系,向北流入琼州海峡,万泉河与陵水河则向东流入南海,昌化江向西流入北部湾,藤桥河向南经三亚藤桥港独立入海。随着DNA序列测定技术和生物信息学的发展,近年来学界陆续涉足了海南岛淡水鱼类生物地理与种群遗传的研究,研究表明海南岛受到地理屏障影响的淡水鱼类在种群遗传上普遍存在一定分化,如细尾白甲鱼(Onychostomalepturum)[5]、东方墨头鱼(Garraorientalis)[6]和拟平鳅(Liniparhomalopteradisparis)[7]等。对于海南岛这一大陆性岛屿进行生物地理学研究,可帮助了解生物类群在大陆与岛屿间分布的演化历史,也可以协助了解该岛屿的地质与气候变化历史。

元江(红河)发源于我国云南省,在越南境内自北部湾入海,近年来的古地磁学研究表明,海南岛原属于广西和越南北部[8- 9],从中生代晚期到新生代早期,海南岛沿红河断裂带向东南方向移动和旋转到现在的位置[10- 11],且在1986年出版的《海南岛淡水及河口鱼类志》[12]中记述了海南岛淡水鱼类与元江及红河水系相同的有62种,所以,元江及其与海南岛西南部水系在地质历史上可能存在某种联系,添加元江样本作为本研究的外类群可帮助了解海南岛淡水鱼类的生物地理过程。

宽额鳢(ChannagachuaHamilton, 1822)亦称南鳢或白边鳢,隶属于鲈形目(Perciformes)鳢科(Channidae)鳢属(Channa),常在溪流和河沟的浅水区域生活。国内主要分布于海南岛各水系[12]、云南各江河及其支流等[13],国外主要分布于印度及东南亚各地[12]。目前,国内外学者主要对宽额鳢的毒理学[14- 15]、生理学[16- 17]、染色体组型分析[18]、同工酶[19]以及宽额鳢与同属鱼类的亲缘关系[20- 21]等方面开展了一些研究。

遗传变异和种群遗传结构的量化对于改善自然种群的管理和保护至关重要[22],但是对于宽额鳢这一物种内的遗传分化和亲缘地理的研究,目前为止国内仅有一篇关于宽额鳢遗传变异的报道[23],该文章主要测定了云南澜沧江和海南岛南渡江、万泉河与昌化江4个水系宽额鳢的线粒体控制区序列,结果表明云南组群和海南组群间高度分化。然而,该研究对于海南岛种群的地理覆盖范围相对缺乏,仅选取了三个水系,且在与海南岛淡水鱼类的亲缘关系上,云南元江相比澜沧江更能说明海南岛的历史生物地理过程。此外,鱼类Cytb基因的长度为1140 bp左右,它的结构与功能在mtDNA的13个蛋白质编码基因中被了解得最清楚,而且进化速率适中,被认为是解决分类和系统进化问题最可信的分子标记之一[24]。因此,本研究旨在通过线粒体Cytb序列的测定和分析,解决如下科学问题：(1)了解海南岛内宽额鳢的谱系分化和遗传多样性现状,试图阐述海南岛内宽额鳢经历的历史生物地理过程。(2)推断可能影响海南岛宽额鳢的遗传多样性和生物地理过程的历史过程或事件。(3)检验海南岛与元江水系宽额鳢之间的亲缘地理关系,重点检验昌化江与元江的关系,以期为未来海南岛淡水鱼类保护区的规划和建设提供理论支持,为我国鱼类的生物多样性保护和管理提供借鉴依据。

1 材料与方法

1.1 样品采集

于2017年7月、10月,在海南岛及云南的12个采样点采集168条宽额鳢个体(样品体长范围为3.9—13.5 cm),将所采集个体按采样的水系进行种群划分,即用水系名称对种群进行命名。各具体采样点及样品信息见表 1,样品分布地图见图1。采集后的样本立即取其新鲜肌肉并用95%的乙醇固定,带回实验室置于-20℃低温冰箱作长期保存。样本存放在华南师范大学生命科学学院鱼类标本室。

表1 宽额鳢样品信息

图1 宽额鳢样品分布图Fig.1 Sampling locations forC.gachua各采样点信息及水系代码见表 1

1.2 DNA的提取及测序

基因组DNA的提取主要使用从上海生工生物工程技术服务有限公司购买的Ezup柱式动物基因组DNA抽提试剂盒所示方法,略有改动。所提取的总DNA溶于TE溶液,-20℃长期保存。使用通用引物L14724(5′-GACTTGAAAAACCACCGTTG- 3′)和H15915(5′-CTCCGATCTCCGGATTACAAGAC- 3′)[25]扩增Cytb片段,引物由上海生工生物工程技术服务有限公司合成。PCR扩增在25 μL反应体系中进行,包括2×Taq PCR Master Mix(诺唯赞)12.5 μL、上下游引物各1.0 μL(浓度为0.01 mmol/L)、DNA模板1.0 μL,ddH2O 9.5 μL,每一次反应都包括阴性对照。具体的反应程序为：94℃预变性4 min,94℃变性30 s,49℃退火45 s,72℃延伸1 min,共35个循环,最后在72℃下延伸8 min。PCR产物通过1% 琼脂糖凝胶电泳检验目的片段的扩增情况后,送往上海生工生物工程技术服务有限公司纯化并进行双向测序。

1.3 序列比对及遗传多样性分析

利用DNASTAR软件包的SeqMan程序对序列进行正反链拼接,再将拼接文件导入Clustal X 2.1程序进行序列比对,比对结果通过MEGA 7.0软件[26]再进一步进行人工校对,并计算序列的碱基组成,可变位点数等。利用DnaSP v5.0软件计算单倍型数及单倍型在各种群的分布情况,并利用单倍型多样性(Haplotype diversity,h)[27]和核苷酸多样性(Current genetic diversity estimates,θπ及Historical genetic diversity estimates,θω)[28]指数评估种群遗传多样性水平。当前遗传多样性(θπ)是基于序列之间的成对变异进行估算,而历史遗传多样性(θω)是基于序列中分离位点的数量进行估算[28],通过比较当前遗传多样性(θπ)和历史遗传多样性(θω),可以深入了解近代进化史的种群动态[29]。

1.4 系统发育分析

利用MEGA7.0软件基于Kimura′s two-parameter模型计算遗传距离,并基于邻接法(Neighbor-joining, NJ)构建NJ树。基于最大似然法(Maximum likelihood, ML)构建的ML树,采用在线软件PhyML 3.0 (http://www.atgc-montpellier.fr/phyml/)构建,并通过SMS[30]采取Akaike信息量准则(Akaike Information Criterion, AIC)得出最佳核苷酸替换模型为GTR+G+I；利用MrBayes 3.1.2进行贝叶斯推断(Bayesian inference, BI)分析,用Treeview显示拓扑结构图。

1.5 种群历史动态及遗传分化分析

为检验各谱系是否发生过扩张,利用DnaSP v5.0软件对各谱系进行Tajima′sD、Fu′Fs检验及核苷酸错配分布(Mismatch-distribution)分析。使用 Network 5软件基于Median-joining法构建单倍型亲缘关系网络图。利用 Arlequin 3.5软件[31]进行分子方差分析(AMOVA),估算遗传变异在群体间和群体内的分布, 并计算群体间遗传分化指数(FST)及其P值。使用DnaSP计算两个遗传分化指数GST和NST,检查系统地理结构的存在。通过BEAUti v1.7.5设置参数,将进化替代模型设为GTR+G+I,将生成的文件导入BEAST v1.7.5中计算节点的分化时间,而后再通过Tracer采用Coalescent：Bayesian Skyline,对每个谱系进行了有效种群大小的历史动态反演,其结果以贝叶斯天际线(Bayesian skyline plots, BSP)方式呈现。

2 结果与分析

本研究中,6个种群的168个序列共获得36个单倍型(序列全长1142bp,可变位点125个),A、C、T、G平均含量分别为30.4%、13.8%、25.3%和30.5%。

2.1 系统发育分析

为更好地了解宽额鳢种群之间的亲缘关系,采用不同方法(NJ,ML和BI)对36个单倍型构建分子系统发育树,由3种方法得出的拓扑结构相似,因此把3种方法构建的系统树合为一致树(图2)。从系统发育树中可以看出,36个单倍型分为两大支系(A支和B支)。A支包含海南岛5个种群的全部单倍型,B支为云南元江种群。其中,A支中可以分为两个亚支,一个是A1亚支,包括海南岛5个种群(CH、LS、TQ、WQ及ND)的大部分单倍型,另一支为A2亚支,包括昌化江的两个采样点即畅好乡和水满乡的全部单倍型。

图2 基于Cytb序列单倍型构建的宽额鳢群体NJ, ML和BI树Fig.2 Neighbor-joining, Maximum likelihood and Bayesian inference tree ofC.gachuabased on Cytbgene haplotypesNJ树：邻接树,Neighbor-joining tree；ML树：最大似然树,Maximum likelihood tree；BI树：贝叶斯推断树,Bayesian inference tree;图中的A与B分别表示宽额鳢海南岛与云南组群的谱系,A1与A2分别为谱系A的亚支；分支节点数字分别表示NJ, ML和BI树的支持率；Hap：单倍型,Haplotype

在单倍型网络图(图3)中,3个主要分支结构与系统发育树结果一致,谱系B与A1、A2亚支至少有40个突变步数。系统发生的结果支持云南元江种群在亲缘关系上与海南岛的5个种群相距较远,且最先分化出来,处于最为古老的分支地位,接着是昌化江的部分个体分化出来,而后其余的海南岛各种群相互聚为一支。

图3 基于Cytb基因构建的宽额鳢单倍型网络图Fig.3 The haplotype network ofC.gachuabased on Cytbgene haplotypes图中 A1、A2、B 分别为宽额鳢系统树中的A1亚支、A2亚支、谱系 B；数字和虚线框代表不同谱系单倍型间的核苷酸替换数目；每个圆圈代表一个单倍型,圆的面积与单倍型频率成正比

2.2 种群遗传多样性

表2显示了宽额鳢各种群的遗传多样性。6个种群的平均单倍型多样性为0.677(范围：0.299—0.872)。除藤桥河种群的单倍型多样性较低(0.299)外,其他种群的单倍型多样性均较高,南渡江种群的单倍型多样性最高(0.872)。对每个种群当前(θπ)和历史(θω)遗传多样性的计算表明各种群的θω几乎均趋近于θπ,说明各种群在历史上维持着相对稳定的状态,几乎无发生种群扩张。

各种群单倍型数目在4—9个之间,每个种群都具有各自的单倍型。特有单倍型共35个,占单倍型总数的 97.22%；共享单倍型仅 1个,分布在昌化江种群与万泉河种群中,占单倍型总数的 2.78%。

2.3 谱系间遗传距离及遗传分化

表3为宽额鳢基于系统发育树的3个谱系之间的遗传距离和遗传分化指数。结果表明,云南宽额鳢的谱系B的遗传距离与海南岛宽额鳢中的两个亚支A1和A2的遗传距离均较大(分别为0.0659,0.0537),表明谱系B与其余两个谱系的关系都较远,但相对A2亚支较近。遗传分化指数FST值也呈现出高度的分化(范围在0.9045—0.9596)。

固定指数NST和GST的比较显示,NST大于GST(分别为0.8184和0.2858)。这一结果表明种群间存在明显的谱系地理结构。

2.4 分子方差分析(AMOVA)

利用分子方差分析方法AMOVA (Analysis of Molecular Variance)可以对研究的种群进行不同层次的归类与划分,估计出组群之间、组群内以及个体之间不同层次表现出的变异占总变异的多少,从而了解变异的分布模式。根据系统树得到的谱系将宽额鳢种群分成3个组群,进行分子方差分析(AMOVA)。其中云南元江水系为一个组群；海南岛诸水系为一个组群,包括昌化江的部分种群、陵水河、藤桥河、万泉河及南渡江；海南岛内昌化江的两个采样点即畅好乡和水满乡为一个组群。结果表明,绝大部分变异发生在组群间(74.73%),而种群内很小(5.87%) (表4),说明宽额鳢种群的变异主要来源于组群间。

2.5 种群动态分析

对各个谱系进行中性检验分析,结果表明,A1亚支和A2亚支的Tajima′sD和Fu′sFs 值均呈现负值,但不显著(P>0.05),表明其在历史上均未发生过明显的种群扩张。单系群B的Tajima′sD和Fu′sFs 值均为正值(0.25163, 2.93079),表明云南元江种群亦未发生过种群扩张(表5)。且核苷酸错配分布分析发现,各谱系动态均呈多峰分布(图4),其中A2亚支的峰形虽接近单峰,但因其顶峰处于边缘临界点,且顶峰过后发生过转折,再加上A2亚支的个体数相对较少,这在核苷酸错配分布分析中会产生一定的影响,结合A2在中性检验中负值不显著的结果,可确定为多峰分布。这都表明了3个谱系在历史进化的过程中处于相对稳定的状态,基本无经历过种群扩张。贝叶斯天际线分析(Bayesian skyline plots, BSP)显示A1亚支与单系群B在约1万年前时有效种群数量趋于减小,A2亚支则未发现种群扩张或减小的趋势(图5)。

表2 基于 Cytb基因的遗传多样性分析

θπ：当前核苷酸多样性指数,即成对序列间核苷酸变异的平均数,the average number of nucleotide differences per site between two sequences；θω：历史核苷酸多样性指数,即序列中分离位点的数目,the number of segregating sites among the sequences

表3 宽额鳢谱系间的遗传距离(对角线下方)和遗传分化指数FST(对角线上方)

表中 A1、A2、B 分别为宽额鳢系统树中的 A1亚支、A2亚支、谱系 B

表4 宽额鳢群体分子方差分析

表5 基于Cytb的宽额鳢Tajima′sD和 Fu′sFs 检验

图4 基于Cytb构建的宽额鳢3个谱系的核苷酸错配分布分析Fig.4 Mismatch-distribution ofC.gachuain three lineages based on Cytbgene实线代表期望值,带有圆环的虚线代表观测值

图5 基于Cytb构建的宽额鳢3个谱系的贝叶斯天际线图Fig.5 Bayesian skyline plots ofC.gachuain three lineages based on CytbgeneNe代表有效种群大小,τ代表生物体的世代时间；黑线表示中值种群大小,彩色线包含区域表示95%置信区间

3 讨论

3.1 宽额鳢种群遗传多样性及遗传结构

种群遗传多样性是生物多样性的重要组成部分,体现了物种对不同环境的适应能力的潜力,遗传多样性水平、形成机制和分布格局可揭示物种起源与进化历程,是物种进化潜能的保证,也是物种保护工作的基础[32]。本研究采集了海南岛与云南不同水系6个种群共计168个野生宽额鳢样本,共检测出36个线粒体单倍型,平均单倍型多样性(h)为0.6773,显示出高单倍型多样性。而核苷酸多样性的积累时间比单倍型多样性的积累时间更漫长,在一定程度上核苷酸多样性比单倍型多样性更能反映群体遗传多样性[33],所以可作为衡量一个种群线粒体DNA遗传变异与种群历史动态的重要指标。本研究中,海南岛与云南各种群的历史遗传多样性(θω)几乎均趋近于当前遗传多样性(θπ),说明各种群在历史上相对稳定,并未因种群扩张或减小而造成θω急剧升高或降低。此外,核苷酸多样性普遍较小这一结果,与东南亚地区关于鳢属鱼类的研究结果[34- 36]十分相似,然而却远低于海南岛其他物种的核苷酸多样性[5- 6, 37- 38]。这就说明,并不是环境压力或瓶颈效应而造成宽额鳢的核苷酸多样性剧烈下降,如果存在瓶颈效应,则相同的地质事件或气候动荡会对相同分布的物种有相似的影响。对于宽额鳢这一鳢属鱼类核苷酸多样性较低的解释是,该物种属于热带、亚热带暖水性淡水鱼类,常栖息于河沟及溪流的浅水区域中[39],生活环境单一,所受外界环境刺激相对较小,故这一物种对于适应环境的要求没有很高,遗传物质发生突变的几率便小了许多,突变没那么容易被保留下来,导致宽额鳢核苷酸多样性水平贫乏。

分子方差分析(AMOVA)中,根据系统树的分支结果进行了分组,确定了种群内5.87%的变异,组群内种群间19.40%的变异以及组群间74.73%的变异(表4),说明宽额鳢种群变异主要来源于组群间,即各谱系间的变异最为明显。谱系间的遗传距离和遗传分化指数FST结果(表3)也表明,谱系B(元江种群)与其他两个谱系A亚支的遗传距离均较大(分别为0.0659, 0.0537),且各谱系间的FST值均大于0.9,说明6条水系宽额鳢Cytb序列的遗传分化主要是由于山脉或水系的隔离导致,海南岛在与大陆分离后,海南岛水系与云南水系不再发生联系,故海南岛和云南的宽额鳢种群发生了较大的遗传变异,而在海南岛的内部水系中,位于昌化江海拔较高两个地点(畅好乡和水满乡)的宽额鳢由于地势原因,与海南岛平原地区种群缺乏基因交流,从而发生了明显的遗传分化。

3.2 宽额鳢的种群历史动态

中国南方较为温和的更新世气候和多样的地形为宽额鳢这一热带物种提供了许多适宜的栖息地。本研究的中性检验及核苷酸错配分布结果证明了3个谱系均不存在种群扩张的现象。BSP结果表明,A2亚支在历史上没有发生过相对剧烈的种群波动,而A1亚支与谱系B都呈现了约1万年前的有效种群数量减小,但其幅度不大。推测这是由于更新世冰期期间,宽额鳢种群不可避免地受到了温度变化和海平面变化的影响而发生一定程度的有效种群数量减小,而A2亚支群体的栖息地(畅好乡和水满乡)由于位于五指山山势较高(1000 m左右)的河流上,受冰河时期温度影响较小,所以可能保持了相对稳定的生态位,并没有经历剧烈的种群波动。

3.3 宽额鳢的亲缘地理过程

亲缘地理学是生物地理学的研究核心[44],其关注生物遗传变异与地理之间的关系。揭示某一地区物种在历史上经历的生物地理过程,有助于我们在制定保护规划和措施的时候做到有的放矢。在本研究中,根据NJ,ML及BI系统树、单倍型分布的特征及单倍型网络图,可证实基于Cytb序列的宽额鳢主要存在两个谱系(A、B),两个谱系均显示出高度的地理结构。谱系A进一步分为两个亚支,其中A1亚支的成员来源于海南岛5个水系的大部分个体,A2亚支的成员仅限于昌化江上游地势较高的2个地区(畅好乡和水满乡),而薛丹等[23]基于宽额鳢群体与Zhou等[5]基于细尾白甲鱼群体的系统地理学研究结果亦均表明昌化江存在2个谱系,这与本研究结果基本一致,证实了昌化江确实存在两个谱系。推测这是由于在海南岛内,宽额鳢的祖先种群最先从岛外定居至昌化江,伴随着五指山和鹦哥岭的隆起,部分祖先种群向东扩散至陵水河、藤桥河,向西扩散至南渡江进而进入万泉河,且由于陵水河、藤桥河、南渡江和万泉河的大部分江河均位于平原,宽额鳢种群间能发生基因交流的机会较大,故这4个种群在系统树上并没有各自形成明显的谱系,即系统树海南组群中的A1亚支,但4个种群依旧在山脉和地理隔离的影响下,均发生了一定程度的分化；而扩散至昌化江上游的宽额鳢祖先种群在山脉隆起过程中,栖息于地势较高(如水满乡、畅好乡等)的河沟等浅水区域中,其中,水满乡是海南岛地理位置最高的乡镇,畅好乡的地势也较高,故推测这一部分祖先种群所在水系由于地势原因而与平原地区的水系发生了隔离,与其他种群之间鲜有交流,从而保留了海南岛中古老的祖先种群的遗传特性,此即为系统树海南组群中,最先分化出来的A2亚支。谱系B的成员均为云南元江种群,且从系统树来看,元江种群处于最为古老的分支地位,为最先分化出来的种群,这也为海南岛种群可能起源于元江(红河)种群提供了证据。

近年来的古地磁学研究表明,海南岛原属于广西和越南北部,而后沿红河断裂带向东南方向移动和旋转到现在的位置,所以海南岛西南部水系与元江(红河)水系在历史上可能存在一定的联系。更新世期间,冰川和冰川消融以及随着海平面的降低和升高已经极大地影响了东南亚的陆地结构和动植物分布[3]。冰期期间,北部湾因海平面下降而暴露,北部湾和海南岛在内的整个地区成为亚洲大陆沿海平原的一部分,并且,越南北部水系(包括元江/红河)通过雷州半岛和海南岛西南部流入南海,这有助于鱼类的迁移扩散[5]。相关地质证据和生物学研究也表明,海南岛和越南北部的种群之间存在密切的进化关系[37,45- 46]。

据推断,在中更新世冰期期间,大陆和岛屿的水系因海平面下降而发生接触,越南北部水系通过雷州半岛和海南岛流入南海,因此元江(红河)和海南岛西南部的昌化江可能存在一条能够联系的古河道,宽额鳢借此进行了由元江(红河)至昌化江的种群扩散。后来的间冰期,北部湾海平面上升,海南岛与大陆隔离,元江(红河)与海南岛西南部水系的基因交流被中断,海南岛和元江(红河)的宽额鳢分成两个分支并随后发生分化。海南岛内部水系亦因五指山鹦哥岭山脉的隆起而产生地理隔离,从而昌化江、南渡江、万泉河、陵水河及藤桥河均发生一定程度的分化。这与Lin等人[47]关于蜡皮蜥(Leiolepisreevesii)的种群遗传结构研究中提出的扩散路线：蜡皮蜥从越南沿海扩散到海南,然后从海南向北扩散至中国大陆的结果十分相似,进一步支持海南岛淡水鱼类与越南北部水系(包括元江/红河)的鱼类关系更为密切。在未来的研究中,需要更多来自中国云南其他水系和越南沿海独立水系的宽额鳢样本,以进一步完整确定其何时何地以及如何迁移至海南岛。

3.4 宽额鳢保护对策的探讨

本研究探讨了宽额鳢海南岛各水系与云南元江的遗传变异与系统地理结构,通过Cytb基因对宽额鳢的研究结果显示,宽额鳢6个种群在Cytb的谱系上呈现了2个主要谱系的关系。其中,云南元江种群无论是遗传分化指数、AMOVA还是遗传距离等都表现出与其他海南地区种群的遗传分化程度最高,存在一定程度的隔离,因此本研究认为云南元江种群很可能已经进化为一个独特的进化显著单元(Evolutionarily Significant Unit, ESU)；在海南组群中,昌化江的部分种群作为独立亚支最先分化了出来,遗传分化指数也较高,故也应作为独特的进化显著单元。而根据系统发育关系、单倍型网络图,并且结合水系地理区域,本文认为6个宽额鳢种群可以划分为3个管理单元(Management Unit, MUs)进行保护,即云南元江种群为一个管理单元；海南组群划分为两个管理单元,一个为昌化江种群,另一个为海南其余种群。

致谢：中国台湾台南第一高级中学的林弘都老师给予数据分析支持,特此致谢。