基于皮肤转录组的半滑舌鳎无眼侧黑化相关基因的筛选鉴定

赵娜 ，何晓旭 ，贾磊，朱春华，张博*

(1.广东海洋大学 水产学院，广东 湛江 524025；2.天津市水产研究所，天津 300457；3.上海海洋大学 水产种质资源发掘与利用教育部重点实验室，上海 201306；4.天津市农学院，天津 300384)

1 引言

动物体色具有迷惑或躲避天敌等自我保护功能，不同体色对鱼类的生存也有着特殊的意义。鱼类在生长发育过程中体色的变化主要受遗传因素的影响，此外，外部环境[1-3]、饲料营养[4]和神经内分泌物[5-7]也会影响鱼类的体色变化[8-9]。多数色素细胞的遗传是保守的[10-11]，如一些与色素相关的基因先后在小鼠（Mus musculus）的皮肤和斑马鱼（Danio rerio）的皮肤中均有发现，而mitf、tyr、dct、sox10和ednr等[12-15]基因在皮肤中的生长发育表达特征也非常相似。此外，在硬骨鱼中发现的slc24a5基因也参与哺乳动物体色或毛发中色素细胞的生长和发育[16]。在鱼类体色生成过程中，黑色素的产生和扩散起主导作用，其中酪氨酸酶基因家族中的酪氨酸酶（TYR）是主效基因[17]。TYR 受多个基因控制，酪氨酸酶受体基因1（tyr1）、性别决定区Y-box 10 基因（sox10）[18]、多巴素异构酶基因（dct）[19]、黑皮质素受体-1（mcir）对酪氨酸酶的产生均具有调节作用[20-22]。

研究表明，饲养密度[23]、水温[24]、光照强度[24]能显著地影响比目鱼（Pleuronectiformes）无眼侧黑化的发生。科学的养殖管理方法如营养成分的摄入控制和亲本筛选[25]也能有效地减少比目鱼体色异常的发生。Nakamura 等[26]认为，缺乏光敏物质会导致黑色素合成不足，如核黄素、类胡萝卜素、维生素A 和维生素D，而核黄素被认为是其中最重要的因素。研究表明[27]，牙鲆（Paralichthysolivaceus）的白化病是视紫红质缺乏的结果，其产生依赖于DHA。而Devresse 等[28]认为，色素沉着似乎与DHA 无关，而是与DHA 和二十碳五烯酸的比例有关。胰蛋白酶、维甲酸相关酶、维甲酸受体及其结合受体（TR、VDR、PPAR 等）的突变及其对仔鱼色素细胞分化的调控，可能影响鱼体两侧色素细胞的命运，进而影响体细胞色素异常[29-30]。不同环境背景色养殖实验中发现，条斑星鲽（Verasper moseri）在浅色养殖环境下mch-R2的mRNA 表达量低于黑色养殖环境，但是两种环境中mch-R1的mRNA 没有变化，表明了mch-R2参与了不同环境背景下体色变化过程，mch-R可能对黑色素聚集激素（MCH）有抑制调节作用[31]。无眼侧黑化消褪的牙鲆垂体MCH mRNA表达水平显著高于无眼侧黑化个体，同时与无眼侧正常个体无显著差异，表明了MCH 及其基因均参与了无眼侧黑化消褪的调控过程[32]。在鱼类中，tyr基因突变导致青鳉（Oryzias latipes）的黑色素细胞异常，皮肤呈现橙色[33]。虹鳟鱼（Oncorhynchusmykiss）中的tyr基因突变显示嵌合白化（伪白化）和完全白化现象[34]。黄颡（Pelteobagrus fulvidraco）的tyr基因突变则显示黑色素缺乏，呈现体色部分白化[35]。综上所述，遗传和环境因素的相互作用影响了鱼类的体色异常。

半滑舌鳎（Cynoglossus semilaevis）的工厂化养殖过程中经常出现高比例体色异常的苗种，即无眼侧出现黑色素的黑化现象[36-37]，主要表现为无眼侧部分覆盖斑块色素，有时甚至完全呈黑底，不同的养殖环境会存在不同程度的黑化现象。体色异常极大地影响半滑舌鳎的销售价格，长期困扰养殖企业，却一直未能得到有效解决[38]。针对半滑舌鳎无眼侧黑化，影响体色变化的新基因有待挖掘和鉴定。转录组是目前研究基因表达的重要手段，是筛选新功能基因的有效方法。本文选取半滑舌鳎体色异常皮肤样本和正常皮肤样本，进行转录组测序的比较分析。以期找出与半滑舌鳎无眼侧黑化相关的表达显著变化的功能基因，为半滑舌鳎体色异常的机制提供分子水平的理论支持。

