Comparative Transcriptome Analysis Reveals the Genetic Basis of Skin Color Variation in Common Carp

Background The common carp is an important aquaculture species that is widely distributed across the world. During the long history of carp domestication, numerous carp strains with diverse skin colors have been established. Skin color is used as a visual criterion to determine the market value of carp. However, the genetic basis of common carp skin color has not been extensively studied. Methodology/Principal Findings In this study, we performed Illumina sequencing on two common carp strains: the reddish Xingguo red carp and the brownish-black Yellow River carp. A total of 435,348,868 reads were generated, resulting in 198,781 assembled contigs that were used as reference sequences. Comparisons of skin transcriptome files revealed 2,012 unigenes with significantly different expression in the two common carp strains, including 874 genes that were up-regulated in Xingguo red carp and 1,138 genes that were up-regulated in Yellow River carp. The expression patterns of 20 randomly selected differentially expressed genes were validated using quantitative RT-PCR. Gene pathway analysis of the differentially expressed genes indicated that melanin biosynthesis, along with the Wnt and MAPK signaling pathways, is highly likely to affect the skin pigmentation process. Several key genes involved in the skin pigmentation process, including TYRP1, SILV, ASIP and xCT, showed significant differences in their expression patterns between the two strains. Conclusions In this study, we conducted a comparative transcriptome analysis of Xingguo red carp and Yellow River carp skins, and we detected key genes involved in the common carp skin pigmentation process. We propose that common carp skin pigmentation depends upon at least three pathways. Understanding fish skin color genetics will facilitate future molecular selection of the fish skin colors with high market values.


Introduction
Coloration is one of the most diverse phenotypic traits in vertebrates, and it exerts multiple adaptive functions, such as species identification, camouflage, warning or threatening of predators, photoprotection, thermoregulation and photoreception [1]. Skin coloration is the result of diverse pigments synthesized by pigment cells or chromatophores, and it is affected by multiple factors, including environmental, nutritional, physiological, or genetic conditions. Among these factors, the most fundamental and important is the genetic basis of skin pigmentation: which genes are likely to be involved and how. Cellular, genetic and genomic approaches have been widely adopted to answer this question.
In humans and mammals, the pigment melanin is the primary determinant of skin color. The availability of the complete human genome sequence and of adequate genomic resources, along with genome-wide association studies (GWAS), has provided insight into the genetic basis of the pigmentation process. In contrast to mammals, which possess only one type of pigment cell (the melanocyte), and amphibians and reptiles, which possess xanthophores, erythrophores and reflecting iridophores, teleost fish possess up to six different pigment cells, including melanophores, xanthophores, erythrophores, iridophores, leucophores and cyanophores [2]. A diversity of pigment cells, associated with a series of cellular, physiological, genetic and environmental factors, makes fish skin pigmentation a complicated biological process. Extensive studies have been conducted on model fish species such as zebrafish and medaka in an effort to unravel the genetic basis of fish skin pigmentation, and dozens of genes have been reported to be involved in the pigmentation process, such as matp, oca4, sox10, kit, ednrb, slc24a5 and many others, through collecting and identifying the pigmentation mutations [3][4][5]. However, few genetic skin color studies have been conducted on non-model fish species.
The common carp (Cyprinus carpio), a freshwater fish that is especially widespread in Europe and Asia, was domesticated more than 2,000 years ago. Over its long history of domestication, the common carp has been introduced into various environments worldwide, resulting in hundreds of strains or varieties that display rich biodiversity, genetic polymorphisms and diverse skin colors. Skin color is an important economic trait for common carp, as it acts as an important criterion for visually determining quality and market value. Several studies have provided initial insights into the genetic basis of skin coloration in common carp. For instance, David et al. reported that the MC1R gene was associated with the development of black pigmentation in the ornamental Koi common carp [6]. A very recent study reported that 80% of identified zebrafish pigmentation genes were also present in the Oujiang color common carp strain [7]. In addition to the direct effects of genes, other regulators, such as microRNAs, play crucial roles in common carp skin coloration by regulating the expression of downstream pigmentation genes [8]. However, some details of the genetic mechanisms underlying common carp skin pigmentation are not well understood, such as the interactions among pigmentation genes and the genetic regulation underlying synthesis of different pigments.
The Xingguo red carp (Cyprinus carpio var. xingguonensis, XGC) and the Yellow River carp (Cyprinus carpio haematopterus Temminck et Schlegel, YRC) are two traditional domestic strains of common carp in China. XGC, a regional strain from Xingguo County in the Jiangxi Province in Southern China, is known for its red color and has approximately 1,300 years of cultural history [9]. YRC, a brownish-black common carp, was originally cultivated along the Yellow River basin thousands of years ago and is now a prominent strain in Northern China. It possesses advantageous traits such as strong cold tolerance, high efficiency of food conversion and good meat quality. The distinct skin colors of these two common carp strains make them suitable models for exploring the genetic basis of skin pigmentation. To better understand skin color genetics, we utilized the powerful approach of comparative transcriptome analysis using next-generation sequencing and examined transcript profiles from the skins of the XGC and YRC common carp strains. We obtained candidate genes that may be involved in the skin pigmentation process and identified gene pathways that may regulate the synthesis of different pigments. Understanding the molecular mechanisms of skin pigmentation in common carp will advance our knowledge of skin color genetics in fish and accelerate the molecular selection of fish species with consumer-favored skin colors.

