高原林蛙线粒体全基因组序列测定及系统发育分析

张湑泽, 董 豹, 哈金强, 赵 英, 魏生楠, 韦迎婷

(1. 青海民族大学生态环境与资源学院, 西宁 810007; 2. 青海省特色经济植物高值化利用重点实验室, 西宁 810007)

高原林蛙(Ranakukunoris),隶属于无尾目(Anura)蛙科(Ranidae)林蛙属(Rana),是青藏高原特有种[1]。高原林蛙曾被认为是中国林蛙的亚种,属于中国林蛙复合体。谢锋等[2]将中国西北地区的林蛙种群分为两个有效种,即中国林蛙和高原林蛙,以Nikol’ skii等于1918年在青海湖采集的高原林蛙(RanaamurensiskukunorisNikol’ skii 1918)标本为模式标本。江建平等[3]采用Cytb基因对高原林蛙不同居群之间的遗传多样性进行分析,基于分子生物学证据确定高原林蛙是林蛙属单独物种,染色体2n=24。Chen等[4]发现高原林蛙存在两性异型现象;Zhou等[5]认为青藏高原隆升和气候环境变化促进中国林蛙复合体物种形成,且中国林蛙部分群体与高原林蛙存在渐渗杂交现象;基于Cytb基因的高原林蛙群体进化分析,发现青藏高原东北部高原林蛙存在多个冰期避难所,并在祁连山南北形成两个高度分化的世系[6];米志平等[7]研究了高原林蛙皮肤结构,发现其皮肤形态特征与其他物种不同,推测可能因为生活在高海拔地区,生存环境气温低,水分不易蒸发;张湑泽等[8]对青藏高原东北部高原林蛙皮肤结构观察发现,高原林蛙对青藏高原环境有明显的适应性;Yang等[9]通过林蛙属3个物种线粒体全基因组序列,分析蛙科的进化关系。Wang等[10]对来自贺兰山的高原林蛙的线粒体全基因组开展测序,进行系统发育分析。但以上研究还未对高原林蛙线粒体基因组具体编码信息、基因结构、遗传多样性和分子进化开展分析。

线粒体DNA不仅具有结构简单、无重复、母系遗传和进化速率快等特点[11-13],还具有很高的专一性和独特性,是动物系统发育分析和物种鉴定的重要工具[14]。探究线粒体结构可以为物种起源、分化和系统发育分析提供良好的参考。如今越来越多的学者通过对线粒体基因组的研究分析物种间的进化关系[15-17]。本研究基于高原林蛙线粒体全基因组,对其组成和结构进行分析,以期对高原林蛙遗传多样性、分类地位和物种保护进行深入研究。

1 材料与方法

1.1 材料

高原林蛙样本于2020年6月采自青藏高原东北部青海湖周边黑马河镇(36°37′26.63 N,100°07′20.82 E,3 194 m)。样本采集后保存于95%乙醇中,采用标准苯酚/氯仿法从高原林蛙肌肉组织中提取基因组DNA,采用 Illumina NovaSeq平台进行测序,测序读长为PE150,测序获得5G数据量,使用fastp(version 0.20.0,https:∥github.com/OpenGene/fastp)软件对原始数据进行过滤,获得Clean Data总碱基数为6 972 984 300 bp,用SPAdes组装成一个完整的高原林蛙线粒体基因组。线粒体基因组全长为21 913 bp,该序列已经提交NCBI基因数据库(GenBank登录号:MZ043820)。另有下载自GenBank的25份其他物种的线粒体全基因组序列,并参照NCBI上所提交的高原林蛙线粒体全基因组序列(MN733918和KU246049)开展分析。实验样本信息见表1。

表1 用于系统进化分析的物种线粒体DNA序列信息Table 1 Mitochondrial DNA sequence of species for phylogenetic analysis

1.2 方法

1.2.1 全序列特征分析

