Transcriptome analysis of skin color variation during and after overwintering of Malaysian red tilapia

The commercial value of red tilapia is hampered by variations in skin color during overwintering. In this study, three types of skin of red tilapia, including the skin remained pink color during and after overwintering (P), the skin changed from pink color to black color during overwintering and remained black color after overwintering (P-B), and the skin changed from pink color to black color during overwintering but recovered to pink color when the temperature rose after overwintering (P-B-P), were used to analyze their molecular mechanisms of color variation. The transcriptome results revealed that the P, P-B, and P-B-P libraries had 43, 42, and 43 million clean reads, respectively. The top 10 abundance mRNAs and specific mRNAs (specificity measure SPM > 0.9) were screened. After comparing intergroup gene expression levels, there were 2528, 1924, and 1939 differentially expressed genes (DEGs) between P-B-P and P-B, P-B-P and P, and P-B and P, respectively. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of color-related mRNAs showed that a number of DEGs, including tyrp1, tyr, pmel, mitf, mc1r, asip, tat, hpdb, and foxd3, might play a potential role in pigmentation. Additionally, the co-expression patterns of genes were detected within the pigment-related pathways by the PPI network from P-B vs. P group. Furthermore, DEGs from the apoptosis and autophagy pathways, such as baxα, beclin1, and atg7, might be involved in the fading of red tilapia melanocytes. The findings will aid in understanding the molecular mechanism underlying skin color variation in red tilapia during and after overwintering as well as lay a foundation for future research aimed at improving red tilapia skin color characteristics.


Introduction
Tilapia is one of the excellent fish species recommended by the Food and Agriculture Organization of the United Nations (Gupta and Acosta 2004). In recent years, tilapia has been widely accepted and has become an export dominant species of aquaculture in China (Pradeep et al. 2014). Red tilapia was obtained by crossing the mutant Oreochromis mossambicus with other tilapia populations such as Oreochromis niloticus and Oreochromis aureus (Li et al. 2003;Jiang et al. 2019). Red tilapia is a valuable fish due to its uniform red skin, the absence of black peritoneum, fast growth, and adaptability to any culture system, and it has a huge market in many parts of the world, such as China, Malaysia, and Thailand (Pradeep et al. 2014). Therefore, most studies on red tilapia have mainly focused on their growth and development in genetic breeding (Wardani et al. 2020;Zhu et al. 2016). However, the key issue restricting the growth of commercial red tilapia cultures is skin color variation during overwintering. Pavlidis et al. (2008) found that water temperature changed the body color by motility of chromatophore in red porgy (Pagrus pagrus), but the molecular mechanism for this change is unclear.
As one of the most diverse phenotypic characteristics in animals, coloration plays numerous adaptive functions such as camouflage, spouse choice, species identification, thermoregulation, and photoreception (Hubbard et al. 2010). Furthermore, in aquaculture species, skin pigmentation patterns can be considered a factor of economic consideration. Therefore, skin colors might play a vital role in quality parameters in certain species. Previous studies have investigated that skin color was affected by many factors, such as genetics, nutrition, physiology, and environmental factors (Jiang et al. 2014;Luo et al. 2021). Water temperature is a major environmental factor for the metabolism, development, and growth of animals (Pavlidis et al. 2008). Many animals are dark under the cold and light under the warm condition (Kats and Van 1986;Sherbrooke and Frost 1989). For example, both the dorsal and ventral skin colors of Rana chiricahuensis in low temperatures (5 ℃) were significantly darker than those exposed to 25 ℃ (Fernandez and Bagnara 1991). Red porgy (Pagrus) showed a darker dorsal skin area at low (15 ℃) water temperatures and lighter skin at 19 ℃ (Pavlidis et al. 2008). The best pigmentation levels were achieved at temperatures from 26 to 30 ℃ in goldfish (Carassius auratus) (Gouveia and Rema 2005). These researches are all focused on the physical or biochemical level. The molecular and cellular mechanisms of regulating skin color variation in fish, especially color variation during and after overwintering in red tilapia, are still unknown.
In our previous study, Illumina RNA-seq and microRNA-seq analysis were conducted on different color varieties of red tilapia (Zhu et al. 2016;Wang et al. 2018). Wang et al. (2018) indicated that the color variation during the overwintering period of red tilapia might be related to the changes in skin melanocytes and tyrosinase (TYR) activity. In our red tilapia breeding procedure, we found three kinds of changes of color, i.e., the skin remained pink color during and after overwintering (P), the skin changed from pink color to black color during overwintering and remained black color after overwintering (P-B), and the skin changed from pink color to black color during overwintering but recovered to pink color when the temperature rose after overwintering (P-B-P) ( Supplementary Fig. S1). In this study, we used RNA-Seq to analyze the transcriptional profiles of P, P-B, and P-B-P skin color of red tilapia during and after overwintering. Particularly, we attempted to screen hundreds of differentially expressed genes (DEGs), which were responsible for skin color variation. Furthermore, the signaling pathways related to color variation during and after overwintering were also examined. Finally, several DEGs were validated by quantitative real-time polymerase chain reaction (qRT-PCR). This study will not only expose the molecular mechanism underlying red tilapia skin color variation during and after overwintering, but also provide valuable genetic information for breeding pure pink color tilapia.