Results and Discussion
Sequencing of short expressed reads from XGC and YRC Diverse skin colors make fish good genetic models for understanding the skin pigmentation process. Various fish colorations are determined by the density and position of different pigment cells, which is believed to be primarily under genetic control. To better understand fish skin color genetics, we conducted a comparative transcriptome analysis between two common carp strains, XGC and YRC, using next-generation sequencing. First, we generated reference sequences for the subsequent analysis of skin genes that are differentially expressed between the two carp strains. Six tissues, including brain, blood, gill, head kidney, muscle and skin were collected and deepsequenced using Illumina HiSeq 2000. As shown in Table 1, a total of 435 million paired-end reads were generated, of which 211 million were from YRC and 224 million were from XGC. The number of reads generated from each tissue ranged from 27.5 million to 36.9 million, with 2 outliers of 50.5 million and 46.3 million from XGC muscle and blood, respectively. After the removal of ambiguous nucleotides, low-quality sequences (Q,20) and short reads (length,15 bp), a total of 422 million clean reads (97%) were selected for further analysis.

Reference sequence assembly and annotation
All of the clean reads were pooled and de novo assembled to generate reference sequences using the Trinity assembler [10]. After the removal of sequence redundancy using CD-hit software [11], a total of 198,781 contigs, with a minimum length of 200 bp, a maximum length of 26,217 bp and an N50 of 1,970 bp, were generated as the reference sequences for subsequent analysis ( Table 2). There were 43,310 contigs longer than 1,000 bp. To assess the quality of the sequencing and de novo assembly, all of the clean reads were mapped back to the assembled sequences. As shown in Table 1, the mapping ratio ranged from 77% to 90% with an average of 83%, indicating a high-quality sequence.
Gene prediction was performed on the assembled contigs using BlastX searches against three protein databases, including the NCBI non-redundant (nr) database, the UniProt database and the Ensembl zebrafish protein database, with an E-value cutoff of 1e 210 . There were 62,343, 57,083, and 69,634 assembled contigs with significant hits against nr, UniProt and zebrafish, respectively (Table 2). Cumulatively, 73,928 assembled contigs had at least one significant hit against at least one of the three databases, allowing for the prediction of 20,028 unique genes.
Gene ontology (GO) annotation was then performed with these 73,928 unique gene-containing contigs using Blast2GO [12]. Of these, 58,841 contigs, corresponding to 16,849 unique genes, were assigned to at least one GO term ( Table 2). As shown in Figure 1, a total of 38 GO terms were assigned, including 10 (26.3%) cellular component terms, 11 (28.9%) molecular function terms and 17 (60.7%) biological process terms. From the GO category of molecular function, binding was the most predominant term, accounting for 63.4% of the sequences annotated in that term, and it was followed by catalytic activity. In the biological processes category, cellular process was the most predominant term (20.2% of the sequences), and it was followed by single-organism process (16.6%) and metabolic process (14.9%).