采用SPAdes v 3.10.1(http:∥cab.spbu.ru/software/spades/)软件组装线粒体基因组,组装不依赖参考基因组。使用Perl脚本进行计算,包括碱基组成与偏移分析、密码子偏好性。使用OGDRAW(https:∥chlorobox.mpimp-golm.mpg.de/OGDraw.html)制作线粒体基因组图谱。使用Mitos2(http:∥mitos2.bioinf.uni-leipzig.de)对组装好的序列进行注释,获得tRNA二级结构。

1.2.2 高原林蛙线粒体基因组编码基因KaKs分析

使用mafft v 7.310(https:∥mafft.cbrc.jp/alignment/software/)软件进行基因序列的比对[18],使用KaKsCalculator v 2.0(https:∥sourceforge.net/projects/kakscalculator2/)软件计算线粒体基因组中编码基因的KaKs值[19]。

1.2.3 遗传距离分析

基于 Kimura 2-parameter 碱基替代模型,使用软件 MEGA11.0[20]计算高原林蛙与其他林蛙属之间的遗传距离。

1.2.4 系统进化分析

基于测序获得的高原林蛙线粒体全基因组序列和GenBank上下载的25条林蛙属物种线粒体全基因组序列以及2条倭蛙属线粒体全基因组序列开展系统发育分析,物种间序列用MAFFT软件(v7.427,auto模式)进行多序列比对,将比对好的数据用BioEdit编辑[21],使用RAxMLv8.2.10[22](https:∥cme.h-its.org/exelixis/software.html)和jModeltest 2.1.7[23]软件进行分析,选用GTR+GAMMA模型, 构建最大似然进化树,Bootstrap=1 000。

2 结果与分析

2.1 高原林蛙线粒体全基因组序列分析

高原林蛙线粒体全基因组序列总长21 913 bp,包括22个tRNA基因、2个rRNA基因、13个蛋白质编码基因和1个D-loop区。其中,蛋白质编码基因序列全长11 287 bp,tRNAs总长1 530 bp,rRNAs总长2 503 bp。

高原林蛙线粒体全基因组序列中A、T、G、C的碱基含量分别为28.41%、30.17%、14.28%和27.13%(图1),A+T的碱基含量为58.58%,G+C的碱基含量为41.42%,表现出明显的A+T偏好性。

图1 高原林蛙线粒体基因组序列结构图谱Figure 1 Sequence map of mitochondrial genome of Rana kukunoris

不同物种在其基因组上各有差异,表2所示为高原林蛙线粒体DNA上的D-loop区、RNA基因、蛋白质编码基因所在位置、蛋白质编码基因的密码子、起始密码子和终止密码子。

表2 高原林蛙线粒体基因组结构Table 2 Mitochondrial genome structure of Rana kukunoris

2.2 高原林蛙蛋白质编码基因分析

高原林蛙线粒体基因组共有13个蛋白质编码基因,包括1个Cytb基因,2个ATP酶(ATP6、ATP8),3个氧化酶亚基(COX1、COX2、COX3)和7个脱氢酶亚基(NAD1、NAD2、NAD3、NAD4、NAD4L、NAD5、NAD6)。在高原林蛙线粒体基因组中,L链和H链都有蛋白质编码基因,但是除了NAD6以外,其他蛋白质编码基因都在H链上。在编码基因中共发现3个重叠基因:COX1/trnS2、ATP6/ATP8、NAD4L/NAD4。重叠基因在表达、调控及进化过程中具有重要作用[24]。

除了NAD1以ATT起始密码子、COX1以ATA起始密码子和NAD4L以GTG起始密码子外,其余10个蛋白质编码基因(NAD2、NAD3、NAD4、NAD5、NAD6、COX2、COX3、ATP6、ATP8)均以ATG为起始密码子。13个蛋白质编码基因总长11 287 bp,A、T、G、C含量分别为25.18%、30.69%、15.08%和29.05%,A+T的总含量为55.87%。所有蛋白质编码基因均呈现AT碱基偏好性(表3)。

表3 高原林蛙线粒体 DNA 蛋白质编码基因核苷酸含量Table 3 Nucleotide content of the protein-coding genes of Rana kukunoris mitochondrial DNA

