Comparative Transcripto me Analysis Depicts Candidate Genes Involved in Skin Color Differentiation in Red Tilapia

This study aimed to investigate the genetic molecular mechanism of body color differentiation and variation of red tilapia, selecting the main genes related to the variation and cultivating the pure and stable red tilapia variety. The effects of different temperature treatments on body color and survival of Guam red tilapia, pearl white red tilapia and Florida red tilapia were compared. Besides, comparative transcriptome analysis was used to screen the candidate genes linked to the skin color differentiation of pearl white red tilapia. Among them, the body color of Guam red tilapia changed when the water temperature dropped to 16 − 14 °C, and continued to drop to about 11 °C, it was discolored in a large area reaching above 90%. According to the differential analysis: Tyrosine Kinase STYK1, HSP70, HSP30 and Transcription factor Sp6 expressions were signicantly increased in the low temperature group, while MC1R, Transcription factor (MafB, jun-D, AP-1, E2F5, ETV6, Sp9, Sp7, E2F1, Sp4) expressions were notably decreased. to different colors: tyr, tyrp1, silv, sox10, slc24a5, CBS and slc7a11 The results showed that the miRNAs associated with red tilapia color mainly included slc7a11, mc1r and asip, predicating that miR-138-5p and miR-722 played an important role in the regulation of pigmentation. However, the molecular mechanism of body color differentiation remains unclear.


Background
Red Tilapia (genus Oreochromis), as an essential tropical freshwater species and one of the most important edible sh in the world, was largely identi ed by its uniform red skin and bright pink peritoneum [1]. Red tilapia is widely cultured because of its rapid growth and tender meat [2].
Fish belongs to poikilothermic animals whose body temperature varies with the environmental temperature [3]. It is well known that water temperature affects the survival, growth, metabolism, reproduction of sh [4]. Tilapia is a warm-water sh with a temperature tolerance of 6-42 °C, but the tolerance to low temperature is in uenced by water quality, sh age and health status [5]. The extreme freezing weather in China always caused great economic loss to the tilapia farming industry [6]. A larger number of red tilapia and wintering seedings that were ready for market were frozen and die under 8-10 °C of water temperature [7]. Additionally, the time of seedling release of the early spring red tilapia with low cold tolerance would be delayed, leading to its market price volatile [8]. Thus, the supply and demand relationship of red tilapia was unbalanced, which had a major impact on the red tilapia breeding industry [9]. Therefore, the physiological study on the cold tolerance of red tilapia can provide data supporting for the breeding of cold-resistant red tilapia in the future.
Differentiation during genetic selection and variation in body color has been a growing limitation to the commercial value of red tilapia [10]. The differentiation of body color in red tilapia is irreversible, while the change of skin color is reversible with the change of ambient temperature [7]. The color of animals' skin is associated with the synthesis of different pigments secreted by pigment cells, heredity, environment, nutrition or physiology [11]. Previous study has shown that the two types of melanin produced by melanocytes extremely in uence on body color phenotype, one of which is eumelanin being responsible for producing black and brown phenotypes, the other is brown melanin contributing to producing the yellow and red phenotypes [12]. The formation of melanin is a complex process requiring the maturation of melanocytes, the synthesis and transportation of melanin, in which regulatory factors and signaling molecules are involved in each stage [13]. The adenylate cyclase pathway, the protein kinase c pathway and the tyrosine pathway are three molecules and transduction pathways in melanocytes that synthesize melanin [14]. The activation of each pathway is conducive to the production of melanin, in which adenylate cyclase plays a key role in the melanin synthesis pathway [15]. The common genetic signaling pathways of two pigments were also discovered [16]. The eumelanin is synthesized from tyrosine in vivo. Under the action of tyrosinase, tyrosine forms a eumelanin through a series of oxidation and catalytic polymerization [17]. It turns to the synthetic pathway of brown melanin after the reduction of tyrosinase activity [17]. Compared with other sh, little is known to the body color differentiation in red tilapia. zhu et al. used comparative transcriptome sequencing to nd the following genes related to skin color in red tilapia of three different colors: tyr, tyrp1, silv, sox10, slc24a5, CBS and slc7a11 [18]. The results showed that the miRNAs associated with red tilapia color mainly included slc7a11, mc1r and asip, predicating that miR-138-5p and miR-722 played an important role in the regulation of pigmentation. However, the molecular mechanism of body color differentiation remains unclear.
Guam strain of red tilapia is beautiful in color which does not change during wintering period with slower growth rate and unstable genetics [7]. The pearly white tilapia is pink-white in color and its growth characteristics are opposite to the Guam strain of red tilapia [19]. The value of red tilapia is determined by the proportion of red area of the sh. And the breeding of red tilapia strains with bright red body and rapid growth needs to be solved for expanding its breeding industry and meeting the market demand [20].
Thus, in this study, transcriptome sequencing was performed on the tissue samples of pearly white tilapia n strips at low temperature and normal temperature to screen out the genes related to body color variation of red tilapia and provide a theoretical basis for the breeding of ne red tilapia strains.

