金鱼Carassius auratus享有中国 “国鱼”之称，是世界上最早被人工养殖的观赏鱼类之一，也是研究遗传育种和种内演化的适宜试验素材[1]。作为观赏鱼，金鱼多姿多彩的体色是观赏的亮点，而一般金鱼在幼苗阶段观赏性不强，因其呈现出野生鲫体色——青灰色或银灰色(五花、透明和紫色金鱼除外)。随着其生长会逐渐呈现出变色过程，最后才表现出具有观赏性的多彩体色，如红、橙、黄、白、青、蓝、紫、墨，或由几种单色组成的复色等。在金鱼发育生长期间，变色过程较为复杂, 其变色程度和速度因品种、养殖水体环境和饵料等因素的不同而呈现出明显的差异[2]。然而，人们对金鱼变色的研究多限于生产经验的表型观察记录[2-4]，尚未有通过转录组水平研究金鱼变色的相关报道。
RNA-Seq技术是以二代高通量测序技术为基础，在特定情况下研究某个物种或特定组织的转录组，并通过转录组分析实现从整体水平对基因结构和功能的研究，进而揭示疾病发生过程和特定生物学过程中的分子机理，这也是研究基础生物学及疾病机理必用的技术方法之一。RNA-Seq主要用于差异表达基因分析、基因表达谱构建、mirco RNA表达谱构建和新基因的发现等方面[5]。目前，采用RNA-seq技术研究鱼类变色的研究主要集中在双冠丽鱼Midas cichlids[6]和红鲫Carassius auratus[7] 等。本研究中，使用RNA-seq技术对不同发育变色期的红望天眼金鱼皮肤(带鳞片)进行高通量转录本测序，得到了转录水平的基因表达信息，通过比较分析这3个不同变色组皮肤组织的差异表达基因，以及差异表达基因的GO功能富集、KEGG通路富集，以期从转录组学水平解析金鱼的变色机制，为揭示金鱼体色发育调控机制提供分子遗传学依据。
1 材料与方法
1.1 材料
试验用金鱼取自北京市农林科学院水产科学研究所小汤山水产良种繁育基地，品种为蛋种绒球红望天眼金鱼，均为全同胞家系并养殖在同一环境下。
1.2 方法
1.2.1 养殖管理与样品采集 于繁殖季节选择一雌一雄亲鱼，人工湿法受精繁殖孵化出苗后，养殖于室内玻璃水族箱(130 cm×60 cm×45 cm)中，采用曝气装置结合XY-380生化棉过滤器，室内自然温度，自然光照下每天用灯光照射约8 h，养殖用水为晾晒超过48 h的地下水。孵化出的幼鱼，15 d内投喂活丰年虾虫卵，15～60 d内投喂由混合红虫和粉料加水制成粒径约4 cm的饼状饲料，60 d后投喂粒径2.0 mm的浮性颗粒饲料。每天用吸便器吸出粪便等脏物，视鱼摄食情况灵活投喂，以保证鱼体健康为主。使用DR900水质检测仪检测水质状况，控制水体的pH为7.0～8.4，溶氧为7.70～9.70 mg/L，亚硝酸盐<0.02 mg/L，氨氮<0.15 mg/L，视水质状况灵活换水。
145日龄时，鱼平均体长约6 cm，这些全同胞个体变色程度不同呈现出3种状态。将这些全同胞个体根据不同变色程度分为3组：未变色组标记为C(caesious，青灰色)，过渡色组标记为V(variant，过渡色即青灰色和橙红色的渐近色)，全变色组标记为Y(orange red，橙红色)(图1)。从每组随机取3尾鱼，用丁子香酚麻醉，解剖剥取每尾鱼的全部皮肤(带鳞片)，随后立即放入5 mL冻存管中，置于液氮中保存待用。
图1 全变色组(Y)、过渡色组(V)和未变色组(C)红望天眼金鱼的体色
Fig.1 Body color of celestial goldfish Carassius auratus in orange-red-color group(Y),variant-color group(V)and caesious-color group(C)
1.2.2 总RNA 提取和质量检测 取3组鱼的皮肤组织样品，采用RNA提取试剂盒(Invitrogen, Carlsbad, CA,USA)，按照说明书操作分别提取每个样品的总RNA，用10 g/L琼脂糖凝胶分析所提取总RNA的完整性及是否存在DNA污染；采用核酸蛋白测定分光光度计(IMPLEN,USA)检测RNA的纯度；采用Agilent 2100 bioanalyzer 精确检测RNA的完整性和总量。
1.2.3 文库构建和测序 采用建库试剂盒(NEBNext UltraTM RNA Library Prep Kit，Illumina)，依说明书操作建库，通过磁珠富集法合成cDNA，经过对cDNA末端修复、加A 尾并连接测序接头和扩增纯化，得到红望天眼皮肤组织cDNA 文库。每个样品建立一个测序文库，对建好的文库进行质检并精准定量，以保证文库有效浓度高于2 nmol/L后，在Illumina HiSeq4000 测序平台进行测序。
1.2.4 测序原始数据质控、序列拼接和转录本注释 通过去除含N、带接头和低质量reads，实现对原始数据过滤，并计算clean data 的GC 含量、Q20和Q30。从NCBI基因组网站(https://www.ncbi.nlm.nih.gov/genome/?term=Carassius+auratus+(goldfish))下载参考基因组和基因模型注释文件,进行转录组测序数据的有参分析。采用HISAT2 2.0.5软件构建参考基因组索引，并比对配对末端clean reads 和参考基因组[8]，利用StringTie软件组装新转录本[9]，然后对新转录本进行Pfam、SUPERFAMILY、GO、KEGG等数据库注释，同时对新基因进行预测。
1.2.5 差异表达基因的比较分析 采用FeatureCounts程序，对映射到每个基因的读数进行计算[10]，并基于基因长度计算每个基因的FPKM(每百万碱基对测序的转录本序列片段的每千个碱基片段的预期数量)[11]。
1.2.6 差异表达基因的富集分析 采用ClusterProfile软件，对筛选所得差异表达基因进行GO功能富集分析[12]和KEGG通路富集分析[13]，显著性富集阈值为padj<0.05，筛选出差异表达基因富集参与到变色的主要生化代谢途径及信号转导途径。
1.2.7 qPCR验证测序结果 以β-actin为内参基因，所得差异表达基因中随机选择10个基因，设计qPCR引物，采用实时荧光定量仪(MA-6000)对转录组数据进行验证。首先采用PrimeScript TM 1st stand cDNA Synthesis逆转录试剂盒，将测序剩余的每个样品总RNA 逆转录成cDNA。实际操作时，每个基因进行3个技术重复和3个生物学重复。PCR反应程序：95 ℃下预变性5 min；95 ℃下循环变性15 s，60 ℃下退火复性30 s，共进行40个循环。以β-actin为内参基因，采用2-△△Ct方法计算每个基因的相对表达量，计算两组间基因表达量倍数的log2值，并与RNA-seq测序结果进行趋势比较。
1.3 数据处理
通过DESeq 21.20.0软件计算基因和转录本表达水平[14]，并计算3组金鱼各3个重复样品基因表达的平均值，对两两组之间的差异表达基因进行比较分析。为了控制错误率，利用Hochberg和Benjamini方法调整P值，筛选标准为padj<0.05。
2 结果与分析
2.1 转录组测序数据的组装和有参比对
通过Illumina HiSeq4000测序平台对9个金鱼皮肤样品分别进行测序，过滤接头、低质量和短序列reads后得到数据，并与参考基因组进行比对(表1)。在去除接头、低质量和短序列后，共获得高质量序列385 389 146条，共产生过滤后碱基量57.8 Gb，每个样品的碱基量均达到5.85 Gb以上。每个样品的测序碱基错误率均小于0.03%，GC含量为45.52%～48.52%，Q30大于93.72%，表明测序数据质量可靠，可继续组装分析。
9个样本与参考基因组总比对率均在91.11%以上，比对到唯一位置率在82.57%以上，比对到基因组外显子区域的reads数及其占clean reads数的比例为89.68%～93.09%(表1)。
表1 红望天眼金鱼样品转录组测序数据统计
Tab.1 Data statistics of transcriptome sequencing of the red celestial goldfish samples
样品sample原始序列数raw reads过滤后序列数clean reads错误率/%error rateQ30/%总比对率/%total map比对到外显子率/%exon比对到内显子率/%intron比对到基因间区率/%intergenic rateY143 557 97242 283 2040.0394.0792.2491.293.764.95Y242 102 19640 750 7060.0393.7291.1189.683.546.78Y346 254 93844 282 8380.0295.0392.3990.853.655.50V146 692 39245 190 1320.0393.8991.6890.494.045.47V246 567 75644 870 1580.0294.2393.3593.092.334.58V345 059 95243 178 2000.0294.2492.5890.734.464.81C145 294 08643 754 8820.0393.9292.8591.513.714.77C243 404 19242 067 9480.0294.1792.1491.223.765.01C340 434 80839 011 0780.0294.2492.2891.784.064.16
2.2 基于金鱼皮肤转录组测序数据的新基因预测
红望天眼金鱼是非模式物种，其基因注释信息通常不是很完善。比对完成后，寻找基因组未被注释的转录区，将比对的序列进行组装，从而挖掘该物种的新转录本或新基因。本研究中，对转录组测序数据进行了新转录本组装，以及Pfam、SUPERFAMILY、GO、KEGG等数据库注释，同时共预测新基因1 629个，并对其进行了GO、KEGG注释。
2.3 不同变色组金鱼皮肤转录组的差异表达基因
通过两两比较不同变色组中存在的差异表达基因，在C与Y组的比较中，筛选出3 955个差异表达基因，其中2 249个基因上调，其余下调；在V与C组的比较中，筛选出2 894个差异表达基因，其中1 287个基因上调，其余下调；在V与Y组的比较中，筛选出927个差异表达基因，其中658个基因上调，其余下调(图2)。Y和V 组比较中筛选出的差异基因最少，Y与C组比较中筛选出的差异基因最多，而V与C组比较中筛选出的差异基因中等。从图1可见，Y和V 组在视觉上体色更接近，Y与C组在视觉上反差最大，V组有少部分与C组体色接近，这与筛选出差异基因数目的趋势一致。
图2 3个变色组皮肤转录组两两比较的差异表达基因
Fig.2 Differential expression genes in the two groups of the skin transcriptome in three body color groups
统计不同变色组间差异基因的重叠情况，结果显示，有25个差异基因属于3组间共有差异基因，964个差异基因属于V和C组独有，有2 021个差异基因属于C和Y组独有，有466个差异基因属于V和Y组独有(图3)。
图3 3个变色组皮肤转录组两两比较的差异表达基因韦恩图
Fig.3 Venn diagram of differential expression genes of the skin transcriptome in the two groups in three body color groups
将3个不同变色组的差异表达基因取并集后作为差异表达基因集进行聚类分析，聚在一起的为表达模式相近的基因。从图4可见,V和Y组差异基因集首先聚在一起，C组与其他两组的差异基因集聚类位置较远。综合图1结果，Y和V 组在视觉上体色更接近，两者与C组相比视觉上反差更大，与聚类热图比较吻合。
图4 3个变色组皮肤转录组差异表达基因聚类热图
Fig.4 Clustering heat map of differential expression genes of the skin transcriptome in three body color groups of the red celestial goldfish
2.4 差异表达基因的GO富集分析
差异基因GO功能富集分析显示，3组共有1 417个差异表达基因得到归类注释。从GO富集分析结果中，按生物过程、细胞组分和分子功能3大类别(每类选最显著的10个Term)及差异基因上调和下调分类绘制柱状图，结果如图5～图7所示。
C和Y组GO功能富集分析显示：有260个差异基因显著富集，生物学过程(biological process，BP)中，基因富集显著的有肽交联、氨基酸跨膜转运、氨基酸转运、有机酸跨膜转运、羧酸跨膜转运、跨膜受体蛋白酪氨酸激酶信号通路、离子跨膜转运、细胞周期阻滞、免疫反应、细胞凋亡调控、细胞凋亡过程的调控、细胞程序性凋亡的调节、生理节律等；在细胞组分(cellular component, CC)中，基因富集显著的有肌球蛋白复合体、肌动蛋白细胞骨架、细胞骨架部分等；分子功能(molecular function, MF)中，基因富集显著的有转移酶活性、谷氨酰胺转移酶蛋白活性、核酸结合调节区、细胞因子间受体结合、DNA结合转录调控、转运受体蛋白酪氨酸激酶活性、有机酸跨膜转运体活性、有机阴离子跨膜转运体活性、氨基酸跨膜转运体活性、羧酸跨膜转运体活性、特异序列的DNA结合和跨膜受体蛋白激酶活性等(图5)。
图5 C和Y组金鱼差异表达基因的GO功能富集
Fig.5 GO enrichment of differentially expressed genes in C group vs Y group
V和C组GO功能富集分析显示，有162个差异基因显著富集，生物学过程中基因富集显著的有阴离子转运和动物器官发育等，细胞组分和分子功能中基因富集的显著性均不高(图6)。
图6 V和C组金鱼差异表达基因的GO功能富集
Fig.6 GO enrichment of differentially expressed genes in V group vs C group
V和Y组GO富集分析显示，有90个差异基因显著富集，生物学过程中基因富集的显著性不高，细胞组分中基因富集显著的有肌球蛋白复合体、肌动蛋白细胞骨架和细胞骨架部分等，分子功能中基因富集显著的有肌动活动、蛋白质谷氨酰胺转移酶活性、转移酶活性、电压门控氯离子通道活性、离子通道活动、氯离子通道活动和电压门控阴离子通道活性等(图7)。
图7 V和Y组金鱼差异表达基因的GO功能富集
Fig.7 GO enrichment of differentially expressed genes in V group vs Y group
2.5 差异表达基因的KEGG通路富集分析
将3个不同变色组金鱼的差异表达基因分别比对到KEGG 数据库中，进行KEGG通路富集分析。选取富集最显著的20 个通路条目分别绘制KEGG 富集散点图，结果如图8～图10所示。C和Y组差异表达基因主要富集在FoxO信号通路、C型凝集素受体信号通路、mTOR信号通路、胰岛素信号通路、细胞因子间受体相互作用、toll样受体信号通路、P53信号通路、细胞衰老和黑色素生成等通路(图8)。V和C组差异表达基因主要富集在C型凝集素受体信号通路、P53信号通路、细胞因子间受体相互作用、toll样受体信号通路、mTOR信号通路、烟酸盐与烟酰胺代谢、果糖与甘露糖代谢、FoxO信号通路、NOD样受体信号通路、黑色素生成和细胞衰老等通路(图9)。V和Y组差异表达基因主要富集在糖酵解/糖异生、酪氨酸代谢和糖代谢信号通路等(图10)。这些结果反映了红望天眼金鱼变色过程中，可能由多种信号通路共同参与。
图8 C和Y组金鱼差异表达基因的KEGG富集
Fig.8 KEGG pathway enrichment of differentially expressed genes in C group vs Y group
图9 V和C组金鱼差异表达基因的KEGG富集
Fig.9 KEGG pathway enrichment of differentially expressed genes in V group vs C group
图10 V和Y组金鱼差异表达基因的KEGG富集
Fig.10 KEGG pathway enrichment of differentially expressed genes in V group vs Y group
2.6 3个变色组的共有表达基因
本研究中，对不同变色组间差异表达基因的重叠情况进行统计，结果有25个差异表达基因属于3组间共有(表2)。查阅这些基因相关研究报道和综述，主要分为3类：1)调控细胞增殖分化，抑制细胞增殖和促进凋亡的基因，主要包括BTG2、生理节奏相关转录因子、丝氨酸/苏氨酸蛋白激酶转录因子、原癌基因c-Fos、P53、真核延伸因子2激酶(e EF2K)、锌指蛋白等；2)与炎症反应、适应性有关的基因，主要包括egl 9同源物3样、磷酸甘露糖变位酶(manB)、类似诱导型一氧化氮合酶、FAM13A蛋白等；3)与色素形成有关的基因，主要包括二羟基吲哚羧酸氧化酶等。从这25个共有差异表达基因可以看出，大部分与细胞增殖分化和凋亡有关，小部分与免疫炎症反应和适应性有关，仅少数与黑色素形成有关。
表2 3个变色组鱼共有差异表达基因
Tab.2 Shared differentially expressed genes in three groups of the goldfish
序号No.基因号gene IDFPKM基因读数 FPKM gene reading变色后Yafter discoloration变色中Vin discoloration变色前Cbefore discoloration基因描述gene description1LOCnovel.4370.345 4921.934 8290—2LOC1130947043.000 03011.043 4100.694 218肌球蛋白重链,快肌样,转录变体X2 myosin heavy chain,fast skeletal muscle-like, transcript variant X23LOC1131176782.755 13414.703 75090.414 690egl 9同源3样,转录变体X1 egl nine homolog 3-like, transcript variant X14LOC1130677960.502 8991.392 2826.922 527磷酸甘露糖变位酶1样,转录变体X2 phosphomannomutase 1-like, transcript variant X25LOC1130607203.996 4929.388 30739.810 790锌指蛋白395样 zinc finger protein 395-like6LOC1130439491.517 0882.931 2938.958 645家族序列相似性13A样蛋白 protein FAM13A-like7LOC1130442000.058 8120.588 2587.013 0045,6-二羟基吲哚-2-羧酸氧化酶样 5,6-dihydroxyindole-2-carboxylic acid oxidase-like8LOC1130996505.332 05412.743 35052.114 410锌指蛋白395样,转录变体X1 zinc finger protein 395-like,transcript variant X19LOC1130609220.509 4452.945 95222.104 700egl 9同源3样,转录变体X1 egl nine homolog 3-like,transcript variant X110LOC11309663059.551 02017.675 0807.905 213BTG2样蛋白 protein BTG2-like11LOC1131065930.095 8231.452 1218.591 9405,6-二羟基吲哚-2-羧酸氧化酶样 5,6-dihydroxyindole-2-carboxylic acid oxidase-like12LOC1131200023.784 6848.643 35823.749 840生理节奏相关的转录抑制因子 circadian-associated transcriptional repressor-like13LOC1130466932.398 3234.070 5488.621 717丝氨酸/苏氨酸蛋白激酶转录因子 serine/threonine-protein kinase WNK4-like, transcript variant X114LOC11309731569.920 90019.736 7505.632 924原癌基因 c-Fos样proto-oncogene c-Fos-like15LOC1131152801.018 9532.856 5406.970 940egl 9同源物2样 egl nine homolog 2-like16LOC1130385625.598 3101.656 7080.474 244BICD家族样运输适配器2 BICD family-like cargo adapter 217LOC1130443759.359 76617.852 01032.669 040胰岛素受体底物2样 insulin receptor substrate 2-like18LOCnovel.8080.829 2161.517 8593.707 884P53四聚化合物底物 P53 tetramerisation motif19LOC1130766170.346 5315.397 9440.930 242肌球蛋白重链,快肌样 myosin heavy chain, fast skeletal muscle-like20LOC1130563138.279 86320.474 16046.301 230真核延伸因子2激酶样,转录变体X6 eukaryotic elongation factor 2 kinase-like, transcript variant X621LOC1130441731.761 5724.186 5058.205 859锚蛋白重复结构域蛋白37样 ankyrin repeat domain-containing protein 37-like22LOC11311484122.110 32044.977 63078.469 800大中性氨基酸转运体小亚基4样 large neutral amino acids transporter small subunit 4-like23LOC1130787553.221 6751.349 2830.494 107类似诱导型一氧化氮合酶 nitric oxide synthase inducible-like24LOC1131099442.371 2325.434 13910.966 180热休克因子蛋白1样 heat shock factor protein 1-like25LOC11305508452.830 290100.095 600164.365 800抑制蛋白结构域蛋白2样,转录变异X2 arrestin domain-containing protein 2-like,transcript variant X2
2.7 测序结果的qPCR验证
从筛选出的差异表达基因中随机选择10个基因，经qPCR分析显示，基本与本研究中RNA-seq转录本测序得到的基因表达高低趋势相吻合，说明本研究转录组测序数据较为可靠，且对不同组间的差异表达基因能较准确地检测并筛选出来(图11)。
图11 不同组间差异表达基因在RNA-seq和q-PCR中的表达趋势比较
Fig.11 Trends of differentially expressed genes in RNA-seq and q-PCR in different groups
3 讨论
3.1 红望天眼金鱼转录组组装
本研究中，采用Illumina HiSeq4000 测序平台对红望天眼金鱼转录组文库进行高通量测序，共得到过滤后碱基量57.8 Gb, 样品Q30碱基百分比均大于93.72%，比多个已发表的鱼类转录组测序质量要高，包括斑马鱼Danic rerio(91.93%)[15]、草金鱼Carassius auratus(93.16%)[16]等。与参考基因组比对，各样品的clean reads总比对率在91.11%以上，且通过qPCR验证基本与测序结果一致，这说明本研究中测序质量较为可靠。
3.2 红望天眼金鱼不同变色组差异表达基因数目与体色的关系
本研究中，红望天眼金鱼变色前(C)为青灰色，变色后(Y)全身为橙红色，视觉差异较大，而过渡色组(V)体色介于两者之间，更偏向于变色后体色(图1)。从筛选出的差异表达基因数目看，金鱼变色过程中体色视觉上差异越大，所筛选出的差异表达基因越多；对其聚类分析显示，V和Y组差异基因集聚在一起，而C组与其他两组差异基因集的聚类位置较远，这也与视觉上体色差异大小趋势一致。本研究结果与对双冠丽鱼[6]不同变色组的转录本差异表达基因聚类分析结果相同。由此可见，总体上在金鱼变色中体色差异越大，相应筛选出的差异表达基因也越多，参与体色变化的调控基因差异也越大。
3.3 红望天眼金鱼不同变色组差异表达基因的GO功能富集和KEGG通路富集分析
3.3.1 GO功能富集 本研究中，通过对3组差异基因的GO功能富集可知，C和Y组比较富集显著的差异表达基因较多，V和Y组比较，富集显著的较少，这与视觉上金鱼变色期差异大小趋势一致。C和Y组比较，GO富集条目在生物学过程、细胞组分和分子功能上均有显著分布，可见变色前和变色后体表色素细胞分布差别极大；V和C组比较，GO富集条目显著集中在生物学过程，细胞组分和分子功能中基因富集差别不大；V和Y组比较，GO富集条目显著集中在细胞组分和分子功能上。本研究中GO富集到的细胞结构与虹彩细胞的细胞结构吻合，这是因为虹彩细胞中色素颗粒的定向运动与微管、肌动蛋白有关[17]。本研究中金鱼差异基因的GO富集结果与刘肖莲等[16]对透明草金鱼RNA-seq差异基因GO分析部分条目较为一致，也富集到肌球蛋白复合体、肌动蛋白细胞骨架、细胞骨架部分中。而Henning等[6]对双冠丽鱼3种不同体色阶段的转录组研究显示，部分差异表达基因富集在免疫早期反应中，与本研究中GO富集到免疫反应的结果一致。本研究中GO富集的跨膜受体蛋白酪氨酸激酶信号通路与黑色素形成有关[18]。
3.3.2 KEGG通路富集 本研究中，3组变色程度不同的红望天眼金鱼差异表达基因KEGG通路富集到的通路主要与细胞增殖分化和凋亡、免疫反应及色素等相关。涉及细胞增殖分化和凋亡的信号通路，如FoxO、BTG2、P53、mTOR、Hedgehog(Hh)和酪氨酸激酶等；涉及免疫反应的信号通路，如FoxO、C型凝集素受体信号通路和toll样受体信号通路等；涉及色素相关的信号通路，如酪氨酸通过一系列氧化还原反应生成的黑色素[18]；涉及糖代谢的通路，如胰岛素信号通路、糖异生及果糖与甘露糖代谢等，而有报道显示，糖代谢与虹彩细胞有关[19]。本研究中，红望天眼金鱼变色过程的差异表达基因KEGG通路富集与红鲫[7]有部分相同，红鲫变色过程差异表达基因也有富集在黑色素信号通路、酪氨酸代谢信号通路、脂肪信号通路、细胞自噬信号通路和细胞凋亡信号通路等。金鱼起源于野生鲫，由红黄色野生鲫演化成金鲫[1]，后经人工饲养繁育而成。有报道显示，红鲫的演变与红色金鱼相似，均经过变色过程，初为银灰色，经60 d左右变为红色[3]。由此推测，调控红望天眼金鱼与红鲫变色的生理和分子机制也有相同之处。本研究中，KEGG通路富集结果与对双冠丽鱼变色研究中富集到黑色素合成信号通路和免疫早期反应的结果一致。这些结果反映了红望天眼金鱼在变色过程中，与红鲫、双冠丽鱼变色过程中所涉及的很多通路相同，推测可能由多种信号通路共同参与完成。
3.4 红望天眼金鱼变色机制分析
3.4.1 鲫属鱼类变色过程观察 鱼类的变色过程比较复杂，其速度和程度可因品种、水体环境和饵料品种等有明显不同[2]。伍惠生[3]报道显示，红鲫的演变与红色金鱼相似，由银灰色经过60 d左右变为红色。也有对红白金鱼变色记载：红白金鱼45 d开始少数变色，60 d变色比例增加，120 d体色基本稳定下来，剩余的青灰色个体在较长一段时间内不会再变色，鱼苗变色均开始于鱼的腹部，逐渐到头部，再逐渐到鱼背部及尾部，最后到达全身[4]。也有研究者发现，有些黑色个体金鱼在2～3年后易褪色成为红色个体[3]。本试验中，金鱼为全同胞的红望天眼金鱼且养殖在同一环境中，一年内变色完成，鱼苗变色也是从腹部开始，慢慢到头、背和尾，最后到全身。可见，品种不同，变色速度和程度也不同，这说明变色现象受遗传影响，其本质上是受基因调控。
3.4.2 饲料对鱼类变色的影响 饲料是鱼类生长和维持正常生命活动所必需的营养来源，尤其是鱼类自身不能合成类胡萝卜素等色素物质，必须从饲料中获取。一方面，饲料中的营养对鱼类体色的影响直接通过色素细胞或色素物质发挥作用；另一方面通过影响鱼类整体生理机能，使鱼类正常生理活动出现障碍，造成鱼类体色发生改变[20]。这可以解释鱼类在营养丰富的情况下，变色速度较快，体色也出现较为明艳的现象。本研究中，红望天眼金鱼KEGG通路部分富集到胰岛素信号通路、糖异生、果糖与甘露糖等糖代谢，以及脂肪酸合成和氨基酸合成通路等，这说明变色过程中糖类等能量营养的重要性。目前，国内外学者对鱼类早期体色褪黑这一复杂生物学过程的研究相对较少，其调控机制尚无明确定论[21]。
3.4.3 鱼类变色过程中黑色素细胞的增殖与凋亡 目前，对观赏鱼体色研究大部分均从色素细胞代谢、色素沉积与形成等方面进行探讨。在鱼类中已鉴定出6 种色素细胞，包括黑色素细胞、黄色素细胞、红色素细胞、虹彩细胞、白色素细胞和蓝色素细胞[22]。对黑色素细胞的相关研究较多，对黄、红色素和虹彩细胞的研究相对较少，对鱼类早期体色褪黑现象也主要从黑色素形成和凋亡通路方面报道较多，而对其他色素细胞增殖和分化通路极少报道。红鲫[7]和双冠丽鱼[6]早期体色发育变色过程，主要由皮肤和鳞片中色素变迁，黑色素细胞的消退，以及红黄色素细胞的增加所造成的。而红色金鱼体色的演变与红鲫相似，均经历变色过程，结合本研究中对红望天眼金鱼变色过程的观察，推测红望天眼金鱼也与以上鱼类的变色过程相同。在哺乳类动物中也有褪色现象，有报道显示，有一种灰色马，刚出生时毛色是灰黑色，6～7年后逐渐褪色成白马，但皮肤终生为灰黑色，同时伴随着老年时易发生黑色素瘤的现象，这是由STX17基因的第6内含子一个4.6 kb序列的加倍所致[23]，但该研究并未明确解释具体原因。鱼类和哺乳类褪色现象有明显的不同，红望天眼金鱼、红鲫和双冠丽鱼变色过程中，明显可观察到鳞片和皮肤均同步褪黑变色，且在红鲫和双冠丽鱼变色过程中，观察到黑色素细胞逐渐凋亡直到消失[6-7]。有学者认为，可能是黑色素颗粒在细胞发育过程中无法继续合成，从而导致黑色素细胞消失，这种现象可能与红鲫的生长发育过程中各种信号通路调控黑色素细胞有关[7]。
3.4.4 红望天眼金鱼变色机制推测 本研究中，对红望天眼金鱼3个不同变色时期的RNA-Seq数据进行3组共有差异表达基因、GO和KEGG富集分析，结果均涉及细胞增殖分化与凋亡、免疫反应和色素形成等相关信号调控。通过在实际生产养殖过程中观察，红望天眼金鱼体色越临近变色，体色发黑越浓，随后黑色逐渐褪去，体色变为橙红，结合红鲫色素细胞变迁规律[7]，笔者初步推测，红望天眼金鱼变色的机制是：红金鱼变色初期黑色素细胞大量增殖，随后开始逐渐凋亡，其他色素细胞则同时逐渐分化增殖，直至完全取代黑色素细胞，期间发生免疫炎症反应，直至变色完成。本研究中，KEGG富集到的通路部分与细胞增殖分化和凋亡有关，如FoxO、BTG2、P53、mTOR、Hh等信号通路，可能这些通路共同参与了调节黑色素细胞增殖凋亡和其他色素细胞的增殖分化过程。在红鲫体色细胞变迁中，发育初期黑色素细胞和其他色素细胞是共存于鱼鳞片和皮肤中的，变色开始后，仅黑色素细胞逐渐发生凋亡[7]。本试验中，红望天眼金鱼发育初期体色也是同红鲫相同为青灰色，不同的是越临近变色，体色发黑越浓，随后黑色逐渐褪去，体色变为橙红。由此推测，红望天眼金鱼极可能也是发育初期黑色素细胞与其他细胞共存，临近变色黑色素细胞增殖更多，随后开始大量凋亡，同时其他色素细胞增殖分化逐渐取代黑色素细胞。
临近红望天眼金鱼变色时是何原因导致了鳞片和皮肤黑色素细胞大量增殖且随后凋亡直至消失，变色过程中为何会伴随着免疫反应，这些问题还需进一步分析。蒋燕玲[24]报道显示，大部分的硬骨鱼类黑色素细胞最先发育，黄色素细胞发育比黑色素细胞晚，红色素细胞发育又晚于黄色素细胞，虹彩色素细胞发育晚于黑色素细胞和黄色素细胞。据此推测，可能是黑色素细胞先发育后，随着鱼类机体整体生长发育的进行，相对于其他色素细胞，营养物质优先供应黑色素细胞生成通路，使得黑色素细胞增殖过快，黑色素含量过高，引发了免疫反应。而鱼类体内适当的黑色素可以保护皮肤免受辐射, 还可抵抗环境污染等因素带来的危害, 但体内过多或过少的黑色素会导致一些疾病如霉病等发生[25]。本研究中，经过RNA-seq测序分析表明，变色前、中、后共有差异表达基因中有5，6-二羟基吲哚羧酸氧化酶，而该酶是黑色素合成的中间体, 可将5，6-二羟基吲哚羧酸(DHICA)氧化脱羧成5，6-二羟基吲哚(DHI)，DHI是一种比DHICA毒性更强的中间体,其浓度过高时可降低体内黑色素细胞数量[26]，从而促进黑色素细胞凋亡。
本研究中，从红望天眼金鱼3个变色组25个共有差异表达基因中挑选参与黑色素形成、细胞增殖与凋亡的相关基因：基因LOC113044200和LOC113106593均涉及二羟基吲哚羧酸氧化酶，该酶在黑色素的形成中具有十分重要的作用[18]；基因LOC113056313与e EF2K有关，e EF2K的表达水平或者活性的改变影响细胞周期、分化、凋亡及自噬等生理进程[27]；基因LOC113120002涉及生理节奏相关的转录抑制因子，而生理节奏相关转录因子参与代谢、内分泌、细胞增殖分化与凋亡及免疫等方面[28]；基因LOC113097315与c-Fos原癌基因有关，而c-Fos作为转录调节因子是核内重要的癌蛋白(Fos)，能够瞬时诱导触发其他特异基因的表达,其在细胞的分化、增殖及凋亡生理过程中起到重要的调节作用[29]。由此可见，基因LOC113044200和LOC113106593参与黑色素形成，基因LOC113056313、LOC113120002和LOC11309-7315参与细胞增殖和凋亡，这5个基因可作为红金鱼变色的候选基因进行后续研究。
4 结论
1)采用RNA-Seq技术对红望天眼金鱼不同变色组转录组文库进行高通量测序，共得到过滤后碱基量57.8 Gb，构建新转录本并预测新基因1 629个，通过差异表达分析筛选出25个共有差异表达基因，通过GO功能富集和KEGG通路富集分析，差异表达基因主要富集到细胞增殖分化与凋亡、免疫反应和黑色素形成等相关通路。
2)初步推测红望天眼金鱼变色的生理机制是红望天眼金鱼变色初期黑色素细胞大量增殖，期间发生免疫炎症反应，随后开始逐渐凋亡，同时其他色素细胞逐渐分化增殖，直至完全取代黑色素细胞，完成变色过程。本试验结果有利于人们进一步理解红金鱼体色变迁的原因，同时为下一步开展调控金鱼变色的功能基因试验和金鱼变色的分子机理研究提供科学参考。
[1] 陈桢.金鱼的家化与变异[M].北京:科学出版社,1959:1-45,133-135.
CHEN Z.The domestication and variation of goldfish[M].Beijing:Science Press,1959:1-45,133-135.(in Chinese)
[2] 李恕,谢柏桐.金鱼的变色培育技术[J].动物学杂志,1993,28(4):48-49.
LI S,XIE B T.The color-changing breeding technology of goldfish[J].Chinese Journal of Zoology,1993,28(4):48-49.(in Chinese)
[3] 伍惠生.中国金鱼[M].天津:天津科学技术出版社,1987.
WU H S.Chinese goldfish[M].Tianjin:Tianjin Science and Tech-nology Press,1987.(in Chinese)
[4] 邱洋洋,陆海燕,宓华明,等.金鱼红白虎头和红白蝶尾及其杂交F1代的体色和体型的分析[J].大连海洋大学学报,2011,26(3):209-214.
QIU Y Y,LU H Y,MI H M,et al.Comparison of pigmentation and morphology in goldfish with tiger-type head or butterfly-type tail and their hybrids[J].Journal of Dalian Ocean University,2011,26(3):209-214.(in Chinese)
[5] 李晓敏.锦鲤红白皮体色表达谱分析及相关基因甲基化表达模式的研究[D].大连:大连海洋大学,2015.
LI X M.Gene expression variation and epigenetic modification associated with skin color in koi carp(Cyprinus carpio var Koi)[D].Dalian:Dalian Ocean University.2015.(in Chinese)
[6] HENNING F,JONES J C,FRANCHINI P,et al.Transcriptomics of morphological color change in polychromatic Midas cichlids[J].BMC Genomics,2013,14:171.
[7] 张永勤.红鲫体色发育的分子调控机制研究[D].长沙:湖南师范大学,2018.
ZHANG Y Q.Molecular regulation mechanism of body color development in red crucian carp(Carassius auratus,red var.)[D].Changsha:Hunan Normal University,2018.(in Chinese)
[8] MORTAZAVI A,WILLIAMS B A,MCCUE K,et al.Mapping and quantifying mammalian transcriptomes by RNA-Seq[J].Nature Methods,2008,5(7):621-628.
[9] PERTEA M,PERTEA G M,ANTONESCU C M,et al.StringTie enables improved reconstruction of atranscriptome from RNA-seq reads[J].Nature Biotechnology,2015,33(3):290-295.
[10] LIAO Y,SMYTH G K,SHI W.FeatureCounts:an efficient general purpose program for assigning sequence reads to genomic features[J].Bioinformatics,2014,30(7):923-930.
[11] BRAY N,PIMENTEL H,MELSTED P,et al.Near-optimal probabilistic RNA-seq quantification[J].Nature Biotechnology,2016,34(5):525-527.
[12] YOUNG M D,WAKEFIELD M J,SMYTH G K,et al.Gene ontology analysis for RNA-seq:accounting for selection bias[J].Genome Biology,2010,11(2):R14.
[13] KANEHISA M,GOTO S.KEGG:kyoto encyclopedia of genes and genomes[J].Nucleic Acids Research,2000,28(1):27-30.
[14] LOVE M I,HUBER W,ANDERS S.Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J].Genome Biology,2014,15(12):550.
[15] 康桂英,董武,杨景峰,等.基于RNA-Seq技术分析对乙酰氨基酚对斑马鱼转录组的影响[J].中国兽医学报,2019,39(10):2033-2039.
KANG G Y,DONG W,YANG J F,et al.Analyzing effect of acetaminophen on zebrafish transcriptome through RNA-Seq technique[J].Chinese Journal of Veterinary Science,2019,39(10):2033-2039.(in Chinese)
[16] 刘肖莲,李春艳,白晓慧,等.透明草金鱼的皮肤转录组分析[J].水产科学,2021,41(1):110-115
LIU X L,LI C Y,BAI X H,et al.Transcriptome analysis of transparent skin in grass goldfish Carassius auratus[J].Fisheries Science,2021,41(1):110-115.(in Chinese)
[17] 刘晓东,陈再忠.鱼类色素细胞及体色调控[J].水产科技情报,2008,35(1):13-18.
LIU X D,CHEN Z Z.Pigment cells and body color regulation in fish[J].Fisheries Science & Technology Information.2008,35(1):13-18.(in Chinese)
[18] SIMON J D,PELES D,WAKAMATSU K,et al.Current challenges in understanding melanogenesis:bridging chemistry,biological control,morphology,and function[J].Pigment Cell & Melanoma Research,2009,22(5):563-579.
[19] HIGDON C W,MITRA R D,JOHNSON S L.Gene expression analysis of zebrafish melanocytes,iridophores,and retinal pigmented epithelium reveals indicators of biological function and developmental origin[J].PLoS One,2013,8(7):e67801.
[20] 叶元土,蔡春芳.鱼类营养与饲料配制[M].北京:化学工业出版社,2013:1-2.
YE Y S,CAI C F.Fish nutrition and feed preparation[M].Beijing:Chemical Industry Press,2013:1-2.(in Chinese)
[21] 陈旭东,邬国强,宋红梅,等.橘色双冠丽鱼体色相关基因mitf的结构及表达调控特性[J].渔业科学进展,2022,43(2):107-118.
CHEN X D,WU G Q,SONG H M,et al.Expression analysis of mitf gene relating to body color variation in Amphilophus citrinellus[J].Progress in Fishery Sciences,2022,43(2):107-118.(in Chinese)
[22] VOLKENING A,SANDSTEDE B.Iridophores as a source of robustness in zebrafish stripes and variability in Danio patterns[J].Nature Communications,2018,9(1):3231.
[23] PIELBERG G R,GOLOVKO A,SUNDSTRÖM E,et al.A cis-acting regulatory mutation causes premature hair graying and susceptibility to melanoma in the horse[J].Nature Genetics,2008,40(8):1004-1009.
[24] 蒋燕玲.橘色双冠丽鱼体色发育变化及体色相关基因TYR的克隆与表达研究[D].上海:上海海洋大学,2016.
JIANG Y L.The development of body color of Amphilophus citrinellus and the cloning and analysis of body color-related gene TYR expression[D].Shanghai:Shanghai Ocean University,2016.(in Chinese)
[25] 孟超,徐皓,黄超,等.鱼类体色相关功能基因的研究进展[J].湖南文理学院学报(自然科学版),2020,32(3):30-35.
MENG C,XU H,HUANG C,et al.Research progress on color-related genes of fish[J].Journal of Hunan University of Arts and Science(Natural Science Edition),2020,32(3):30-35.(in Chinese)
[26] POPA I L,MILAC A L,SIMA L E,et al.Cross-talk between dopachrome tautomerase and caveolin-1 is melanoma cell phenotype-specific and potentially involved in tumor progression[J].Journal of Biological Chemistry,2016,291(24):12481-12500.
[27] 王根柱,戚欣,李静.真核延伸因子2激酶与肿瘤[J].药学学报,2015,50(7):808-813.
WANG G Z,QI X,LI J.Eukaryotic elongation factor 2 kinase and cancer[J].Acta Pharmaceutica Sinica,2015,50(7):808-813.(in Chinese)
[28] 杨柳青,石汉平.昼夜节律与肿瘤[J].肿瘤代谢与营养电子杂志,2017,4(2):338-342.
YANG L Q,SHI H P.Circadian rhythm and cancer[J].Electronic Journal of Metabolism and Nutrition of Cancer,2017,4(2):338-342.(in Chinese)
[29] 张果忠,熊杰,华晓宁.癫痫与原癌基因c-fos的研究近况[J].现代中西医结合杂志,2002,11(10):985-986,988.
ZHANG G Z,XIONG J,HUA X N.Recent research progress on epilepsy and proto-oncogene c-fos[J].Modern Journal of Integrated Chinese Traditional and Western Medicine,2002,11(10):985-986,988.(in Chinese)