军曹鱼响应低氧胁迫转录组SNP位点鉴定及其功能注释分析

杨二军，杨林桐，王维政，黄建盛,2,3*，张健东,2,3，王忠良,2,3，陈刚,2,3*

( 1. 广东海洋大学 水产学院，广东 湛江 524088；2. 南方海洋科学与工程广东省实验室（湛江），广东 湛江 524025；3. 广东省水产经济动物病原生物学及流行病学重点实验室，广东 湛江 524088)

1 引言

溶解氧（Dissolved Oxygen, DO）对水生动物的发育和健康至关重要[1]，同时也是水生动物生存不可缺少的环境因子之一。溶解氧不足时，机体的生长发育、呼吸、新陈代谢、繁殖及免疫功能等生化反应和生理功能会受到影响，严重时会发生死亡[2]。随着全球气候变暖和水体富营养化引发的一系列生态问题，水体低氧化现象越来越严重。调查发现，我国近岸海域出现的低氧层（DO浓度约为3 mg/L）水域范围越来越大，对生活在其中的海洋生物造成了极大的威胁，低氧现象已成为一个被高度重视和广泛关注的问题[3-4]。

军曹鱼(Rachycentron canadum)，隶属于鲈形目(Periformes)，军曹鱼科(Rachycentridae)，军曹鱼属(Rachycentron)，又名海鲡、海龙鱼等，是热带和亚热带海域的广盐性、肉食性洄游鱼类。因其具有生长快、抗病能力强、营养丰富、肉质鲜美等特点，使其成为我国南方沿海地区深水网箱养殖的重要品种之一[5]。近年来，随着高密度集约化养殖模式的推广，军曹鱼近海浮动式网箱养殖和深海网箱养殖发展迅速，然而受养殖密度大、管理不善、养殖生态环境恶化、天气、温度等因素的影响，军曹鱼在养殖过程中容易出现缺氧、抗病能力下降等情况，导致病害日趋频繁，严重制约了军曹鱼的健康养殖发展[6-7]。针对军曹鱼养殖过程中出现的缺氧问题，本实验室前期开展了低氧胁迫对军曹鱼幼鱼生长和血清生化指标、氧化应激、能量代谢以及肠道菌群影响的研究，结果表明，低氧胁迫对军曹鱼幼鱼的生理状况有显著影响[7-12]。因此，开展军曹鱼的遗传改良和优良性状品种培育的工作，对促进我国军曹鱼养殖业健康快速的发展具有现实意义。

单核苷酸多态性（Single Nucleotide Polymorphism，SNP）是指因单个核苷酸变异引起的DNA序列多态性，其形式包括单个碱基的转换、颠换、插入和缺失，是一种常见的基因突变[13]。SNP作为第三代分子标记方法，是基因型和表型之间最好的研究载体，相比于前几代的标记方法，它具有位点丰富、分布广泛、遗传稳定性强、检测快速、自动化分析便利和成本低等特点，已成为最具前景的遗传标记方法之一，被广泛应用于生物、农学、医学等众多领域[14-15]。近年来，SNP分子标记技术在鱼类[16-17]、甲壳类[18-19]和贝类[20-21]等水产动物的研究中也取得了重要进展。

随着高通量测序技术的快速发展以及测序成本的下降，利用转录组测序筛选目标性状SNP位点已成为研究遗传多样性、分子育种等分子标记的常用手段[22]。研究发现，利用转录组测序技术和物种基因组信息，更易获得与目标性状相关的SNP[23]。目前，通过转录组测序来筛选大量SNP位点的研究已在多种水产动物中被报道，如马氏珠母贝(Pinctada fucata)[20]、方斑东风螺(Babylonia areolata)[23]、罗氏沼虾(Macrobrachium rosenbergii)[24]、团头鲂(Megalobrama amblyocephala)[25]、鳙鱼(Hypophthalmichthys nobilis)[26]等。关于军曹鱼群体遗传的研究资料较为匮乏，曹丹煜[27]曾利用转录组测序结果筛选与军曹鱼盐度适应相关SNP位点，但未对其响应低氧胁迫SNP位点挖掘和功能注释进行研究。本研究基于低氧胁迫后的军曹鱼幼鱼肠道转录组测序结果，筛选出大量SNP位点，对SNP所在基因的SNP-Unigene进行功能注释，侧重对免疫防御及低氧通路中的相关基因进行SNP的挖掘与筛选，以期为今后军曹鱼耐低氧和抗病品系的遗传选育、分子辅助育种研究提供基础资料。

