瑞氏红鲂(Satyrichthys rieffeli)基因组survey分析及线粒体基因组注释*

廖贤晖 王乙婷 瞿印权 刘 启 高天翔①

(1.浙江海洋大学水产学院 浙江舟山 316022; 2.武汉万摩科技有限公司 湖北武汉 430076)

高通量测序已成为一种广泛用于动植物基因组测序技术(孟刚等, 2021; 王子寅等, 2023)。通过高通量测序技术可以进行基因组survey 研究, 包括基因组大小估计、重复序列比例和GC 含量计算, 以及对遗传多样性、基因组结构和遗传改良的了解(Barchiet al, 2011; Roweet al, 2011; Xuet al, 2014; Shiet al,2018), 其结果可用于调查基因组图谱和提取线粒体基因组。对于缺乏基因组数据的非模式物种来说, 针对基因组的survey 分析是分子机理研究和基因资源开发的前提(Kimet al, 2011)。

线粒体是半自主细胞性细胞器, 普遍存在于真核细胞内(Boore, 1999)。脊椎动物的线粒体 DNA(Mitochondrial DNA, mtDNA)基因结构紧凑, 长度一般在15～20 kb 之间, 并且编码的基因库非常保守(张方等, 1998)。脊椎动物的线粒体基因组中一般包括了22 个tRNA 基因、2 个rRNA 基因、13 个蛋白质编码区(PCGs)、1 个非编码区包括D-loop 区和1 个轻链复制起始区(OL) (Satohet al, 2016)。线粒体基因组具有结构简单、编码区域高度保守、母系遗传、进化速率快、拷贝数高等鲜明特征(Ruanet al, 2020), 被广泛应用于群体遗传学、进化生物学和谱系遗传学等方面(Koet al, 2018; Zenget al, 2018; 匡卫民等, 2019;Chanet al, 2019)。线粒体物种特异性DNA 片段, 如核糖体RNA (12S rRNA 和16S rRNA)、细胞色素b(Cytb)和细胞色素c氧化酶I (COI)通常用于鱼类物种鉴定(Cutarelliet al, 2018)。目前有较多的学者用线粒体基因组全序列构建系统发育树来比较进化关系,相比于其他的COI、12S rRNA 和16S rRNA 等单一基因序列构建的系统发育树可以减少因信息分布的不均一性和序列长度太短而造成的进化树置信度降低带来的影响, 使得研究结果更精确、可靠(肖家光, 2015)。

1 材料与方法

1.1 样品采集与保存

1.2 实验方法

1.2.1 基因组survey 测序 用标准的苯酚-氯仿法提取肌肉组织的总基因组DNA (Maniatiset al, 1982)。检测合格的DNA 样本通过超声波破碎仪随机打断成长度为300～350 bp 的片段, 经末端修复、加A 尾、加测序接头、纯化、PCR 扩增等步骤完成整个文库制备。构建好的文库通过Illumina nova 进行PE150 测序。使用fastp 软件对测序数据进行数据过滤与去冗余操作(Kajitaniet al, 2014), 过滤后得到的高质量数据将用于后续的杂合度、GC 含量以及基因组大小等信息分析。

1.2.2 基因组大小预估和K-mer 分析 利用过滤后的数据, 取K=17, 参照Liu 等(2013)的K-mer 分析方法对瑞氏红鲂基因组大小、杂合率和重复率等基因特征进行进一步分析。基于K-mer 分析的结果, 获得了关于峰的深度和预测的最佳K-mer 数量的信息,并用于估计基因组的大小。

1.2.3 基因组组装 使用SOAPdenovo 软件(Luoet al, 2012)对过滤去除低质量数据之后的高质量测序数据进行基因组预组装分析。采用K=48 构建Contig和Scaffold, 利用高质量数据进行全基因组装, 使用SOAPdenovo 软件将过滤后的reads 比对拼接到该组装序列上, 获得原始基因组序列及碱基深度。

1.2.4 线粒体基因组组装和注释 借助NOVOPlasty 2.6.3 软件从原始的基因组测序数据中提取线粒体基因组序列(Dierckxsenset al, 2017), 参数采用默认设置, 参考序列为GenBank 瑞氏红鲂COI 序列(GenBank No.MK777566.1)。通过MitoFish将得到的数据进行注释和可视化(Iwasakiet al, 2013)。利用MEGA 11 计算基因组碱基组成､ 氨基酸使用频率及同义密码子相对使用频率(匡卫民等, 2019)。