高原林蛙线粒体基因组中各蛋白质编码基因氨基酸使用情况如图2所示,结果表明,13个蛋白质编码的基因中氨基酸使用最频繁的是Leu和Ser。密码子使用最频繁的为GCC、GCU、CAA、ATA、TCC、ACA,而使用较少的密码子为GCG、CAG、CCG、ACG、TCG。

图2 高原林蛙线粒体基因组蛋白质编码基因的密码子偏好性Figure 2 RSCU of protein coding gene in mitochondrial genome of Rana kukunoris

2.3 高原林蛙tRNA基因结构分析

高原林蛙线粒体基因组上共有22个tRNA基因,总长度1 530 bp,范围从65 bp(trnC)～73 bp(trnL2,trnN),A、T、G、C的含量分别为29.35%、28.37%、22.16%和20.13%,其中,8个tRNA基因(trnP、trnQ、trnA、trnN、trnC、trnY、trnS2、trnE)位于L链上,其余14个tRNA基因在H链上。

通过生物信息学分析对高原林蛙线粒体基因组编码22个tRNA的二级结构进行预测,结果表明:22个tRNA基因都有经典的三叶草结构(图3);在22个tRNA基因的二级结构中出现了碱基错配现象,其中最多的为A-C错配。

图3 高原林蛙线粒体基因组编码的tRNA二级结构Figure 3 Secondary structure of tRNA encoded by mitochondrial genome of Rana kukunoris

2.4 高原林蛙rRNA基因结构分析

高原林蛙的mtDNA中2个rRNA的长度分别为928 bp(rrnS)和1 575 bp(rrnL),两者相距69 bp,被trnV隔开,A、T、G、C的含量分别为32.60%、23.89%、18.58%和24.93%,A+T的偏倚度为0.154。

2.5 高原林蛙线粒体基因组编码基因KaKs分析

高原林蛙线粒体基因组13个编码基因与林蛙属其他物种线粒体基因组序列进行KaKs分析(表4),Ka/Ks值均小于1(P<0.05),表明高原林蛙线粒体基因组编码基因未受到正选择作用[25-26]。

表4 高原林蛙线粒体基因组编码基因KaKs分析Table 4 Analysis of KaKs encoding gene in mitochondrial genome of Rana kukunoris

2.6 林蛙属遗传距离分析

对包括高原林蛙在内的林蛙属24个物种进行基于线粒体全基因组的遗传距离分析。结果表明(表5):对不同地区高原林蛙,青海湖高原林蛙种群与若尔盖高原林蛙种群遗传距离最近,与贺兰山高原林蛙种群遗传距离较远,遗传距离为0.003;对不同物种来说,高原林蛙与中国林蛙、桓仁林蛙遗传距离最近,遗传距离分别为0.032和0.046;与北美洲木蛙、佛罗里达沼泽蛙、牛蛙遗传距离最远,遗传距离分别为0.188、0.183和0.183。

表5 基于线粒体全基因组的林蛙属种内遗传距离分析Table 5 Intraspecific genetic distance of Rana genus

2.7 高原林蛙与林蛙属其他物种的系统进化关系

根据高原林蛙和林蛙属其他物种线粒体全基因组序列,以高山倭蛙和倭蛙为外群,采用最大似然法(maximum likelihood method,ML)构建系统发育树,结果如图4所示,青海湖高原林蛙与若尔盖高原林蛙聚为一支,再与贺兰山高原林蛙聚为一支,与中国林蛙、桓仁林蛙形成姊妹群,东北林蛙、昭觉林蛙等聚为一支;韩国林蛙、昆嵛林蛙和黑龙江林蛙聚为一支,这表明三者亲缘关系最近,可能是因为三者分布的地理位置较近。欧洲林蛙与比利斯牛蛙和中亚林蛙聚为一支,越南趾沟蛙[27]和武夷林蛙聚为一支。牛蛙与佛罗里达沼泽蛙聚为一支,与北美洲木蛙形成姊妹群。

高山倭蛙和倭蛙设置为外群。图4 基于林蛙属26条线粒体全基因组序列构建的最大似然树Figure 4 The maximum likelihood tree constructed based on 26 mitochondrial genome sequences of the genus Rana

3 讨论与结论