2 材料和方法

2.1 实验用鱼及饲养

军曹鱼幼鱼由广东海洋大学湛江东海岛生物研究基地提供。实验在广东恒兴饲料实业股份有限公司863基地进行，养殖设施为室内24 h持续充气的流水养殖系统，水体为500 L。选取大小规格整齐、活力好、体格健壮的军曹鱼210尾（体重为（50.44±2.78）g），先在养殖水槽中暂养1周使其适应新环境后再进行低氧胁迫实验，每天投喂2次，及时用虹吸管清理粪便和剩余饵料，暂养期间水体溶解氧在6 mg/L以上，水温为（29±1）℃，pH为7.8～8.0，盐度为28～30，氨氮浓度小于0.02 mg/L，自然光照。

2.2 实验设计及采样

将实验用鱼分为对照组（CT组）和低氧组（HT组），每组设置3个重复实验，每个重复实验有35尾鱼。HT组溶解氧浓度为（3.15±0.21）mg/L，通过覆盖塑料薄膜、调整充气量和水流速度的大小来降低水体溶解氧浓度。水体溶解氧采用碘量法(GB 7489-87)进行测定，CT组无任何处理，溶解氧浓度为 （6.18±0.23）mg/L。实验处理周期为28 d。军曹鱼采样前24 h禁止饲喂，采样时先经丁香酚麻醉，测定鱼体的生物学指标，然后置于冰上解剖，取出整个肠道并剔除脂肪、结缔组织和内容物，用无菌生理盐水冲洗干净后，放置于2.0 mL的冻存管中，投入液氮速冻，运回实验室后转入-80℃冰箱保存进行后续实验。

2.3 转录组测序数据

实验样品送至广州基迪奥生物科技有限公司进行转录组测序工作，使用TRIzol法提取总RNA，纯度（NanoDrop 2000）和浓度（Agilent 2100）检测合格后，每组建3个转录组文库，每个文库以2尾鱼的肠道RNA等效混合进行建库，利用Illumina Hiseq-2000进行测序[23]。转录组测序共获得302 588 246条原始读数（Raw Reads），经去除含有接头、重复及测序质量较低的原始读数后共获得293 736 220条干净读数（Clean Reads）（数据未发表）。

2.4 SNP位点检测及所在Unigene的注释与功能分析

使用hisat2（http://ccb.jhu.edu/software/hisat2/index.shtml）全局比对方式，将不同长度的干净读数比对到本实验室军曹鱼参考基因组（NCBI Bioproject ID:PRJNA634421）获得比对文件，使用samtools sort对其进行排序[28]；将CT组和HT组中所有样本排序后的比对文件和相应的参考基因组序列输入到bcftools（https://github.com/samtools/bcftools）软件中，使用bcftools mpileup命令将两组军曹鱼比对后的bam文件转换成vcf文件，进行SNP分型；使用bcftools filter对vcf文件进行过滤，将质量值小于20的低质量突变位点过滤掉，保留质量值不小于20的结果；使用annovar（http://www.openbioinformatics.org/annovar/）软件对过滤后的vcf文件进行注释；使用Excel统计SNP位点的数量及分布情况[29-31]。再将包含SNP的基因序列比对到GO、KOG和KEGG数据库进行功能注释。同时，从转录组数据中筛选免疫和低氧相关的KEGG通路，再由其通路筛选出差异表达的免疫防御和低氧相关基因，进行SNP位点分析[23]。

3 结果与分析

