基于转录组测序探究乌鳢皮肤白化的分子机制

李 岩，周 燕，雷 骆，苏 建，樊 威，罗 煜，李俊婷，高 贺，周朝伟

( 1.西南大学 水产学院，重庆 402460；
2.中国水产科学研究院 珠江水产研究所，广东 广州 510380；
3.四川省内江市农业科学院，四川 内江 641000 )

体色是脊椎动物最重要的表型特征之一，具有择偶、物种识别、隐藏、警告等重要的功能[1]。色素沉着失调会引起体色的变异，其中最典型的表现是皮肤白化，通常被认为是由先天性代谢紊乱所致[2]。目前，公认的导致动物白化的原因主要有两种：黑色素细胞缺乏和黑色素合成障碍。在鱼类上，黑色素细胞合成载有大量黑色素的黑素体，黑素体的发育涉及到细胞器的转运，经过4个不同的阶段后，含有黑色素的成熟黑素体被转运到角质层，最终使鱼体出现不同的体色[3-4]。若生理、环境或遗传上的任何改变干扰了动物黑色素细胞发育与分化或黑色素合成过程中的任一环节，均会导致动物白化现象的发生，但是影响鱼类皮肤色素沉着最基本的因素是遗传[5]。鱼类频繁出现白化现象，如白化的斑马鱼(Daniorerio)[6]、骨颊海鲇(Osteogeneiosusmilitaris)[7]和叉牙鲷(Sarpasalpa)[8]等。然而，大量的研究集中在分离和鉴定与鱼类白化性状相关的基因方面，分子机制相关研究较为缺乏，有关色素沉着相关基因和其调控网络之间的联系仍需要进一步探究。

乌鳢(Channaargus)属鲈形目鳢属，俗称乌鱼和黑鱼，白乌鳢(C.argusvar.)为乌鳢白化变异种，俗称白乌鱼、白乌棒，主要分布于四川省嘉陵江中下游流域、沱江流域，是我国重要的肉食性经济鱼类[9]。与普通乌鳢相比，白乌鳢由于其营养价值高、体色特殊受到了养殖户及观赏鱼从业者的关注。截至目前，乌鳢白化形成的分子机制尚未查明。转录组学技术可广泛应用于分析生物不同性状、不同生理状态和不同组织中基因表达变化情况，有利于挖掘与特定性状密切相关的基因。笔者利用Illumina HiSeq技术对乌鳢和白乌鳢的体色差异进行研究，探究乌鳢皮肤白化的分子调控机制，以期对其种质资源保护和遗传改良提供理论基础。

1.1 试验材料及样品采集

试验所用白乌鳢、乌鳢均由内江市农业科学院水产研究所提供，出苗时间、生长环境、喂养条件均一致。白乌鳢为乌鳢白化变异种，体色多为灰白色，纵观背部到腹部，呈现出从青灰色到灰白色逐渐过渡的趋势(图1)；
乌鳢身体背部较暗，呈黑褐色，腹部色浅，体侧具多角形不规则深色斑纹(图1)。随机挑选健康状况良好的白乌鳢和乌鳢各3尾，体长(38.00±2.09) cm，体质量(513±23.78) g。取3尾白乌鳢的体侧皮肤组织(图1a)，标记为白乌鳢1、白乌鳢2、白乌鳢3；
取3尾乌鳢体侧花色均匀的皮肤组织(图1b)，标记为乌鳢1、乌鳢2、乌鳢3；
共6个组织样品。将试验鱼用丁香酚麻醉后，快速取下其皮肤组织，置于 1.5 mL无酶EP管，迅速放入液氮中速冻后转入-80 ℃中保存，用于随后的转录组测定。另取15尾白乌鳢的体侧皮肤组织和15尾乌鳢体侧花色均匀的皮肤组织，每尾鱼取3块皮肤组织，冻存于-80 ℃，用于随后的实时荧光定量PCR(qRT-PCR)分析。

图1 白乌鳢与乌鳢皮肤采样部位Fig.1 Skin sampling location of albino northern snakehead C. argus var. and northern snakehead C. argus

1.2 RNA-seq文库制备及测序