2 材料和方法

2.1 鱼类取样准备

从天津乾海源水产养殖有限公司获得3 尾无眼侧黑化的工厂化养殖半滑舌鳎，平均体重为（0.61±0.02）kg，平均体长为（28.2±0.1）cm。用浓度为120 mg/L的间氨基苯甲酸乙酯甲磺酸盐（MS-222）（北京格林恒兴生物科技有限公司）将鱼麻醉后，使用灭菌的手术剪和镊子剪开实验鱼的皮肤并采集样本，获取无眼侧色素沉着的皮肤组织和非色素沉着的正常皮肤组织（图1），并在-80℃下保存，直至提取RNA。

图1 半滑舌鳎无眼侧黑化鱼体色特征Fig.1 The body color feature of the Cynoglossus semilaevis with no eye side blackened fish

2.2 RNA 测序、数据过滤及生物信息学分析

使用Trizol 试剂（日本Takara）提取总RNA。安捷伦2100 生物分析仪（安捷伦科技公司，美国圣克拉拉）用于评估RNA 质量和浓度。无眼侧黑化皮肤组织中的样品标记为csc1、csc2 和csc3；正常皮肤样品标记为csh1、csh2 和csh3。测序平台为Illumina HiSeq™2500（Illumina 公司，圣地亚哥，加利福尼亚州，美国）。用带Oligo（dT）的磁珠富集mRNA，再以片段化后的mRNA 为模板，用六碱基随机引物合成cDNA 第一链，并加入缓冲液、dNTPs（A、U、G、C）、核糖核酸酶H 和DNA 聚合酶I 合成cDNA 第二链；再经过磁珠纯化并加EB 缓冲液洗脱，对洗脱之后的双链cDNA进行末端修复、加碱基A、加测序接头处理，然后用磁珠进行片段大小选择、降解含U 链，并进行PCR 扩增，从而完成整个文库制备工作。根据物种的参考基因组和基因信息（GCA_000523025.1），将转录组测序所产生的数据比对到参考基因组上。获得过滤后数据后，将其与参考基因组进行序列比对，获取测得条目的定位信息。

2.3 差异表达分析和样品相关性分析

用bowtie2 软件将过滤后数据比对到参考转录本序列，而后利用RSEM 软件，调用bowtie2 的比对结果进行统计，得到每个样品比对到每个转录本上的测得条目数目，并对其进行FPKM 转换。FPKM 是指每百万测得片段中来自某一基因平均每1 000 碱基长度的测得片段数目，其同时考虑了测序深度和基因长度对测得片段计数的影响，是目前最为常用的基因表达水平估算方法。用R 语言包edgeR 进行差异分析，筛选阈值为错误发现率（False Discovery Rate，FDR）小于0.05，log2FC（FC 为实验组与对照组可测量比值）大于1 或log2FC 小于-1。利用皮尔森相关系数进行相关性分析，相关系数越接近1，表明样品之间表达模式的相似度越高。相关性系数是介于-1～1 之间的实数，当相关性系数介于-1～0 时，表明变量之间存在负相关关系；当相关性系数介于0～1 时，表明变量之间存在正相关关系；当相关性系数为0 时，二者之间不存在相关性。

2.4 差异基因GO 和KEGG 富集及功能筛选

根据分析目的筛选出差异基因后，研究差异基因在Gene Ontology（简称GO，http://www.geneontology.org/）中的分布状况。KEGG（Kyoto Encyclopedia of Genes and Genomes）是有关通路的主要公共数据库。通路显著性富集分析以KEGG 通路为单位，应用超几何检验，找出差异基因相对于所有有注释的基因显著富集的通路。GO 富集分析方法为Hyper-geometric 分布，与KEGG 富集分析方法一样，我们选取FDR≤0.05的GO 条目作为显著富集的GO 条目。FDR≤0.05 的通路定义为在差异表达基因中显著富集的通路。我们挑选了富集最显著的20 条通路，绘制差异表达基因KEGG 富集散点图。KEGG 富集程度通过富集因子、Qvalue 和富集到此通路上的基因个数来衡量。其中富集因子指该通路中富集到的差异基因个数与注释基因个数的比值。富集因子越大，表示富集的程度越大。Qvalue 是做过多重假设检验校正之后的Pvalue，Qvalue 的取值范围为0～1，越接近0，表示富集越显著。差异基因功能筛选集中于比目鱼体色异常的相关6 条代谢通路，这6 条KEGG 代谢通路如下：花生四烯酸代谢通路、甲状腺激素合成通路、光转导通路、黑色素生成通路、视黄醇代谢通路和PPAR 信号通路。