Identification of differentially expressed genes in XGC skin compared with YRC skin
Based on the criteria that |fold-change| §2 and p-value!0.05, a total of 4,367 of the 198,781 assembled contigs showed significantly different expression in XGC skin compared with YRC skin. Differentially expressed contigs represented 2,012 unique genes, of which 874 genes were up-regulated in XGC and 1,138 genes were down-regulated ( Figure 2, Table S1). To validate the differentially expressed genes identified by comparative transcriptome analysis, we randomly selected 20 representative genes for qRT-PCR confirmation of differential expression. The gene expression patterns in the result of qRT-PCR were compared with the data obtained from the comparative transcriptome analysis. As shown in Figure 3, the qRT-PCR expression patterns of 19 out of the 20 randomly selected differentially expressed genes were in agreement with the results of the comparative transcriptome analysis with only slight differences in expression levels, indicating that there was no consistent bias in the expression patterns (i.e., in the direction of the differential expression or in the degree of fold change) for either method. Melting-curve analysis showed that a single product was amplified for all tested genes, indicating that the reference assembly was largely accurate and that it did not contain a large number of chimeric transcripts.
After carefully investigating the differentially expressed gene list from the skin tissues of XGC and YRC, we identified many putative pigmentation genes that may be associated with melanin synthesis, such as those encoding tyrosinase-related protein 1 (TYRP1), premelanosome protein (SILV/PMEL), cysteine/glutamate transporter (xCT/SLC7a11) and agouti signaling protein (ASIP). Melanin is the major pigment present in vertebrates [13]. There are two types of melanin, eumelanin and pheomelanin. Eumelanin corresponds to a brown/black color, while pheomelanin corresponds to a red/yellow color. As their names suggest, the skin color of XGC is red, while YRC is brownish-black, indicating that the amount and density of eumelanin are much lower in XGC skin than in YRC skin. Conversely, we hypothesized that the amount and density of pheomelanin might be much higher in XGC skin than in YRC skin. The observed differential expression patterns of melanin synthesis pathway genes endorsed this hypothesis. As shown in Table 3, the expression level of TYRP1 was 5.7-fold up-regulated in YRC skin compared with XGC, suggesting that TYRP1 is one key gene that contributes to brownish-black coloring in common carp by acting in the eumelanin synthesis pathway. TYRP1 is expressed specifically in melanocytes, and it plays a crucial role in pigmentation and color patterning by affecting melanin synthesis, stabilizing tyrosinase protein, modulating tyrosinase catalytic activity, maintaining melanosome structure, and affecting melanocyte proliferation and cell death [14][15][16]. TYRP1 was the first cloned pigmentation gene, and it was later mapped to the brown locus in mouse [17]. Mutations in TYRP1 result in failure to form a black coat color, as observed in mouse [18], dog [19], cat [20], cow [21], and sheep [22]. In human cells, TYRP1 is only detectable in cells containing eumelanin [23], and it is involved in the pigmentation differences among human populations [24]. These studies suggest essential roles for the TYRP1 gene in eumelanogenesis. Similar roles for TYRP1 have been reported in teleost fish. Knockdown of the TYRP1 gene in zebrafish and medaka revealed that black eumelanin formation essentially relies on the presence of TYRP1 [25]. Our results confirmed a role for TYRP1 in the skin coloration of the common carp. Another significant eumelanin synthesis pathway gene we identified was SILV/PMEL. SILV encodes premelanosome protein, which is extensively expressed in pigment cells [26] and can catalyze the conversion of indole-5,6-quinone carboxylic acid into eumelanin [27]. Mutations in SILV causing pigmentation phenotypes have been reported in a number of vertebrate species. Inactivating SILV in mice led to a substantial reduction in eumelanin, indicating that SILV plays a critical role in maintaining efficient epidermal pigmentation [26]. In mice, SILV is required for normal melanosome development in skin melanocytes, choroid melanocytes and retinal pigment epithelium cells [26]. In the domesticated chicken, SILV mutations inhibited the production of all eumelanin in plumage and skin [28]. Similarly, in horses, SILV mutations caused a dilution of black pigment in the  [29]. Consistent with all these studies, our results showed that the expression level of SILV in YRC was significantly higher than in XGC, exhibiting a 5.2-fold up-regulation (Table 3). This result further confirmed the role of SILV in common carp eumelanogenesis.
Several other genes in our differentially expressed gene list are involved in the production of pheomelanin, the most significant of which are ASIP and xCT. In our results (Table 3), both ASIP and xCT were significantly up-regulated in XGC (by 3.2-fold and 3.5fold, respectively), suggesting a significant effect on the production of pheomelanin in common carp. The functions of ASIP and xCT in pheomelanin synthesis have been widely studied in other species. ASIP can regulate pigmentation by antagonizing the binding of a-MSH to Mc1r, thus switching melanin synthesis from eumelanin to pheomelanin [30]. The ASIP protein also downregulates genes such as TYRP1, DCT, and TYR to inhibit the synthesis of eumelanin [31], which is in agreement with our result that ASIP gene expression is up-regulated while TYRP1 gene expression is down-regulated in XGC skin. The important role of ASIP in determining skin or coat color has been demonstrated in human [32], mouse [33], fox [34], rat [35], pig [36] and other species. xCT (SLC7a11), a cysteine/glutamate exchanger, mediates the cellular uptake of cysteine [37]. Cysteine is an important component of pheomelanin; therefore, xCT directly affects the pheomelanin synthesis pathway [38]. Loss of xCT expression caused marked inhibition of pheomelanogenesis [38]. Conversely, high-level expression of xCT contributed to the production of pheomelanin, as indicated by our results.