1.2.5 系统进化分析 GenBank 中下载18 种 鲉形目鱼类的线粒体基因组, 以豹鲂(Dactylopterus volitans)为外群, 基于12 种蛋白质编码基因(除ND6)进行系统发育分析, 确定瑞氏红鲂的分类地位。利用MEGA X 和贝叶斯法(Bayesian)构建最大似然树(ML)和贝叶斯树(BI) (Ronquistet al, 2003), 采用MEGA X 方法估计了最优核苷酸替换模型(Kumaret al, 2018), 基于BIC 标准得到最优GTR+G+I 模型, 确定构建系统发育树的最佳模型。利用MEGA X 软件构建了1 000 个重复的最大似然系统发育树, 探究系统进化关系｡

2 结果

2.1 基因组survey 结果

2.1.1 测序数据统计及质量评估 使用Illumina nova 测序平台进行双末端测序, 使用瑞氏红鲂基因组总DNA 构建300 bp 的文库, 测序得到raw reads 57.59 Gb, 质控后得到约49.13 Gb 的clean reads。瑞氏红鲂质控后其Q20值为96.73%,Q30值为91.72%(表1), 显示本次瑞氏红鲂的高通量测序结果有较高的准确性, 结果可用于后续的分析。

表1 瑞氏红鲂survey 测序结果Tab.1 Results of the survey sequencing of S. rieffeli

表1 瑞氏红鲂survey 测序结果Tab.1 Results of the survey sequencing of S. rieffeli

注: Q20: 质量百分比≥20; Q30: 质量百分比≥30

名字 碱基数 碱基大小/Gb Q20/%Q30/%GC 含量/%原始读数 383 941 058 57.59 96.57 91.64 41.29有效读数 329 460 024 49.13 96.73 91.72 41.24

表2 瑞氏红鲂高通量测序结果与NT 库比对Tab.2 The results of S. rieffeli comparision with the NT library

表2 瑞氏红鲂高通量测序结果与NT 库比对Tab.2 The results of S. rieffeli comparision with the NT library

物种 比对数量/bp 比对总数/bp 占比%Epinephelus 石斑鱼属 1 185 5 945 19.93 Cottoperca 杜父鲈images/BZ_285_912_795_944_826.png属 672 5 945 11.3 Lateolabrax 花鲈属 549 5 945 9.23 Plectropomus 鳃棘鲈属 530 5 945 8.92 Sparus 鲷属 276 5 945 4.64 Scophthalmus 菱鲆属 267 5 945 4.49 Sebastes 平鲉 属 160 5 945 2.69 Nibea 黄姑鱼属 144 5 945 2.42 Anarrhichthys 狼鳗属 134 5 945 2.25 Trachurus 竹荚鱼属 122 5 945 2.05 Myripristis 锯鳞鱼属 121 5 945 2.04 Sphaeramia 圆竺鲷属 120 5 945 2.02 Larimichthys 黄鱼属 120 5 945 2.02 Pseudochaenichthys 拟冰鱼属 108 5 945 1.82 Acanthopagrus 棘鲷属 108 5 945 1.82 Cyprinus 鲤属 107 5 945 1.80 Perca 鲈属 100 5 945 1.68

2.1.3 K-mer 分析 使用K-mer 分析预测了瑞氏红鲂的基因组特征及杂合率和重复比例。选取K=17, 预测瑞氏红鲂基因组大小为827 Mb, 修正后的基因组大小为813 Mb; 基因组杂合率为0.92%,重复序列比例为45.97% (表3)。其K-mer 深度在48时出现一个期望深度, 且没有出现其他峰值(图1),说明瑞氏红鲂基因组杂合度较低。

图1 瑞氏红鲂K-mer 深度分布图Fig.1 Distribution of the K-mer depth of S. rieffeli

表3 瑞氏红鲂基因组特征统计Tab.3 Statistics of S. rieffeli genome characteristics

表3 瑞氏红鲂基因组特征统计Tab.3 Statistics of S. rieffeli genome characteristics

样品 K-mer 数 K-mer 深度 基因组大小/Mb 修正后基因组大小/Mb 杂合率/% 重复率/%S. rieffeli 43 850 343 050 48 827 813 0.92 45.97