2.5 候选基因的qPCR 定量验证

用qPCR 方法检测黑化鱼正常皮肤样本（csh1、csh2 和csh3）和黑化皮肤样本（csc1、csc2 和csc3）中的候选mRNA 的表达。采用北京百奥创新科技有限公司核酸提取试剂盒和PowerUp-SYBR-Green-Master-Mix 进行mRNA 提取和qPCR。在42℃加热2 min，在冰浴中孵育2 min，加入提取RNA 的试剂，在42℃加热15 min，提取总RNA，然后在85℃下保持5 min。将RNA 储存在-80℃下。cDNA 合成按照PrimeScript™RT reagent Kit 试剂盒操作方法进行，cDNA 储存在-20℃下备用。使用QuantStudio6 Flex 实时PCR 系统（美国加利福尼亚州赛默飞世尔科学公司）在96 孔板中进行qPCR，反应体积为20 μL：AceQ®qPCR SYBR®Green Master Mix (Without ROX) 10.0 μL，Primer 1(10 μmol/L) 0.4 μL，Primer 2 (10 μmol/L) 0.4 μL，cDNA 1.0 μL，ddH2O 8.2 μL，配制过程在冰上完成。反应程序包括：50℃预热2 min，95℃预热2 min，随后在95°C 下，持续30 s，58°C 持续30 s，95°C 持续15 s，进行40 个循环。选择β-肌动蛋白的表达水平作为正常化的内源性对照。每个样品设置3 个技术重复，并使用2-ΔΔCt方法计算mRNA 的相对表达量。

3 结果

3.1 RNA-seq 数据披露

通过RNA-seq，6 个样本共产生416 464 144 个过滤后数据，平均数据量在10 G 左右（9 208～11 897 M），Q20 和Q30 分别在92%和97%左右。数据结果表明，6 个样品的测序质量良好，数据量充足，详细数据如表1 所示，转录组数据已经上传至美国国立生物技术信息中心，数据获取编号为（PRJNA665385）。

表1 6 个半滑舌鳎皮肤转录组测序样品过滤后数据产出及质控分析Table 1 Clean reads output and quality control analysis of six Cynoglossus semilaevis skin transcriptome sequencing samples

对两组6 个样品的RNA-seq 的过滤后数据进行基因组比对统计，6 个样品的过滤后数据数量分布在30 772 761～39 743 985 之间，唯一比对到参考基因组上的测得条目百分比分布在76.47%～81.36%之间，不能比对到参考基因组的测得条目数目（双端统计）普遍低于20%（表2）。以上结果进一步说明6 个样品的转录组测序数据质量过关，可作进一步分析。

表2 6 个半滑舌鳎皮肤转录组测序样品基因组比对统计表Table 2 Genomic comparison statistics of six Cynoglossus semilaevis skin transcriptome sequencing samples

3.2 差异表达分析及样品相关性分析

csc1、csc2 和csc3 样本之间的相关系数处于0.97～1之间，接近1，表明变量之间存在正相关关系，样品之间表达模式为强相关。csh1、csh2 和csch3 样本之间的相关系数处于0.76～1 之间，表明变量之间存在正相关关系，样品之间表达模式为强相关，但是相关性弱于csc 样本。csc1、csc2 和csc3 样本和csh1、csh2和csch3 样本的相关系数处于0.18～0.45 之间，表明变量之间存在正相关关系，样品之间表达模式处于弱相关和中度相关（图2）。

图2 测序样品相关性分析热图Fig.2 Correlation analysis heat map of sequencing samples

图3a 显示，样本无眼侧黑化皮肤csc1、csc2、csc3和无眼侧正常皮肤csh1、csh2、csch3 之间差异基因总数1 665 个，其中上调基因497 个，下调基因1 168 个。通过转录组测序得到黑化半滑舌鳎的无眼侧正常和无眼侧黑化皮肤的转录组序列信息，对转录组数据分析共得到837 540 个位点，其中A/G、C/G、C/T、G/A之间转换数量最多。通过图3b 也可以观察到，csc与csh 两组之间部分基因的表达量存在明显差异（图3）。

图3 半滑舌鳎正常皮肤与黑化皮肤样品间差异表达基因火山图和差异表达基因丰度热图Fig.3 Volcano map of the number of differentially expressed genes and heat map of abundance of differentially expressed gene between normal and hypermelanotic skin samples of Cynoglossus semilaevis

3.3 差异表达基因的GO 及 KEGG 富集分析及体色功能差异基因筛选

