2022, Vol. 29 
2. 中国水产科学研究院黑龙江水产研究所，淡水鱼类育种国家地方联合工程实验室，农业农村部淡水水产生物技术与遗传育种重点实验室，黑龙江省特殊生境鱼类种质特性与抗逆育种重点实验室，黑龙江 哈尔滨 150070
2. National & Local United Engineering Laboratory of Freshwater Fish Breeding; Key Laboratory of Freshwater Aquatic Biotechnology and Genetic Breeding, Ministry of Agriculture and Rural Affairs; Heilongjiang Province’s Key Laboratory of Fish Stress Resistance Breeding and Germplasm Characteristics on Special Habitats, Heilongjiang River Fisheries Research Institute, Chinese Academy of Fishery Sciences, Harbin 150070, China
瓦氏雅罗鱼(Leuciscus waleckii)隶属于鲤形目(Cypriniformes)、鲤科(Cyprinidae)、雅罗鱼亚科(Leuciscinae)、雅罗鱼属，是我国北方知名的土著鱼类，主要分布于黑龙江流域、辽河、黄河、内陆淡水及盐碱湖泊。瓦氏雅罗鱼具有极强的盐碱耐受能力，能够耐受内蒙古达里湖碱度53.57 mmol/L, pH 9.6的极端不良环境，并形成优势种[1]。据地质学和地理学考证，达里湖在更新世早期(11500~7600 CAL.BP)是一个淡水湖泊，和古辽河水系相连；到更新世中晚期(3450 CAL. BP至今)，由于气候干燥，降水减少，地质下陷等原因，演变为封闭的内陆盐碱湖泊[2]。而栖息于此的瓦氏雅罗鱼(碱水种)，被逐渐升高的碱度胁迫和选择，在不足1万年的时间里实现了对达里湖高碱度高pH环境的快速适应[3]。室内碳酸盐碱度急性毒理研究发现，瓦氏雅罗鱼碱水种的耐碱能力是淡水种的1.5倍[4]。群体遗传学研究显示，瓦氏雅罗鱼碱水种和淡水种出现了明显的群体遗传结构分化，推测地理隔离和局域适应(local adaptation)可能是造成二者分化的主要原因[1,5-6]。因此，瓦氏雅罗鱼独特的地质和演化历史为揭示硬骨鱼类在古环境中适应极端环境的分子机制研究提供了理想的研究模型。
全基因组高通量测序技术的发展，使从基因组尺度解析瓦氏雅罗鱼盐碱适应的遗传基础成为可能。Xu等[6]通过对瓦氏雅罗鱼达里湖和黑龙江种群基因组SNP多态性进行扫描和比较分析，发现了一系列与硬骨鱼适应极端碱性环境相关的候选基因，功能分析发现这些基因参与了酸碱调节、离子转运、含氮废物排泄等多种调控通路。笔者近期通过消除地理隔离和空间背景差异，聚焦碱水和淡水两种异质环境，基于全基因组差异的SNP位点筛选鉴定了21个受盐碱正向选择的候选基因[5]。这些研究表明达里湖种群已进化出独特的基因组适应性机制以应对极端盐碱环境。
相比于SNP，插入缺失(insertion-deletion, InDel)作为第二种最常见的基因组变异类型，在稳定性和多态性上的表现更加突出[7]。目前在动植物中这种分子标记已经开始应用于种群遗传结构分析以及通过全基因组关联研究(genome-wide association studies, GWAS)开展的分子辅助育种[8-10]，但在鱼类基因组中相关研究报道相对较少。
本研究选取瓦氏雅罗鱼达里湖和松花江种群为研究对象，利用全基因组重测序技术，筛选碱水种和淡水种基因组差异InDel突变位点，结合基因组选择清除分析和基因选择压力分析筛选受盐碱正向选择的InDel位点及其关联基因，并对候选InDel关联基因进行选择压力及蛋白质结构预测分析。旨在通过全基因重测序挖掘InDel变异位点，为瓦氏雅罗鱼适应极端盐碱环境提供全面的、潜在的遗传基础证据，也为鱼类育种新型分子标记的开发提供参考和依据。
1 材料与方法 1.1 样本来源及DNA提取瓦氏雅罗鱼分别采自内蒙古达里湖和松花江流域绥滨段，两个采样点的水质参数参照Chang等[1]和常玉梅等[11]的报道。每个群体各采集5尾样本，达里湖种群(碱水种，DL)平均体重和体长分别为(139.18±10.96) g和(20.90±0.68) cm; 松花江种群(淡水种，SH)平均体重和体长分别为(125.85± 3.67) g和(20.10±0.24) cm。鱼体经MS-222 (Sigma- Aldrich，美国)麻醉后用经肝素钠浸润过的1.5 mL注射器从尾静脉采集血液，放入液氮保存备用。血液基因组DNA使用DNeasy Blood and Tissue Kit (Qiagen，德国)试剂盒进行提取，1%的琼脂糖凝胶电泳检测DNA完整性，并使用Nanodrop 8000 (Thermo Scientific，美国)检测DNA浓度和纯度。检测合格的DNA置于−20 ℃保存备用。
1.2 DNA文库构建及高通量测序检验合格的DNA样品(>1 μg)通过CovarisS2随机打断成长度为350 bp的片段。采用TruSeq Library Construction Kit (Illumina，美国)进行建库，严格使用说明书推荐的试剂和耗材。DNA片段经末端修复、加ployA尾、加测序接头、纯化、PCR扩增等步骤完成整个文库制备过程。构建好的文库经质检后在Illumina Hiseq4000平台上进行PE150双末端测序。文库构建及高通量测序由北京诺禾致源科技股份有限公司合作完成。
1.3 测序数据质量评估原始测序数据经过FastX-Toolkit软件(v0.0.14) (http://hannonlab.cshl.edu/fastx_toolkit/)进行过滤，过滤流程如下：(1) 过滤掉总体质量偏低(碱基质量分数Q<20)的reads序列；(2) 过滤去除Q>20但有效碱基数含量小于50%的reads序列；(3) 过滤掉3ʹ端Q<10的碱基序列，确保测序碱基错误率小于0.1; (4) 过滤掉所有reads序列中的接头序列；(5) 过滤掉测序深度不足而无法识别的碱基占比大于10%的reads序列；(6) 过滤掉长度<20 bp的reads序列。
1.4 InDel calling和统计分析使用Burrows-Wheeler Alignment软件(BWA, v0.7.17)[12]以默认参数将过滤后的数据与瓦氏雅罗鱼参考基因组(NCBI: GCA_900092035.1)[6]进行比对，使用NGS QC工具包[13]对过滤reads进行进一步的质量控制；采用Samtools软件(v1.9)[14]中的rmdump去除重复序列，使用基因组GATK工具包(v.3.5)[15]对InDel周围的reads进行比对，其中RealignerTargetCreator查找需要进行精细比对的位置信息；InDelRealigner对初始比对文件进行精细比对；最后使用BaseRecalibrator进行碱基质量分数的重新校准。
1.5 差异InDel位点挖掘使用Samtools软件[14]mpileup鉴定InDel变异位点，过滤标准如下：(1) Q≥20; (2) 每个个体测序覆盖深度(QualByDepth)>3; (3) 序列比对质量分数均方根(RMSMappingQuality)>20。
使用SnpEff软件(v4.1)[16]根据瓦氏雅罗鱼参考基因组和注释文件创建本地数据库，对获得的高质量InDel突变位点进行位置和效应信息注释。完成注释后，首先收集达里湖和松花江种群基因组上所有位置坐标相同的InDel区域位点，其次与松花江种群个体间共享位点(同一种群个体间InDel位点变化一致)比较，筛选出达里湖种群差异共享位点及其对应基因。此外，笔者前期基于遗传分化指数(fixation index, Fst)和种群核苷酸多样性比值π ratio (SH/DL)进行了达里湖种群和松花江种群基因组上的选择信号检测，将Fst和π ratio (SH/DL)占前5%的部分(Fst≥0.375, π ratio≥4.160)作为基因组具有强选择性扫描信号的区域[5]。在本研究中，笔者以强选择扫描信号区域的Fst和π ratio (SH/DL)值为参考，过滤得到位于该区域的达里湖种群InDel差异共享位点及其候选基因。
1.6 基因功能富集分析将候选基因递交至KOBAS 3.0[17] (http://kobas. cbi.pku.edu.cn/)进行GO和KEGG富集分析[18-19]。在Fisher精确检验的基础上，采用Benjamini- Hochberg多重测试对P值进行调整，设定0.05为显著性阈值。
1.7 进化选择压力分析根据富集分析结果，选择参与离子转运和酸碱平衡关键调控过程的候选基因mRNA序列进行选择压力分析。首先从雅罗鱼基因组数据库(http://www.fishbrowser.org/database/amur-ide/)中调取达里湖种群候选基因mRNA序列，然后检索GenBank数据库(www.ncbi.nlm.nih.gov/genbank)收集底鳉(Fundulus heteroclitus)、尼罗罗非鱼(Oreochromis niloticus)、大西洋鲑(Salmo salar)、大黄鱼(Larimichthys crocea)、斑马鱼(Danio rerio)、鲤(Cyprinus carpio)、鲫(Carassius auratus)等12种硬骨鱼对应基因mRNA序列(表1)。将每个物种所得序列与达里湖种群序列进行Blast比对以确保序列相似性。通过MEGA X软件ClustalW程序[20]以默认比对参数进行多序列比对，比对结果利用TBtools软件[21]进行修剪后通过最大似然法(maxmium likelihood, ML)构建系统发育树确定各物种间的系统发育关系，自验值设定为1000以验证进化树的可信程度。使用EasyCodeml软件(v4.8)[22]中的CodeML算法估算候选基因的非同义/同义替换率比率dN/dS (ω)，系统发育树(.tree)和转换为.phy的多序列比对文件作为输入数据。基于该软件中的位点模型检验候选基因的正向选择假设，该模型允许选择压力随位点的不同而变化，即假设系统发育树的不同分支所受选择压力相同，但不同氨基酸位点经历的选择压力不同。本研究考虑了4对替换模型，M0 (one ratio) vs M3 (discrete), M1a (Nearly Neutral) vs M2a (Positive Selection), M7 (β) vs M8 (β and ω>1)以及M8 (β and ω>1) vs M8a (β and ω=1)，与模型M1a、M0、M7相对应的M2a、M3、M8s假定条件均允许ω>1。使用似然比检验(LRT)比较了模型的拟合度，以模型中的Bayes Empirical Bayes (BEB)标准来确定正向选择压力位点，若后验概率≥0.95，表明该位点受到显著的正向选择压力。此外，为避免假阳性正向选择位点的出现，本研究还考虑了同义率变化和重组的其他测试，包括混合效应模型(MEME)和固定效应模型(FEL)。两种测试方法通过Datamonkey在线服务器(http:// www.datamonkey.org)完成，满足以下条件的位点鉴定为正向选择位点：(1) β+>α; (2) 似然比检验P<0.1。通过EasyCodeML、MEME和FEL其中两种或两种以上测试的位点被认为是强正向选择压力位点。
1.8 正向选择位点功能预测为了确定受正向选择氨基酸位点的功能，将瓦氏雅罗鱼候选基因氨基酸序列使用HMMER[23]与斑马鱼序列进行同源搜索；使用MEGAX软件的MUSCLE程序进行比对以确定瓦氏雅罗鱼正向选择位点在斑马鱼上的等效位置；使用Pfam在线分析系统[24]确定斑马鱼同源序列上的等效位点是否处于功能结构域或功能结构域位点周围。此外，使用蛋白空间结构分析服务器I-TASSER[25]通过折叠辨识模拟的方法预测候选基因的蛋白质三维结构，进一步了解基因编码蛋白的功能意义，模型评估通过C-score和TM-score完成。对于C-score，可信度区间为[−5,2]; 对于TM-score，大于0.5表明拓扑结构具有较高的准确性，而小于 0.17表明拓扑结构准确性较差[25-26]。最终的蛋白质三维结构模型通过PyMOL软件(https://www. schrodinger.com/pymol/)进行展示。
|
表1 不同物种mRNA序列登录号 Tab. 1 mRNA accession number of different species |
样品经过测序和质控共产生了3243686532个高质量clean reads，包含486.57 G的测序数据。其中DL组有1659962056 reads和249 G碱基数，SH组有1583724476 reads和237.57 G碱基数(表2)。与初始reads相比，保留下来的高质量clean reads占96.7%，其中，phred score>20的reads占97.3%, phred score>30的reads数占92.6%。GC含量较稳定，在39.30~39.88之间浮动。reads比对至参考基因组后，统计结果显示，DL组总的Mapping reads数为1364105229，平均比对率为82.17%，平均测序深度约为54.4 X; SH组总的Mapping reads数为1287038455，平均比对率为81.27%，平均测序深度约为51.34 X (表2)。
2.2 InDel检测与分布统计DL和SH两组所有个体共有983528个InDels被保留下来，在各个染色体上，分布的InDels数量为21244~67765 (图1)。总体上各个染色体上分布的InDels数量有所差异或偏倚，平均每条染色体上每Mb分布2.16~6.89个(图2)。在保留的983528个InDels中，DL组有420095个，SH组有563433个。InDels长度为1~29 bp，插入多于缺失，主要以1~8 bp的小片段为主。单碱基InDels最多，其中 DL 组插入数和缺失数分别为294851个和30515个，SH组分别为353277个和26539个。其次为双碱基InDels，其中DL组插入数和缺失数分别为18584个和13097个，SH组分别为30870个和21302个。随着InDels碱基数的增加，InDels数量整体呈下降趋势(图3)。
|
|
表2 瓦氏雅罗鱼达里湖和松花江种群个体测序数据及比对结果统计 Tab. 2 Summary of sequencing and mapping statistics of Leuciscus waleckii individuals from Lake Dali Nor (DL) and Songhua River (SH) populations |
02-0184-16-F001.jpg "点击查看原图") |
图1 瓦氏雅罗鱼每条染色体上的插入/缺失位点(InDels)数量分布 Fig. 1 The number of insertion-deletion sites (InDels) on each chromosome of Leuciscus waleckii |
DL组总InDels数量为420095个。大多数InDels位于基因间区域(328682个，78.24%)，少部分位于基因区域(76289个，18.16%)、基因上游(8906个，2.12%)和下游(6218个，1.48%)，其中基因区域包括内含子(66173, 86.74%)、外显子(6340, 8.31%)和非翻译区(3′-UTR, 2609 个，3.42%; 5′-UTR, 1169个，1.53%)。编码区内移码插入总数(frameshift insertion)为1818个(28.68%)，编码区内移码缺失总数(frameshift deletion)为2191个(34.56%)，编码区内InDels突变获得终止子总数(stop-gain)为186个(2.94%)，编码区内 InDels突变丢失终止子总数(stop-lost)为280个(4.41%); 非编码区内移码插入总数(non-frameshift insertion)和非编码区内移码缺失总数(non-frameshift deletion)分别为793个(12.50%)和1072个(16.91%) (图4a)。
02-0184-16-F002.jpg "点击查看原图") |
图2 插入/缺失位点(InDels)在瓦氏雅罗鱼各染色体的分布情况 Fig. 2 Insertion-deletion sites (InDels) distribution on the chromosomes of Leuciscus waleckii |
02-0184-16-F003.jpg "点击查看原图") |
图3 不同长度插入/缺失位点(InDels)在瓦氏雅罗鱼基因组的分布a. 瓦氏雅罗鱼达里湖种群不同长度InDels基因组分布；b. 瓦氏雅罗鱼松花江种群不同长度InDels基因组分布. Fig. 3 The distribution of different size insertion-deletion sites (InDels) in Leuciscus waleckii genomea. The distribution of InDels length in Leuciscus waleckii from Lake Dali Nor (DL) population; b. The distribution of InDels length in Leuciscus waleckii from Songhua River (SH) population. |
SH组总InDel数量为563433个。与DL组分布相似，大多数InDels位于基因间区域(430068个，76.33%)，少部分位于基因区域(120743个，21.43%)、基因上游(7043个，1.25%)和下游(5579个，0.99%)，其中基因区域包括内含子(93709个，77.61%)、外显子(20164个，16.70%)和非翻译区(3′-UTR, 3513个，2.91%; 5′-UTR, 3357个，2.78%)。编码区内移码插入总数为5767个(28.60%)，编码区内移码缺失总数为6769个(33.57%)，编码区内 InDels突变获得终止子总数为645个(3.20%)，编码区内 InDel 突变丢失终止子总数为250个(1.24%); 非编码区内移码插入总数和非编码区内移码缺失总数分别为3438个(17.05%)和3295个(16.34%)(图4b)。
02-0184-16-F004.jpg "点击查看原图") |
图4 瓦氏雅罗鱼不同种群插入/缺失位点(InDels)区域功能分类a. 瓦氏雅罗鱼达里湖种群InDels区域功能分类；b. 瓦氏雅罗鱼松花江种群InDels区域功能分类. Fig. 4 Functional classification of the detected InDels for different Leuciscus waleckii population groupsa. Functional classification of the detected InDels for Leuciscus waleckii from Lake Dali Nor (DL); b. Functional classification of the detected InDels for Leuciscus waleckii from Songhua River (SH). |
DL组与SH组共筛选得到8176个差异共享InDels位点，其中基因间区域为5755个(70.39%)，基因上游和下游分别为60个(0.74%)和55个(0.68%)，基因区域中非翻译区3′-UTR为23个(0.29%), 5′-UTR为22个(0.27%)，内含子为2074个(25.37%)，外显子为187个(2.26%)。通过计算Fst和π ratio (SH/DL)在两个种群中共筛选到3259个强选择信号区域[5]，将DL差异共享位点与强选择信号区域取交集，共得到325个差异共享位点关联到176个基因，将这些基因确定为最终的功能候选基因。
2.5 基因功能富集分析176个候选基因共富集到1613个GO条目和181个KEGG 通路中，其中显著富集的GO条目有193个，KEGG通路有27个(P<0.05)。GO分析显示生物过程(biological process, BP)相关的候选基因富集程度较高的为免疫系统，包括凋亡正调控(GO: 0043065)、免疫应答(GO: 0045087)、白细胞介素-2合成正调控(GO: 0032743)等过程；其次为离子跨膜转运(GO: 0034220)和碱性氨基酸转运体活性(GO: 0015174)相关过程的候选基因也存在显著富集。细胞组分(cellular component, CC)中与能量代谢有关的候选基因存在显著富集，包括线粒体(GO: 0005739)、线粒体内膜(GO: 000743)、呼吸链复合物IV (GO: 0045277)等过程。分子功能(molecular function, MF)中转运调节基因富集程度较高，包括Na+跨膜转运活性(GO: 0015081)和ABC型转运蛋白活性(GO: 0140359)等(图5)。KEGG分析显示参与细胞凋亡与免疫应答相关的功能区域显著富集，包括NOD样受体通路(04621)、Toll样受体通路(04620)、自噬(04215)、TNF信号通路(04668)和p53信号通路(04115)等；此外，可能与逆境生境适应相关的功能区域也有显著富集，包括Ca2+信号通路(04020)、心肌收缩通路(04260)、肾素-血管紧张素系统(04614)以及应激激素区域，如催乳素通路(04917)等(图6)。
02-0184-16-F005.jpg "点击查看原图") |
图5 达里湖瓦氏雅罗鱼种群盐碱适应候选基因GO分析 Fig. 5 GO annotation of candidate genes for alkali-saline adaptation in Leuciscus waleckii from Lake Dali Nor population |
02-0184-16-F006.jpg "点击查看原图") |
图6 达里湖瓦氏雅罗鱼种群盐碱适应候选基因KEGG分析 Fig. 6 KEGG pathway analysis of candidate genes for alkali-saline adaptation in Leuciscus waleckii from Lake Dali Nor populations |
基于种群强选择信号区域、InDels差异共享位点、富集的GO条目类型(图5)以及基因生物学功能，abcc1、atp2b1、slc4a4、slc7a2、aqp4被确定为与DL组离子转运和酸碱平衡等过程相关的关键调控基因，这5个基因都处于强选择信号区域中，并且都至少包含有一个差异共享InDel突变位点(表3)。将这些关键调控基因首先利用EasyCodeML软件进行选择压力分析。对于aqp4，在M0模型下，ω=0.39436, ω值远小于1，没有直接的证据表明存在正向选择(表4)，比较M1a和M2a，似然比检验2ΔlnL=0, P值为1，说明M2a并不优于M1a，分别比较M7 vs M8和M8a vs M8，似然比检验分别为2ΔlnL=9.44 (P<0.01)和2ΔlnL=0.015 (P>0.5)，说明M8a模型并不优于M8，综合表明aqp4基因主要受到较强的净化选择约束，具有很高的保守性。对于基因atp2b1，在模型M0中，进化速率ω0=1.10315，说明atp2b1在进化过程中受到了一定的正向选择，对4个位点模型进行似然比检验发现，M3、M2a和M8模型明显优于其相应的假设模型M0、M1a、M7与M8a (P< 0.01)(表4)，表明各点之间的选择压力存在差异。其余基因slc4a4、slc7a2、abcc1的M0模型表明它们在进化过程中都受到了不同程度的净化选择(ω<1)，但M3、M2a和M8模型明显优于其相应的假设模型M0、M1a、M7与M8a (P<0.01)，因此尽管在进化过程中受到不同程度的净化选择，但在这些基因M8模型中也存在着发生正向选择的氨基酸位点。在模型M8中，分别以EasyCodeML、MEME和FEL计算检测了abcc1、atp2b1、slc4a4、slc7a2和aqp4 5个候选基因的正向选择位点，分别检测到2、16、16、8和0个强正向选择位点，其中abcc1的两个位点达到极显著水平，其余的位点均为显著水平(表5)。
|
|
表3 达里湖瓦氏雅罗鱼种群盐碱适应关键调控基因插入/缺失位点信息 Tab. 3 The detailed information of the key regulatory genes and the related InDel loci for alkali-saline adaptation in Leuciscus waleckii from Lake Dali Nor (DL) populations |
|
|
表4 达里湖瓦氏雅罗鱼种群盐碱适应关键调控基因适应性进化分析 Tab. 4 Adaptive evolution analysis of key regulatory genes for alkali-saline adaptation in Leuciscus waleckii from Lake Dali Nor (DL) populations |
将上述正向选择位点所在氨基酸序列与斑马鱼同源序列比对后确定了瓦氏雅罗鱼在斑马鱼正向选择位点的等效位置。abcc1的两个氨基酸位点虽然达到了极显著水平，但未发现与斑马鱼相关功能结构域及其附近位点的重合区域。atp2b1受选择位点主要集中在阳离子转运ATP酶功能域的两端(N端和C端)以及E1-E2ATP酶功能域。slc4a4受选择位点与斑马鱼功能重叠区域主要集中在带3蛋白胞质段和HCO3–共转运结构域。slc7a2受选择位点主要集中在斑马鱼氨基酸通透酶结构域(图7)。由于slc4a4和slc7a2的部分受选择位点处于共转运和跨膜载体蛋白等关键调控结构域中，因此本研究将这两个基因的氨基酸序列提交至I-TASSER进行三维蛋白结构预测。在SLC4A4蛋白结构中，检测到的正向选择氨基酸位点660th (Phe)与配体结合位点重合，144th (Phe)位点非常接近配体结合位点。在SLC7A2蛋白结构中，正向选择位点45th (Leu)接近配体结合位点，这些位点可能在配体结合中起关键作用，此外，在两个蛋白结构预测中均未发现酶激活位点(图8)。
|
|
表5 正向选择氨基酸位点 Tab. 5 Amino acid sites undergoing positive selection |
02-0184-16-F007.jpg "点击查看原图") |
图7 基于瓦氏雅罗鱼正向选择位点预测斑马鱼功能结构域等效位置 Fig. 7 Prediction of equivalent sites of Danio rerio’s functional domains based on positive selection sites of Leuciscus waleckii |
02-0184-16-F008.jpg "点击查看原图") |
图8 SLC4A4和SLC7A2三维蛋白结构及正向选择位点所处位置 Fig. 8 Positive selection sites and three dimensional structure of SLC4A4 and SLC7A2 |
本研究借助于高通量测序技术，筛选鉴定了瓦氏雅罗鱼两个生态种群的InDels位点，得到的InDels数量与类型基本一致，说明两个种群基因组上的InDels分布均匀，没有因为起源背景或生活环境差异呈现偏态。总体来看，InDels插入数大于缺失数，数量集中在1~8 bp，其中1~3 bp最多。此外，两个种群InDels注释变体大部分分布于基因间区和内含子等非编码区域，编码区中分布的位点最少，这与人类[27]、荷斯坦牛[28]、藏酋猴[29]基因组数据上观察到的结果类似。
在InDel注释结果的基础上，通过每个种群个体间共享InDels位点信息筛选得到瓦氏雅罗鱼达里湖种群差异共享InDels位点，通过与前期基于SNP位点筛选得到的差异基因组区域整合[5]，确定了在基因组上一些可能与盐碱适应相关的候选基因。生物体对非生物胁迫的应答通常涉及多种机制[30]，筛选得到的候选基因在功能富集分析中显著富集于几个不同的功能区域，这与瓦氏雅罗鱼盐碱适应机制复杂、多样、系统性强的特点相符[11]。
通路富集分析发现与离子转运相关的基因显著富集，其编码蛋白在细胞膜的流动性和渗透性方面发挥重要作用。如受正向选择压力较强的slc7a2和slc4a4。研究显示，SLC7A2作为渗透酶参与跨质膜转运和吸收阳离子氨基酸，这种不依赖于 Na+的渗透调节系统被证实可通过氨基酸的有效转运使哺乳动物能够在胁迫环境下维持细胞的正常功能[31]。NBC1是基因slc4a4的表达产物，属于Na+/HCO3–共转运蛋白，是重要的HCO3–转运家族成员之一，对维持鱼类血液酸碱平衡发挥重要作用[32-33]。相比于淡水种，达里湖种群受到了包括总碱度、pH值以及氨氮代谢变化产生的协同毒性[11]，使得血液中HCO3–浓度增加，破坏体内酸碱平衡。课题组前期研究发现，在相同碱度胁迫下，随着胁迫时间的延长，NBC1 mRNA在达里湖种群中显著高表达，而在松花江种群表达量显著降低[34]，表明达里湖瓦氏雅罗鱼通过积极主动地排碱保持体内离子和酸碱平衡。关键调控基因atp2b1编码质膜钙泵蛋白，能够将Ca2+从细胞质转运到细胞外来维持跨质膜的钙浓度梯度，更重要的是，质膜钙泵与渗透调节过程密切相关[35]，一项尼罗罗非鱼盐度耐受实验表明，在0~32的盐度递增条件下尼罗罗非鱼鳃和脑组织atp2b1和Na+/K+ ATPase (NKA) mRNA均显著高表达，从而实现在高渗环境下Na+/Cl−转运和Ca2+内流[36]。课题组前期关于达里湖瓦氏雅罗鱼的碱度耐受实验同样发现NKA随着碱度耐受时间的延长呈现显著高表达[34]，因此ATP2B1参与的Ca2+泵激活可能与NKA的表达变化存在密切关系。此外，水体pH的升高将导致H+浓度的降低，使得鱼体排氨(NH3+H+=NH4+)受阻，最终导致“氨中毒”[37]，本研究鉴定得到的关键调控基因abcc1的表达产物是一种转运谷胱甘肽结合物的膜转运蛋白[38]，在欧洲舌齿鲈(Dicentrarchus labrax)和尼罗罗非鱼中，ABCC1定位在细胞底端，将内容物运输至血液，对排出细胞内积累的环境毒物以及组织防御等过程至关重要[39-41]。达里湖瓦氏雅罗鱼能够对抗PNH3梯度继续主动排氨[11]，推测ABCC1对于防止细胞内氨中毒发挥重要作用，但细胞解毒机制以及与氨氮代谢之间的关系仍需进一步探索。
通路分析还筛选到与免疫过程相关的基因显著富集，如xiap、ube2l、ube2d2、psme4b属于泛素-蛋白酶体系统，这一系统提供了控制蛋白质质量的基本方式，即增加蛋白的水解作用，加快蛋白质的更新速率，从而有选择性地清除在极端环境下由于生物合成错误、氧自由基或变性而引起的异常折叠或损坏[42]，表明修复细胞和清除坏死的细胞器也是瓦氏雅罗鱼盐碱适应的重要机制。
除了离子交换、转运及免疫的基因显著富集外，本研究还发现部分候选基因富集到了钙离子信号通路，包括编码钙黏蛋白(cdh4, cdh7, cdh11, cdhr1)，紧密连接蛋白(cldn1)和细胞骨架蛋白(camsap1)相关的基因。其中cldn1在生物体处于压力环境时能够通过维持上皮细胞的完整性来充当物理屏障，防止溶质和水自由通过上皮细胞和内皮细胞之间的细胞间隙。camsap1编码细胞骨架蛋白，当细胞暴露于应激环境时，细胞可通过快速重组各种细胞骨架网络从而做出应答[43]。
适应性进化(正向选择)是通过环境变化对蛋白质施加选择压力，以提高生物体在局域环境中的适应性[44]。本研究推测达里湖瓦氏雅罗鱼在盐碱适应的进化过程中，可能通过调控基因的SNP或InDel等序列多态性，促使基因编码或调控序列的功能发生变化，以此达到对极端环境的适应。因此，除了对盐碱适应相关基因的功能研究外，本研究还重点关注了这些抗性相关基因在进化过程中受到的自然选择压力，对参与离子转运和酸碱平衡过程的关键调控基因进化速率进行了分析。笔者在选择压力检测分析中首先采用了EasyCodeML 中的位点模型，因为在自然界的进化过程中，选择压力可能只是作用于某些位点，而不是作用于整个基因中[45]。研究结果显示，slc4a4、slc7a2、abcc1虽然主要经历了净化选择，但仍检测到小部分的正向选择氨基酸位点，并且slc4a4和slc7a2部分正向选择位点处于HCO3–共转运蛋白和氨基酸渗透关键功能结构域中。推测虽然整条基因在进化过程中受到较强的净化选择，但某些功能结构域的正向选择位点可能与盐碱胁迫环境的适应性进化、基因结构的改变有着重要的联系，盐碱环境胁迫可能是进化的主要驱动力。
4 结论本研究对瓦氏雅罗鱼达里湖种群和松花江种群各5 尾个体进行了基因组重测序，共鉴定得到983528个InDels位点。通过比较两个种群之间的差异InDels位点并整合已知的基因组强选择信号区域，筛选得到与盐碱适应相关的325个InDels位点及其关联候选基因176个。富集分析显示候选基因主要在离子转运、酸碱平衡及炎症免疫反应等生物学过程发挥重要作用。进一步的选择压力分析显示，关键调控基因abcc1、atp2b1、slc4a4、slc7a2和aqp4在进化过程中净化选择均起主导作用，但处于关键功能结构域的正向选择状态位点可能有助于瓦氏雅罗鱼达里湖种群快速适应极端盐碱环境。本研究结果将为瓦氏雅罗鱼盐碱适应相关基因或标记的开发提供重要参考。
| [1] |
Chang Y M, Tang R, Sun X W, et al. Genetic analysis of population differentiation and adaptation in Leuciscus waleckii[J]. Genetica, 2013, 141(10-12): 417-429..》Google Scholar
|
| [2] |
Geng K, Zhang Z C. The geomorphic characteristics and evolution of the lakes in Dalairuoer area of Neimenggu Plateau during the Holocene[J]. Journal of Beijing Normal University (Natural Science), 1988, 24(4): 94-101. [耿侃，张振春. 内蒙古达来诺尔地区全新世湖群地貌特征及其演化[J]. 北京师范大学学报(自然科学版), 1988, 24(4): 94-101.].》Google Scholar
|
| [3] |
Chi B J, Chang Y M, Yan X C, et al. Genetic variability and genetic structure of Leuciscus waleckii Dybowski in Wusuli River and Dali Lake[J]. Journal of Fishery Sciences of China, 2010, 17(2): 228-235. [池炳杰，常玉梅，闫学春，等. 瓦氏雅罗鱼达里湖群体和乌苏里江群体的遗传多样性和遗传结构分析[J]. 中国水产科学，2010, 17(2): 228-235.].》Google Scholar
|
| [4] |
Chi B J. Genetic variability of Leuciscus waleckii (Dybowski) and construction of the cDNA substractive library of saline-alkaline tolerance[D]. Shanghai: Shanghai Ocean University, 2010: 12-17. [池炳杰. 瓦氏雅罗鱼遗传多样性分析和耐盐碱cDNA消减文库的构建[D]. 上海：上海海洋大学，2010: 12-17.].》Google Scholar
|
| [5] |
Wang S Y, Kuang Y Y, Liang L Q, et al. Resequencing and SNP discovery of Amur ide (Leuciscus waleckii) provides insights into local adaptations to extreme environments[J]. Scientific Reports, 2021, 11: 5064..》Google Scholar
|
| [6] |
Xu J, Li J T, Jiang Y L, et al. Genomic basis of adaptive evolution: The survival of Amur ide (Leuciscus waleckii) in an extremely alkaline environment[J]. Molecular Biology and Evolution, 2017, 34(1): 145-159..》Google Scholar
|
| [7] |
Adedze Y M N, Lu X, Xia Y C, et al. Agarose-resolvable InDel markers based on whole genome re-sequencing in cucumber[J]. Scientific Reports, 2021, 11: 3872..》Google Scholar
|
| [8] |
Misra G, Badoni S, Anacleto R, et al. Whole genome sequencing-based association study to unravel genetic architecture of cooked grain width and length traits in rice[J]. Scientific Reports, 2017, 7: 12478..》Google Scholar
|
| [9] |
Grobet L, Royo Martin L J, Poncelet D, et al. A deletion in the bovine myostatin gene causes the double-muscled phenotype in cattle[J]. Nature Genetics, 1997, 17(1): 71-74..》Google Scholar
|
| [10] |
Chen M Y, Yan H L, Wang K, et al. Goat SPEF2: Expression profile, indel variants identification and association analysis with litter size[J]. Theriogenology, 2019, 139: 147-155..》Google Scholar
|
| [11] |
Chang Y M, Liang L Q. Advances of research of physiological and molecular mechanisms related to alkali-saline adaptation for fish species inhabiting alkali-saline water[J]. Journal of Fisheries of China, 2021, 45(5): 798-812. [常玉梅，梁利群. 耐盐碱鱼类的生理和分子机制研究进展[J]. 水产学报，2021, 45(5): 798-812.].》Google Scholar
|
| [12] |
Langmead B, Trapnell C, Pop M, et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[J]. Genome Biology, 2009, 10(3): R25..》Google Scholar
|
| [13] |
Patel R K, Jain M. NGS QC Toolkit: A toolkit for quality control of next generation sequencing data[J]. PLoS One, 2012, 7(2): e30619..》Google Scholar
|
| [14] |
Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data[J]. Bioinformatics, 2011, 27(21): 2987-2993..》Google Scholar
|
| [15] |
McKenna A, Hanna M, Banks E, et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next- generation DNA sequencing data[J]. Genome Research, 2010, 20(9): 1297-1303..》Google Scholar
|
| [16] |
Cingolani P, Platts A, Wang L L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff[J]. Fly, 2012, 6(2): 80-92..》Google Scholar
|
| [17] |
Xie C, Mao X Z, Huang J J, et al. KOBAS 2.0: A web server for annotation and identification of enriched pathways and diseases[J]. Nucleic Acids Research, 2011, 39(suppl_2): W316-W322..》Google Scholar
|
| [18] |
Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes[J]. Nucleic Acids Research, 2000, 28(1): 27-30..》Google Scholar
|
| [19] |
Kanehisa M. Toward understanding the origin and evolution of cellular organisms[J]. Protein Science, 2019, 28(11): 1947-1951..》Google Scholar
|
| [20] |
Kumar S, Stecher G, Li M, et al. MEGA X: Molecular evolutionary genetics analysis across computing platforms[J]. Molecular Biology and Evolution, 2018, 35(6): 1547-1549..》Google Scholar
|
| [21] |
Chen C J, Chen H, Zhang Y, et al. TBtools: An integrative toolkit developed for interactive analyses of big biological data[J]. Molecular Plant, 2020, 13(8): 1194-1202..》Google Scholar
|
| [22] |
Gao F L, Chen C J, Arab D A, et al. EasyCodeML: A visual tool for analysis of selection using CodeML[J]. Ecology and Evolution, 2019, 9(7): 3891-3898..》Google Scholar
|
| [23] |
Finn R D, Clements J, Eddy S R. HMMER web server: Interactive sequence similarity searching[J]. Nucleic Acids Research, 2011, 39(suppl_2): W29-W37..》Google Scholar
|
| [24] |
Mistry J, Chuguransky S, Williams L, et al. Pfam: The protein families database in 2021[J]. Nucleic Acids Research, 2021, 49(D1): D412-D419..》Google Scholar
|
| [25] |
Yang J Y, Zhang Y. I-TASSER server: New development for protein structure and function predictions[J]. Nucleic Acids Research, 2015, 43(W1): W174-W181..》Google Scholar
|
| [26] |
Zhang Y. I-TASSER server for protein 3D structure prediction[J]. BMC Bioinformatics, 2008, 9: 40..》Google Scholar
|
| [27] |
Fujimoto A, Nakagawa H, Hosono N, et al. Whole-genome sequencing and comprehensive variant analysis of a Japanese individual using massively parallel sequencing[J]. Nature Genetics, 2010, 42(11): 931-936..》Google Scholar
|
| [28] |
Jiang J P, Gao Y H, Hou Y L, et al. Whole-genome resequencing of Holstein bulls for indel discovery and identification of genes associated with milk composition traits in dairy cattle[J]. PLoS One, 2016, 11(12): e0168946..》Google Scholar
|
| [29] |
Fan Z X, Zhao G, Li P, et al. Whole-genome sequencing of Tibetan macaque (Macaca thibetana) provides new insight into the macaque evolutionary history[J]. Molecular Biology and Evolution, 2014, 31(6): 1475-1489..》Google Scholar
|
| [30] |
White C R, Alton L A, Frappell P B. Metabolic cold adaptation in fishes occurs at the level of whole animal, mitochondria and enzyme[J]. Proceedings of the Royal Society B: Biological Sciences, 2012, 279(1734): 1740-1747..》Google Scholar
|
| [31] |
Ellingsen S, Narawane S, Fjose A, et al. Sequence analysis and spatiotemporal developmental distribution of the Cat-1-type transporter slc7a1a in zebrafish (Danio rerio)[J]. Fish Physiology and Biochemistry, 2020, 46(6): 2281-2298..》Google Scholar
|
| [32] |
Esbaugh A J, Heuer R, Grosell M. Impacts of ocean acidification on respiratory gas exchange and acid-base balance in a marine teleost, Opsanus beta[J]. Journal of Comparative Physiology B, 2012, 182(7): 921-934..》Google Scholar
|
| [33] |
Lee Y C, Yan J J, Cruz S A, et al. Anion exchanger 1b, but not sodium-bicarbonate cotransporter 1b, plays a role in transport functions of zebrafish H+-ATPase-rich cells[J]. American Journal of Physiology Cell Physiology, 2011, 300(2): C295-C307..》Google Scholar
|
| [34] |
Chang Y M, Zhao X F, Liew H J, et al. Effects of Bicarbonate Stress on Serum Ions and Gill Transporters in Alkali and Freshwater Forms of Amur Ide (Leuciscus waleckii)[J]. Frontiers in physiology, 2021, 12: 676096..》Google Scholar
|
| [35] |
Verbost P M, Schoenmakers T J, Flik G, et al. Kinetics of ATP- and Na+-gradient driven Ca2+ transport in basolateral membranes from gills of freshwater- and seawater-adapted tilapia[J]. The Journal of Experimental Biology, 1994, 186: 95-108..》Google Scholar
|
| [36] |
Rengmark A H, Slettan A, Lee W J, et al. Identification and mapping of genes associated with salt tolerance in tilapia[J]. Journal of Fish Biology, 2007, 71(sc): 409-422..》Google Scholar
|
| [37] |
Chen S G, Hou J L, Yao N, et al. Comparative transcriptome analysis of Triplophysa yarkandensis in response to salinity and alkalinity stress[J]. Comparative Biochemistry and Physiology Part D: Genomics and Proteomics, 2020, 33: 100629..》Google Scholar
|
| [38] |
Keppler D, Cui Y, König J, et al. Export pumps for anionic conjugates encoded by MRP genes[J]. Advances in Enzyme Regulation, 1999, 39: 237-246..》Google Scholar
|
| [39] |
Ferreira M, Costa J, Reis-Henriques M A. ABC transporters in fish species: A review[J]. Frontiers in Physiology, 2014, 5: 266..》Google Scholar
|
| [40] |
Costa J, Reis-Henriques M A, Wilson J M, et al. P-glycoprotein and CYP1A protein expression patterns in Nile tilapia (Oreochromis niloticus) tissues after waterborne exposure to benzo(a)pyrene (BaP)[J]. Environmental Toxicology and Pharmacology, 2013, 36(2): 611-625..》Google Scholar
|
| [41] |
Ferreira M, Santos P, Rey-Salgueiro L, et al. The first demonstration of CYP1A and the ABC protein(s) gene expression and activity in European seabass (Dicentrarchus labrax) primary hepatocytes[J]. Chemosphere, 2014, 100: 152-159..》Google Scholar
|
| [42] |
Benedito-Palos L, Ballester-Lozano G F, Simó P, et al. Lasting effects of butyrate and low FM/FO diets on growth performance, blood haematology/biochemistry and molecular growth-related markers in gilthead sea bream (Sparus aurata)[J]. Aquaculture, 2016, 454: 8-18..》Google Scholar
|
| [43] |
Kim S, Coulombe P A. Emerging role for the cytoskeleton as an organizer and regulator of translation[J]. Nature Reviews Molecular Cell Biology, 2010, 11(1): 75-81..》Google Scholar
|
| [44] |
Morgan C C, Loughran N B, Walsh T A, et al. Positive selection neighboring functionally essential sites and disease-implicated regions of mammalian reproductive proteins[J]. BMC Evolutionary Biology, 2010, 10: 39..》Google Scholar
|
| [45] |
Kamath P L, Getz W M. Adaptive molecular evolution of the Major Histocompatibility Complex genes, DRA and DQA, in the genus Equus[J]. BMC Evolutionary Biology, 2011, 11: 128..》Google Scholar
|