3.1 SNP位点分析

3.1.1 SNP位点数量

利用软件SOAPsnp将CT组和HT组的转录组数据进行分析，26 120条Unigene中共存在431 845个SNP位点，其中CT组有215 953个SNP位点，HT组有215 892个SNP位点。所有SNP位点中，纯合SNP位点有152 642个（CT组76 550个，HT组76 092个）（表1）。SNP发生频率约为1/171 bp（约171 bp有1个SNP位点）。SNP密度频数分析显示，每1 000 bp含有0～5个SNP位点的Unigene出现的频率最高（图1）。

表1 SNP位点数量概况Table 1 The number of the SNP sites

图1 SNP密度频数分布Fig. 1 The density frequency distribution of SNP

3.1.2 SNP的分布及转换、颠换信息

对于Unigene中SNP位点数目统计得出，CT组和HT组中分别含有1个SNP位点的Unigene有4 337和4 384条，含有2～10个SNP位点的Unigene有8 299和8 258条；含有10个以上SNP位点的Unigene有400条和402条（图2）。

图2 SNP的分布统计Fig. 2 The statistics distribution of SNP

对军曹鱼SNP突变类型的数量分析表明，在纯合SNP位点中，转换数量明显高于颠换数量，即CT组SNP位点中，转换类型131 054个（60.69%），颠换类型84 899个（39.31%）；HT组SNP位点中，转换类型131 587个（60.95%），颠换类型84 305个（39.05%）。在核苷酸的变异类型中，C-T发生转换的频率最高，分别占纯合SNP数的44.60%（CT组34 145个）和45.21%（HT组34 400个）；T-A发生颠换的频率最高，分别占纯合SNP数的16.45%（CT组12 591个）和16.37%（HT组12 458个）（图3）。

图3 SNP突变类型数量统计Fig. 3 The numbers of different mutation types statistics of SNP

3.1.3 SNP位置分析及突变类型分析

对军曹鱼肠道转录组SNP序列位置变异类型分析结果显示（图4），CT组和HT组中SNP位点总数分别为215 953和215 892个，位于内含子和基因间区的SNP位点较多。在CT组和HT组中，位于内含子的SNP位点分别占45.13%和44.62%；位于基因间区的SNP位点分别占24.66%和25.17%；位于外显子的SNP位点分别占20.53%和20.43%；位于基因下游的SNP位点分别占13.54%和13.71%；位于3′非编码区的SNP位点分别占8.74%和8.76%；位于基因上游的SNP位点分别占4.45%和4.44%；位于5′非编码区的SNP位点分别占2.20%和2.15%；位于基因上游和下游的SNP位点分别为2 031个和2 028个; 位于剪切位点的SNP位点分别为523个和518个。

图4 SNP位置鉴定Fig. 4 Identification of SNP location

对SNP的突变类型分析结果显示（图5），突变的SNP位点数量较少，其中同义突变SNP位点数量占比最高，CT组为10.10%，HT组为10.04%；非同义突变SNP位点数量在CT组中占比为9.15%，HT组占比为9.11%。其他少量突变类型包括移码缺失、无义突变和非移码缺失等。

图5 SNP突变类型鉴定Fig. 5 Identification of SNP variants type

3.2 SNP-Unigene的注释与功能

将26 120条SNP-Unigene与GO数据库比对，结果显示，分布于生物学过程（Biological Process，BP）、分子功能（Molecular Function，MF）和细胞组分（Cellular Component，CC）3大功能分类的GO条目分别为8 043条、7 537条和7 854条，其中参与生物学过程的SNP-Unigene数量最多（图6）。将SNP-Unigene序列比对到KOG数据库进行注释和分类，结果显示，共有14 334条SNP-Unigene匹配到相应的注释信息，根据功能分为25类，其中以“信号转导机制”和“一般功能预测”类较多，分别包含2 822条和2 548条SNPUnigene（图7）。将26 120条SNP-Unigene序列进行KEGG富集分析，发现有5 160条SNP-Unigene富集到349条KEGG代谢通路中，其中以“信号转导” “全球概览地图”和“传染病”通路中富集的SNP-Unigene较多，分别为1 541条、1 123条和1 153条（图8）。