Results
Effect of temperature on apparent body color of red tilapia When the temperature dropped to 15 °C, body color of the PR and FR with the various body weight began to change; body color of the PR and FR changed over large areas when the temperature was decreased to 11 °C, and the discoloration ratio is all over 90%. When the temperature dropped to 10.2 °C, GR (10 g) began to change color, as temperature continued to reduce, the total body color change ratio was about 6%. However, for GR (80 g and 120 g), the decrease of temperature has no effect on the change of body color (Table 1).

Transcriptome sequencing and assembly
In the gradient cooling culture process, the apparent body color of WR changed obviously. Thus, 3 n tissues of WR in the normal water temperature (20 °C) and cooling process (after 24 h of body color change) were selected for transcriptome sequencing. In total, a mean of 8, 194, 153 paired-end (PE) clean reads for each library were obtained after data ltering and the sequencing quality was high with Q30 ratio larger than 96% for all samples ( Table 2). Read_num: Number of sequences sequenced; Base_num (bp): Number of bases sequenced; Length (bp): Sequencing read length; GC%: Ratio of G base and C base to total base number; Q20: Ratio of the number of bases having a mass value greater than 20 (error rate less than 1%) to the total number of bases; Q30: Ratio of the number of bases whose mass value is greater than 30 (error rate less than 0.1%) to the total number of bases.
We pooled all the clean reads from six libraries and de novo assembled them using Trinity. A total of 117,609 contigs, corresponding to unique genes of which length ranged from 201 bp to 9,868 bp. There were 27,273 (23.19%) contigs longer than 1,000 bp. The detailed length distribution of transcripts from the red tilapia was in Fig. 1.

GO annotation and analysis of DEGs
Gene ontology (GO) annotation was then performed based on the Nr annotation and 34,542 contigs (29.37%) were assigned to GO term (Fig. 2). As shown in Fig. 2, a total of 49 terms were assigned, including 23 (46.94%) biological process category, 13 (26.53%) cellular component terms and 13 (26.53%) molecular function terms. In the biological process category, cellular process was the most abundant term, followed by metabolic process biological regulation and pigmentation. For the cellular component category, cell and cell part were the predominant terms which were followed by organelle. For the molecular function category, binding was the main term, and it was followed by catalytic molecular transducer and transcription regulator (Fig. 2).

Functional annotation
All contigs were compared with four data bases including EMBL database, the UniProt-SwissProt database, KEGG database and GO database for functional annotation. There were 27997, 35370, 34542 assembled contigs that had signi cant hits against EMBL database, the UniProt-SwissProt database and GO database, respectively (Table 3, Fig. 2). To identify the biological pathway in the red tilapia transcriptome, all contigs were mapped to the KEGG database and were associated with 216 pathways (Supplementary Table 1). Identifying candidate genes involved in skin color and pigmentation To show the differences in the skin color of red tilapia, we made a comparative analysis of two skin transcriptomes. Based on the standard that |logFC| ≥ 1 and FDR ≤ 0.5, 3, 228 differentially expressed transcripts were identi ed between HLF7 and HLF3, which include 1,327 up-regulated transcripts and 1,901 down-regulated transcripts in HLF7. Candidate gene enrichment analysis identi ed 49 candidate genes involved in the processes of skin color and pigmentation from the two transcriptome data sets.
Interestingly, six genes were revealed to be speci cally expressed in HLF7, while 16 genes were only found expressed in HLF3 (Table 4).