采用动物组织总RNA提取试剂盒(天根生化科技有限公司，北京)提取总RNA，再通过1%的琼脂糖凝胶电泳分析28S RNA和18S RNA降解程度以及是否有污染，使用NanoDrop 2000核酸检测仪(Thermo公司,美国)检测RNA纯度,使用Qubit©2.0 Flurometer (Life Technologies，美国)测定总RNA浓度，用Agilent 2100(安捷伦科技公司，美国)精确检测RNA的完整性。样品检测合格后，用带有Oligo(dT)的磁珠富集真核生物mRNA。随后加入裂解缓冲液将mRNA打断成短片段，以这些mRNA为模板，用六碱基随机引物random hexamers合成一链cDNA，加入缓冲液、dNTPs、DNA polymerase Ⅰ和RNase H，合成二链cDNA，用AMPure XP beads纯化双链cDNA。纯化的双链cDNA先进行末端修复、加poly(A)尾并连接测序接头，再用AMPure XP beads进行片段大小选择。最后进行PCR富集得到最终的cDNA文库。库检合格后，将不同文库按照有效浓度及目标下机数据量的需求集中后进行HiSeq双端测序，每个样本测6 Gb数据量，称为原始测序数据(Raw reads)，并对其原始测序数据进行数据过滤，去除含adapter和含N比例大于10%的数据，去除低质量数据(质量值Q≤10的碱基数占整条数据的50%以上)，得到的数据称为过滤后数据(Clean reads)，用于后续的信息分析。试验用6个样品均由武汉菲沙基因信息有限公司进行转录组测序及分析。

1.3 差异表达基因的组装筛选、功能注释和富集分析