Sample collection
The red tilapia used in this study was obtained from the pilot experimental station of Freshwater Fisheries Research Center (FFRC), affiliated with the Chinese Academy of Fishery Sciences. The red tilapia (initial weight: 500 ± 20 g) were cultivated in a 330m 2 plastic film shed pond at water temperatures of 18 ± 1 ℃ and fed twice a day (morning and evening) during the winter. When overwintering ended in April of the next year, a few whole pink tilapias changed from pink color to black color. Then, we classified the red tilapia with body color variation and tagged them with passive integrated transponder (PIT) tags. All these fish continue to be cultured in the same pond without a plastic film shed, and the water temperature gradually rose with the ambient temperature (final water temperature: 24 ± 1 ℃). During the whole period, the fish was cultured under natural light and the range of water pH value was 7.33-7.87. The dissolved oxygen (DO) was more than 6 mg/L, and NH 4 -N was less than 0.5 mg/L in the water. Some red tilapia reversed body color from black to pink after the water temperature rose.
Skin tissues were collected from four P-B (pink changed to black) red tilapia, four P (pink unchanged) red tilapia, and four P-B-P (black return to pink) red tilapia individuals, respectively. All fresh tissue samples were immediately snap-frozen in liquid nitrogen and then stored at − 80 ℃ until use.
RNA extraction, cDNA library construction, and sequencing Total RNA was obtained from red tilapia samples using RNA TRIzol (Invitrogen, UK) according to the manufacturer's protocol, and genomic DNA was removed using DNase (New England Biolabs). RNA purity was assessed using the Nanodrop-2000 (Thermo Scientific, USA). The ratio of A260:A280 in all RNA samples was higher than 1.9, and that of A260:A230 was higher than 1.8. Total RNA integrity was then checked using a Bioanalyzer RNA 6000 Pico Kit (Agilent Technologies). Samples with an RNA Integrity Number (RIN) > 8 were retained for subsequent analysis.
A total of twelve RNA samples of three different skin colors (four samples per skin color) were prepared and used for library construction. The libraries were constructed by TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Firstly, mRNA was purified from total RNA by Poly-T oligoattached magnetic beads and then fragmented under elevated temperature. Then first-strand and second-strand cDNA were subsequently synthesized. Secondly, the doublestranded cDNA was purified for end repair, dA tailing, adaptor ligation, and DNA fragment enrichment. Finally, PCR was performed, and aimed products were purified. The final product was assessed for its size distribution using Bioanalyzer DNA High Sensitivity Kit (Agilent Technologies). Each library was conducted on the Illumina X-Ten for 2 × 150 bp pairend (PE) sequencing.
Quality control and mapping to the reference genome Quality of all raw reads was conducted by FastQC (Andrews 2014) software. At the initial filtering step, SOAPnuke (Chen et al. 2018b) was used to discard poor quality reads, including adaptor reads and lowquality reads (reads more than 50% bases with a quality value less than 5). Then the clean reads were mapped onto the reference genome independently by HISAT2 version 2.1.0 (Kim et al. 2015) with default values. And RSeQC packages (version 2.6.4)  were used to make a comprehensive evaluation of RNA-seq data after alignment, including sequencing saturation, mapped reads distribution, coverage uniformity, strand specificity, and transcript level RNA integrity.