Discussion
We employed the transcriptome analysis to investigate the different gene expression of n tissues of pearl white red tilapia. The results showed that the tyrosine protein kinase STYK1 HLF7, HSP 70, HSP 30, HSP 90 and Transcription factor Sp6 expression increased signi cantly, while MC1R, Transcription factor (MafB, Jun-D, AP − 1, E2F5, ETV6, Sp9, Sp7, E2F1, Sp4) expression signi cantly decreased in low temperature group, indicating that the color change of red tilapia at low temperature resulted from the coregulation of multiple genes.
Heat shock protein (HSP70) is one of the most conservative and important proteins in HSPP family which is synthesized by the reaction of organism to the physical, chemical and biological stress agents in the environment [21]. HSP70 gene can produce abundant heat stress proteins, remove the denatured proteins, protect the cells from the damage of denatured proteins, and reduce the impact of external stress environment on the body of sh when the cells are under stress [22]. Studies have found that the HSP70 mRNA expression in the liver was signi cantly increased after 12 hrs of low-temperature stress on the sh [23]. It has been demonstrated that both low temperature and high temperature caused the rapid expression of HSP70 gene in the body cells of platy bream and the synthesis of heat stress protein [23]. In our study, the expression of HSP70 in HLF7 cultured at low temperature was signi cantly higher than that in normal temperature group, suggesting that the low temperature environment outside led to the changes in the physiological function of red tilapia. Moreover, HSP70 in the body of red tilapia would be highly expressed to protect the normal physiological function of sh, thus enhancing the adaptability of red tilapia to the environment.
The biological basis of body color change in sh is the number, morphology and distribution of different pigment cells on skin and scales [24]. Compared with only a kind of pigment cell in the mammalian (melanocytes), bony shes were found to have 6 kinds of pigment cells, including the melanocytes, red Methods pigment cells, yellow pigment cells, rainbow, white and blue pigment cells, making the formation and changes of the sh body color more complicated and different red tilapia showing different body colors [25]. It has also been reported that MC1R gene, ASIP gene and TYRP1 gene are the major candidate genes of sh body color formation [26]. MC1R gene, as a key gene to produce melanin, expressed slightly in HLF7 of the low temperature group, and it was not the main reason for the formation of black spots [27]. It is speculated that the formation of melanin by MC1R gene was affected by the albino gene. Meanwhile, the MC1R gene, which was related to the synthesis of melanin, represents a critical role in the formation of sh body color. However, MC1R expression was low in HLF7 in black-spotted red tilapia of the lowtemperature group in our study, probably because it was an acute regulator and made pigments signi cantly deposited in the short term. And hence the surface of the skin was colored with streaks, and then black markings or red markings, but the expression level was stable.
It has been revealed that the body color of Mozambique tilapia became black in the dark after 25 hrs, but there was no signi cant change in MC1R expression in the epidermis and the α-MSH level in the serum. Therefore, the MC1R gene could also act as an acute regulator in the regulation of sh pigmentation, but the pigment is signi cantly deposited in a short period sine this activity expression cannot be maintained for a long time [28]. So even if the red tilapia epidermis exhibited a certain color, there was no signi cant change in MC1R mRNA expression in sh tissues.
Tyrosinase is the rate-limiting enzyme of melanin synthesis vital for the rate and quantity of eumelanin and pheomelanin [29]. Under the action of tyrosinase, tyrosine forms eumelanin through a series of oxidation and catalytic polymerization. It turns to the synthetic pathway of brown melanin when the activity of tyrosine is reduced [30]. Among them, eumelanin is responsible for producing black and brown phenotypes, and brown melanin for producing yellow and red phenotypes [31]. During the formation of melanin, hormones, inorganic ions which are in vivo and external light, UV light physicochemical factors may affect the activity of TYR. Besides, the of the eumelanin to pheomelanin by ASIP gene was linked to the expression of the tyrosine family gene (TYR, TYRP1, DCT) and reduced its expression [32]. Currently, the expression of gene related to the tyrosinase synthesis becomes an important indirect indicator of the energy of melanin formation in sh due to the absence of an effective method for determining the melanin content and body color index of sh [33].