对无眼侧黑化组和无眼侧正常组皮肤进行GO富集分析时发现，在生物过程水平上差异基因主要富集在细胞过程、信号组织过程、代谢过程、生物调控、生物学过程这5 个GO 通路；在细胞组分上主要富集在膜、细胞、细胞组分、膜组分、细胞器这5 个GO 通路；在细胞功能上主要富集在连接、催化活性这2 个GO 通路（图4a）。由KEGG 富集散点图（图4b）可见，无眼侧黑化和无眼侧正常皮肤表达差异显著的基因参与的通路主要有：丙酮酸代谢、C5-支链二元酸代谢、甲烷代谢及苯丙氨酸、酪氨酸和色氨酸的生物合成等。

图4 半滑舌鳎正常与黑化皮肤转录组测序差异基因富集注释的GO 通路和KEGG 通路Fig.4 GO pathways and KEGG pathways analysis of differential expressed genes enrichment annotations of transcriptome sequencing of normal and hypermelanotic skin of Cynoglossus semilaevis

3.4 定量PCR 验证结果

从转录组富集结果中选取的6 条KEGG 代谢通路如下：花生四烯酸代谢通路、甲状腺激素合成通路、光转导通路、黑色素生成通路、视黄醇代谢通路和过氧化物酶体增殖物激活受体信号通路。在以上6 条通路中选取差异最显著的前10 个基因进行定量表达验证。对黑化半滑舌鳎皮肤中差异表达基因进行相对表达水平定量分析发现，黑化半滑舌鳎皮肤差异表达基因的表达趋势与RNA-seq 定量结果一致，表明转录组测序结果可靠。针对6 个测序样本，对6 个候选体色相关mRNA 定量qPCR 表达检测结果绘制箱线图（图5）。其中5 个基因txndc、alox15b、ptgs2、ptgis、atp1a2a在两组样本中均有显著差异（p＜0.05），且有4 个基因在黑化组中表达均高于对照组，只有atp1a2a表达在黑化组表达较低，而acbd7基因则在两组中无显著差异（p＞0.05）。

图5 6 个候选体色相关mRNA 在半滑舌鳎6 个测序样品中的qPCR 表达检测结果箱线图Fig.5 Box plot of qPCR expression results of 6 candidate body color-related mRNAs in 6 sequencing samples of Cynoglossus semilaevis

4 讨论

色素的异常沉着是比目鱼养殖过程中的一个主要问题，体色异常显著降低了比目鱼的市场价值[39]。近年来，关于半滑舌鳎体色异常的研究多有文献报道[40-41]，其主要集中于已知体色调控基因的表达模式、着色剂饲料增色效果、理化因子对体色的影响等[42-43]，而从组学角度挖掘半滑舌鳎体色调控新基因的研究则比较缺乏。本研究运用二代测序技术对黑化半滑舌鳎的皮肤转录组进行测序分析，得到了半滑舌鳎无眼侧黑化皮肤组织和非色素沉着的正常皮肤组织的转录组数据，从6 条代谢通路（花生四烯酸代谢通路、甲状腺激素合成通路、光转导通路、黑色素生成通路、视黄醇代谢通路和过氧化物酶体增殖物激活受体信号通路）中筛选出半滑舌鳎体色异常相关的差异基因，为今后半滑舌鳎进行体色异常排查提供了参考。选取其中差异最显著的10 个基因进行验证后发现有5 个差异基因与测序结果趋势相一致，分别为txndc、alox15b、ptgs2、ptgis、atp1a2a，且前4 个基因均在黑化组中上调，只有atp1a2a在黑化组中下调表达，暗示了这5 个差异基因可能参与了无眼侧皮肤黑化的调控。

TXNDC 是含硫氧还原蛋白，此基因编码内质网蛋白质的二硫键异构酶家族成员，该酶催化蛋白质折叠和巯基-二硫键互换反应。硫氧还原蛋白5（TXNDCS），是近年发现的PDI 家族成员之一，具有抗凋亡、抗氧化损伤和促进细胞增殖、促进血管形成、参与细胞炎症、能量代谢等多种生物功能[44-45]。ALOX15B 是花生四烯酸15 脂氧合酶，编码该酶的基因编码参与脂肪酸氢过氧化物生产的结构相关的非血红素铁双加氧酶的脂氧合酶家族成员，其所编码的蛋白质仅将花生四烯酸转化为15S-氢过氧二十碳四烯酸[46-47]。PTGS2 是前列腺素内过氧化物合酶2，又称环加氧酶2，是花生四烯酸（AA）转化为前列腺素（PGs）的关键酶。在前列腺素合成初期起催化作用，可将AA 代谢成各种PGs 类产物[48]，这些PGs 参与维持机体内多种生理过程，是胚胎移植初期时胚胎和子宫信息交流的重要介质之一，在调节发情周期、妊娠维持和分娩中起着关键作用[49-50]。PTGIS 是前列环素合酶，是细胞色素P450（CYP450）超家族的8 族（CYP8）成员，也是AA 的主要产物[51]。研究发现，PTGIS 在许多生理和病理过程中起着重要作用，该基因编码前列腺素12 合酶，并催化前列腺素I2（PGI2）的合成[52]。atp1a2是ATP 酶A2 亚单位，该基因编码的蛋白质属于P 型阳离子，转运ATP 酶家族和Na+/K+-ATP 酶亚家族[53]。