图6 SNP-Unigene GO功能注释分析Fig. 6 GO annotation analysis of SNP-Unigene

图7 SNP-Unigene KOG功能注释分析Fig. 7 KOG function annotation analysis of SNP-Unigene

图8 SNP-Unigene的KEGG功能注释分析Fig. 8 KEGG classification of SNP-Unigene

3.3 免疫防御相关SNP-Unigene的富集及特异SNP位点

根据KEGG信号通路的富集分析，发现共有3 417条免疫防御相关的SNP-Unigene被注释到“MAPK信号通路”等35条与免疫相关的通路中，其中以“MAPK信号通路”中注释的SNP-Unigene最多（262条），其次分别为“神经活性配体-受体相互作用”（228条）、“Rap1信号通路”（190条）、“Ras信号通路”（164条）等信号通路（表2）。

表2 免疫防御相关SNP-Unigene的KEGG富集分析Table 2 KEGG enrichment analysis of immune-related SNP-Unigene

3.4 差异表达免疫防御及低氧相关基因SNP位点分析

基于军曹鱼低氧胁迫后肠道转录组分析，得到大量的差异基因（筛选条件为错误发生率（False Discovery Rate，FDR）小于0.05且|log2FC|大于1，其中FC为差异倍数）。对“MAPK信号通路”等7个免疫通路及“HIF-1信号通路”中差异基因进行了SNP位点分析（表3）。

表3 CT和HT转录组中差异表达免疫防御及低氧相关基因SNP位点分析Table 3 SNP identified in differentially expressed immune and hypoxia-related genes in CT and HT transcriptomes

4 讨论

4.1 低氧胁迫下军曹鱼转录组数据中SNP位点鉴定及功能注释

SNP作为一种分子遗传标记，在水生动物的遗传选育中发挥着重要作用[25]。本研究从CT转录组、HT转录组数据中分别获得大量的SNP位点，SNP发生频率依次为0.00 586和0.005 858（1/171 bp），高于盐度胁迫对军曹鱼（1/730 bp）[27]、半滑舌鳎（Cynoglossus semilaevis）（1/491 bp）[32]、大 黄 鱼（Larimichthy crocea）（1/506 bp）[33]、罗氏沼虾（Macrobrachium rosenbergii）（1/1 152 bp）[24]、大西洋鲑（Salmo salar）（1/641 bp）[34]SNP发生频率，但低于马氏珠母贝（1/100 bp）[20]。有研究表明，不同物种中SNP分布频率差异很大，可能与研究物种、遗传背景和生存环境有关[19]，也可能是因为在表达序列标签中筛选SNP的要求和条件不同[21]。SNP碱基替换是产生基因突变和推动进化的原始动力[35]。碱基替换类型分为转换(A-G和T-C的相互置换)和颠换(A-T、C-G、A-C、T-G的相互置换)。理论上，生物中发生颠换与转换比例为1∶2[23]。本研究中，SNP转换类型的比例约为60%，颠换类型的比例约为40%，转换大于颠换，符合“转换偏差”原理。在核苷酸变异类型中，C-T发生转换的比例最高，与大多数水产动物的研究结果相一致[19,36]，原因可能是CpG二核苷酸上甲基化的胞嘧啶残基容易脱去氨基形成胸腺嘧啶[37]，或受物种进化过程中承受的压力影响[38]。

