Transcriptome Analysis of Skin Color Variation During and after Overwintering of Malaysian Red Tilapia

Transcriptome Analysis of Skin Color Variation During and after Overwintering of Malaysian Red Tilapia | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Transcriptome Analysis of Skin Color Variation During and after Overwintering of Malaysian Red Tilapia bingjie Jiang, lanmei Wang, mingkun Luo, jianjun Fu, wenbin Zhu, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-730510/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 13 Apr, 2022 Read the published version in Fish Physiology and Biochemistry → Version 1 posted 4 You are reading this latest preprint version Abstract The commercial value of red tilapia is hampered by variations in skin color during overwintering. In this study, three types skin of red tilapia, including the skin remained pink color during and after overwintering (WP), the skin changed from pink color to black color during overwintering and remained black color after overwintering (PB), and the skin changed from pink color to black color during overwintering but recovered to pink color when the temperature rose after overwintering (BP), were used to analyze their molecular mechanisms of color variation. The transcriptome results revealed that the PB, WP, and BP libraries had 42, 43, 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 BP and PB, BP and WP, and PB and WP, 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 PPI network from PB_WP 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. Physiology General Biochemistry red tilapia skin color transcriptome overwintering differentially expressed transcripts Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Introduction Tilapia is one of the excellent fish species recommended by the Food and Agriculture Organization of the United Nations (Gupta et al. , 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 ). Red tilapia is a valuable fish due to its uniform red skin, the absence of black peritoneum, very 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 pattern 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 metabolism, development and growth of animals (Pavlidis et al., 2008 ). Many animals are dark under the cold and light under the warm condition (Kats et al. , 1986, Sherbrooke et al. , 1989). For example, both the dorsal and ventral skin colors of Rana chiricahuensis in low temperatures (5°C) were significantly darker than those exposed to 25°C (Fernandez et al. , 1991). Pagrus showed a darker dorsal skin area at low (15°C) water temperatures and lighter skin at 19°C (Pavlidis et al., 2008 ). The best pigmentation levels were achieved at temperatures from 26 to 30°C in goldfish ( Carassius auratus ) (Gouveia et al. , 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, an 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 overwintering period of red tilapia might be related to the changes of skin melanocytes and tyrosinase (TYR) activity. In our red tilapia breeding procedure, we found three kind changes of color, i.e. the skin remained pink color during and after overwintering (WP), the skin changed from pink color to black color during overwintering and remained black color after overwintering (PB), and the skin changed from pink color to black color during overwintering but recovered to pink color when the temperature rose after overwintering (BP). In this study, we used RNA-Seq to analyze the transcriptional profiles of WP, PB and BP 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. Materials And Methods 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. Whole pink red tilapia (initial weight: 500 ± 20 g) were cultivated in 2 mu plastic shed pond at the water temperatures of 18 ± 1℃ and fed twice a day (morning and evening) during the winter. In April of the next year, a few whole pink tilapias changed from pink color to black color. Then, the red tilapia with body color variation was cultured in the same environment without plastic shed, and the water temperature gradually rose with the ambient temperature. The red tilapia with reversible body color (black to pink) was selected. Skin tissues were collected from four PB (pink changed to black) red tilapia, four WP (pink unchanged) red tilapia, and four BP (black return to pink) red tilapia individuals, respectively. All fresh tissue samples were immediately snap-frozen in liquid nitrogen and then stored at − 80°C 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 were higher than 1.9, and that of A260:A230 were higher than 1.8. Total RNA integrity was then subsequently checked using an 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 were fragmented under elevated temperature. Then first strand and second strand cDNA were subsequently synthesized. Secondly, the double stranded 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 pair-end (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., 2018 ) was used to discard poor quality reads, including adaptor reads and low-quality reads (reads with more than 50% bases with quality value less than 5). Then the clean reads were mapped onto reference genome independently by HISAT2 version 2.1.0 (Kim et al., 2015 ) with default values. And RSeQC packages (version 2.6.4) (Wang et al., 2012 ) were used to make a comprehensive evaluation on RNA-seq data after alignment, including sequencing saturation, mapped reads distribution, coverage uniformity, strand specificity, transcript level RNA integrity etc . Differential expression analysis Based on the HISAT2 alignment BAM file, featureCounts 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 genes. 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. GO and KEGG enrichment analysis of differentially expressed genes Gene Ontology (GO) term ( http://www.geneontology.org/ ) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp/ ) enrichment analysis of DEGs were performed using clusterProfiler package (Yu et al., 2012 ), KOBAS v3.0 (Xie et al., 2011 ), respectively. The calculated p -value goes through Bonferroni Correction, taking P .adjust value ≤ 0.05 as a threshold. Quantitative real-time PCR analysis Total RNAs were extracted as described above. Each RNA sample was treated by 5× PrimeScript™ RT Master Mix (Takara) to remove residual genomic DNA and reverse transcribed into cDNA. All primer pairs (Additional file 1: 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 Real-time 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 nuclease-free water. Amplification was performed with an initial denaturation at 95°C for 5 min, followed by 40 cycles of 95°C for 10 s, 60°C for 30 s and 72°C 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 β-actin gene as the reference control. Data were analyzed statistically with SPSS 20 (IBM, Chicago, IL, USA) by t-test. Thresholds for statistical significance were set at P <0.05 (significant) and P <0.01 (highly significant). Co-expression analysis of protein-protein interaction (PPI) 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 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 PPI network and count the number of nodes in DEGs (Franz et al., 2016 ). The node genes with node degree above 10 were selected as the key target genes. Results Overview of the RNA-Seq data To better understand skin color variation of red tilapia during and after overwintering, the mRNA libraries of PB, WP and BP were determined and analyzed by Illumina sequencing technology. In Table 1 , we presented the Q30 percentage, GC percentage and the other indexes to describe the libraries. In total, an average of 42,721,496, 43,921,378 and 43,386,254 raw reads were obtained from the BP, PB and WP libraries, respectively. After filtering the low-quality reads and removing adaptor sequences, the average of 42,664,181, 43,888,258 and 43,347,484 clean reads were retrieved for further analysis. The percentage of G + C content and Q30 ratio was an average of 48.23% and 94.27%, indicated a high-quality sequence. All raw transcriptome data were submitted to the NCBI Short Read Archive (SRA) database with the accession numbers PRJNA690595. Table 1 The specific statistics for each library sequencing and quality control Sample Raw reads Raw base(G) Clean reads Clean base(G) Q30(%) GC(%) BP-1 40,313,772 6.05 40,254,280 6.00 95.23 48.31 BP-2 44,376,218 6.66 44,317,948 6.61 95.01 47.34 BP-3 40,589,108 6.09 40,530,756 6.05 95.01 47.84 BP-4 45,606,884 6.84 45,553,738 6.80 95.12 48.52 PB-1 46,324,834 6.95 46,289,448 6.88 93.32 47.41 PB-2 43,749,772 6.56 43,722,494 6.49 93.97 48.04 PB-3 40,425,896 6.06 40,402,734 5.98 93.81 48.37 PB-4 45,185,008 6.78 45,138,356 6.71 94.34 47.72 WP-1 43,218,496 6.48 43,170,434 6.42 93.80 49.02 WP-2 43,187,756 6.48 43,150,158 6.42 93.63 49.12 WP-3 43,362,090 6.50 43,324,748 6.46 93.92 48.15 WP-4 43,776,672 6.57 43,744,596 6.51 94.04 48.86 Average 43,343,042 6.50 43,299,974 6.44 94.27 48.23 To assess the quality of sequencing and reassembly, all clean reads were mapped to 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 (Additional file 1: 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% (Additional file 1: 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 were shown in Additional file 1: 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 were the most, while the number of genes with RPKM ≥ 100 were less than 1%. Expression profiling of mRNAs To further understand mRNA expression of differences 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 the further analysis (SPM > 0.9, Additional file 2: Table S1). In detail, there were 119, 294, and 53 specific mRNAs in BP, PB, and WP skins respectively. The KEGG results of specific mRNAs showed that metabolic pathways, ribosome biogenesis in eukaryotes and oxidative phosphorylation were dominant pathway in BP skin and regulation of actin cytoskeleton, melanogenesis, tight junction and tyrosine metabolism were dominant in PB skin (Additional file 2: Table S2). Furthermore, we also analyzed the top 10 abundance mRNAs of differences skin color in red tilapia. As shown in Table 2 , Granulin and tat gene were abundant from WP skin and other pigment-related genes were abundantly expressed in them, including oca2 and slc45a4 . Tyrp1b gene was the most abundant in PB skin and melanin synthesis gene accounts for the largest proportion of abundance expressed genes, such as tyrp1b , pmelb , tyr , pmela , and tyrp1a . In addition, baxa gene showed dominantly expression in BP skin. Table 2 Ten most abundant genes of three difference colors in red tilapia during and after overwintering WP PB BP granulin tyrp1b baxα tat cavin2a si:ch73-86n18.1 ENSONIG00000007142 zgc:101810 crtac1a oca2 pmelb ENSONIG00000040753 ENSONIG00000040391 tyr sdhb oni-mir-24a-4 myh9a ppdpfa col10a1a ENSONIG00000019137 loxa ENSONIG00000039537 aqp3 ENSONIG00000009288 slc45a4 pmela ENSONIG00000018504 map3k7cl tyrp1a atg7 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 (BP_PB, BP_WP, and PB_WP) revealed the expression trend of each pair (Fig. 1 a). We also constructed a histogram of DEGs in the three skin tissues (Fig. 1 b). Compared with the PB skin, there were 2,528 DEGs in BP skin, of which 1,420 were up-regulated and 1,108 were down-regulated. A total of 1,091 DEGs were up-regulated in BP skin compared with WP skin, while 833 DEGs were down-regulated. There were 1,939 DEGs displaying greater abundance in PB skin compared with WP skin, of which 1,387 DEGs were up-regulated and 552 genes were down-regulated. Among these genes, 22 DEGs were detected as shared genes in each comparison group, in which 11 were known DEGs and seemed to play a key role in the color variation process (Fig. 2 c). 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. The results showed that the reliability of the sequencing result was high. Functional enrichment analysis of DEGs The top GO function enrichment terms for the pairwise comparisons among three samples were shown in Additional file 2: Table S3. After GO annotation, all DEGs were classified into different biological processes, molecular function and cellular component. In each comparison, the top 50 of GO categories were selected in three different categories. The detailed annotations of each category were depicted in Additional file 1: Fig. S1a-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 terms 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 a further study. KEGG analysis of the pathways To further explore the biological functions of the DEGs, an enrichment analysis based on KEGG database was performed. A total of 141 KEGG pathways were listed in this study (Additional file 2: Table S4). The DEGs between the BP and PB skins were 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 BP and WP 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, taurine and hypotaurine metabolism were significantly enriched between the PB and WP skins ( P < 0.05). Since fish skin color was mainly correlated with the synthesis of different pigments, we were interested in 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, etc. 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 pigmentation-related pathways, such as melanogenesis, tyrosine metabolism, Wnt signaling pathway and MAPK signaling pathway. These genes may play a potential role in 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 up-regulated in WP and BP group compared with the PB group, while the rest genes including tyrosinase ( tyr ), tyrosinase-related protein 1 ( tyrp1 ), telanocortin receptor 1 ( mc1r ), microphthalmia-associated transcription factor ( mitf ), premelanosome protein ( pemel ) etc . were down-regulated in WP and BP group compared with PB group. In addition, some DEGs were involved in autophagy and apoptotic pathways such as baxα , beclin1 and atg7 . The expression of these genes in BP skin of red tilapia was significantly higher than that of PB skin, which play 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 PB_WP groups, we identified the mutual correlation of their protein products using the STRING online tool (Fig. 5 a). The genes (or proteins) from the pigmentation-related pathway were integrated together, and the relationship of them were extensive and strong. Among them, asip , tat , hpdb , 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 degree were shown in Fig. 5 b. These hub genes included tyr, mc1r , oca2 , mitfb , slc45a2 , tyrp1b , dct , asip , kit , kitlg , pmela , pmelb , mitfa , and egfra . Among these genes, tyr gene showed the highest node degree, which was 18. Discussion To explore the different expression patterns among the three types skins, we performed differential gene expression analysis. When comparing PB skin to BP skin and WP skin, the results showed that more DEGs were up-regulated in PB 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, there genes rich in melanin synthesis were abundantly expressed in PB skin, while a few of pigment genes were abundantly expressed in WP. 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 BP_PB, BP_WP and PB_WP comparison groups, of which st3gal1 , plxna4 , fabp11a and aqp3 plays 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 st3gal1 gene suppresses melanoma invasion and significantly reduces the survive ability of aggressive melanoma cells of human in metastatic environment (Pietrobono et al., 2020 ). Plexins family can functionally activate tyrosine kinase receptors in mammalian, such as MET, RON, HER2, and KDR (Swiercz et al., 2008 ). Overexpression of aqp3 gene can promoted the proliferation and migration of human hepatocytes (Chen et al., 2018 ). Aqp3 can reduce the differentiation and inhibit the apoptosis of stem cells in human through reducing the expressions of related genes in 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 signalling 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 (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 BP_PB group and PB_WP group. Several studies have shown that high expression of ribosomal protein related 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 participate in oxidative phosphorylation, cardiac muscle contraction signal pathways in Pristella maxillaris (Bian et al ., 2013) 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 BP skin. Similar results were found in comparative analysis of BP_PB skin and PB_WP skin. It is 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 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 PB skin. It showed that PB skin required more melanin synthesis than BP and WP skins. Meanwhile, we found that some DEGs between the BP and PB skin expressed abundantly in the process of apoptosis and autophagy pathways. The identification of genes in these pigmentation-related term and pathways are informative and are worthy of further study. Studies have suggested that the mRNA expression levels of genes including tyr , tyrp1 , mc1r , mitf , pemel were higher in PB 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 WP skin samples, the expression levels of tyr and tyrp1 were the highest in PB 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 WP colors, PB 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, 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, guppy, zebrafish and koi carp (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 three skin colors used in this study. 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 PB skin when compared to WP and BP skin samples. Similarly, we observed the top 10 abundance mRNAs in PB skin, including tyrp1b , tyr , pmelb , pmela , tyrp1a, etc. 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 PB_WP group. It was speculated that those proteins (or genes) could be co-regulated 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 , etc , 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, zebrafish and medaka (Ceinos et al., 2015 , Guillot et al., 2012 , Jose´ et al., 2005 ). In our study, we observed higher expression of asip gene in BP skin transcripts and lower expression of mc1r , 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. Higher level of tat and hpda gene would directly reduce tyrosine level, thereby inhibited the synthesis of melanin (Zhang et al., 2008 ). In addition, tat gene was the most abundant in WP skin, suggested 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 et al. , 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 BP skin samples compared to the PB skin, indicated 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 apoptosis gene, such as baxα , was significantly increased in the body color transformation from the BP skin to the PB and WP skins. Meanwhile, the mRNA levels of beclin1 and autophagy-related genes 7 ( atg7 ), as the autophagy genes, were all upregulated in the BP skin compared with PB 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 interacts 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 24h after birth (Komatsu et al., 2005 ). Baxα and atg7 gene were among the top 10 abundance mRNAs in BP skin, further confirmed that 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 consist with 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. Declarations Funding projects This study was supported by the National Natural Science Foundation-Youth Fund Project (31802290). Conflicts of interest/Competing interests The authors have no relevant financial or non-financial interests to disclose. Code availability Not applicable Availability of data and material The authors declare that all other data supporting the findings of this study are available within the article and its supplementary information files. Authors' contributions Zaijie Dong conceived the study; Wenbin Zhu provided the experimental materials; Lanmei provided the funding for the experiment; Bingjie Jiang performed the experiments and wrote the paper; Jianjun Fu provided technical assistance in experiments; Mingkun Luo revised the manuscript; Zaijie Dong reviewed the manuscript; All authors read and approved the manuscript. Ethics approval The sampling scheme and experimental protocols were subject to approval by the Bioethical Committee of Freshwater Fisheries Research Center (FFRC) of the Chinese Academy of Fishery Sciences (CAFS) (BC 2013863, 9/2013). The methods of samples handled and experimental procedures carried out in accordance with the guidelines for the care and use of animals for scientific purposes issued by the Ministry of Science and Technology, Beijing China (No.398, 2006). Availability of data and material/ Data availability All data generated or analysed during this study are included in this published article (and its supplementary information files). Consent to participate Not applicable Consent for publication Not applicable References Andrews S (2014) FastQC: a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ Bian FF, Yang XF, Ou ZJ et al (2019) Morphological characteristics and comparative transcriptome analysis of three different phenotypes of Pristella maxillaris . Front Genet 10:698 Braasch I, Liedtke D, Volff JN, Schartl M (2010) Pigmentary function and evolution of tyrp1 gene duplicates in fish. Pigment Cell Melanoma Research 22(6):839–850 Ceinos RM, Guillot R, Kelsh RN et al (2015) Pigment patterns in adult fish result from superimposition of two largely independent pigmentation mechanisms. Pigment Cell Melanoma Research 28(2):196–209 Cheli Y, Ohanna M, Ballotti R, Bertolotto C (2010) Fifteen-year quest for microphthalmia-associated transcription factor target genes. Pigment Cell Melanoma Research 23(1):27–40 Chen G, Shi Y, Liu M, Sun J et al (2018) circHIPK3 regulates cell proliferation and migration by sponging miR-124 and regulating AQP3 expression in hepatocellular carcinoma. Cell Death Dis 9(2):175 Chen YX, Chen YS, Shi CM et al (2018) Soapnuke: a mapreduce acceleration supported software for integrated quality control and preprocessing of high-throughput sequencing data. Oxford Open 7(1):gix120 Dong ZJ, Luo MK, Wang LM et al (2020) MicroRNA-206 regulation of skin pigmentation in koi carp ( Cyprinus carpio L.). Front Genet 11:47 Francisco S, María ME, Celia JC, Simon PH et al (2000) New insights on the structure of the mouse silver locus and on the function of the silver protein. Pigment Cell Melanoma Research 13(8):118–124 Franz M, Lopes CT, Huck G, Dong Y, Sumer O, Bader GD (2016) Cytoscape.Js a graph theory library for visualisation and analysis. Bioinformatics 32(2):309–311 Fujimura N, Taketo MM, Mori M et al (2009) Spatial and temporal regulation of Wnt/β-catenin signaling is essential for development of the retinal pigment epithelium. Dev Biol 334(1):31–45 Gouveia L, Rema P (2005) Effect of microalgal biomass concentration and temperature on ornamental goldfish ( Carassius auratus ) skin pigmentation. Aquac Nutr 11(1):19–23 Gross JB, Borowsky R, Tabin CJ (2009) A novel role for mc1r in the parallel evolution of depigmentation in independent populations of the cavefish Astyanax Mexicanus . Plos Genetics 5(1):1000326 Guillot R, Ceinos RM, Cal R et al (2012) Transient ectopic overexpression of agouti-signalling protein 1 ( asip1 ) induces pigment anomalies in flatfish. Plos One 7(12):e48526 Gupta MV, Acosta BO (2004) A review of global tilapia farming practices. Aquaculture Asia 9(1):7–12 Gutierrez GB, Wiener P, Williams JL (2007) Genetic effects on coat color in cattle: dilution of eumelanin and´ phaeomelanin pigments in an F2-backcross Charolais × Holstein population. BMC Genet 8(1):56 Henning F, Jones JC, Franchini P et al (2013) Transcriptomics of morphological color change in polychromatic Midas cichlids. BMC Genom 14(1):171–171 Higdon CW, Mitra RD, Johnson SL et al (2013) Gene expression analysis of zebrafish melanocytes, iridophores, and retinal pigmented epithelium reveals indicators of biological function and developmental origin. Plos One 8(7):e67801 Huang R, Kitajama C, Ke C (1997) A Study on Induced factor of White Skin of Flat Fish Paralichthys Olivaceus (Temminck et Schlegel). Modern Fisheries Information 4(9):21–23 Hubbard JK, Uy JA, Hauber ME, Hoekstra HE, Safran RJ (2010) Vertebrate pigmentation: from underlying genes to adaptive function. Trends Genet 26:231–239 Jennifer R, Pia RL, Natalie LR, Julia RD et al (2008) Mc1r Pathway regulation of zebrafish melanosome dispersion. Zebrafish 5(4):289–295 Jiang YL, Zhang SH, Xu J et al (2014) Comparative transcriptome analysis reveals the genetic basis of skin color variation in common carp. Plos One 9(9):e108200 Jose´ MCR, Tatjana H, Helgi BS, Richard EP (2005) Gene structure of the goldfish agouti-signaling protein: a putative role in the dorsal-ventral pigment pattern of fish. Endocrinology 146(3):1597–1610 Kim D, Langmead B, Salzberg SL (2015) Hisat: a fast spliced aligner with low memory requirements. Nature methods 12(4):357–360 Komatsu M, Waguri S, Ueno T et al (2005) Impairment of starvation-induced and constitutive autophagy in atg7 deficient mice. J Cell Biol 169(3):425–434 Krauss J, Geiger-Rudolph S, Koch I, Nusslein-Volhard C, Irion U (2015) A dominant mutation in tyrp1a leads to melanophore death in zebrafish. Pigm Cell Melanoma Research 27(5):827–830 Li XJ, Li SF, Feng JH et al (2003) Preliminary study on salinity tolerance of Israel red tilapia. Journal of Shanghai Fisheries University 12(3):205–208 Li XM, Song YN, Xiao GB, Zhu BH, Xu GC, Sun MY et al (2015) Gene expression variations of red-white skin coloration in common carp ( Cyprinus carpio ). Int J Mol Sci 16:21310–21329 Liao Y, Gordon K, Smyth, Shi W (2014) Featurecounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7):923–930 Lister JA, Close J, Raible DW (2001) Duplicate mitf genes in zebrafish: complementary expression and conservation of melanogenic potential. Dev Biol 237(2):333–344 Liu CS, Liu LY, Zhang Y, Jing HY (2020) Molecular mechanism of aqp3 in regulating differentiation and apoptosis of lung cancer stem cells through Wnt/GSK-3β/β-Catenin pathway. Jbuon 25(4):1714–1720 Luo MK, Lu GQ, Yin HR et al (2021) Fish pigmentation and coloration: molecular mechanisms and aquaculture perspectives. Reviews in Aquaculture 00:1–18 Minvielle F, Bed'Hom B, Coville JL et al (2010) The "silver" Japanese quail and the mitf gene: causal mutation, associated traits and homology with the "blue" chicken plumage. BMC Genet 11(1):1–7 Nakano T, Ono K, Takeuchi M (1992) Levels of zinc, iron, and cupper in the skin of abnormal pigmented Japanese flounder. Nippon Suisan Gakkaishi 58(11):2207 Oltvai ZN, Milliman CL, Korsmeyer SJ (1993) Bcl-2 heterodimerizes in vivo with a conserved homolog, bax , that accelerates programed cell death. Cell 74(4):609–619 Pan JB, Hu SC, Wang H et al (2012) PaGeFinder: quantitative identification of spatiotemporal pattern genes. Bioinformatics 28(11):1544–1545 Parmar MB, Venkatachalam AB, Wright JM (2012) The evolutionary relationship of the transcriptionally active fabp11a (intronless) and fabp11b genes of medaka with fabp11 genes of other teleost fishes. Febs Journal 279(13):2310–2321 Pavlidis M, Karkana M, Fanouraki E et al (2008) Environmental control of skin color in the red porgy, Pagrus pagrus . Aquac Res 39(8):837–849 Pietrobono S, Anichini G, Sala C et al (2020) ST3GAL1 is a target of the SOX2-GLI1 transcriptional complex and promotes melanoma metastasis through AXL. Nat Commun 11(1):5865 Pradeep PJ, Srijaya TC, Hassan A et al (2014) Optimal conditions for cold-shock induction of triploidy in red tilapia. Aquacult Int 22(3):1163–1174 Robinson MD, Mccarthy DJ, Smyth GK (2010) Edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139–140 Simon DJ, Peles D, Wakamatsu K, Ito S (2009) Current challenges in understanding melanogenesis: bridging chemistry, biological control, morphology and function. Pigment Cell Melanoma Research 22(5):563–579 Skarnes WC, Rosen B, West AP et al (2011) A conditional knockout resource for the genome-wide study of mouse gene function. Nature 474(7351):337–342 Swiercz JM, Worzfeld T, Offermanns S (2008) ErbB-2 and met reciprocally regulate cellular signaling via plexin-B1. J Biol Chem 283(4):1893–1901 Szklarczyk D, Franceschini A, Wyder S et al (2015) STRING v10: protein-protein networks, integrated over the tree of life. Nucleic Acids Res 43:D447–D452 Takacs-Vellai K, Vellai T, Puoti A, Passannante M et al (2005) Inactivation of the autophagy gene bec-1 triggers apoptotic cell death in C. elegans . Curr Biol 15(16):1513–1517 Tezuka A, Yamamoto H, Yokoyama J et al (2011) The mc1r gene in the guppy ( Poecilia reticulata ): Genotypic and phenotypic polymorphisms. BMC Research Notes 4(1):31 Thomas AJ, Erickson CA (2009) Foxd3 regulates the lineage switch between neural crest-derived glial cells and pigment cells by repressing mitf through a non-canonical mechanism. Development 136(11):1849–1858 Voisey J, Box NF, Daal AV (2001) A polymorphism study of the human agouti gene and its association with mc1r . Pigment Cell Res 14(4):264–267 Wang LG, Wang SQ, Li W (2012) RSeQC: quality control of RNA-seq experiments. Bioinformatics 28(16):2184–2185 Wang LM, Song FB, Zhu WB et al (2018) Effects of temperature on body color of Malaysian red tilapia during overwintering period. Journal of Fisheries of China 42(1):72–79 Wang XY, Yang J, Yao Y, Shi X, Yang GS, Li X (2020) Aqp3 facilitates proliferation and adipogenic differentiation of porcine intramuscular adipocytes. Genes 11(4):453 Wardani WW, Alimuddin A, Junior MZ et al (2020) Growth performance, robustness against stress, serum insulin, igf-1 and glut4 gene expression of red tilapia ( Oreochromis sp. ) fed diet containing graded levels of creatine. Aquaculture Nutrition 27(1):274–286 Wu X, Zhao JD, Ruan YY et al (2018) Sialyltransferase st3gal1 promotes cell migration, invasion, and TGF-β1-induced EMT and confers paclitaxel resistance in ovarian cancer. Cell Death Dis 9(11):1102 Xie C, Mao XZ, Huang JJ, Yang D et al (2011) KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res 39(suppl_2):W316–W322 Yu GC, Wang LG, Han YY, He QY (2012) Clusterprofiler: an r package for comparing biological themes among gene clusters. Omics-a Journal of Integrative Biology 16(5):284–287 Zhang Y, Wang L, Zhang S, Yang H, Tan H (2008) HmgA, transcriptionally activated by hpda , influences the biosynthesis of actinorhodin in Streptomyces Coelicolor . FEMS Microbiol Lett 280(2):219–225 Zhang YP, Wang ZD, Guo YS et al (2015) Morphological characters and transcriptome profiles associated with black skin and red skin in crimson snapper ( Lutjanus erythropterus ). Int J Mol Sci 16(11):26991–27004 Zhang YQ, Jin H et al (2017) Comparative transcriptome analysis of molecular mechanism underlying gray-to-red body color formation in red crucian carp ( Carassius auratus , red var.). Fish Physiol Biochem 43(6):1387–1398 Zhu WB, Wang LM, Dong ZJ et al (2016) Comparative transcriptome analysis identifies candidate genes related to skin color differentiation in red tilapia. Sci Rep 6(1):31347 Cite Share Download PDF Status: Published Journal Publication published 13 Apr, 2022 Read the published version in Fish Physiology and Biochemistry → Version 1 posted Reviews received at journal 01 Aug, 2021 Reviewers invited by journal 21 Jul, 2021 Editor assigned by journal 20 Jul, 2021 First submitted to journal 18 Jul, 2021 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-730510","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":40949061,"identity":"455c946a-f292-4997-890f-d0471735fbfa","order_by":0,"name":"bingjie Jiang","email":"","orcid":"","institution":"Nanjing Agricultural University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"bingjie","middleName":"","lastName":"Jiang","suffix":""},{"id":40949062,"identity":"47a6a784-3c91-4abc-a55f-32a455546401","order_by":1,"name":"lanmei Wang","email":"","orcid":"","institution":"Chinese Academy of Fishery Sciences Freshwater Fisheries Research Center","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"lanmei","middleName":"","lastName":"Wang","suffix":""},{"id":40949063,"identity":"edb1c6f3-7dda-48e8-a991-ae0a4f7b0755","order_by":2,"name":"mingkun Luo","email":"","orcid":"","institution":"Chinese Academy of Fishery Sciences Freshwater Fisheries Research Center","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"mingkun","middleName":"","lastName":"Luo","suffix":""},{"id":40949064,"identity":"7d2b3dcc-fda2-48a7-a3d1-6eeb9483dcdd","order_by":3,"name":"jianjun Fu","email":"","orcid":"","institution":"Chinese Academy of Fishery Sciences Freshwater Fisheries Research Center","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"jianjun","middleName":"","lastName":"Fu","suffix":""},{"id":40949065,"identity":"2c4ad243-10d8-4307-be99-1db94cc875e1","order_by":4,"name":"wenbin Zhu","email":"","orcid":"","institution":"Chinese Academy of Fishery Sciences Freshwater Fisheries Research Center","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"wenbin","middleName":"","lastName":"Zhu","suffix":""},{"id":40949066,"identity":"5400d0fe-a654-4b1f-b84a-020d3fd761dc","order_by":5,"name":"zaijie dong","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAq0lEQVRIiWNgGAWjYDCCAyCCx4afvYFELWmSPQdI08JwmAQtfMd7D7/8IXNegkf68AOGnzuI0CJ55lyahQTPbQkevjQDxt4zRGgxuJFjZmDAc7vOnoeHgZmxjVgtCTznJHhI0WL84ADPARK0SJ45Y8bYwJMM1MJmcLCXGC18x3uMP/7ssQNqYX744CcxWoCATYKxB8I6QJwGBgbmDww/iFU7CkbBKBgFIxIAANJKMjkbJh4jAAAAAElFTkSuQmCC","orcid":"","institution":"Nanjing Agricultural University","correspondingAuthor":true,"submittingAuthor":false,"prefix":"","firstName":"zaijie","middleName":"","lastName":"dong","suffix":""}],"badges":[],"createdAt":"2021-07-18 16:22:20","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-730510/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-730510/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1007/s10695-022-01073-5","type":"published","date":"2022-04-13T20:49:07+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":11861796,"identity":"cc4e2110-c649-4a7b-b9d0-05b830b1e1be","added_by":"auto","created_at":"2021-07-27 22:22:10","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":406182,"visible":true,"origin":"","legend":"DEGs in BP, WP and PB skins. a, Volcano plot of differential mRNA expression levels among the three pairwise comparisons. The gray, red, and blue dots represent non-significant, up-regulated and down-regulated transcripts, respectively; b, Number of DEGs among the three pairwise comparisons. The red and blue color stand for up-regulated and down-regulated expression, respectively; c, DEGs number and Venn diagram of the overlap of the different groups","description":"","filename":"1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-730510/v1/0e5c750f67395c28feb9d26d.jpg"},{"id":11861795,"identity":"fe86090d-57c5-4b75-b5aa-0516edbc1ecb","added_by":"auto","created_at":"2021-07-27 22:22:09","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":199314,"visible":true,"origin":"","legend":"Comparison of mRNA expression levels among the 25 DEGs obtained using qRT-PCR validation and RNA sequencing. Log-fold changes are expressed as the ratio of gene expression after normalization to β-actin","description":"","filename":"2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-730510/v1/65cd565c05b43c20a243af60.jpg"},{"id":11862043,"identity":"86be069e-0b79-4a5b-bf94-99438c80a70e","added_by":"auto","created_at":"2021-07-27 22:25:10","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":403557,"visible":true,"origin":"","legend":"Pigmentation-related pathways based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Gene number: number of genes in each pathway; Rich factor: ratio of the number of target genes divided by the total number of genes in each pathway","description":"","filename":"3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-730510/v1/721d87a1f4e02f9bedcb6e2e.jpg"},{"id":11861797,"identity":"67621440-702e-48fa-84d8-7e4529d13c36","added_by":"auto","created_at":"2021-07-27 22:22:10","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":508282,"visible":true,"origin":"","legend":"Heatmap of showing the expressions of selected DEGs in WP, PB and BP skins. Note: Each row in the map represents a DEG and column represents condition used; Log10 normalized expression value is used for constructing heat-map","description":"","filename":"4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-730510/v1/ea93115f01859b22b1501497.jpg"},{"id":11861799,"identity":"6aa35446-76aa-4fba-855a-88163332db6f","added_by":"auto","created_at":"2021-07-27 22:22:10","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":573113,"visible":true,"origin":"","legend":"The mutual protein-protein interactions of candidate genes in PB_WP comparison. a, PPI network; b, The genes of the PPI according to the node degree over 10","description":"","filename":"5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-730510/v1/4fdf868f1f51181a668ad48f.jpg"},{"id":20316632,"identity":"1b01653c-1395-4a12-bc24-de3d5ad1dd3a","added_by":"auto","created_at":"2022-04-13 20:49:10","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":888699,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-730510/v1/966a4eb0-a34d-4ff1-b7d9-5eca44f9c234.pdf"}],"financialInterests":"","formattedTitle":"\u003cp\u003eTranscriptome Analysis of Skin Color Variation During and after Overwintering of Malaysian Red Tilapia\u003c/p\u003e","fulltext":[{"header":"Introduction","content":"\u003cp\u003eTilapia is one of the excellent fish species recommended by the Food and Agriculture Organization of the United Nations (Gupta \u003cem\u003eet al.\u003c/em\u003e, 2004). In recent years, tilapia has been widely accepted and has become an export dominant species of aquaculture in China (Pradeep et al., \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). Red tilapia was obtained by crossing the mutant \u003cem\u003eOreochromis mossambicus\u003c/em\u003e with other tilapia populations such as \u003cem\u003eOreochromis niloticus\u003c/em\u003e and \u003cem\u003eOreochromis aureus\u003c/em\u003e (Li et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2003\u003c/span\u003e). Red tilapia is a valuable fish due to its uniform red skin, the absence of black peritoneum, very 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., \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). Therefore, most studies on red tilapia have mainly focused on their growth and development in genetic breeding (Wardani et al., \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Zhu et al., \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). However, the key issue restricting the growth of commercial red tilapia cultures is skin color variation during overwintering. Pavlidis et al. (\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2008\u003c/span\u003e) found that water temperature changed the body color by motility of chromatophore in red porgy (\u003cem\u003ePagrus pagrus\u003c/em\u003e), but the molecular mechanism for this change is unclear.\u003c/p\u003e \u003cp\u003eAs 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., \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). Furthermore, in aquaculture species, skin pigmentation pattern 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., \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2014\u003c/span\u003e; Luo et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Water temperature is a major environmental factor for metabolism, development and growth of animals (Pavlidis et al., \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). Many animals are dark under the cold and light under the warm condition (Kats \u003cem\u003eet al.\u003c/em\u003e, 1986, Sherbrooke \u003cem\u003eet al.\u003c/em\u003e, 1989). For example, both the dorsal and ventral skin colors of \u003cem\u003eRana chiricahuensis\u003c/em\u003e in low temperatures (5\u0026deg;C) were significantly darker than those exposed to 25\u0026deg;C (Fernandez \u003cem\u003eet al.\u003c/em\u003e, 1991). \u003cem\u003ePagrus\u003c/em\u003e showed a darker dorsal skin area at low (15\u0026deg;C) water temperatures and lighter skin at 19\u0026deg;C (Pavlidis et al., \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). The best pigmentation levels were achieved at temperatures from 26 to 30\u0026deg;C in goldfish (\u003cem\u003eCarassius auratus\u003c/em\u003e) (Gouveia \u003cem\u003eet al.\u003c/em\u003e, 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.\u003c/p\u003e \u003cp\u003eIn our previous study, an Illumina RNA-seq and microRNA-seq analysis were conducted on different color varieties of red tilapia (Zhu et al., \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2016\u003c/span\u003e; Wang et al., \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). Wang et al. (\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) indicated that the color variation during overwintering period of red tilapia might be related to the changes of skin melanocytes and tyrosinase (TYR) activity. In our red tilapia breeding procedure, we found three kind changes of color, i.e. the skin remained pink color during and after overwintering (WP), the skin changed from pink color to black color during overwintering and remained black color after overwintering (PB), and the skin changed from pink color to black color during overwintering but recovered to pink color when the temperature rose after overwintering (BP). In this study, we used RNA-Seq to analyze the transcriptional profiles of WP, PB and BP 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.\u003c/p\u003e"},{"header":"Materials And Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eSample collection\u003c/h2\u003e \u003cp\u003eThe 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. Whole pink red tilapia (initial weight: 500\u0026thinsp;\u0026plusmn;\u0026thinsp;20 g) were cultivated in 2 mu plastic shed pond at the water temperatures of 18\u0026thinsp;\u0026plusmn;\u0026thinsp;1℃ and fed twice a day (morning and evening) during the winter. In April of the next year, a few whole pink tilapias changed from pink color to black color. Then, the red tilapia with body color variation was cultured in the same environment without plastic shed, and the water temperature gradually rose with the ambient temperature. The red tilapia with reversible body color (black to pink) was selected.\u003c/p\u003e \u003cp\u003eSkin tissues were collected from four PB (pink changed to black) red tilapia, four WP (pink unchanged) red tilapia, and four BP (black return to pink) red tilapia individuals, respectively. All fresh tissue samples were immediately snap-frozen in liquid nitrogen and then stored at \u0026minus;\u0026thinsp;80\u0026deg;C until use.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003eRNA extraction, cDNA library construction, and sequencing\u003c/h2\u003e \u003cp\u003eTotal RNA was obtained from red tilapia samples using RNA TRIzol (Invitrogen, UK) according to the manufacturer\u0026rsquo;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 were higher than 1.9, and that of A260:A230 were higher than 1.8. Total RNA integrity was then subsequently checked using an Bioanalyzer RNA 6000 Pico Kit (Agilent Technologies). Samples with an RNA Integrity Number (RIN)\u0026thinsp;\u0026gt;\u0026thinsp;8 were retained for subsequent analysis.\u003c/p\u003e \u003cp\u003eA 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\u0026rsquo;s instructions. Firstly, mRNA was purified from total RNA by Poly-T oligoattached magnetic beads, and then were fragmented under elevated temperature. Then first strand and second strand cDNA were subsequently synthesized. Secondly, the double stranded 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\u0026times;150 bp pair-end (PE) sequencing.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eQuality control and Mapping to the reference genome\u003c/h2\u003e \u003cp\u003eQuality of all raw reads was conducted by FastQC (Andrews \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) software. At the initial filtering step, SOAPnuke (Chen et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) was used to discard poor quality reads, including adaptor reads and low-quality reads (reads with more than 50% bases with quality value less than 5). Then the clean reads were mapped onto reference genome independently by HISAT2 version 2.1.0 (Kim et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2015\u003c/span\u003e) with default values. And RSeQC packages (version 2.6.4) (Wang et al., \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2012\u003c/span\u003e) were used to make a comprehensive evaluation on RNA-seq data after alignment, including sequencing saturation, mapped reads distribution, coverage uniformity, strand specificity, transcript level RNA integrity \u003cem\u003eetc\u003c/em\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eDifferential expression analysis\u003c/h2\u003e \u003cp\u003eBased on the HISAT2 alignment BAM file, featureCounts v1.6.2 (Liao et al., \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) was used to estimate and quantify gene expression with default parameters, so as to generate the raw read count of each RNA genes. Gene expression was normalized by reads per kilobase of exon per million reads mapped (RPKM). Finally, edgeR (Robinson et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2010\u003c/span\u003e) was used to identify the DEGs by pairwise comparisons. The difference was considered significant if the |logFC| \u0026ge;1 and FDR (False Discovery Rate)\u0026thinsp;\u0026le;\u0026thinsp;0.05.\u003c/p\u003e \u003cp\u003eTo further understand the mRNA expression of each sample in red tilapia, the specificity measure (SPM) was introduced to analyze all screened mRNA (FPKM value\u0026thinsp;\u0026ge;\u0026thinsp;1, at least 3 samples) by PaGeFinder algorithms (Pan et al., \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). 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.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003eGO and KEGG enrichment analysis of differentially expressed genes\u003c/h2\u003e \u003cp\u003eGene Ontology (GO) term (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.geneontology.org/\u003c/span\u003e\u003c/span\u003e) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.kegg.jp/\u003c/span\u003e\u003c/span\u003e) enrichment analysis of DEGs were performed using clusterProfiler package (Yu et al., \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2012\u003c/span\u003e), KOBAS v3.0 (Xie et al., \u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e2011\u003c/span\u003e), respectively. The calculated \u003cem\u003ep\u003c/em\u003e-value goes through Bonferroni Correction, taking \u003cem\u003eP\u003c/em\u003e.adjust value\u0026thinsp;\u0026le;\u0026thinsp;0.05 as a threshold.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eQuantitative real-time PCR analysis\u003c/h2\u003e \u003cp\u003eTotal RNAs were extracted as described above. Each RNA sample was treated by 5\u0026times; PrimeScript\u0026trade; RT Master Mix (Takara) to remove residual genomic DNA and reverse transcribed into cDNA. All primer pairs (Additional file 1: 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 Real-time PCR System (Bio-Rad, CA, USA) in 25 \u0026micro;L reactions containing 12.5 \u0026micro;L SYBR Premix Ex Taq II (2\u0026times;) (Takara Bio), 1 \u0026micro;L each primer (10 \u0026micro;M), 2 \u0026micro;L PCR template (cDNA) and 8.5 \u0026micro;L of nuclease-free water. Amplification was performed with an initial denaturation at 95\u0026deg;C for 5 min, followed by 40 cycles of 95\u0026deg;C for 10 s, 60\u0026deg;C for 30 s and 72\u0026deg;C 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\u003csup\u003e\u0026minus;(ΔΔCt)\u003c/sup\u003e method with \u003cem\u003eβ-actin\u003c/em\u003e gene as the reference control. Data were analyzed statistically with SPSS 20 (IBM, Chicago, IL, USA) by t-test. Thresholds for statistical significance were set at \u003cem\u003eP\u003c/em\u003e \u003c0.05 (significant) and \u003cem\u003eP\u003c/em\u003e \u003c0.01 (highly significant).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003eCo-expression analysis of protein-protein interaction (PPI) network analysis\u003c/h2\u003e \u003cp\u003eThe DEGs were imported into the Search Tool for the Retrieval of Interacting Genes (STRING, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://string-db.org\u003c/span\u003e\u003c/span\u003e) database to obtain the PPI information (Szklarczyk et al., \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). Only validated interactions with a composite score greater than 0.4 were considered significant. Cytoscape 3.6.0 software was used to construct PPI network and count the number of nodes in DEGs (Franz et al., \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). The node genes with node degree above 10 were selected as the key target genes.\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cdiv class=\"Section2\" id=\"Sec11\"\u003e\n \u003ch2\u003eOverview of the RNA-Seq data\u003c/h2\u003e\n \u003cp\u003eTo better understand skin color variation of red tilapia during and after overwintering, the mRNA libraries of PB, WP and BP were determined and analyzed by Illumina sequencing technology. In Table \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e, we presented the Q30 percentage, GC percentage and the other indexes to describe the libraries. In total, an average of 42,721,496, 43,921,378 and 43,386,254 raw reads were obtained from the BP, PB and WP libraries, respectively. After filtering the low-quality reads and removing adaptor sequences, the average of 42,664,181, 43,888,258 and 43,347,484 clean reads were retrieved for further analysis. The percentage of G\u0026thinsp;+\u0026thinsp;C content and Q30 ratio was an average of 48.23% and 94.27%, indicated a high-quality sequence. All raw transcriptome data were submitted to the NCBI Short Read Archive (SRA) database with the accession numbers PRJNA690595.\u003c/p\u003e\n \u003cdiv class=\"gridtable\"\u003e\u0026nbsp;\u003ctable border=\"1\" id=\"Tab1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003eThe specific statistics for each library sequencing and quality control\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003ccolgroup cols=\"7\"\u003e\u003c/colgroup\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eSample\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eRaw reads\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eRaw base(G)\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eClean reads\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eClean base(G)\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eQ30(%)\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eGC(%)\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBP-1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40,313,772\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.05\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40,254,280\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.00\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e95.23\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48.31\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBP-2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e44,376,218\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.66\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e44,317,948\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.61\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e95.01\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e47.34\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBP-3\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40,589,108\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.09\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40,530,756\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.05\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e95.01\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e47.84\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBP-4\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e45,606,884\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.84\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e45,553,738\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.80\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e95.12\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48.52\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ePB-1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e46,324,834\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.95\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e46,289,448\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.88\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e93.32\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e47.41\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ePB-2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,749,772\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.56\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,722,494\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.49\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e93.97\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48.04\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ePB-3\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40,425,896\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.06\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e40,402,734\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e5.98\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e93.81\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48.37\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ePB-4\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e45,185,008\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.78\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e45,138,356\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.71\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e94.34\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e47.72\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWP-1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,218,496\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.48\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,170,434\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.42\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e93.80\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e49.02\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWP-2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,187,756\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.48\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,150,158\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.42\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e93.63\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e49.12\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWP-3\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,362,090\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.50\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,324,748\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.46\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e93.92\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48.15\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eWP-4\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,776,672\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.57\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,744,596\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.51\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e94.04\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48.86\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eAverage\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,343,042\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.50\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e43,299,974\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e6.44\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e94.27\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e48.23\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n \u003c/div\u003e\n \u003cp\u003eTo assess the quality of sequencing and reassembly, all clean reads were mapped to Nile tilapia (\u003cem\u003eOreochromis niloticus\u003c/em\u003e) 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 (Additional file 1: 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% (Additional file 1: Table S3).\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv class=\"Section2\" id=\"Sec12\"\u003e\n \u003ch2\u003eAnalysis of gene expression level of the red tilapia transcriptome\u003c/h2\u003e\n \u003cp\u003eThe RPKM method was used to estimate gene expression. The distribution of RPKM values for each sample were shown in Additional file 1: 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\u0026thinsp;\u0026le;\u0026thinsp;RPKM\u0026thinsp;\u0026le;\u0026thinsp;1 were the most, while the number of genes with RPKM\u0026thinsp;\u0026ge;\u0026thinsp;100 were less than 1%.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv class=\"Section2\" id=\"Sec13\"\u003e\n \u003ch2\u003eExpression profiling of mRNAs\u003c/h2\u003e\n \u003cp\u003eTo further understand mRNA expression of differences 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 the further analysis (SPM\u0026thinsp;\u0026gt;\u0026thinsp;0.9, Additional file 2: Table S1). In detail, there were 119, 294, and 53 specific mRNAs in BP, PB, and WP skins respectively. The KEGG results of specific mRNAs showed that metabolic pathways, ribosome biogenesis in eukaryotes and oxidative phosphorylation were dominant pathway in BP skin and regulation of actin cytoskeleton, melanogenesis, tight junction and tyrosine metabolism were dominant in PB skin (Additional file 2: Table S2). Furthermore, we also analyzed the top 10 abundance mRNAs of differences skin color in red tilapia. As shown in Table \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e, \u003cem\u003eGranulin\u003c/em\u003e and \u003cem\u003etat\u003c/em\u003e gene were abundant from WP skin and other pigment-related genes were abundantly expressed in them, including \u003cem\u003eoca2\u003c/em\u003e and \u003cem\u003eslc45a4\u003c/em\u003e. \u003cem\u003eTyrp1b\u003c/em\u003e gene was the most abundant in PB skin and melanin synthesis gene accounts for the largest proportion of abundance expressed genes, such as \u003cem\u003etyrp1b\u003c/em\u003e, \u003cem\u003epmelb\u003c/em\u003e, \u003cem\u003etyr\u003c/em\u003e, \u003cem\u003epmela\u003c/em\u003e, and \u003cem\u003etyrp1a\u003c/em\u003e. In addition, \u003cem\u003ebaxa\u003c/em\u003e gene showed dominantly expression in BP skin.\u003c/p\u003e\n \u003cdiv class=\"gridtable\"\u003e\u0026nbsp;\u003ctable border=\"1\" id=\"Tab2\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 2\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003eTen most abundant genes of three difference colors in red tilapia during and after overwintering\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003ccolgroup cols=\"3\"\u003e\u003c/colgroup\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eWP\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003ePB\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eBP\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003egranulin\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003etyrp1b\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ebax\u0026alpha;\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003etat\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ecavin2a\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003esi:ch73-86n18.1\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eENSONIG00000007142\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ezgc:101810\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ecrtac1a\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eoca2\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003epmelb\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eENSONIG00000040753\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eENSONIG00000040391\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003etyr\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003esdhb\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eoni-mir-24a-4\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003emyh9a\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eppdpfa\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ecol10a1a\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eENSONIG00000019137\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eloxa\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eENSONIG00000039537\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eaqp3\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eENSONIG00000009288\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eslc45a4\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003epmela\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eENSONIG00000018504\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003emap3k7cl\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003etyrp1a\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eatg7\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n \u003c/div\u003e\n \u003cp\u003e\u003cbr\u003e\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv class=\"Section2\" id=\"Sec14\"\u003e\n \u003ch2\u003eDifferential gene expression (DEGs) identified in different skin patterns\u003c/h2\u003e\n \u003cp\u003eIn comparative transcriptome analysis, many genes showed different expression levels in three skin color samples. Under the criteria of FDR\u0026thinsp;\u0026le;\u0026thinsp;0.05 and |logFC| \u0026ge;1, the volcano plots of three pairwise comparisons (BP_PB, BP_WP, and PB_WP) revealed the expression trend of each pair (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ea). We also constructed a histogram of DEGs in the three skin tissues (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eb). Compared with the PB skin, there were 2,528 DEGs in BP skin, of which 1,420 were up-regulated and 1,108 were down-regulated. A total of 1,091 DEGs were up-regulated in BP skin compared with WP skin, while 833 DEGs were down-regulated. There were 1,939 DEGs displaying greater abundance in PB skin compared with WP skin, of which 1,387 DEGs were up-regulated and 552 genes were down-regulated. Among these genes, 22 DEGs were detected as shared genes in each comparison group, in which 11 were known DEGs and seemed to play a key role in the color variation process (Fig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003ec). They were \u003cem\u003est3gal1\u003c/em\u003e, \u003cem\u003esi:ch73-86n18.1\u003c/em\u003e, \u003cem\u003eplxna4\u003c/em\u003e, \u003cem\u003emalb\u003c/em\u003e, \u003cem\u003efabp11a\u003c/em\u003e, \u003cem\u003eplcb4\u003c/em\u003e, \u003cem\u003esdhb\u003c/em\u003e, \u003cem\u003esi:dkey-65b12.6\u003c/em\u003e, \u003cem\u003esi:ch211-157c3.4\u003c/em\u003e, \u003cem\u003eagr1\u003c/em\u003e, and \u003cem\u003eaqp3\u003c/em\u003e, respectively.\u003c/p\u003e\n \u003cp\u003eTo 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. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e, 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. The results showed that the reliability of the sequencing result was high.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv class=\"Section2\" id=\"Sec15\"\u003e\n \u003ch2\u003eFunctional enrichment analysis of DEGs\u003c/h2\u003e\n \u003cp\u003eThe top GO function enrichment terms for the pairwise comparisons among three samples were shown in Additional file 2: Table S3. After GO annotation, all DEGs were classified into different biological processes, molecular function and cellular component. In each comparison, the top 50 of GO categories were selected in three different categories. The detailed annotations of each category were depicted in Additional file 1: Fig. S1a-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 terms 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 a further study.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv class=\"Section2\" id=\"Sec16\"\u003e\n \u003ch2\u003eKEGG analysis of the pathways\u003c/h2\u003e\n \u003cp\u003eTo further explore the biological functions of the DEGs, an enrichment analysis based on KEGG database was performed. A total of 141 KEGG pathways were listed in this study (Additional file 2: Table S4). The DEGs between the BP and PB skins were involved in ribosome, oxidative phosphorylation, ribosome biogenesis in eukaryotes, cardiac muscle contraction, RNA degradation, RNA polymerase and DNA replication were significantly enriched (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05). The DEGs were significantly enriched in some genetic information processing between the BP and WP 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, taurine and hypotaurine metabolism were significantly enriched between the PB and WP skins (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05). Since fish skin color was mainly correlated with the synthesis of different pigments, we were interested in 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, \u003cem\u003eetc.\u003c/em\u003e were identified, which were related to the skin color regulation and pigmentation (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003e).\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv class=\"Section2\" id=\"Sec17\"\u003e\n \u003ch2\u003eCandidate genes related to melanophore\u003c/h2\u003e\n \u003cp\u003eThe KEGG pathway analysis results showed that 32 candidate genes were involved in pigmentation-related pathways, such as melanogenesis, tyrosine metabolism, Wnt signaling pathway and MAPK signaling pathway. These genes may play a potential role in skin color variation of red tilapia during and after overwintering. The heatmap of these genes (Fig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003e) indicated that four genes including agouti signaling protein (\u003cem\u003easip\u003c/em\u003e), tyrosine amino transferase (\u003cem\u003etat\u003c/em\u003e), hydroxyphenylpyruvate hydroxylase \u003cem\u003e(hpdb)\u003c/em\u003e and forkhead transcription factor 3 (\u003cem\u003efoxd3\u003c/em\u003e) were up-regulated in WP and BP group compared with the PB group, while the rest genes including tyrosinase (\u003cem\u003etyr\u003c/em\u003e), tyrosinase-related protein 1 (\u003cem\u003etyrp1\u003c/em\u003e), telanocortin receptor 1 (\u003cem\u003emc1r\u003c/em\u003e), microphthalmia-associated transcription factor (\u003cem\u003emitf\u003c/em\u003e), premelanosome protein (\u003cem\u003epemel\u003c/em\u003e) \u003cem\u003eetc\u003c/em\u003e. were down-regulated in WP and BP group compared with PB group. In addition, some DEGs were involved in autophagy and apoptotic pathways such as \u003cem\u003ebax\u0026alpha;\u003c/em\u003e, \u003cem\u003ebeclin1\u003c/em\u003e and \u003cem\u003eatg7\u003c/em\u003e. The expression of these genes in BP skin of red tilapia was significantly higher than that of PB skin, which play an important role in black-to-pink skin color transformation in red tilapia.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv class=\"Section2\" id=\"Sec18\"\u003e\n \u003ch2\u003eCorrelation of candidate genes at protein levels\u003c/h2\u003e\n \u003cp\u003eBased on the candidate genes in PB_WP groups, we identified the mutual correlation of their protein products using the STRING online tool (Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003ea). The genes (or proteins) from the pigmentation-related pathway were integrated together, and the relationship of them were extensive and strong. Among them, \u003cem\u003easip\u003c/em\u003e, \u003cem\u003etat\u003c/em\u003e, \u003cem\u003ehpdb\u003c/em\u003e, \u003cem\u003efox3\u003c/em\u003e 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 degree were shown in Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eb. These hub genes included \u003cem\u003etyr, mc1r\u003c/em\u003e, \u003cem\u003eoca2\u003c/em\u003e, \u003cem\u003emitfb\u003c/em\u003e, \u003cem\u003eslc45a2\u003c/em\u003e, \u003cem\u003etyrp1b\u003c/em\u003e, \u003cem\u003edct\u003c/em\u003e, \u003cem\u003easip\u003c/em\u003e, \u003cem\u003ekit\u003c/em\u003e, \u003cem\u003ekitlg\u003c/em\u003e, \u003cem\u003epmela\u003c/em\u003e, \u003cem\u003epmelb\u003c/em\u003e, \u003cem\u003emitfa\u003c/em\u003e, \u003cem\u003eand egfra\u003c/em\u003e. Among these genes, \u003cem\u003etyr\u003c/em\u003e gene showed the highest node degree, which was 18.\u003c/p\u003e\n\u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eTo explore the different expression patterns among the three types skins, we performed differential gene expression analysis. When comparing PB skin to BP skin and WP skin, the results showed that more DEGs were up-regulated in PB 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, there genes rich in melanin synthesis were abundantly expressed in PB skin, while a few of pigment genes were abundantly expressed in WP. 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 BP_PB, BP_WP and PB_WP comparison groups, of which \u003cem\u003est3gal1\u003c/em\u003e, \u003cem\u003eplxna4\u003c/em\u003e, \u003cem\u003efabp11a\u003c/em\u003e and \u003cem\u003eaqp3\u003c/em\u003e plays a vital role in regulating cell proliferation, migration and invasion (Wu et al., \u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Wang et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). It has been reported that silencing of \u003cem\u003est3gal1\u003c/em\u003e gene suppresses melanoma invasion and significantly reduces the survive ability of aggressive melanoma cells of human in metastatic environment (Pietrobono et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Plexins family can functionally activate tyrosine kinase receptors in mammalian, such as MET, RON, HER2, and KDR (Swiercz et al., \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). Overexpression of \u003cem\u003eaqp3\u003c/em\u003e gene can promoted the proliferation and migration of human hepatocytes (Chen et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). \u003cem\u003eAqp3\u003c/em\u003e can reduce the differentiation and inhibit the apoptosis of stem cells in human through reducing the expressions of related genes in Wnt/GSK-3 β/β-catenin pathway (Liu et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). \u003cem\u003eFabp11a\u003c/em\u003e were probably involved in cellular uptake and transport of fatty acids, targeting of fatty acids to transport systems and several signalling pathways in \u003cem\u003eOryzias latipes\u003c/em\u003e (Parmar et al., \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). 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.\u003c/p\u003e \u003cp\u003eGO 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., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2013\u003c/span\u003e), Midas cichlids (\u003cem\u003eAmphilophus\u003c/em\u003e) (Henning et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2013\u003c/span\u003e) and common carp (Li et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). KEGG pathway analysis showed that many DEGs were significantly enriched in oxidative phosphorylation, ribosome, ribosome biogenesis in eukaryotes and cardiac muscle contraction in the BP_PB group and PB_WP group. Several studies have shown that high expression of ribosomal protein related genes was associated with the black color in mice (Skarnes et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). Four of the five highly expressed genes in pigment cells of zebrafish were ribosomal protein (Higdon et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2013\u003c/span\u003e). There were many DEGs participate in oxidative phosphorylation, cardiac muscle contraction signal pathways in \u003cem\u003ePristella maxillaris\u003c/em\u003e (Bian \u003cem\u003eet al\u003c/em\u003e., 2013) and \u003cem\u003eLutjanus erythropterus\u003c/em\u003e (Zhang et al., \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). In addition, the KEGG results of specific mRNAs showed that ribosome biogenesis in eukaryotes, oxidative phosphorylation were dominant pathways in BP skin. Similar results were found in comparative analysis of BP_PB skin and PB_WP skin. It is suggested that these pathways may play an important role in color variation in red tilapia.\u003c/p\u003e \u003cp\u003eWe also found that the DEGs of three skin colors of tilapia were mainly enriched in 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., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2009\u003c/span\u003e, Zhang et al., \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). For specific mRNA, the tyrosinase metabolism and melanogenesis pathways were dominant in PB skin. It showed that PB skin required more melanin synthesis than BP and WP skins. Meanwhile, we found that some DEGs between the BP and PB skin expressed abundantly in the process of apoptosis and autophagy pathways. The identification of genes in these pigmentation-related term and pathways are informative and are worthy of further study.\u003c/p\u003e \u003cp\u003eStudies have suggested that the mRNA expression levels of genes including \u003cem\u003etyr\u003c/em\u003e, \u003cem\u003etyrp1\u003c/em\u003e, \u003cem\u003emc1r\u003c/em\u003e, \u003cem\u003emitf\u003c/em\u003e, \u003cem\u003epemel\u003c/em\u003e were higher in PB 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. \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2010\u003c/span\u003e, Simon et al., \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2009\u003c/span\u003e). Therefore, TYR, TYRP1 and DCT are critical enzymes for the formation of melanin. Mutations or dysfunction of \u003cem\u003etyr\u003c/em\u003e or \u003cem\u003etyrp1\u003c/em\u003e genes lead to melanocyte death or extensive hypopigmentation in zebrafish (Krauss et al, \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). In our study, compared to WP skin samples, the expression levels of \u003cem\u003etyr\u003c/em\u003e and \u003cem\u003etyrp1\u003c/em\u003e were the highest in PB skin. This was also consistent with the pigmentation of red crucian carp (Zhang et al, \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). \u003cem\u003eMitf\u003c/em\u003e is a member transcription factor involved in the development of melanocytes, retinal cells, osteoclasts and mast cells (Minvielle et al., \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). It has been reported that \u003cem\u003emitf\u003c/em\u003e could directly regulate the expression of multiple genes necessary for the development of melanophores, including \u003cem\u003etyr\u003c/em\u003e, \u003cem\u003etyrp1\u003c/em\u003e, and \u003cem\u003edct\u003c/em\u003e (Cheli et al., \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). Compared with WP colors, PB skin color was caused by the increase of melanin content, suggesting that \u003cem\u003emitf\u003c/em\u003e may play a potential role in regulating the differentiation and development of melanocytes. \u003cem\u003eMc1r\u003c/em\u003e gene is a key gene in melanogenesis in animals. Alpha-melanocyte stimulating hormone (α-MSH) binds to \u003cem\u003emc1r\u003c/em\u003e, resulting in the decrease of cAMP level. Consequently, melanin biosynthesis process was triggered (Voisey et al., \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2001\u003c/span\u003e). Previous studies have shown that \u003cem\u003emc1r\u003c/em\u003e mutations were associated with skin color variation in many fish species, such as cavefish, guppy, zebrafish and koi carp (Gross et al., \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2009\u003c/span\u003e, Tezuka et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2011\u003c/span\u003e, Richardson \u003cem\u003eet al.\u003c/em\u003e, 2008, Dong et al., \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Similarly, we observed a significant difference in \u003cem\u003emc1r\u003c/em\u003e expression between the red tilapia of three skin colors used in this study. \u003cem\u003ePmel\u003c/em\u003e gene acts as a scaffold in the melanosome by creating a proteolytic fibrillary matrix where melanin is deposited (Solano \u003cem\u003eet al.\u003c/em\u003e, 2000). \u003cem\u003ePmel\u003c/em\u003e mutations promoted pigment dilution in many animals (Gutierrez et al., \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2007\u003c/span\u003e). Here \u003cem\u003epmel\u003c/em\u003e was significantly up-regulated in PB skin when compared to WP and BP skin samples. Similarly, we observed the top 10 abundance mRNAs in PB skin, including \u003cem\u003etyrp1b\u003c/em\u003e, \u003cem\u003etyr\u003c/em\u003e, \u003cem\u003epmelb\u003c/em\u003e, \u003cem\u003epmela\u003c/em\u003e, \u003cem\u003etyrp1a, etc.\u003c/em\u003e 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 PB_WP group. It was speculated that those proteins (or genes) could be co-regulated 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 \u003cem\u003etyr, mc1r\u003c/em\u003e, \u003cem\u003emitfb\u003c/em\u003e, \u003cem\u003etyrp1b\u003c/em\u003e, \u003cem\u003edct\u003c/em\u003e, \u003cem\u003epmela\u003c/em\u003e, \u003cem\u003epmelb\u003c/em\u003e, and \u003cem\u003emitfa\u003c/em\u003e, \u003cem\u003eetc\u003c/em\u003e, were consistent with the results of mRNA expression levels.\u003c/p\u003e \u003cp\u003eRegarding the black to pink stage, the mRNA expression levels of \u003cem\u003easip\u003c/em\u003e, \u003cem\u003etat\u003c/em\u003e, \u003cem\u003ehpdb\u003c/em\u003e and \u003cem\u003efoxd3\u003c/em\u003e were all up-regulated. \u003cem\u003eAsip\u003c/em\u003e gene product blocks melanogenesis by antagonizing the binding of α-MSH to \u003cem\u003emc1r. Asip\u003c/em\u003e mutations were associated with skin color variation in Psetta maxima, zebrafish and medaka (Ceinos et al., \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2015\u003c/span\u003e, Guillot et al., \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2012\u003c/span\u003e, Jose\u0026acute; et al., \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2005\u003c/span\u003e). In our study, we observed higher expression of \u003cem\u003easip\u003c/em\u003e gene in BP skin transcripts and lower expression of \u003cem\u003emc1r\u003c/em\u003e, which further establishes the role of \u003cem\u003easip\u003c/em\u003e as an antagonistic of the \u003cem\u003emc1r\u003c/em\u003e gene. TAT and HPDB catalyze the substrate tyrosine to form homogeneous acid (HGA). Homogentisate1, 2-dioxygenase catalysis HGA to produce melatonin. Higher level of \u003cem\u003etat\u003c/em\u003e and \u003cem\u003ehpda\u003c/em\u003e gene would directly reduce tyrosine level, thereby inhibited the synthesis of melanin (Zhang et al., \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). In addition, \u003cem\u003etat\u003c/em\u003e gene was the most abundant in WP skin, suggested that it affected the skin variation in red tilapia. \u003cem\u003eFoxd3\u003c/em\u003e 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 \u003cem\u003emitf\u003c/em\u003e at the early phase of neural crest migration (Thomas \u003cem\u003eet al.\u003c/em\u003e, 2009). In addition, overexpression of \u003cem\u003efoxd3\u003c/em\u003e in melb-a mouse melanoblasts blocked the expression of \u003cem\u003emitf\u003c/em\u003e (Lister et al., \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2001\u003c/span\u003e). \u003cem\u003eFoxd3\u003c/em\u003e was significantly up-regulated in BP skin samples compared to the PB skin, indicated that \u003cem\u003efoxd3\u003c/em\u003e might play a significant role in the black-to-pink color transformation in red tilapia.\u003c/p\u003e \u003cp\u003eIn addition, autophagy and apoptotic pathways were able to control the transition from black to pink in red tilapia. In detail, the mRNA level of apoptosis gene, such as \u003cem\u003ebaxα\u003c/em\u003e, was significantly increased in the body color transformation from the BP skin to the PB and WP skins. Meanwhile, the mRNA levels of \u003cem\u003ebeclin1\u003c/em\u003e and autophagy-related genes 7 (\u003cem\u003eatg7\u003c/em\u003e), as the autophagy genes, were all upregulated in the BP skin compared with PB skin. \u003cem\u003eBaxα\u003c/em\u003e, as one of the homologous proteins of BCL-2, could determine survival or death by an apoptotic stimulus. Overexpression of \u003cem\u003ebaxa\u003c/em\u003e may accelerate cell death (Oltvai et al., \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e1993\u003c/span\u003e). \u003cem\u003eBeclin1\u003c/em\u003e plays a key role in regulating autophagy and cell death by interacts with either BCL-2 or PI3k class III (Takacs-Vellai et al., \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2005\u003c/span\u003e). \u003cem\u003eAtg7\u003c/em\u003e activates the ubiquitin-like protein ATGL2, which binds to \u003cem\u003eatg5\u003c/em\u003e and extends the autophagic vesicle membrane. Whole body knock-out of \u003cem\u003eatg7\u003c/em\u003e in mice led to death within 24h after birth (Komatsu et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2005\u003c/span\u003e). \u003cem\u003eBaxα\u003c/em\u003e and \u003cem\u003eatg7\u003c/em\u003e gene were among the top 10 abundance mRNAs in BP skin, further confirmed that appearance of autophagy may lead to melanocyte reduction.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eIn 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 consist with 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.\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003eFunding projects\u003c/h2\u003e\n\u003cp\u003eThis study was supported by the National Natural Science Foundation-Youth Fund Project (31802290).\u003c/p\u003e\n\u003ch2\u003eConflicts of interest/Competing interests\u003c/h2\u003e\n\u003cp\u003eThe authors have no relevant financial or non-financial interests to disclose.\u003c/p\u003e\n\u003ch2\u003eCode availability\u003c/h2\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003ch2\u003eAvailability of data and material\u003c/h2\u003e\n\u003cp\u003eThe authors declare that all other data supporting the findings of this study are available within the article and its supplementary information files.\u003c/p\u003e\n\u003ch2\u003eAuthors\u0026apos; contributions\u003c/h2\u003e\n\u003cp\u003eZaijie Dong conceived the study; Wenbin Zhu provided the experimental materials; Lanmei provided the funding for the experiment; Bingjie Jiang performed the experiments and wrote the paper; Jianjun Fu provided technical assistance in experiments; Mingkun Luo revised the manuscript; Zaijie Dong reviewed the manuscript; All authors read and approved the manuscript.\u003c/p\u003e\n\u003ch2\u003eEthics approval\u003c/h2\u003e\n\u003cp\u003eThe sampling scheme and experimental protocols were subject to approval by the Bioethical Committee of Freshwater Fisheries Research Center (FFRC) of the Chinese Academy of Fishery Sciences (CAFS) (BC 2013863, 9/2013). The methods of samples handled and experimental procedures carried out in accordance with the guidelines for the care and use of animals for scientific purposes issued by the Ministry of Science and Technology, Beijing China (No.398, 2006).\u003c/p\u003e\n\u003ch2\u003eAvailability of data and material/ Data availability\u003c/h2\u003e\n\u003cp\u003eAll data generated or analysed during this study are included in this published article (and its supplementary information files).\u003c/p\u003e\n\u003ch2\u003eConsent to participate\u003c/h2\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003ch2\u003eConsent for publication\u003c/h2\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eAndrews S (2014) FastQC: a quality control tool for high throughput sequence data. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.bioinformatics.babraham.ac.uk/projects/fastqc/\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBian FF, Yang XF, Ou ZJ et al (2019) Morphological characteristics and comparative transcriptome analysis of three different phenotypes of \u003cem\u003ePristella maxillaris\u003c/em\u003e. Front Genet 10:698\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBraasch I, Liedtke D, Volff JN, Schartl M (2010) Pigmentary function and evolution of \u003cem\u003etyrp1\u003c/em\u003e gene duplicates in fish. Pigment Cell Melanoma Research 22(6):839\u0026ndash;850\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCeinos RM, Guillot R, Kelsh RN et al (2015) Pigment patterns in adult fish result from superimposition of two largely independent pigmentation mechanisms. Pigment Cell Melanoma Research 28(2):196\u0026ndash;209\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCheli Y, Ohanna M, Ballotti R, Bertolotto C (2010) Fifteen-year quest for microphthalmia-associated transcription factor target genes. Pigment Cell Melanoma Research 23(1):27\u0026ndash;40\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen G, Shi Y, Liu M, Sun J et al (2018) circHIPK3 regulates cell proliferation and migration by sponging miR-124 and regulating AQP3 expression in hepatocellular carcinoma. Cell Death Dis 9(2):175\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen YX, Chen YS, Shi CM et al (2018) Soapnuke: a mapreduce acceleration supported software for integrated quality control and preprocessing of high-throughput sequencing data. Oxford Open 7(1):gix120\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDong ZJ, Luo MK, Wang LM et al (2020) MicroRNA-206 regulation of skin pigmentation in koi carp (\u003cem\u003eCyprinus carpio\u003c/em\u003e L.). Front Genet 11:47\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFrancisco S, Mar\u0026iacute;a ME, Celia JC, Simon PH et al (2000) New insights on the structure of the mouse silver locus and on the function of the silver protein. Pigment Cell Melanoma Research 13(8):118\u0026ndash;124\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFranz M, Lopes CT, Huck G, Dong Y, Sumer O, Bader GD (2016) Cytoscape.Js a graph theory library for visualisation and analysis. Bioinformatics 32(2):309\u0026ndash;311\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFujimura N, Taketo MM, Mori M et al (2009) Spatial and temporal regulation of Wnt/β-catenin signaling is essential for development of the retinal pigment epithelium. Dev Biol 334(1):31\u0026ndash;45\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGouveia L, Rema P (2005) Effect of microalgal biomass concentration and temperature on ornamental goldfish (\u003cem\u003eCarassius auratus\u003c/em\u003e) skin pigmentation. Aquac Nutr 11(1):19\u0026ndash;23\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGross JB, Borowsky R, Tabin CJ (2009) A novel role for \u003cem\u003emc1r\u003c/em\u003e in the parallel evolution of depigmentation in independent populations of the cavefish \u003cem\u003eAstyanax Mexicanus\u003c/em\u003e. Plos Genetics 5(1):1000326\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGuillot R, Ceinos RM, Cal R et al (2012) Transient ectopic overexpression of agouti-signalling protein 1 (\u003cem\u003easip1\u003c/em\u003e) induces pigment anomalies in flatfish. Plos One 7(12):e48526\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGupta MV, Acosta BO (2004) A review of global tilapia farming practices. Aquaculture Asia 9(1):7\u0026ndash;12\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGutierrez GB, Wiener P, Williams JL (2007) Genetic effects on coat color in cattle: dilution of eumelanin and\u0026acute; phaeomelanin pigments in an F2-backcross \u003cem\u003eCharolais\u003c/em\u003e \u0026times;\u003cem\u003eHolstein\u003c/em\u003e population. BMC Genet 8(1):56\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHenning F, Jones JC, Franchini P et al (2013) Transcriptomics of morphological color change in polychromatic Midas cichlids. BMC Genom 14(1):171\u0026ndash;171\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHigdon CW, Mitra RD, Johnson SL et al (2013) Gene expression analysis of zebrafish melanocytes, iridophores, and retinal pigmented epithelium reveals indicators of biological function and developmental origin. Plos One 8(7):e67801\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHuang R, Kitajama C, Ke C (1997) A Study on Induced factor of White Skin of Flat Fish \u003cem\u003eParalichthys Olivaceus\u003c/em\u003e (Temminck et Schlegel). Modern Fisheries Information 4(9):21\u0026ndash;23\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHubbard JK, Uy JA, Hauber ME, Hoekstra HE, Safran RJ (2010) Vertebrate pigmentation: from underlying genes to adaptive function. Trends Genet 26:231\u0026ndash;239\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJennifer R, Pia RL, Natalie LR, Julia RD et al (2008) \u003cem\u003eMc1r\u003c/em\u003e Pathway regulation of zebrafish melanosome dispersion. Zebrafish 5(4):289\u0026ndash;295\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJiang YL, Zhang SH, Xu J et al (2014) Comparative transcriptome analysis reveals the genetic basis of skin color variation in common carp. Plos One 9(9):e108200\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJose\u0026acute; MCR, Tatjana H, Helgi BS, Richard EP (2005) Gene structure of the goldfish agouti-signaling protein: a putative role in the dorsal-ventral pigment pattern of fish. Endocrinology 146(3):1597\u0026ndash;1610\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim D, Langmead B, Salzberg SL (2015) Hisat: a fast spliced aligner with low memory requirements. Nature methods 12(4):357\u0026ndash;360\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKomatsu M, Waguri S, Ueno T et al (2005) Impairment of starvation-induced and constitutive autophagy in \u003cem\u003eatg7\u003c/em\u003e deficient mice. J Cell Biol 169(3):425\u0026ndash;434\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKrauss J, Geiger-Rudolph S, Koch I, Nusslein-Volhard C, Irion U (2015) A dominant mutation in \u003cem\u003etyrp1a\u003c/em\u003e leads to melanophore death in zebrafish. Pigm Cell Melanoma Research 27(5):827\u0026ndash;830\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi XJ, Li SF, Feng JH et al (2003) Preliminary study on salinity tolerance of Israel red tilapia. Journal of Shanghai Fisheries University 12(3):205\u0026ndash;208\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi XM, Song YN, Xiao GB, Zhu BH, Xu GC, Sun MY et al (2015) Gene expression variations of red-white skin coloration in common carp (\u003cem\u003eCyprinus carpio\u003c/em\u003e). Int J Mol Sci 16:21310\u0026ndash;21329\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiao Y, Gordon K, Smyth, Shi W (2014) Featurecounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7):923\u0026ndash;930\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLister JA, Close J, Raible DW (2001) Duplicate \u003cem\u003emitf\u003c/em\u003e genes in zebrafish: complementary expression and conservation of melanogenic potential. Dev Biol 237(2):333\u0026ndash;344\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu CS, Liu LY, Zhang Y, Jing HY (2020) Molecular mechanism of \u003cem\u003eaqp3\u003c/em\u003e in regulating differentiation and apoptosis of lung cancer stem cells through Wnt/GSK-3β/β-Catenin pathway. Jbuon 25(4):1714\u0026ndash;1720\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLuo MK, Lu GQ, Yin HR et al (2021) Fish pigmentation and coloration: molecular mechanisms and aquaculture perspectives. Reviews in Aquaculture 00:1\u0026ndash;18\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMinvielle F, Bed'Hom B, Coville JL et al (2010) The \"silver\" Japanese quail and the \u003cem\u003emitf\u003c/em\u003e gene: causal mutation, associated traits and homology with the \"blue\" chicken plumage. BMC Genet 11(1):1\u0026ndash;7\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNakano T, Ono K, Takeuchi M (1992) Levels of zinc, iron, and cupper in the skin of abnormal pigmented Japanese flounder. Nippon Suisan Gakkaishi 58(11):2207\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eOltvai ZN, Milliman CL, Korsmeyer SJ (1993) Bcl-2 heterodimerizes in vivo with a conserved homolog, \u003cem\u003ebax\u003c/em\u003e, that accelerates programed cell death. Cell 74(4):609\u0026ndash;619\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePan JB, Hu SC, Wang H et al (2012) PaGeFinder: quantitative identification of spatiotemporal pattern genes. Bioinformatics 28(11):1544\u0026ndash;1545\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eParmar MB, Venkatachalam AB, Wright JM (2012) The evolutionary relationship of the transcriptionally active \u003cem\u003efabp11a\u003c/em\u003e (intronless) and \u003cem\u003efabp11b\u003c/em\u003e genes of medaka with \u003cem\u003efabp11\u003c/em\u003e genes of other teleost fishes. Febs Journal 279(13):2310\u0026ndash;2321\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePavlidis M, Karkana M, Fanouraki E et al (2008) Environmental control of skin color in the red porgy, \u003cem\u003ePagrus pagrus\u003c/em\u003e. Aquac Res 39(8):837\u0026ndash;849\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePietrobono S, Anichini G, Sala C et al (2020) ST3GAL1 is a target of the SOX2-GLI1 transcriptional complex and promotes melanoma metastasis through AXL. Nat Commun 11(1):5865\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePradeep PJ, Srijaya TC, Hassan A et al (2014) Optimal conditions for cold-shock induction of triploidy in red tilapia. Aquacult Int 22(3):1163\u0026ndash;1174\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRobinson MD, Mccarthy DJ, Smyth GK (2010) Edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139\u0026ndash;140\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSimon DJ, Peles D, Wakamatsu K, Ito S (2009) Current challenges in understanding melanogenesis: bridging chemistry, biological control, morphology and function. Pigment Cell Melanoma Research 22(5):563\u0026ndash;579\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSkarnes WC, Rosen B, West AP et al (2011) A conditional knockout resource for the genome-wide study of mouse gene function. Nature 474(7351):337\u0026ndash;342\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSwiercz JM, Worzfeld T, Offermanns S (2008) ErbB-2 and met reciprocally regulate cellular signaling via plexin-B1. J Biol Chem 283(4):1893\u0026ndash;1901\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSzklarczyk D, Franceschini A, Wyder S et al (2015) STRING v10: protein-protein networks, integrated over the tree of life. Nucleic Acids Res 43:D447\u0026ndash;D452\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTakacs-Vellai K, Vellai T, Puoti A, Passannante M et al (2005) Inactivation of the autophagy gene \u003cem\u003ebec-1\u003c/em\u003e triggers apoptotic cell death in \u003cem\u003eC. elegans\u003c/em\u003e. Curr Biol 15(16):1513\u0026ndash;1517\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTezuka A, Yamamoto H, Yokoyama J et al (2011) The \u003cem\u003emc1r\u003c/em\u003e gene in the guppy (\u003cem\u003ePoecilia reticulata\u003c/em\u003e): Genotypic and phenotypic polymorphisms. BMC Research Notes 4(1):31\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eThomas AJ, Erickson CA (2009) \u003cem\u003eFoxd3\u003c/em\u003e regulates the lineage switch between neural crest-derived glial cells and pigment cells by repressing \u003cem\u003emitf\u003c/em\u003e through a non-canonical mechanism. Development 136(11):1849\u0026ndash;1858\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVoisey J, Box NF, Daal AV (2001) A polymorphism study of the human agouti gene and its association with \u003cem\u003emc1r\u003c/em\u003e. Pigment Cell Res 14(4):264\u0026ndash;267\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang LG, Wang SQ, Li W (2012) RSeQC: quality control of RNA-seq experiments. Bioinformatics 28(16):2184\u0026ndash;2185\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang LM, Song FB, Zhu WB et al (2018) Effects of temperature on body color of Malaysian red tilapia during overwintering period. Journal of Fisheries of China 42(1):72\u0026ndash;79\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang XY, Yang J, Yao Y, Shi X, Yang GS, Li X (2020) \u003cem\u003eAqp3\u003c/em\u003e facilitates proliferation and adipogenic differentiation of porcine intramuscular adipocytes. Genes 11(4):453\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWardani WW, Alimuddin A, Junior MZ et al (2020) Growth performance, robustness against stress, serum insulin, \u003cem\u003eigf-1\u003c/em\u003e and \u003cem\u003eglut4\u003c/em\u003e gene expression of red tilapia (\u003cem\u003eOreochromis sp.\u003c/em\u003e) fed diet containing graded levels of creatine. Aquaculture Nutrition 27(1):274\u0026ndash;286\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWu X, Zhao JD, Ruan YY et al (2018) Sialyltransferase \u003cem\u003est3gal1\u003c/em\u003e promotes cell migration, invasion, and TGF-β1-induced EMT and confers paclitaxel resistance in ovarian cancer. Cell Death Dis 9(11):1102\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXie C, Mao XZ, Huang JJ, Yang D et al (2011) KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res 39(suppl_2):W316\u0026ndash;W322\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYu GC, Wang LG, Han YY, He QY (2012) Clusterprofiler: an r package for comparing biological themes among gene clusters. Omics-a Journal of Integrative Biology 16(5):284\u0026ndash;287\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang Y, Wang L, Zhang S, Yang H, Tan H (2008) HmgA, transcriptionally activated by \u003cem\u003ehpda\u003c/em\u003e, influences the biosynthesis of actinorhodin in \u003cem\u003eStreptomyces Coelicolor\u003c/em\u003e. FEMS Microbiol Lett 280(2):219\u0026ndash;225\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang YP, Wang ZD, Guo YS et al (2015) Morphological characters and transcriptome profiles associated with black skin and red skin in crimson snapper (\u003cem\u003eLutjanus erythropterus\u003c/em\u003e). Int J Mol Sci 16(11):26991\u0026ndash;27004\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang YQ, Jin H et al (2017) Comparative transcriptome analysis of molecular mechanism underlying gray-to-red body color formation in red crucian carp (\u003cem\u003eCarassius auratus\u003c/em\u003e, red var.). Fish Physiol Biochem 43(6):1387\u0026ndash;1398\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhu WB, Wang LM, Dong ZJ et al (2016) Comparative transcriptome analysis identifies candidate genes related to skin color differentiation in red tilapia. Sci Rep 6(1):31347\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":true,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"fish-physiology-and-biochemistry","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"fish","sideBox":"Learn more about [Fish Physiology and Biochemistry](https://www.springer.com/journal/10695)","snPcode":"10695","submissionUrl":"https://submission.nature.com/new-submission/10695/3","title":"Fish Physiology and Biochemistry","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"red tilapia, skin color, transcriptome, overwintering, differentially expressed transcripts","lastPublishedDoi":"10.21203/rs.3.rs-730510/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-730510/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe commercial value of red tilapia is hampered by variations in skin color during overwintering. In this study, three types skin of red tilapia, including the skin remained pink color during and after overwintering (WP), the skin changed from pink color to black color during overwintering and remained black color after overwintering (PB), and the skin changed from pink color to black color during overwintering but recovered to pink color when the temperature rose after overwintering (BP), were used to analyze their molecular mechanisms of color variation. The transcriptome results revealed that the PB, WP, and BP libraries had 42, 43, and 43\u0026nbsp;million clean reads, respectively. The top 10 abundance mRNAs and specific mRNAs (specificity measure SPM\u0026thinsp;\u0026gt;\u0026thinsp;0.9) were screened. After comparing intergroup gene expression levels, there were 2528, 1924 and 1939 differentially expressed genes (DEGs) between BP and PB, BP and WP, and PB and WP, 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 \u003cem\u003etyrp1\u003c/em\u003e, \u003cem\u003etyr\u003c/em\u003e, \u003cem\u003epmel\u003c/em\u003e, \u003cem\u003emitf\u003c/em\u003e, \u003cem\u003emc1r\u003c/em\u003e, \u003cem\u003easip\u003c/em\u003e, \u003cem\u003etat\u003c/em\u003e, \u003cem\u003ehpdb\u003c/em\u003e and \u003cem\u003efoxd3\u003c/em\u003e, might play a potential role in pigmentation. Additionally, the co-expression patterns of genes were detected within the pigment-related pathways by PPI network from PB_WP group. Furthermore, DEGs from the apoptosis and autophagy pathways, such as \u003cem\u003ebaxα\u003c/em\u003e, \u003cem\u003ebeclin1\u003c/em\u003e, and \u003cem\u003eatg7\u003c/em\u003e, 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.\u003c/p\u003e","manuscriptTitle":"Transcriptome Analysis of Skin Color Variation During and after Overwintering of Malaysian Red Tilapia","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2021-07-27 22:22:08","doi":"10.21203/rs.3.rs-730510/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"editorInvitedReview","content":"","date":"2021-08-01T09:01:39+00:00","index":0,"fulltext":""},{"type":"reviewersInvited","content":"","date":"2021-07-21T13:41:05+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2021-07-21T02:27:10+00:00","index":"","fulltext":""},{"type":"submitted","content":"Fish Physiology and Biochemistry","date":"2021-07-18T12:21:56+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"fish-physiology-and-biochemistry","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"fish","sideBox":"Learn more about [Fish Physiology and Biochemistry](https://www.springer.com/journal/10695)","snPcode":"10695","submissionUrl":"https://submission.nature.com/new-submission/10695/3","title":"Fish Physiology and Biochemistry","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"9b9d75ae-2a6c-4e5c-8263-03a298d6b41e","owner":[],"postedDate":"July 27th, 2021","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[{"id":6016805,"name":"Physiology"},{"id":6016806,"name":"General Biochemistry"}],"tags":[],"updatedAt":"2022-04-13T20:49:07+00:00","versionOfRecord":{"articleIdentity":"rs-730510","link":"https://doi.org/10.1007/s10695-022-01073-5","journal":{"identity":"fish-physiology-and-biochemistry","isVorOnly":false,"title":"Fish Physiology and Biochemistry"},"publishedOn":"2022-04-13 20:49:07","publishedOnDateReadable":"April 13th, 2022"},"versionCreatedAt":"2021-07-27 22:22:08","video":"","vorDoi":"10.1007/s10695-022-01073-5","vorDoiUrl":"https://doi.org/10.1007/s10695-022-01073-5","workflowStages":[]},"version":"v1","identity":"rs-730510","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-730510","identity":"rs-730510","version":["v1"]},"buildId":"-HB7Z8yhvgn0wM9Nzuekk","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}