这5 个差异基因中有3 个基因均与AA 有一定的相关性，即ALOX15B 是AA15 脂氧合酶，PTGS2 是AA 转化为PGs 的关键酶，PTGIS 是AA 的主要产物。AA 是一种不饱和脂肪酸，由磷脂酶 A2（Phospholipase A2，PLA2）催化而来，可通过环氧合酶（Cyclooxygenase，COX）、脂氧合酶（Lipoxygenase，LOX）与CYP450 途径生成多种代谢产物。在海水鱼类中，AA 主要在生长、存活、调节免疫、抗应激及繁育等方面发挥重要作用[54-56]。目前，关于AA 在鱼类方面的研究主要集中在饲料中改变AA 含量，分析其对鱼类的生长性能、脂质代谢、脂肪酸营养、肠道组织、色素沉积等方面的影响。Sargent 等[55]指出，在比目鱼的开口饵料中添加一定量的AA，将对背部皮肤色素沉积产生较大影响。Estévez 等[57]认为，饮食中增加的AA 含量对色素沉积有负面影响。McEvoy 等[58]、Copeman 等[59]、Bell 等[60]的研究均表明，高AA 水平与各种比目鱼的色素沉着有关[58-60]。Willey 等[61]用富含AA 的轮虫(Rotifer)和卤虫（Artemia）喂养大西洋牙鲆（Paralichthys dentatus），发现增加AA 饲料和色素沉着发生率存在剂量依赖效应。Villalta 等[62]改变了塞内加尔舌鳎（Cynoglossus senegalensis）日粮中AA 的含量，发现随着AA 的增多伴随着PG 含量的升高，并且PG 含量与色素沉积存在着紧密的相关性。Lund 等[63]在欧洲鳎（Solea solea）中发现，体色异常只与AA 的绝对量相关，而与AA 和其他多不饱和脂肪酸（Polyunsaturated Fatty Acid，PUFA）的比例无关。Boglino 等[64]研究发现，PGE2 在塞内加尔舌鳎中也引起了体色异常的结果。以上研究都表明AA 的变化会影响比目鱼类体表色素沉着。导致AA 诱导比目鱼色素沉着的确切机制尚不清楚，但已存在3 种假说：第一，改变饲料中AA 含量可能导致神经组织膜的次优生化组成，而神经组织参与控制变态、黑色素合成和色素团分化[57]；第二，AA 衍生的二十碳五烯酸的过量产生会导致鱼类经历非生物化学诱导的应激，从而导致色素失调[55]；第三，酪氨酸酶的合成和降解是调节色素沉着所必需的，饮食中n-6多不饱和脂肪酸可以改变酪氨酸酶的合成和降解[65]。有研究将正常体色个体与由AA 修饰引起的变色日本牙鲆进行对比研究，发现黑素皮质激素（a-黑素细胞刺激激素a-MSH 和肾上腺皮质激素ACTH）在色素形成过程中起重要作用[66]。而MSH 的作用主要为激活酪氨酸酶，并促进酪氨酸酶合成，从而促进黑色素合成，使皮肤颜色加深。

本研究筛选出半滑舌鳎体色异常相关的5 个差异基因，其中alox15b、ptgs2可以将AA 转变或者使其代谢为其他物质，从而改变半滑舌鳎中AA 的含量。而ptgis不仅是AA 的主要产物，还是CYP450 中的成员，饲料中营养素对鱼类体色的影响主要就是通过鱼体色素细胞或色素体起作用的。因此，本文推测alox15b、ptgs2、ptgis基因是通过影响AA 进而影响了半滑舌鳎无眼侧体色的色素沉积。本研究为半滑舌鳎无眼侧体色异常的进一步研究提供了分子线索，关于AA影响鱼类体色异常的分子机制仍有待进一步研究。