根据SNP在基因组的分布位置可分为3类：基因编码区SNPs（cSNPs）、基因间SNPs（iSNPs）和基因周边SNPs（pSNPs）[39]。其中cSNP和pSNP会对基因的表达产生影响[22]。根据对生物遗传性状的影响，cSNPs又分为同义突变和非同义突变，前者编码序列的变化不会影响翻译后的氨基酸序列，故不会导致蛋白质的功能变化；后者编码序列的变化会引起翻译后氨基酸序列的变化，改变蛋白质的生物学功能[40]。由于受到较强的选择压力，非同义SNP位点的突变经常会影响生物的表型，而关于非同义SNP位点的研究已在团头鲂[25]、大黄鱼[41]有相关报道。对军曹鱼肠道转录组序列SNP位置和突变类型分析发现，位于基因编码区的SNP位点以及非同义突变SNP位点占比较多。因此，研究基因编码区非同义SNP可能对军曹鱼具有重要的生物学意义。此外，GO、KOG及KEGG功能注释结果表明，大部分SNP-Unigene涉及信号转导、传染病、代谢、癌症、内分泌系统和免疫系统等多种途径，表明军曹鱼幼鱼遭受低氧胁迫后，这些SNP-Unigene参与各种生命活动，并影响机体的各种生物学性状。

4.2 低氧胁迫下军曹鱼免疫相关基因SNP位点分析

免疫相关基因的SNP在水产遗传育种中的应用对鱼类抗逆性性状的选育具有重要意义[41]。在本研究中，为了解军曹鱼免疫防御机制，对转录组测序数据进一步分析，结果发现，3 417条SNP-Unigene被注释到35条与免疫防御相关的通路中，大部分SNP-

Unigene参与“MAPK信号通路” “神经活性配体-受体相互作用” “Rap1信号通路” “Ras信号通路”等免疫通路。我们从相关免疫通路中筛选出了大量差异表达的免疫相关基因及其SNP位点。因此，这里对部分关键基因进行后续的研究分析。

IL1R1是细胞表面重要的免疫基因，是白细胞介素-1β（IL-1β）的受体，在机体炎症反应和免疫信号转导过程中非常重要[42]。Lam等[43]研究报道，慢性缺氧时，大鼠颈动脉体细胞中促炎细胞因子及其受体（IL-1β和IL1R1）的表达量显著增加，表明其在颈动脉体局部炎症中起着重要功能。在本研究中，军曹鱼受低氧胁迫后，肠道IL1R1基因的表达量显著上调，提示IL1R1基因可能在军曹鱼免疫调节和抵御低氧应激中发挥着重要作用。此外，在基因多态性方面，已发现IL1R1基因的SNP与多种人类疾病相关[44]。胡亮[45]对岗巴绵羊的候选特异性SNP位点进行基因注释和富集分析，发现IL1R1和IL1R2基因参与炎症反应及自身免疫性调节。因此，暗示IL1R1基因可作为军曹鱼免疫防御的候选基因，这为今后进行IL1R1基因的多态性与疾病关联研究提供了参考。

Wnt信号通路是一个高度保守的系统，通过Wnt蛋白与细胞表面受体相互作用，调控着细胞增殖、分化和黏附等多种复杂的生物学过程[23,46]。Wnt配体蛋白与跨膜受体卷曲蛋白(FZD)和低密度脂蛋白受体相关蛋白(LRP5/6)形成的细胞表面受体复合物结合，从而激活细胞内的Wnt信号通路[47]。在本研究中，军曹鱼受低氧胁迫后，Wnt信号通路中FZD3和LRP5基因表达量显著下调，且存在多个SNP位点。目前有关鱼类中FZDs家族基因的研究主要涉及基因克隆、组织表达及在鱼类发育过程中的生理作用[48]。张丽晗等[48]研究了铜胁迫下黄颡鱼（Pelteobagrus fulvidraco）FZD家 族 基 因（FZD2、FZD3A、FZD4和FZD10）的mRNA水平。FZD3在小鼠和人的毛囊发育、细胞迁移等方面发挥着关键作用[49]，已有研究者对该基因多态性位点进行分析，筛选影响性状的SNP位点。Zhao等[49]对中国美利奴羊FZD3基因的多态性进行检测，发现SNP1和SNP2与羊毛性状有显著的相关性，可作为绵羊育种的潜在遗传标记。LRP5基因在Wnt信号通路中扮演着重要角色，其突变与骨质疏松、癌症等疾病的发生密切相关[50]。已有研究表明[50]，LRP5基因多态性与2型糖尿病患者的β细胞功能和脂质代谢有关。由此推测，FZD3和LRP5基因可能参与调控军曹鱼低氧应激过程中的免疫反应，具体机制有待进一步探讨，而研究这些免疫基因的SNP对低氧胁迫下军曹鱼免疫防御相关SNP位点的开发也具有潜在意义。