线粒体基因组在真核生物中具有高度保守性,是遗传进化分析的重要分子标记[28-29]。线粒体基因组测序及分析是两栖动物遗传多样性研究的重要方式[30]。林蛙属线粒体基因组序列长度大约为18 kb,通常包括13个蛋白质编码基因,22个tRNA,2个rRNA和1个控制区(CR)。高原林蛙线粒体基因组长度约为22 kb,与昆嵛林蛙接近[31]。研究采集的高原林蛙样品来自高原林蛙模式标本产地青海湖,具有更好的代表性。测定的青海湖高原林蛙线粒体基因组的组成结构与林蛙属其他物种相似,全长21 913 bp,其中,控制区长度为6 323 bp,AT含量58.58%,呈现AT碱基偏好性[32]。贺兰山高原林蛙线粒体全基因组由13个蛋白质编码、22个tRNA、2个rRNA基因和1个控制区(CR)组成,总长度为16 644 bp,控制区长度为837 bp,核苷酸成分为A 27.49%、T 29.06%、G 14.77%、C 28.68%,AT含量为56.55%[10]。若尔盖高原林蛙线粒体全基因组组成基因与以上两者相同,其总长度为18 863 bp,控制区长度为3 272 bp,种群核苷酸含量与其接近[9]。三者的线粒体基因组除控制区外序列长度分别为15 590、15 807和15 591 bp,无明显差异。因此,造成青海湖高原林蛙线粒体基因组全长较长的原因主要是其控制区长度较长。

在13个蛋白质编码基因中发现3个重叠:ATP8/ATP6、NAD4L/NAD4和NAD5/NAD6,这种重叠在其他两个高原林蛙线粒体中也都存在。在高原林蛙线粒体基因组中,L链和H链都有蛋白质编码基因,除了NAD6外,其余蛋白质编码基因都在H链上,高原林蛙3个种群结构基本相似。高原林蛙13个蛋白质编码基因中氨基酸使用最频繁的是Leu和Ser,密码子使用最频繁的有GCC、GCU、CAA、ATA、TCC和ACA。高原林蛙22个tRNA的二级结构都具有经典的三叶草结构,出现碱基错配现象,其中最多的为A-C错配。在林蛙属线粒体全基因组序列中,这种序列长度和碱基组成差异主要是因为D-loop区重复以及串联重复元件数目变化[9,33]。

基于高原林蛙线粒体基因组13个蛋白编码基因与林蛙属其他物种线粒体基因组的KaKs分析表明,高原林蛙线粒体基因组编码基因中没有检测到正选择信号,可能与两栖动物基因组相对较慢的进化速率有关[34]。青藏高原东北部地区分布的高原林蛙居群之间存在河流孤岛效应,基于Cytb基因的高原林蛙谱系地理分析发现,高原林蛙在祁连山地区具有明显的冰期避难所,形成两个世系:一个世系在末次冰期后经历了种群数量的突然扩张,向南迁徙到达横断山区;另一世系在祁连山北麓不同河流周边形成了斑块状生境[6]。高原林蛙贺兰山种群与中国林蛙黄土高原种群分布存在重叠,这种地理重叠可能导致线粒体基因组交流与渐渗,并与中国林蛙产生渐渗杂交现象,导致遗传结构产生了一定变化[5]。这表明mtDNA序列数据可以作为基因渗入研究的遗传标记[35]。研究发现,青海湖高原林蛙与若尔盖高原林蛙聚为一支,再与贺兰山高原林蛙形成姊妹群,进一步支持了高原林蛙南北世系的有效性[6],也支持高原林蛙与中国林蛙存在最近共同祖先。

两栖动物由于对气候变化敏感,且生存繁衍离不开水,分布的范围往往表现出边缘效应和距离隔离[36],其种群衰退和灭绝速度更快[37]。遗传多样性是物种适应环境变化的基础,对物种起源、分化、多样性分析具有重要意义[38]。青藏高原东部地区高原林蛙种群分布广泛,遗传多样性高,应通过对栖息地保护,保护其遗传多样性。研究结果将为高原林蛙的起源、分化和遗传多样性研究提供数据支持。