Enrichment and pathway analysis
We attempted to categorize the 2,012 differentially expressed genes based on their likely functions by using gene annotation, GO term enrichment analysis, and KEGG pathway analysis. Using the Blast2GO software, all the differentially expressed genes were  classified into different cellular, biological and functional gene ontologies. Together with the GO annotation of the reference sequences, the web-based program WEGO [39] revealed 17 significantly over-represented GO terms at the 2 nd GO level that were enriched in XGC skin relative to YRC skin. The top 10 of these are shown in Table 4, and they include pigmentation (GO:0043473), binding (GO:0005488), enzyme regulator activity (GO:0030234), biological regulation (GO:0065007), multicellular organism process (GO:0032501), multi-organism process (GO: 0051704), developmental process (GO:0032502), cellular process (GO:0009987), localization (GO:0051179), and extracellular region part (GO:0044421). The primary interest of this study was to detect expression signatures indicative of fish skin pigmentation; therefore, the 363 genes in the pigmentation GO category were considered to be most informative for further pathway analysis.
As described in a pathway analysis of non-model fish species [40], downstream pathway analysis was conducted using KEGG pathway analysis combined with manual literature searches. We mainly focused on pigmentation-related pathways to reveal the mechanism of skin color variation in the common carp. We included 1) melanin biogenesis; 2) the Wnt signaling pathway; and 3) the MAPK signaling pathway. Table 3 lists the key differentially expressed genes involved in each pathway. Melanin biogenesis is under complex controls by multiple agents via pathways. Slominski et al. [41] elucidated that myriad hormonal factors are involved in regulating melanogenesis, such as melanocortins (MSH), ß-endorphin, endothelins, histamines, eicosanoids, catecholamines, c-kit ligand estrogens, androgens, vitamin D, serotonin, melatonin, dopamine, acetylcholine, agouti proteins and their receptors. Among those hormonal regulators, MSH with its receptor Mc1r is the most important positive regulators, whereas agouti proteins stand out among the negative regulators [41]. Besides hormonal regulators, nutritional regulators such as aromatic amino acid L-tyrosine and L-DOPA play critical roles in melanogenesis as well [42]. The molecular and genetic basis of pigmentation including the structure of tyrosinase, the development and differentiation of pigment cell, is largely conserved between mammals and teleost [41,43]. The putative gene pathways involved in the common carp skin pigmentation process are shown in Figure 4, indicating the critical roles of the hormonal and nutritional regulators in melanogenetic pathways in common carp ( Figure 4). Briefly, MSH binds to the G-proteincoupled receptor Mc1r, resulting in up-regulated cAMP levels, which in turn trigger the eumelanin biosynthesis process. The synthesis of eumelanin is then catalyzed successively by Tyr, DCT, TYRP1 and SILV [44] (Figure 4). The switch from eumelanogenesis to pheomelanogenesis during melanin biosynthesis mainly depends on the presence of ASIP, an inverse agonist of Mc1r, which consequently make melanocytes start producing pheomelanin and stop producing eumelanin.
Signaling pathway such as Wnt and MAPK is implicated in numerous development and physiological process. Several studies have been reported that both Wnt and MAPK signaling pathways take part in the development of pigment cells [45][46][47]. Intriguingly, our results indicated that both pathways are highly likely to be involved in melanin biosynthesis in common carp. A group of identified differentially expressed genes involved in the Wnt or MAPK signaling pathway have the potential affecting the melanin biosynthesis (Table 3). For instance, WNT5B, a member of the Wnt protein family, is expressed in pigment cells located at the extended posterior edges of the caudal fin in zebrafish [48]. In our results, WNT5B was 2.5-fold up-regulated in YRC, suggesting that WNT5B may be involved in the skin pigmentation process. SFRP5 (secreted frizzled-related protein 5), a member of the SFRP family that modulates Wnt signaling transduction, is highly expressed in the retinal pigment epithelium [49]. This finding, combined with our result that the expression level of SFRP5 is higher in the skin of brownish-black YRC than in red XGC (Figure 3), raises the likelihood that SFRP5 plays an important role in the Wnt signaling pathways involved in eumelanin synthesis.
The melanogenetic pathway in the skin could interact with other regulatory pathways such as nervous system, immune system, or circulatory system, via various organic compounds [50]. For instance, catecholamine, a compound produced in the skin, can be oxidized to form melanin, and also can affects the skin immune functions and systemic functions through communicating with specific cell surface receptors [50]. However, the details of interactions of regulatory pathways in common carp skin remain to be further investigated.