表4 Contigs 组装结果统计Tab.4 Statistics of the contigs assembly

2.2 线粒体基因组注释

2.2.1 线粒体基因组组成及其特征 从survey 结果里面提取并组装线粒体基因组, 结果显示瑞氏红鲂的线粒体全序列长16 527 bp, 其中有38 个基因获得注释, 分别由22 个tRNA、13 个蛋白质编码基因(PCGs:APT6、ATP8、Cyt-b、COXI～COXIII、ND1～ND6和ND4L)、2 个rRNA (12S rRNA和16S rRNA)和1 个D-loop 区组成(图2, 表5)。

图2 瑞氏红鲂线粒体基因组的圆形图谱Fig.2 The circular map of the mitogenome of S. rieffeli

表5 瑞氏红鲂的线粒体基因组结构特征Tab.5 Structural characteristics of the mitochondrial genome of S. rieffeli

表5 瑞氏红鲂的线粒体基因组结构特征Tab.5 Structural characteristics of the mitochondrial genome of S. rieffeli

点 终止位点 长度/bp 氨基酸 起始密码子 终止密码子 反密码子 基基因 起始位 因间隔/bp 链tRNA-Phe 1 68 68 GAA 0 H 12S rRNA 69 1 014 946 0 H tRNA-Val 1 015 1 086 72 TAC 0 H 16S rRNA 1 087 2 785 1 699 0 H tRNA-Leu 2 786 2 859 74 TAA 0 H ND1 2 860 3 834 975 324 ATG TAG 0 H tRNA-Ile 3 839 3 908 70 GAT 4 H tRNA-Gln 3 908 3 978 71 TTG –1 L tRNA-Met 3 978 4 046 69 GAT –1 H ND2 4 047 5 092 1 046 348 ATG TA 0 H tRNA-Trp 5 093 5 163 71 TCA 0 H tRNA-Ala 5 165 5 233 69 TGC 1 L tRNA-Asn 5 235 5 307 73 GTT 1 L tRNA-Cys 5 347 5 412 66 GCA 39 L tRNA-Tyr 5 413 5 482 70 GTA 0 L COI 5 484 7 034 1 551 516 GTG TAA 1 H tRNA-Ser 7 035 7 105 71 GCT 0 L tRNA-Asp 7 109 7 181 73 GTC 3 H COII 7 189 7 879 691 230 ATG T 7 H tRNA-Lys 7 880 7 953 74 TTT 0 H ATPase 8 7 955 8 122 168 55 ATG TAA 1 H ATPase 6 8 113 8 795 683 227 ATG TA 0 H COIII 8 796 9 580 785 261 ATG TA 0 H tRNA-Gly 9 581 9 652 72 TCC 0 H ND3 9 653 10 001 349 116 ATG T 0 H tRNA-Arg 10 002 10 070 69 TCG 0 H ND4L 10 071 10 367 297 98 ATG TAA 0 H ND4 10 361 11 741 1 381 460 ATG T –5 H tRNA-His 11 742 11 810 69 GTG 0 H tRNA-Ser 11 811 11 878 68 TGA 0 H tRNA-Leu 11 884 11 956 73 TAG 5 H ND5 11 957 13 759 1 803 600 ATG TAA 0 H ND6 13 792 14 313 522 173 ATG TAG 32 L tRNA-Glu 14 314 14 382 69 TTC 0 L Cyt b 14 388 15 528 1 141 380 ATG T 5 H tRNA-Thr 15 529 15 600 72 TGT 0 H tRNA-Pro 15 600 15 669 70 TGG –1 L D-loop 15 670 16 527 858 0 H

除8 个tRNA 基因(tRNA-Gln、Ala、Asn、Cys、Tyr、Ser、Glu、Pro)和1 个PCGs 基因(ND6)在线粒体基因组的L 链上, 其他29 个基因和D-loop 区均在H 链上, 基因间都未发生重排。其中11 个基因发生了间隔, 共99 bp,tRNA-Cys基因间隔最长, 共39 bp,其次是基因ND6, 共32 bp; 有4 个基因发生了重叠,共8 bp, 重叠的基因是tRNA-Gln、tRNA-Met、ND4和tRNA-Pro, 其中基因ND4 重叠最多, 有5 bp。瑞氏红鲂线粒体全基因碱基组成: A (27.4%)、T(26.9%)、G (16.8%)、C (28.9%) (表6), 各基因的GC含量在44.4%～51.4%之间, 其中A+T 的含量(54.3%)大于G+C 的含量(45.7%)。