Conclusions
In this study, the expressive quantity of tyrosine protein kinase in the low temperature group HLF7 was signi cantly higher than that in the normal temperature group HLF3, which indicated that the melanin synthesis ability was enhanced in the low temperature group HLF7. The result was consistent with the tendency of red tilapia to become darker as the temperature gradient decreased. The pathway of melanin synthesis is a complex process involving the participation and regulation of polygene. Nevertheless, the overall melanin content, the ratio of eumelanin to brown melanin and the molecular metabolic mechanism of red tilapia body color variation was still needed to be further studied.

Sample collection
Each 30 Guam red tilapia (GR), pearl white red tilapia (WR) and Florida red tilapia (FR), which were bred in the same year and weighed 10 g, 80 g and 120 g, respectively were selected for the experiment. The animal experiment is complying with the ARRIVE guidelines and carried out in accordance with the U.K. Animals (Scienti c Procedures) Act, 1986 and associated guidelines, EU Directive 2010/63/EU for animal experiments. Different red tilapia strains were placed in diverse aquarium for 5 days at 20 °C to adapt to the living environment of the aquarium for subsequent temperature culture test. Red tilapia strains were fasted for 24 hrs before the trial. The circulating water ltration system was adopted in this experiment, and feeding was conducted at 9:00 am and 4:00 pm, respectively per day. The laboratory temperature was arti cially cooled, from 20 °C to 6 °C whose gradient was 0.2 °C/2 h. During this period, water temperature was recorded every 2 hrs during the day and every 4 hrs at night. Meanwhile, data such as the deaths, the temperature at which the body color began to change, the location and degree of discoloration of red tilapia in each group were recorded. In the normal water temperature (20 °C) and cooling process (after 24 h of body color change), 3 n tissues of pearl white red tilapia were extracted for transcriptome sequencing respectively, stored in a -80 °C freezer before sequencing.

RNA isolation and cDNA library construction
Total RNA was extracted from tissue samples using MagZol Reagent. Then concentration and mass of each RNA sample was detected by nucleic acid protein detector (OD260: OD280) (Eppendorf, Germany) and agarose gel electrophoresis, respectively. The quali ed RNA samples were enriched with magnetic beads with Oligo (dT). cDNA strand was synthesized with six base primers using enriched mRNA as template, and cDNA library was obtained by a series of modi cations and PCR enrichment. The cDNA library was then quantitatively analyzed by Qbuit3.0, and the library insert size was detected by Agilent 2100 Bioanalyzer to ensure that the library effective concentration and quality meet the requirements of the RNA sequencing.
Transcriptome sequencing and De novo assembly of transcriptome Each quali ed cDNA library was sequenced using Illumina HiSeq2500 for 2 × 125 bp pair-end (PE) sequencing. Quality control of all raw reads was conducted by Fastqc (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) software. Raw Data obtained from sequencing were ltered, and connector sequences and low-quality Reads were removed to obtain high quality Clean Data. Gene assembly analysis was performed on high quality clean data after quality control. Trinity (v2.2.0) analysis software based on DE Bruijn graphs was used for subsequent analysis to assemble the transcriptome sequencing data from scratch. The sequences output from Trinity software can be used for further analysis.

Assembled sequence annotation Functional annotation
Functional annotation of the assembled sequences was performed by homology searches against the UniProt-SwissProt (The Universal Protein Resource) database, the GO (Gene Ontology) database, the KEGG (Kyoto Encyclopedia of Genes and Genomes) database and the TrEMBL protein database, and searches were conducted by the BLAST program.

RT-PCR analysis
Total RNA was isolated from samples by TRIzol method (Invitrogen, CA, USA), and cDNA was obtained by reverse transcription after DNA elimination. The results were calculated by ABI Prism 7500 RT-PCR system (Applied Biosystem, Foster City, CA, USA). Ampli cation was performed in a mixture (10 µl

Availability of data and materials
All data generated or analyzed during this study are included in this manuscript.
Ethics approval and consent to participate Page 13/17 Not applicable.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Figure 1
Transcriptome sequencing and assembly The detailed length distribution of transcripts from the red tilapia.