Differential expression analysis
Based on the HISAT2 alignment BAM file, feature-Counts v1.6.2 (Liao et al. 2014) was used to estimate and quantify gene expression with default parameters so as to generate the raw read count of each RNA gene. Gene expression was normalized by reads per kilobase of exon per million reads mapped (RPKM). Finally, edgeR (Robinson et al. 2010) was used to identify the DEGs by pairwise comparisons. The difference was considered significant if the |logFC|≥ 1 and FDR (false discovery rate) ≤ 0.05.
To further understand the mRNA expression of each sample in red tilapia, the specificity measure (SPM) was introduced to analyze all screened mRNA (FPKM value ≥ 1, at least 3 samples) by PaGeFinder algorithms (Pan et al. 2012). SPM values greater than 0.9 were used as the selection criterion for specific genes. The higher the SPM value, the more the specific gene expression in the sample.

Quantitative real-time PCR analysis
Total RNAs were extracted as described above. Each RNA sample was treated by 5 × PrimeScript™ RT Master Mix (Takara, Japan) to remove residual genomic DNA and reverse transcribed into cDNA. All primer pairs (Supplementary Table S1) were designed based on the unigene sequences and then synthesized by Sangon Biotech. (Shanghai, China). Real-time PCR was performed on a CFX-96 Realtime PCR System (Bio-Rad, CA, USA) in 25 μL reactions containing 12.5 μL SYBR Premix Ex Taq II (2 ×) (Takara Bio), 1 μL each primer (10 μM), 2 μL PCR template (cDNA), and 8.5 μL of nucleasefree water. Amplification was performed with an initial denaturation at 95 ℃ for 5 min, followed by 40 cycles of 95 ℃ for 10 s, 60 ℃ for 30 s, and 72 ℃ for 30 s. All the reactions were conducted in triplicate, which six biological replicates. At the end of the PCR cycle, the relative expression was calculated using the 2 −(ΔΔCt) method with the β-actin gene as the reference control. The F-test and Shapiro-Wilk test in SPSS 25 (IBM, Chicago, IL, USA) were used for homogeneity and normality tests, respectively. Differences in gene expression between P vs. P-B skins were assayed with SPSS 25 (IBM, Chicago, IL, USA) by t-test (one-tail). Thresholds for statistical significance were set at P < 0.05 (significant) and P < 0.01 (highly significant). The values of log2 (fold change) were calculated statistically with Excel sheet software. Positive values indicated up-regulation of genes expression, and negative values represented downregulation of genes expression. Finally, the expression patterns of genes were compared with the RNA-Seq sequencing results.