表6 线粒体基因组核苷酸组成比例Tab.6 The composition ratio of mitochondrial genome nucleotide

图3 瑞氏红鲂13 个蛋白质编码基因的AT-skew 和GC-skew 值Fig.3 The AT-skew and GC-skew values of 13 protein-coding genes in S. rieffeli

图4 瑞氏红鲂tRNA-Ser (GCT)二级结构预测图Fig.4 Prediction of secondary structure of tRNA-Ser (GCT) in S.rieffeli

图5 瑞氏红鲂线粒体基因密码子使用频率Fig.5 Frequency of usage of mitochondrial gene codon of S. rieffeli

图6 瑞氏红鲂线粒体基因相对同义密码子使用频率Fig.6 Frequency of relative synonymous codon usage (RSCU) in S. rieffeli

2.2.5 非编码区 非编码区又称之为 D-环区(displacement-loop region), 位 于tRNA-Pro和tRNA-Phe基因之间, 是整个线粒体基因组序列和长度变异最大的区域, 但其中也包含保守片段。瑞氏红鲂的非编码区各碱基含量分别为T (30.0%)、C (21.4%)、A (32.3%)、G (16.3%), 表现出明显的AT 偏好性, A+T的含量达到了62.3%, G+C 的含量仅为37.7%。

2.3 系统发育关系

3 讨论

近年来, 高通量测序技术的迅速发展为解决非模式海洋鱼类的基因组问题提供了有效的技术手段。利用K-mer 方法可以估计非模型物种的基因组大小,无须任何先验知识(Liet al, 2019b), 其方法已被应用于绒杜父鱼(Hemitripterus villosus) (赵蕊蕊等, 2022)、褐菖 鲉 (Sebastiscus marmoratus) (Xuet al, 2020)、菖鲉属(Sebastiscus) (Jiaet al, 2021) 等 鲉形目鱼类的研究。目前瑞氏红鲂的分子信息还不完善, 本研究survey 结果显示GC 含量在41%左右出现峰值(双端测序结果都出现在41%左右), 且没有出现明显的偏差, 表明测序结果没有偏向性; 并获得49.13 Gb 的clean reads,Q20大于90%, 基因组大小为827 Mb, 修正后的基因组大小为813 Mb, 其杂合率为0.92%, 重复率为45.97%。瑞氏红鲂与大多数的海洋鱼类相似, 其基因组大小稍大于绒杜父鱼713.18 Mb (赵蕊蕊等, 2022) 、褐菖 鲉812.86 Mb (Xuet al, 2020)、艾氏蛇鳗(Ophichthus lithinus) 755.24 Mb (宁子君等,2022), 稍小于许氏平鲉 (Sebastes schlegelii) 846.36 Mb、朝鲜平 鲉(Sebastes koreanus) 832.53 Mb 、金斑平鲉(Sebastes nudus) 813.12 Mb (Xuet al, 2019)、双带缟虎(Tridentiger bifasciatus) 887.60 Mb (Zhaoet al,2022), 导致物种间基因组大小的差异可能是由于重复序列含量差异所导致。瑞氏红鲂的重复率为45.97%, 较绒杜父鱼38.61%、 褐菖 鲉39.65%、艾氏蛇鳗43.30%、 许氏平 鲉44.10%、 朝鲜平 鲉43.65%、金斑平 鲉41.40%、双带缟虎32.60%均较大。此外,瑞氏红鲂的杂合率为0.92%, 较绒杜父鱼0.26%、褐菖 鲉0.17%、艾氏蛇鳗0.70%、 许氏平 鲉0.22%、朝鲜平 鲉0.20%、 金斑平 鲉0.31%、双带缟虎0.47%均较大, 一般来说, 杂合率大于0.8%, 重复率大于60%就称为复杂基因组(高胜寒等, 2018), 由于瑞氏红鲂的杂合率高于0.8%, 重复率低于60%, 其基因组不属于复杂基因组类型。在本研究中, 组装的初步结果显示Contigs 的N50 为712 bp, N50 数量为254 317, 其结果满足基因组组装要求。通过过滤后的clean reads 组装瑞氏红鲂的基因组, 这是目前第一个基于二代测序技术组装瑞氏红鲂的基因组, 该基因组为瑞氏红鲂的进化生物学研究提供了基础数据, 也为进一步探索该物种的基因组特征提供了参考。由于二代测序技术存在读长短、组装结果的连续性无法保证以及由于基因组的杂合导致过度组装或组装不彻底等问题, 因此本研究缺少对GC-depth和含量的相关分析, 后续开展三代测序以获得瑞氏黄鲂高质量的基因组。