本研究结果表明，参与mTOR信号通路的差异基因LIPIN1表达水平上调，也存在多个SNP位点。Lipin1是由LIPIN1基因编码的磷脂酸磷脂酶，具有双向调控脂肪代谢的功能，既能催化磷脂酸去磷酸化得到二脂酰甘油，从而促进甘油三酯和磷脂的合成；还具有转录激活功能，参与调控脂肪代谢相关基因表达[51-52]。关于LIPIN蛋白家族的研究主要集中在蛋白功能和基因表达方面，多态性方面以小鼠、猪、牛等哺乳动物为主，水产经济动物较为匮乏。有报道称[53]，LIPIN1基因受HIF-1调控，其表达水平上调与非脂肪细胞中的甘油三酯和脂滴的积累有关，提示LIPIN1是细胞适应低氧应激的重要中介物。He等[54]在猪的LIPIN1基因的编码区和UTR检测到5个SNP位点与猪的脂肪沉积性状存在显著相关。这些研究表明，LIPIN1基因可能对军曹鱼低氧胁迫过程中脂质代谢有一定程度的影响，可为LIPIN1基因多态性与目标性状关联研究提供了新思路。

4.3 低氧胁迫下军曹鱼HIF-1信号通路中差异基因SNP位点分析

除了免疫防御相关信号通路所涉及的差异基因外，还有多种基因是通过其他信号通路来共同应答鱼类的低氧胁迫。HIF-1信号通路是鱼类遭受低氧应激时的关键通路之一，通路中的相关基因不仅参与机体代谢的调控，而且还会通过抑制血红细胞的增殖和凋亡，提高血氧亲和力，以此来适应低氧应激[55]。在本研究结果中发现，HIF-1信号通路中存在PIK3CA、HK1、ANGPT1、EPAS1、GAPDH、ANGPT2、EDN1和ENO3等8个差异基因也同样含有多个SNP位点。

PIK3CA基因是编码I类磷脂酰肌醇3激酶（PI3Ks）的催化亚单位，其突变能激活PI3K/AKT信号通路，进而参与细胞增殖、凋亡和蛋白质合成等多种生理过程[56-57]。Zhang等[56]研究发现，SmPIK3CA基因参与大菱鲆的免疫反应，并在该基因中筛选出4个抗病相关的SNPs。EPAS1作为激活低氧调控基因表达的关键转录因子[58]，在细胞代谢和氧化还原反应中起着重要作用。关于EPAS1基因SNP位点的研究多集中在藏族人和高原哺乳动物，在鱼类中的研究较少。在藏族人和高原动物低氧适应研究中发现，EPAS1基因都受到选择且信号显著。在牦牛中发现EPAS1基因多态性与血红蛋白(HGB)关系紧密[58]。在藏绵羊EPAS1基因与低氧适应相关性研究中发现，该基因第8内含子C871T 变异位点的等位基因(T)突变频率与海拔高度呈现正相关，且TT基因型能显著增加血红蛋白浓度[59]。由此推测，HIF-1信号通路中筛选出来的PIK3CA等8个显著差异基因可能对军曹鱼低氧相关SNP位点的挖掘具有重要的参考价值。

综上所述，本研究基于低氧胁迫下军曹鱼肠道转录组测序数据，获得了大量与免疫、低氧相关差异基因的SNP位点，丰富了军曹鱼分子标记的数据，这些候选标记将为今后深入了解军曹鱼响应低氧胁迫的分子机制、SNP分子标记开发和遗传连锁图谱构建等研究提供重要的理论依据。