Ethics statement
This study was approved by the Animal Care and Use committee of the Centre for Applied Aquatic Genomics at the Chinese Academy of Fishery Sciences.

Fish sampling, RNA isolation and sequencing
As described in our previous study [51], the two strains of common carp used in this study were sampled from distinct breeding populations. XGC individuals were collected from the National Fish Hatchery of Xingguo red carp, Xingguo, Jiangxi province, China. YRC individuals were collected from the Henan Academy of Fishery Sciences, Henan province, China. Six tissues from 18 individuals of each strain were collected; the tissues included brain, blood, gill, head kidney, skin and muscle. All tissues were placed in 2 ml of RNAlater (Qiagen, Hilden, Germany), and kept at 220uC until RNA extraction. RNA isolation was performed using TRIzol (Invitrogen, Carlsbad, CA, USA) with DNase I in accordance with the manufacturer's instructions. RNA quality was verified using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Equal amounts of high-quality RNA from each tissue were pooled and sent to HudsonAlpha Genomic Services Laboratory (Huntsville, AL, USA) for sequencing on an Illumina HiSeq2000 platform. All De novo assembly of reference sequences All raw sequencing reads were trimmed by removing adaptor sequences, ambiguous nucleotides, low quality sequences (Q,20) and short reads (length below 30 bp) using CLC Genomics Workbench (CLC bio, Aarhus, Denmark). The Trinity software was then used to assemble all of the cleaned reads using default parameters. The assembled sequences were then filtered using the CD-hit program to reduce redundancy. The resulting contigs that were larger than 200 bp were considered as the final nonredundant transcripts and were used as the reference sequences.