本研究组装出来的线粒体基因组大小为16 527 bp,线粒体结构未出现基因重排现象, GC 含量为45.7%,出现明显的AT 偏倚, 其结果符合硬骨鱼类线粒体碱基组成偏好于碱基A 和T 的特点(Broughtonet al,2001; 蒙子宁等, 2004; Consuegraet al, 2015), 碱基G的含量为16.8%, 表现出明显的反G 偏倚(孟刚等,2021), 与大多数鱼类研究结果一致(丁少雄等,2006)。13 个PCGs 中,COI是以GTG 作为起始密码子, 其余均为ATG 起始密码子, 其中密码子ATG 是翻译效率最高的一个(Consuegraet al, 2015)。PCGs终止密码子中出现T--和TA-这类不完整的终止密码子, 是因为这些基因之后是一个编码在同一条链上的基因, 允许转录在没有完整密码子的情况下终止(Hechtet al, 2017)。不完全终止密码子的存在在鱼类的线粒体基因组中很常见, 可以通过在mRNA 加工过程中添加一个poly A 尾巴来完成(Ojalaet al, 1981;Liuet al, 2009)。PCGs 的AT-skew 和GC-skew 为负值(表6), 一般来说, AT 倾斜的幅度小于GC 倾斜, 而且在许多情况下, 没有统计学意义, 在这里, AT-skew 低于GC-skews (绝对值), 这符合传统的偏好(Yuet al,2019)。AT-skew 的最大值和GC-skew 的最小值均出现在ND6 中, 且该基因的AT/GC-skew 值波动较大,核苷酸偏斜可能是由于在复制和转录过程中突变压力和选择压力之间的平衡, 为基因复制提供了一个潜在的方向(McLeanet al, 1998; Touchonet al, 2008)。tRNA 基因结构中仅tRNA-Ser(TCG)基因缺失DHU臂的结构, 这是鱼类线粒体基因组的共同特征, DHU臂的缺失改变了tRNA-Ser(TCG)在线粒体基因中的识别潜力, 使其更容易被识别(Hardtet al, 1993)。瑞氏红鲂线粒体基因组的12S rRNA和16S rRNA基因定位于H 链tRNA-Phe和tRNA-Leu(UUR)基因之间,中间以tRNA-Val基因为间隔,12S rRNA基因比16S rRNA基因更保守(Satohet al, 2016)。当RSCU 值=1时, 表示密码子的使用频率与其他简并密码子没有差异; 当氨基酸的RSCU 值均不等于1, 说明每个氨基酸的使用都有不同程度的偏倚(周丹等, 2013; 孟乾等, 2020)。D-loop 区在tRNA-Pro和tRNA-Phe之间, 与其他大多数脊椎动物的排列一样(Weiet al, 2010; Liet al, 2019a)。

目前通过线粒体基因组构建系统发育树来确定物种的进化关系已经成为系统发育分析的常用方法,比如周志雄等(2020)通过线粒体基因组构建系统发育树得出东方 鲀属内系统发育关系较为复杂, 且存在野生环境下种间自然杂交现象; Miya 等(2015)通过过去15 年线粒体基因组发展应用做出的总结, 认为线粒体基因组分析和高通量测序技术会成为物种分类重要方法。本研究基于12 个PCGs 构建系统发育树,以豹鲂作为外群, 分析瑞氏红鲂的系统发育关系。ML 树和BI 树结果表明鲉科在ML 树中与外群遗传距离较近, 聚在一起, 与在BI 树中结果不一致,但是两个树结果都显示 鲉科与其他科遗传距离都较远。瑞氏红鲂与须叉吻鲂聚为一支, 作为黄鲂科鱼类, 与形态学研究结果一致(Kawai, 2008)。

4 结论

致谢 感谢赵宸枫在采集样品和鉴定方面的帮助及杨天燕老师的修改意见, 谨致谢忱。