Co-expression analysis of protein-protein interaction network analysis
The DEGs were imported into the Search Tool for the Retrieval of Interacting Genes (STRING, https:// string-db. org) database to obtain the protein-protein interaction (PPI) information (Szklarczyk et al. 2015). Only validated interactions with a composite score greater than 0.4 were considered significant. Cytoscape 3.6.0 software was used to construct a PPI network and count the number of nodes in DEGs (Franz et al. 2016). The node genes with node degrees above 10 were selected as the key target genes.

Overview of the RNA-Seq data
To better understand the skin color variation of red tilapia during and after overwintering, the mRNA libraries of P, P-B, and P-B-P were determined and analyzed by Illumina sequencing technology. In Table 1, we presented the Q30 percentage, To assess the quality of sequencing and reassembly, all clean reads were mapped to the Nile tilapia (Oreochromis niloticus) genome within the range of known gene annotations. We found that 90.48-94.62% of the clean reads could be mapped to the Nile tilapia reference genome (Supplementary Table S2). In particular, the percentage of multiple mapped reads and unique mapped reads for all libraries averaged 5.34% and 92.54%. In addition, to consider the total mapping rate of sequencing reads and genomes for transcriptome sequencing, we also need to understand the distribution of mapped reads. The proportion of all reads operation in the CDs area exceeded 68.14%, and the ratio of matched with the intron area was the lowest, less than 8.71% (Supplementary Table S3).

Analysis of gene expression level of the red tilapia transcriptome
The RPKM method was used to estimate gene expression. The distribution of RPKM values for each sample was shown in Supplementary Table S4. A total of 33,437 genes were identified in the skin of red tilapia, and the expressed genes accounted for more than 63.26% of the total. The number of genes with 0 ≤ RPKM ≤ 1 was the most, while the number of genes with RPKM ≥ 100 was less than 1%.

Expression profiling of mRNAs
To further understand mRNA expression of different skin colors in red tilapia, SPM analysis was conducted for each sample, in which the expressed mRNAs were filtered. 9755 mRNAs participated in SPM analysis with their mean RPKM value in each group, and 465 specific mRNAs were screened for further analysis (SPM > 0.9; Supplementary Table S1). In detail, there were 119, 294, and 53 specific mRNAs in P-B-P, P-B, and P skins, respectively. The KEGG results of specific mRNAs showed that metabolic pathways, ribosome biogenesis in eukaryotes, and oxidative phosphorylation were dominant pathways in P-B-P skin, and regulation of actin cytoskeleton, melanogenesis, tight junction, and tyrosine metabolism was dominant in P-B skin (Supplementary Table S2). Furthermore, we also analyzed the top 10 abundance mRNAs of different skin colors in red tilapia. As shown in Table 2, granulin and tat genes were abundant in P skin and other pigment-related genes were highly expressed, including oca2 and slc45a4. The melanin synthesis genes accounted for the largest proportion of abundantly expressed genes, such as tyrp1b, pmelb, tyr, pmela, and tyrp1a, while the tyrp1b gene was the most abundant in P-B skin. In addition, the baxa gene showed dominantly expression in P-B-P skin. Differential gene expression (DEGs) identified in different skin patterns In comparative transcriptome analysis, many genes showed different expression levels in three skin color samples. Under the criteria of FDR ≤ 0.05 and |logFC|≥ 1, the volcano plots of three pairwise comparisons (P-B-P vs. P-B, P-B-P vs. P, and P-B vs. P) revealed the expression trend of each pair (Fig. 1a).
We also constructed a histogram of DEGs in the three skin tissues (Fig. 1b). Compared with the P-B skin, there were 2528 DEGs in P-B-P skin, of which 1420 were up-regulated and 1108 were down-regulated. A total of 1091 DEGs were up-regulated in P-B-P skin compared with P skin, while 833 DEGs were downregulated. There were 1939 DEGs displaying greater abundance in P-B skin compared with P skin, of which 1387 DEGs were up-regulated and 552 genes were down-regulated. Among these genes, 22 DEGs were detected as shared genes in each comparison group, of which 11 were known DEGs and seemed to play a key role in the color variation process (Fig. 1c). They were st3gal1, si: ch73-86n18.1, plxna4, malb, fabp11a, plcb4, sdhb, si:dkey-65b12.6, si:ch211-157c3.4, agr1, and aqp3, respectively. To verify the credibility of the sequencing results, we randomly selected 25 DEGs related to pigment biosynthesis for qRT-PCR, including 13 up-regulated genes and 12 down-regulated genes. As shown in Fig. 2, the expression patterns of all down-regulated genes were consistent with the sequencing result, and 12 of the 13 up-regulated genes' expression patterns were consistent with the sequencing results. Of which, 24 of 25 genes were significantly different by both RNA-Seq and qRT-PCR methods while dctd genes only significantly differed by RNA-Seq method. The results showed that the reliability of the sequencing result was high. The top GO function enrichment terms for the pairwise comparisons among three samples were shown in Supplementary Table S3. After GO annotation, all DEGs were classified into different biological processes, molecular functions, and cellular components. In each comparison, the top 50 of GO categories were selected in three different categories. The detailed annotations of each category were depicted in Supplementary Fig. S2a-c. In the molecular function category, binding and catalytic activity were the most mapped terms. In the biological process category, cellular process, metabolic process, biological regulation, and regulation of biological process were the most mapped terms. In the cellular component category, cell, cell part, and membrane were the main mapped terms. Furthermore, a few DEGs were mapped to pigmentation-related terms such as developmental pigmentation (GO:0048066), melanocyte differentiation (GO:0030318), melanosome transport (GO: 0032402), retinal pigment epithelium development (GO:0003406), and pigmentation (GO:0043473). These genes enriched in pigmentation-related processes are informative and worthy of further study.

KEGG analysis of the pathways
To further explore the biological functions of the DEGs, an enrichment analysis based on the KEGG database was performed. A total of 141 KEGG pathways were listed in this study (Supplementary Table S4). The DEGs between the P-B-P and P-B skins involved in ribosome, oxidative phosphorylation, ribosome biogenesis in eukaryotes, cardiac muscle contraction, RNA degradation, RNA polymerase, and DNA replication were significantly enriched (P < 0.05). The DEGs were significantly enriched in some genetic information processing between the P-B-P and P skins, including ribosome biogenesis in eukaryotes, spliceosome, RNA degradation, and RNA polymerase. Ten pathways of oxidative phosphorylation, ribosome biogenesis in eukaryotes, tight junction, adrenergic signaling in cardiomyocytes, cardiac muscle contraction, GNRH signaling pathway, mucin-type o-glycan biosynthesis, sphingolipid metabolism, glycosphingolipid biosynthesis, and taurine and hypotaurine metabolism were significantly enriched between the P-B and P skins (P < 0.05). Since fish skin color was mainly correlated with the synthesis of different pigments, we were interested in the pigments biosynthesis pathway. Several pathways including oxidative phosphorylation, ribosome, Wnt (wingless-type MMTV integration site family) signaling pathway, MAPK (mitogen-activated protein kinase) signaling pathway, cell cycle, melanogenesis, tyrosine metabolism, autophagy pathway, and apoptosis pathway were identified, which were related to the skin color regulation and pigmentation (Fig. 3). Candidate genes related to melanophore The KEGG pathway analysis results showed that 32 candidate genes were involved in pigmentationrelated pathways, such as melanogenesis, tyrosine metabolism, Wnt signaling pathway, and MAPK signaling pathway. These genes may play a potential role in the skin color variation of red tilapia during and after overwintering. The heatmap of these genes (Fig. 4) indicated that four genes including agouti signaling protein (asip), tyrosine amino transferase (tat), hydroxyphenylpyruvate hydroxylase (hpdb), and forkhead transcription factor 3 (foxd3) were upregulated in the P and P-B-P groups compared with the P-B group, while the rest genes including tyrosinase (tyr), tyrosinase-related protein 1 (tyrp1), telanocortin receptor 1 (mc1r), microphthalmia-associated transcription factor (mitf), and premelanosome protein (pemel) were down-regulated in the P and P-B-P groups compared with the P-B group. In addition, some DEGs were involved in autophagy and apoptotic pathways such as baxα, beclin1, and atg7. The expression of these genes in P-B-P skin of red tilapia was significantly higher than that of P-B skin, which plays an important role in black-to-pink skin color transformation in red tilapia.

Correlation of candidate genes at protein levels
Based on the candidate genes in P-B vs. P groups, we identified the mutual correlation of their protein products using the STRING online tool (Fig. 5a). The genes (or proteins) from the pigmentation-related pathway were integrated together, and the relationship between them was extensive and strong. Among them, asip, tat, hpdb, and fox3 genes were significantly down-regulated, and other genes were significantly up-regulated. In addition, fourteen hub nodes in a PPI network with more than 10 nodes degrees were shown in Fig. 5b. These hub genes included tyr, mc1r, oca2, mitfb, slc45a2, tyrp1b, dct, asip, kit, kitlg, pmela, pmelb, mitfa, and egfra. Among these genes, the tyr gene showed the highest node degree, which was 18. Fig. 4 Heatmap showing the expressions of selected DEGs in P, P-B, and P-B-P skins. Each row in the map represents a DEG, and the column represents the condition used; Log10 normalized expression value is used for constructing heatmap

Discussion
To explore the different expression patterns among the three types of skins, we performed differential gene expression analysis. When comparing P-B skin to P-B-P skin and P skin, the results showed that more DEGs were up-regulated in P-B skin, indicating that the formation of black skin is complex and that more genes are needed to participate in the process. Combined with the 10 abundance mRNAs result, the genes rich in melanin synthesis were abundantly expressed in P-B skin, while a few of pigment genes were abundantly expressed in P. It was further suggested that melanin genes were involved in the body color variation of red tilapia during overwintering. Eleven known DEGs were shared by P-B-P vs. P-B, P-B-P vs. P, and P-B vs. P comparison groups, of which st3gal1, plxna4, fabp11a, and aqp3 played a vital role in regulating cell proliferation, migration, and invasion (Wu et al. 2018;Wang et al. 2020). It has been reported that silencing of the st3gal1 gene suppresses melanoma invasion and significantly reduces the survival ability of aggressive melanoma cells in humans in a metastatic environment (Pietrobono et al. 2020). Plexins family can functionally activate tyrosine kinase receptors in mammals, such as MET, RON, HER2, and KDR (Swiercz et al. 2008). Overexpression of the aqp3 gene can promote the proliferation and migration of human hepatocytes (Chen et al. 2018a). Aqp3 can reduce the differentiation and inhibit the apoptosis of stem cells in humans through reducing the expressions of related genes in the Wnt/GSK-3 β/β-catenin pathway (Liu et al. 2020).
Fabp11a were probably involved in cellular uptake and transport of fatty acids, targeting of fatty acids to transport systems, and several signaling pathways in Oryzias latipes (Parmar et al. 2012). In this study, all of these genes expression suggested that skin color variation of red tilapia during overwintering might be related to the proliferation, migration, and differentiation of melanocytes. GO enrichment analysis of DEGs revealed that variations in pigmentation were related to cellular components and biological processes. Most of the DEGs clusters were consistent with previous works on fish, such as zebrafish (Higdon et al. 2013), Midas cichlids (Amphilophus) (Henning et al. 2013), and common carp (Cyprinus carpio) (Li et al. 2015). KEGG pathway analysis showed that many DEGs were significantly enriched in oxidative phosphorylation, ribosome, ribosome biogenesis in eukaryotes, and cardiac muscle contraction in the P-B-P vs. P-B group and P-B vs. P group. Several studies have shown that high expression of ribosomal proteinrelated genes was associated with the black color in mice (Skarnes et al. 2011). Four of the five highly expressed genes in pigment cells of zebrafish were ribosomal protein (Higdon et al. 2013). There were many DEGs that participated in oxidative phosphorylation, cardiac muscle contraction signal pathways in Pristella maxillaris (Bian et al. 2019) and Lutjanus erythropterus (Zhang et al. 2015). In addition, the KEGG results of specific mRNAs showed that ribosome biogenesis in eukaryotes, oxidative phosphorylation were dominant pathways in P-B-P skin. Similar results were found in comparative analysis of P-B-P vs. P-B skin and P-B vs. P skin. It was suggested that these pathways may play an important role in color variation in red tilapia.
We also found that the DEGs of three skin colors of tilapia were mainly enriched in the MAPK signaling pathway, Wnt signaling pathway, tyrosine metabolism, and melanogenesis. Tyrosinase metabolism and melanogenesis pathways have been reported in mammals. And both the Wnt and MAPK signaling pathways are involved in melanophore development in vertebrates (Fujimura et al. 2009;Zhang et al. 2017). For specific mRNA, the tyrosinase metabolism and melanogenesis pathways were dominant in P-B skin. It showed that P-B skin required more melanin synthesis than P-B-P and P skins. Meanwhile, we found that some DEGs between the P-B-P and P-B skin expressed abundantly in the process of apoptosis and autophagy pathways. The identification of genes in these pigmentation-related terms and pathways is informative and is worthy of further study.
Studies have suggested that the mRNA expression levels of genes including tyr, tyrp1, mc1r, mitf, and pemel were higher in P-B skin. TYR carries out tyrosine hydroxylation to L-DOPA, which is the first step in the biosynthetic pathway of melanin. Under the action of dopachrome tautomerase (DCT) and TYRP1, the dopaquinone (DOPA) chrome was rapidly oxidized and polymerized to form melanin (Braasch et al. 2010;Simon et al. 2009). Therefore, TYR, TYRP1, and DCT are critical enzymes for the formation of melanin. Mutations or dysfunction of tyr or tyrp1 genes lead to melanocyte death or extensive hypopigmentation in zebrafish (Krauss et al. 2015). In our study, compared to P skin samples, the expression levels of tyr and tyrp1 were the highest in P-B skin. This was also consistent with the pigmentation of red crucian carp (Zhang et al. 2017). Mitf is a member transcription factor involved in the development of melanocytes, retinal cells, osteoclasts, and mast cells (Minvielle et al. 2010). It has been reported that mitf could directly regulate the expression of multiple genes necessary for the development of melanophores, including tyr, tyrp1, and dct (Cheli et al. 2010). Compared with P skin colors, P-B skin color was caused by the increase of melanin content, suggesting that mitf may play a potential role in regulating the differentiation and development of melanocytes. Mc1r gene is a key gene in melanogenesis in animals. Alpha-melanocyte stimulating hormone (α-MSH) binds to mc1r, resulting in the decrease of cAMP level. Consequently, the melanin biosynthesis process was triggered (Voisey et al. 2001). Previous studies have shown that mc1r mutations were associated with skin color variation in many fish species, such as cavefish (Astyanax Mexicanus), guppy (Poecilia reticulata), zebrafish, and koi carp (Cyprinus carpio) (Gross et al. 2009;Tezuka et al. 2011;Richardson et al. 2008;Dong et al. 2020). Similarly, we observed a significant difference in mc1r expression between the red tilapia of the three skin colors used in this study. The Pmel gene acts as a scaffold in the melanosome by creating a proteolytic fibrillary matrix where melanin is deposited (Solano et al. 2000). Pmel mutations promoted pigment dilution in many animals (Gutierrez et al. 2007). Here, pmel was significantly up-regulated in P-B skin when compared to P and P-B-P skin samples. Similarly, we observed the top 10 abundance mRNAs in P-B skin, including tyrp1b, tyr, pmelb, pmela, and tyrp1a. All genes involved in melanin production, transport, and structural proteins for melanin have been verified in red tilapia. In addition, we noticed that the proteins from pigment-related pathways were distinctly integrated together in PPI networks in P-B vs. P group. It was speculated that those proteins (or genes) could be coregulated in skin color variation in red tilapia during overwintering. Among them, the most important top 10 genes based on the key nodes in the PPI network, including tyr, mc1r, mitfb, tyrp1b, dct, pmela, pmelb, and mitfa, were consistent with the results of mRNA expression levels.
Regarding the black to pink stage, the mRNA expression levels of asip, tat, hpdb, and foxd3 were all up-regulated. Asip gene product blocks melanogenesis by antagonizing the binding of α-MSH to mc1r. Asip mutations were associated with skin color variation in Psetta maxima (Scophthalmus maximus), zebrafish, and medaka (Oryzias latipes) (Ceinos et al. 2015, Guillot et al. 2012, José et al. 2005. In our study, we observed higher expression of the asip gene in P-B-P skin transcripts and lower expression of the mc1r gene, which further establishes the role of asip as an antagonistic of the mc1r gene. TAT and HPDB catalyze the substrate tyrosine to form homogeneous acid (HGA). Homogentisate1, 2-dioxygenase catalysis HGA to produce melatonin. A higher level of tat and hpda gene would directly reduce tyrosine level, thereby inhibiting the synthesis of melanin (Zhang et al. 2008). In addition, the tat gene was the most abundant in P skin, suggesting that it affected the skin variation in red tilapia. Foxd3 is a good candidate for the negative regulator of melanophore development. It can affect the lineage between neural or glial and pigment cells by repressing mitf at the early phase of neural crest migration (Thomas and Erickson 2009). In addition, overexpression of foxd3 in melb-a mouse melanoblasts blocked the expression of mitf (Lister et al. 2001). Foxd3 was significantly up-regulated in P-B-P skin samples compared to the P-B skin, indicating that foxd3 might play a significant role in the black-to-pink color transformation in red tilapia.
In addition, autophagy and apoptotic pathways were able to control the transition from black to pink in red tilapia. In detail, the mRNA level of the apoptosis gene, such as baxα, was significantly increased in the body color transformation from the P-B-P skin to the P-B and P skins. Meanwhile, the mRNA levels of beclin1 and autophagy-related genes 7 (atg7), as the autophagy genes, were all up-regulated in the P-B-P skin compared with P-B skin. Baxα, as one of the homologous proteins of BCL-2, could determine survival or death by an apoptotic stimulus. Overexpression of baxa may accelerate cell death (Oltvai et al. 1993). Beclin1 plays a key role in regulating autophagy and cell death by interacting with either BCL-2 or PI3k class III (Takacs-Vellai et al. 2005). Atg7 activates the ubiquitin-like protein ATGL2, which binds to atg5 and extends the autophagic vesicle membrane. Whole-body knock-out of atg7 in mice led to death within 24 h after birth (Komatsu et al. 2005). Baxα and atg7 gene were among the top 10 abundance mRNAs in P-B-P skin, further confirming that the appearance of autophagy may lead to melanocyte reduction.

Conclusions
In conclusion, we performed a transcriptome study of various skin colors during and after overwintering in red tilapia. We screened the top 10 abundance mRNAs, specific mRNAs and identified significant DEGs by pairwise comparison. These specifically expressed mRNAs provide the basis for further studies to clarify the role in skin variation of red tilapia. GO and KEGG analysis of specific mRNAs and DEGs identified numerous signaling pathways. We elucidated 32 candidate genes involved in skin color variation of red tilapia and constructed a PPI network consisting of these genes for revealing the mechanisms of color variation. These findings will help us learn more about the molecular mechanism of skin pigmentation in red tilapia. More specially, it provides valuable genetic data for breeding improved red tilapia strains with consistent skin color.