Transcriptome annotation and ontology
The assembled reference sequences were used as query sequences to search against the NCBI nr database, the UniProt database and the Ensembl zebrafish protein database. Searches were conducted using the BlastX program with an E-value cutoff of 1e 210 . Gene ontology was performed by importing the zebrafish Blast results into Blast2GO software. GO terms were then assigned to each sequence automatically. The annotation output was categorized by cellular component, molecular function and biological process.

Differential gene expression analysis
Clean reads from the skin tissues of each common carp strain were first aligned to their respective reference sequences using Bowtie 2 [52]. As described previously [53], RSEM was then used to calculate and estimate gene or isoform abundances. The expression level of each transcript in each sample was then normalized using edgeR [54]. Transcripts with fold change values larger than 2 and p values lower than 0.05 were included in subsequent analyses as the differentially expressed genes.

qRT-PCR analysis
A total of 20 genes with significantly different expression in the two carp strains were selected for validation using qRT-PCR analysis. The house-keeping gene b-actin was used as an internal reference. All primer sequences are listed in Table S2. The firststrand cDNA was synthesized using the SuperScript III RT kit (Invitrogen) according to the manufacturer's instructions. All cDNA samples were diluted to 100 ng/ml before use in qRT-PCR reactions on an ABI PRISM 7500 Real-Time Detection System (Life Technologies). The amplifications were performed in a total volume of 15 ml and included 7.5 ml of 2X SYBR Green Master Mix reagent, 1 ml of cDNA and 0.3 ml of each primer (10 mM). The thermal cycling profile consisted of an initial denaturation at 95uC for 10 min followed by 40 cycles of denaturation at 95uC for 15 s and annealing/extension at 60uC for 1 min. An additional temperature-ramping step from 95uC to 65uC was used to produce the melting curve. All reactions were conducted in triplicate and included negative controls with no template. Expression differences between XGC and YRC skin were assessed for statistical significance using a randomization test in REST software. The expression levels of genes were normalized to the levels of b-actin in the same sample.

GO enrichment analysis
The web-based program WEGO was used with default parameters for statistical analysis of GO term overrepresentation among the genes with differential expression in the skin of the two carp strains. The study set represented the frequency of GO terms in the differentially expressed genes, while the population set corresponded to the whole reference gene set.

Conclusions
To better understand the genetic mechanisms of skin color variation in common carp, we conducted a comparative analysis of the skin transcriptomes from two common carp strains with distinct skin colors, the red XGC and the brownish-black YRC. We detected 2,012 unique genes that were differentially expressed in the two strains. Further annotation, GO enrichment and pathway analysis indicated that key pigmentation-related genes were involved in at least three pathways that regulate pigment synthesis in the two strains. These results provide us with a valuable basis for understanding of the molecular mechanisms of skin pigmentation in teleost fish. Such understanding will facilitate the genetic selection and breeding of common carp with marketfavored colors.

Supporting Information
Table S1 Detailed information of the differentially expressed genes. The fold change cutoff is 2, and p-value cutoff is 0.05. (XLSX)