以乌鳢的genome为参考基因(http://gigadb.org/dataset/100279)，对高质量过滤后数据进行注释。使用FPKM法[10]估算基因表达量，即基因的长度和比对到该基因的次数。用R语言包DEseq2进行差异分析[11]，以|log2(fold change)|≥2且P≤0.05作为阈值。错误发现率(FDR)是通过对差异显著性P值进行校正得到的。采用公认的Benjamini-Hochberg校正方法对原有假设检验得到的显著性P值进行校正，并最终采用错误发现率作为差异表达基因筛选的关键指标。利用Hyper-geometric distribution分析法，将差异基因、上调基因、下调基因分别进行GO功能注释和KEGG通路富集分析，错误发现率≤0.05，表示差异基因显著富集。

1.4 单核苷酸多态性、插入缺失标记分析方法

用软件Samtools进行去重复(Mark duplicates)[12]和软件GATK进行局部重比对(Local realignment)、碱基质量值校正(Base recalibration)[13]等处理，对乌鳢和白乌鳢共6个样本进行单核苷酸多态性和插入缺失标记检测，取两款软件得到结果的交集，作为最终的高质量突变。随后，采用SnpEff(4.3p)软件对各样本进行单核苷酸多态性和插入缺失标记注释，以获得单核苷酸多态性和插入缺失标记的位置，明确其位于哪个基因/基因间区，是否为错义突变等。对单核苷酸多态性、插入缺失标记的检测结果进行注释，从而实现DNA水平差异基因挖掘和差异基因功能注释等。

1.5 定量表达分析

为验证转录组测序结果的准确性，利用qRT-PCR分析，挑选出转录组测序结果中前10个具有差异表达的基因(干扰素诱导的GTP酶5样基因、XLOC_000762、NLR家族CARD域基因、嘌呤核苷磷酸化酶样异结构X1基因、补体C7、角鲨烯单氧酶基因、XLOC_022395、C-C基序趋化因子17样的基因、唾液酸结合凝集素、超长链酰基辅酶A合成酶样的基因)和酪氨酸代谢通路中的4个差异表达基因[酪氨酸酶相关蛋白1(tyrp1)基因、CF_GLEAN_10017395、乙醇脱氢酶5(adh5)基因、黑尿酸酶(hgd)基因]进行验证。使用与转录组测序同批次的RNA进行反转录(白乌鳢皮肤组织样本，n=15；
乌鳢皮肤组织样本，n=15)，cDNA的合成参照PrimerScript©RT reagent Kit (宝生物，中国)试剂盒说明书进行。使用软件Primer 6设计引物，以18S RNA和β-actin为内参进行qRT-PCR验证(表1)。使用SYBR©Premix Ex TaqTM (宝生物，中国)荧光定量试剂盒，利用qTOWER 3.0荧光定量Roche Light Cycler 96 (Roche，瑞士)荧光定量仪进行qRT-PCR定量验证。扩增程序：95 ℃预变性2 min；
95 ℃ 5 s，60 ℃ 30 s，进行40个循环。采用2-△△Ct法对目的基因进行相对定量分析；
使用SPSS 22软件进行差异分析，P<0.01为差异极显著，P<0.05为差异显著。

表1 qRT-PCR验证引物序列Tab.1 Sequences of qRT-PCR amplification

2.1 测序结果统计及质量评估

利用3个白乌鳢皮肤组织和3个乌鳢皮肤组织，构建了6个皮肤组织转录组文库，得到209 332 480条过滤后数据，平均过滤后数据为34 888 747，Q20碱基百分比平均为 95.78%，Q30碱基百分比平均为90.19%，碱基G和C的数量总和占总碱基数量的百分比平均为47.31%，比对率平均为58.48%(表2)。本次测序结果中Q30碱基百分比平均在90%以上，Q20碱基百分比平均在95%以上(表2)。试验结果表明，本次测序结果质量较好。此外，为提高本试验数据的准确性和可靠性，笔者分析了3个白乌鳢和3个乌鳢皮肤组织之间的相关性，结果表明，白乌鳢和乌鳢皮肤的RNA文库中3个生物学重复之间的相关性很高(图2)。这说明笔者利用转录组技术测定获得的数据准确性和可靠性较高，可用于后续分析。

图2 6个文库相关性分析Fig.2 Clustering figure of correlation analysis of 6 libraries

表2 转录组测序数据评估统计Tab.2 Summary of RNA-seq data assessment

2.2 差异表达基因分析

通过R语言中的edgeR程序，计算得到694个差异表达基因(|log2(fold change)|≥2，P<0.01)，相对于乌鳢，白乌鳢皮肤中有223个表达上调基因和471个表达下调基因(图3)。

图3 乌鳢与白乌鳢皮肤差异表达基因火山图(a)和柱形图(b)Fig.3 Volcano plots (a) and column chart (b) of differentially expressed genes between the albino norther snakehead C. argus var. and norther snakehead C. argus

2.3 差异表达基因功能分析

对694个差异表达基因进行GO功能注释，共注释到1241个GO条目上(图4)，其中：594条注释到生物过程，包含了细胞进程和代谢过程等；
387条注释到细胞组分，包括细胞器和细胞联接等；
260条注释到分子功能，包括信号传感器活性和转运活性等。注释的3个一级功能富集最多的GO条目分别为细胞过程(n=97)、膜(n=84)和结合功能(n=124)。

图4 差异表达基因的GO富集分析结果Fig.4 The results of GO functional enrichment of the differentially expressed genes

对白乌鳢上调差异表达基因(223个)进行KEGG通路富集，最显著的前20条通路见图5，包括吞噬体、细胞间受体相互作用、蛋白质消化吸收、氨基糖和核苷酸糖代谢及HIF-1信号通路等。白乌鳢下调的471个差异表达基因富集最显著的前20条KEGG通路见图6，包括细胞黏附分子通路、凝血与补体级联反应信号通路、视黄醇代谢通路、细胞色素P450外源物代谢通路和酪氨酸代谢等。

图5 白乌鳢上调差异基因KEGG富集结果Fig.5 The results of KEGG enrichment of differentially up-regulated genes in albino northern snakehead C. argus var.

图6 白乌鳢下调差异基因KEGG富集结果Fig.6 The results of KEGG enrichment of differentially down-regulated genes in albino northern snakehead C. argus var.

此外，转录组测序结果显示，酪氨酸代谢通路(Ko00350)为显著的富集通路。此通路上共检测到37个基因，比较白乌鳢和乌鳢皮肤组织基因的转录水平发现，这37个基因里只有tyrp1、CF_GLEAN_10017395、adh5、hgd基因这4个基因为差异表达基因。

2.4 单核苷酸多态性、插入缺失标记分析

利用Samtools和GATK检测到纯合型单核苷酸多态性位点最大值为12 152，最小值为5443，杂合型单核苷酸多态性位点最大值为54 382，最小值为5588；
相较于乌鳢，白乌鳢纯合型单核苷酸多态性位点数目多，而杂合型单核苷酸多态性位点数少(表3)。此外，纯合型插入缺失标记位点最大值为2003，最小值为916，杂合型插入缺失标记位点最大值为1858，最小值为1072；
相较于乌鳢，白乌鳢纯合型插入缺失标记位点数目多，杂合型插入缺失标记位点数少(表4)。

表3 乌鳢及白乌鳢转录组测序数据中单核苷酸多态性位点数量统计Tab.3 Summary of single nucleotide polymorphism (SNP) identified in the RNA-seq data of norther snakehead C. argus and albino northern snakehead C. argus var.

表4 乌鳢及白乌鳢转录组测序数据中插入缺失标记数量统计Tab.4 Summary of INDEL identified in the RNA-seq data of northern snakehead C. argus and albino northern snakehead C. argus var.

基于转录组单核苷酸多态性和插入缺失标记的结果，在经过筛选后的6578个基因里共找到18 767个单核苷酸多态性位点，其中10.03%属于非同义编码位点，33.31%为下游位点，16.89%为上游位点，2.88%的基因间位点，24.41%为同义编码位点，12.32%为内含子位点(图7a)。在差异表达基因中找到3个基因符合筛选标准，共发现5个突变位点(表5)，其中hgd、adh5基因属于酪氨酸代谢通路里的差异表达基因，表明这2个单核苷酸多态性可能会造成黑色素合成受阻，从而导致白乌鳢皮肤白化。另外，在符合筛选标准的804个基因中找到917个插入缺失标记位点，其中0.87%属于移码位点，15.38%为内含子位点，3.49%为基因间位点，51.69%为下游位点，28.57%为上游位点(图7b)。将804个基因与色素沉淀及酪氨酸代谢相关基因比对，未找到符合条件的插入缺失标记位点。

图7 单核苷酸多态性(a)和插入缺失标记(b)的分布Fig.7 The distribution of SNPs (a) and INDELs (b)

表5 突变差异基因Tab.5 Mutation of DEGs

2.5 转录组数据验证

qRT-PCR结果显示，选取的转录组测序结果中前10个具有差异表达的基因以及在酪氨酸代谢通路中4个差异表达基因，在乌鳢(n=15)与白乌鳢(n=15)皮肤组织中表达水平趋势与转录组数据相吻合(图8)，表明转录组结果是真实可信的。

图8 转录组分析获得的14个差异表达基因的qRT-PCR验证结果Fig.8 Results of the qRT-PCR verification from 14 differentially expressed genes obtained by transcriptome analysis1.干扰素诱导的GTP酶5样基因；
2.XLOC_000762；
3.NLR家族CARD域基因；
4.嘌呤核苷磷酸化酶样异结构X1基因；
5.补体C7基因；
6.角鲨烯单氧酶基因；
7.XLOC_022395；
8.C-C基序趋化因子17样的基因；
9.唾液酸结合凝集素基因；
10.超长链酰基辅酶A合成酶样的基因；
11.酪氨酸酶相关蛋白1基因；
12.CF_GLEAN_10017395；
13.乙醇脱氢酶5基因；
14.黑尿酸酶基因.1.interferon-inducible GTPase 5-like gene;2.XLOC_000762;3.NACHT gene;4.purine nucleosidephosphorylase-like isoform X1 gene;5.complement component C7 gene;6.squalene monooxygenase gene;7.XLOC_022395;8.C-C motif chemokine 17-like gene;9.sialic acid-binding lectin-like gene;10.very long-chain acyl-CoA synthetase-like gene;11.tyrp1 gene;12.CF_GLEAN_10017395;13.adh5 gene;14.hgd gene.

白乌鳢作为我国名贵的经济鱼类，皮肤颜色呈白色或灰白色，这种漂亮的体色常常成为其市场价值的判断标准。为提供更多遗传信息推动白乌鳢遗传育种工作，笔者利用转录组测序技术从分子角度探究乌鳢和白乌鳢的皮肤颜色差异。受外界环境(光照、温度和摄食)和自身遗传学(基因)等多种因素的影响，体色变异在鱼类发育过程中较为常见[14]，已有研究证明，遗传因素是鱼类体色变化的主要决定性因素，由许多基因构成的调控网络共同控制色素细胞的发育和形成，最终影响鱼类体色变化[5]。随着测序技术的不断进步，学者们利用转录组测序技术获得大量测序数据，从而可以更深入地从遗传学的角度揭示鱼类皮肤颜色调控机制[12]。

3.1 差异表达基因分析

刘肖莲等[15]对透明草金鱼(Carassiusauratus)的研究表明，透明草金鱼与不透明草金鱼皮肤转录组共存在181条差异表达基因，主要富集在MAPK通路、嘌呤代谢通路、内质网的蛋白质加工通路，推测草金鱼皮肤透明机制的调节很可能与抑制黑色素形成以及虹彩细胞中鸟嘌呤合成、色素颗粒转移过程有关。而史东杰等[16]对三色锦鲤(Cyprinuscarpiokoi)肤色调控基因的研究表明，采取墨色、红色、白色部位皮肤样品测序并两两比较，获得的差异表达基因主要富集在细胞过程和环境信息处理通路，说明三色锦鲤的皮肤生长代谢途径可能与细胞增殖分化相关[16]。此外，普通鲤(C.carpio)和金边鲤(C.carpiovar.jinbian)皮肤转录组数据显示，2475个差异基因主要富集在半胱氨酸代谢途径、cGMP-PKG信号通路和心肌细胞途径，表明这些代谢途径可能与金边鲤皮肤颜色变异有关[17]。然而，通过对同属鲤属的瓯江彩鲤(C.carpiovar.color)进行研究，发现影响其皮肤不同颜色的主要代谢通路为嘌呤代谢通路、糖酵解与糖代谢合成通路和氧化磷酸化通路[18]。由此可见，鱼类体色形成的机制十分复杂，关于鱼类体色沉积的调控机制尚不清楚，相关研究仍然十分匮乏。白乌鳢是乌鳢的白化变异种[19]，因其体色特殊极具观赏价值。本研究中，对乌鳢和白乌鳢皮肤组织进行转录组测序得到694个差异表达基因，结合上述研究，可以推测这694个差异表达基因可能参与皮肤颜色的调控。

3.2 差异表达基因功能分析

将得到的差异表达基因进行功能富集，发现其主要与能量代谢、细胞联接及酪氨酸代谢通路有关，尤其是酪氨酸代谢通路被显著富集。在脊椎动物中，酪氨酸代谢通路已被证实与黑色素合成相关[20]。酪氨酸酶(Tyr)、酪氨酸酶相关蛋白1(Tyrp-1)、酪氨酸酶相关蛋白2(Tyrp-2)这3种限制酶是催化真黑色素形成的关键酶，三者协同调控合成黑色素[21]。在本研究中，酪氨酸代谢通路上的tyrp1基因表达量显著上调，tyrp1能抑制黑色素细胞死亡，与鱼类体表色素合成紧密相关。在褐牙鲆(Paralichthysolivaceus)中也发现了与本研究结果类似的结果，正常褐牙鲆皮肤中tyrp1a基因表达量远高于白化褐牙鲆[22]；
tyrp1a、tyrp1b基因为tyrp1基因的同源基因，主要存在于绝大多数的硬骨鱼类中[23]。湖栖鳍虾虎鱼(Gobiopteruslacustris)皮肤转录组测序结果分析显示，酪氨酸代谢通路上的tyr、tyrp1、tyrp2基因表达量均显著下降，这很可能是造成湖栖鳍虾虎鱼皮肤呈透明状的原因[24]。Zhou等[9]研究指出，乌鳢幼鱼皮肤中的与黑色素合成紧密相关的tyr基因的表达量均高于白乌鳢幼鱼；
但在本研究中，成年白乌鳢和成年乌鳢的tyr、tyrp2基因表达量无显著变化。tyr基因家族作为调节黑色素合成的关键基因，其表达量在鱼的不同生长阶段存在显著差异，这可能是幼鱼的色素模式逐渐向成鱼的色素模式过渡所导致的[25]。类似的研究结果在虹鳟(Oncorhynchusmykiss)和金鳟(黄色突变型虹鳟)中也有报道，虹鳟的tyr基因的表达量在早期发育阶段显著高于金鳟，然而在12月龄时虹鳟和金鳟tyr基因的表达量无显著差异[26]。tyr基因是影响鱼类皮肤白化的关键基因，但对鱼类皮肤白化调控的分子机制存在种间特异性[27]。

因此，本研究中乌鳢皮肤白化机制中很可能存在其他调节色素沉着的因素，而tyr、tyrp2基因可能不是影响其色素沉着的关键基因。已有研究表明，皮肤中的黑色素沉着还与黑色素小体形成、黑色素小体转运等因素密切相关，大量载有黑色素的黑色素小体被转运至邻近的皮肤角质细胞，使黑色素最终在皮肤角质细胞内沉积[28]。在对虹鳟皮肤细胞转录组的研究中发现，虹鳟皮肤上的黑色斑点区域差异表达基因主要富集在细胞联接进程上，如细胞黏附分子通路[29]。这一研究结果与本研究结果相同，表明在乌鳢中细胞黏附分子通路上的差异表达基因表达量上调也可能是影响真皮层黑色素细胞沉积的关键因素之一。然而，在本研究中，酪氨酸代谢通路上的adh5和hgd基因表达量同样也出现显著上调。乙醇脱氢酶是参与机体乙醇代谢的重要酶，已在人体上鉴定出5种，adh5基因编码乙醇脱氢酶5[30]。哺乳动物相关研究表明，乙醇脱氢酶不仅参与乙醇代谢，还参与肾上腺素、多巴胺、脂肪酸等物质的代谢[31]，而其对黑色素合成的影响还未在其他物种中报道。hgd基因编码尿黑酸酶，尿黑酸酶可以降解苯丙氨酸和酪氨酸分子[32]。

3.3 单核苷酸多态性、插入缺失标记分析

在对乌鳢和白乌鳢皮肤转录组数据突变位点分析中，共在6578个基因里找到18 767个单核苷酸多态性位点，针对酪氨酸代谢通路上的差异表达基因，笔者通过与乌鳢全基因组比对，经过单核苷酸多态性位点预测，在hgd和adh5基因上发现了错义突变，而与黑色素合成紧密相关的tyr基因在两种鱼中无显著差异。有研究报道，鱼类的白化可能受一个或多个基因影响，不同物种间存在差异[33]。研究表明：hps4基因上内含子2和外显子3交界处99 bp片段缺失可能导致斑点叉尾(Ictaluruspunctatus)白化[34]；
不仅tyrp1a基因发生突变可导致斑马鱼(Daniorerio)黑色素细胞死亡[35]，oca2基因突变也可造成斑马鱼白化[36]；
对斑马鱼胚胎的研究表明，单独敲除tyrp1的1个等位基因并不会造成色素缺失，而同时敲除2个等位基因，导致严重黑色素缺失[23]；
oca2和mclr基因的突变都可能改变丽脂鲤(Astyanax)的体色[37-38]；
有学者利用转录组测序发现，与体色紧密相关的差异基因在野生型黄颡鱼(Pelteobagrusfulvidraco)和白化型黄颡鱼皮肤中的表达量差异显著，而单核苷酸多态性、插入缺失标记分析未发现可能导致白化的突变位点[39]。因此结合本研究结果和已有报道可以推测，乌鳢白化并不是由某一个基因突变或某一个基因的表达量减少造成的，而主要由两种机制共同进行调控：(1)hgd和adh5基因上单核苷酸多态性产生错义突变导致黑色素沉着缺乏；
(2)酪氨酸代谢通路中的差异基因在白乌鳢皮肤中表达量下降，导致黑色素沉着量降低。

笔者从分子角度对乌鳢的白化机制进行解释，通过高通量测序技术对乌鳢与白乌鳢皮肤组织进行转录组测序，筛选出694个差异表达基因，通过功能注释发现，差异表达基因主要富集在酪氨酸通路上。此外，分析单核苷酸多态性发现，位于酪氨酸通路上的hgd和adh5基因产生错义突变，这可为乌鳢后续的多态性检测、群体遗传多样性分析等奠定基础，但有关遗传机制的细节还有待进一步的研究。