Skip to main content
  • Research article
  • Open access
  • Published:

Comparative Transcriptome and DNA methylation analyses of the molecular mechanisms underlying skin color variations in Crucian carp (Carassius carassius L.)

Abstract

Background

Crucian carp is a popular ornamental strain in Asia with variants in body color. To further explore the genetic mechanisms underlying gray and red body color formation in crucian carp, the skin transcriptomes and partial DNA methylation sites were obtained from red crucian carp (RCC) and white crucian carp (WCC). Here, we show significant differences in mRNA expression and DNA methylation sites between skin tissues of RCC and WCC.

Results

Totals of 3434 and 3683 unigenes had significantly lower and higher expression in WCC, respectively, compared with unigenes expressed in RCC. Some potential genes for body color development were further identified by quantitative polymerase chain reaction, such as mitfa, tyr, tyrp1, and dct, which were down-regulated, and foxd3, hpda, ptps, and gch1, which were up-regulated. A KEGG pathway analysis indicated that the differentially expressed genes were mainly related to mitogen activated protein kinase (MAPK), Wnt, cell cycle, and endocytosis signaling pathways, as well as variations in melanogenesis in crucian carp. In addition, some differentially expressed DNA methylation site genes were related to pigmentation, including mitfa, tyr, dct, foxd3, and hpda. The differentially expressed DNA methylation sites were mainly involved in signaling pathways, including MAPK, cAMP, endocytosis, melanogenesis, and Hippo.

Conclusions

Our study provides the results of comparative transcriptome and DNA methylation analyses between RCC and WCC skin tissues and reveals that the molecular mechanism of body color variation in crucian carp is strongly related to disruptions in gene expression and DNA methylation during pigmentation.

Background

Fish display numerous colors and color combinations. Some species have been used as ornamental fish, in which body color and pattern are important characteristics. The crucian carp (Carassius auratus L.) is a freshwater fish in the family Cyprinidae (order Cypriniformes) and is one of the most commonly kept aquarium fish [1]. As a colored fish with variants in body color, the crucian carp is a popular ornamental strain, as they are produced in a variety of colors including red, orange/gold, yellow, white, and brown, or black. Red crucian carp (RCC; C. auratus red var.), or the common goldfish, is one of the earliest domesticated fish and has become a popular ornamental fish owing to its red/orange body color. White crucian carp (WCC; C. auratus cuvieri Temminck et Schlegel) is back-gray in body color [2, 3].

Pigment cells on the body surface provide body color pigment patterns and are usually classified according to pigment composition. Fish body color is controlled by the distribution of pigment cells, such as melanocytes, xanthophores, erythrophores, and iridocytes [4,5,6]. Melanocytes, contain a large number of melanin granules and are able to absorb the specific wavelengths of incident light that make fish appear grey-black. Xanthophores and erythrophores hold carotenoids and pteridines and provide the yellow, orange and red colors of fish.

The genetic control rules for body color formation in RCC may be differ from those in WCC. Actually, melanocytes appear during the RCC embryo and larval stages. Then, the body color of RCC changes from dark-gray to red-orange as the melanocytes disappear and body color begins to form [7]. It has been suggested that body color formation in RCC is subject to complex controls by multiple agents not only via the pigment biosynthetic pathway and chromatophore differentiation pathway but also the autophagy and apoptosis pathways [8]. DNA methylation is an important epigenetic modification of the eukaryotic genome and plays an important role maintaining the biological functions of higher organisms, such as normal cellular functions, genetic markers, embryonic development, aging, and human tumorigenesis [9, 10]. Li (2015) performed the DNA methylation analyses for red skin and white skin of koi carp, there showed that the DNA methylation levels of two selected DEGs inversely correlated with gene expression, indicating the participation of DNA methylation in the coloration [11]. In this study, we chose the methylation-RAD detection method to further explore the mechanisms underlying gray and red body color formation in crucian carp. We compared the skin tissues of RCC and WCC using transcriptome sequencing and methylation-RAD sequencing analyses (simple genomic methylation site detection method) [12, 13]. We analyzed the differentially expressed genes (DEGs) and different DNA methylation levels using the transcriptome and MethylRAD data. This study attempts to reveal the molecular and genetic mechanisms as well as epigenetic modification of coloration formation in two kinds of crucian carp.

Results

Transcriptome assembly and annotation

After filtering low quality and short sequences, we obtained about 54 and 52 million clean reads in the skin tissues of RCC and WCC, respectively. The complete clean reads for these libraries have been uploaded onto the NCBI Sequence Read Archive site SRS2441209, SRX3105778, SRR5947250, and SRS2441217. After eliminating redundant sequences and filtering short sequence (bp ≤ 400), 56,564 and 56,612 assembled transcripts were obtained from the skin tissues of RCC and WCC, respectively. Average read size, Q20 percentage, and other Parameters are presented in Table 1.

Table 1 Overview of sequencing and assembly

Functional enrichment analysis of annotated Unigenes

Twenty four thousand nine hundred seventy one unigenes were obtained from the crucian carp skin, which annotated with Gene Ontology (GO). The GO term classifications were similar between RCC and WCC, with little variation in molecular function (Fig. 1a). Using COG annotation, 16,804 unigenes in the crucian carp group were assigned into 25 functional categories. A total of 34,193 annotated sequences were enriched into 258 predicted KEGG pathways in the transcriptome from the skin of crucian carp. Detailed annotated information is given in Additional file 1: Table S1.

Fig. 1
figure 1

Summary of GO and KEGG pathway enrichment analysis. a Gene ontology classification of unigenes between the RCC and the WCC. “Red” represents the unigenes that are expressed in the RCC; “White” represents the unigenes that have expression inthe WCC. b Differentially gene expression in the skin of the RCC and the WCC. c Gene ontology classification of differentially expressed unigenes between the RCC and the WCC. “Up” represents the unigenes that are over-expressed in the RCC; “Down” represents the unigenes that have reduced expression in the WCC. d KEGG enrichment top30 bubble chart. The x-axis Enrichment Score is an enrichment score. The more the number of different genes is the larger bubble. The bubble color is changed by red-green, and its enrichment p-value is large

Comparative and enrichment analysis of differentially expressed genes

Using the MARS method [14], we defined an expression level of >2-fold changes and false discovery rates (FDR) of ≤0.001 to filter genes with significant differential expression between the skin tissues of RCC and WCC. A total of 7117 genes were differentially expressed between RCC and WCC. Overall, 3683 were up-regulated in RCC and 3434 were down-regulated compared with their levels in WCC (Fig. 1b). The DEGs of RCC compared to those of WCC were mainly in the cellular component, molecular function, and biological process categories (Fig. 1c). And the analysis of KEGG pathway revealed that multiple pathways, including the “MAPK signaling pathway”, “Wnt signaling pathway”, “cell cycle”, “VEGF signaling pathway”, “endocytosis”, “melanogenesis” and “tyrosine metabolism” were clearly enriched in the crucian carp DEGs (Fig. 1d).

Analysis of some pigmentation related expressed genes and RNA-Seq data validation

To examine the reliability of the RNA-seq results, eight DEGs (tyr, mitfa, tyrp1, dct, foxd3, hpda, gch1, and ptps) involved in the development of body color were selected for validation using qPCR. As shown in Fig. 2, the qRT-PCR expression patterns of the eight DEGs were in agreement with the RNA-seq data. A number of DEGs related to the occurrence of pigment cells were selected.

Fig. 2
figure 2

Comparison of relative fold change of Illumina sequencing with qRT-PCR results between the RCC skin and the WCC skin

As shown in Table 2 and Fig. 3, the mRNA expression levels of key genes for tyrosine synthesis, such as those encoding MITFa, TYR, tyrosinase-related protein 1 (TYRP1) and DCT, were all down-regulated in RCC skin compared to that in WCC skin. Furthermore, the mRNA expression levels of foxd3, hpda, ptps, and gch1 were up-regulated in skin tissues from RCC compared to those in WCC.

Table 2 Analysis of differentially expressed genes related to pigmentation
Fig. 3
figure 3

Real-time PCR analysis for a part of the differentially expressed genes. The lowercase letters show the significant differences (*P < 0.05) (means ± SD of relative expression; n = 3 for each group)

Analysis of DNA methylation site distribution and DNA methylation levels

The MethylRAD analysis was used to characterize the cytosine methylation pattern in the skin tissues of RCC and WCC. The high-quality methylation tag libraries of the two samples were 18.08–21.45% comparable to unique positions in the genome (Additional file 1: Table S2). The sequencing depths of the DNA methylation sites (CCGG sites and CCWGG sites) in each sample are shown in Table 3. Totals of 96,122 and 7961, CCGG and CCWGG DNA methylation sites were found in RCC, with average methylation coverage of 30.71 and 9.18, respectively. Totals of 88,536 and 9534 CCGG and CCWGG DNA methylation sites were found in WCC with average methylation coverage of 28.35 and 9.02, respectively. These results show that the DNA methylation level of RCC was slightly higher than that of WCC (Table 3). The RCC and WCC differential methylation sites (CCGG/CCWGG) were identified on the chromosomes of common carp (Fig. 4a). We obtained the DNA methylation site annotation of the common carp genome, and found that the DNA methylation site distribution curve had TSS representing an upstream sequence centered on the transcription initiation site, and TTS representing a downstream sequence centered on the transcription termination site (Fig. 4b). We also discovered that these DNA methylation sites were mainly distributed in the gene body, intron, and intergenic areas, and that there were slightly more DNA methylation sites in RCC than that in WCC (Fig. 4c).

Table 3 DNA methylation site coverage depth in two coloration crucian carp skins
Fig. 4
figure 4

DNA methylation site on gene DNA methylation levels and function components distribution. a The distribution of differential DNA methylation sites on chromosomes. b DNA methylation level change line chart. The red curve represents the methylation level of the RCC, and the green curve represents the methylation level of the WCC. c DNA methylation site on different gene function components distribution histogram. 1stexon: first exon; Body: gene area; TSS200: gene starting sites within 200 bp upstream area; TSS1500: gene starting sites within 1500 bp upstream area; Exon: intron area; Intergenic: intergenic region. d KEGG enrichment analysis DNA methylation gene on top30 bubble chart. The x-axis Enrichment Score is an enrichment score. The more the number of different genes is the larger bubble. The bubble color is changed by red-blue-green-yellow, and its enrichment p-value is large

Enrichment analysis of differential methylation site-related genes

The function of a gene was described by a GO function enrichment analysis of the gene where the differential DNA methylation sites were located (Note Download link: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/951/615/GCF_000951615.1_common_carp_genome/GCF_000951615.1_common_carp_genome_genomic.gff.gz). The GO enrichment analysis top30 bar graph is shown in Additional file 1: Fig. S1. The “MAPK signaling pathway”, “cAMP signaling pathway”, “Hippo signaling pathway-fly”, “Endocytosis”, “Melanogenesis” and “Hippo signaling pathway” were significantly enriched in crucian carp. Some of these pathways are associated with pigmentation synthesis (Fig. 4d).

Partial differential expression DNA methylation site gene analysis

Genes with the same trends of change in expression levels were determined between the RNA-seq and MethylRAD data (Table 4). We obtained the DNA methylation site distribution of the TYR gene that is closely related to melanin synthesis (Fig. 5a). Figure 5b shows a comparison of the DNA methylation site coverage level for each tag in the TYR gene between RCC and WCC.

Table 4 Analysis of differentially expressed genes and their DNA methylation levels
Fig. 5
figure 5

The tyrosine genomic DNA methylation site analysis. a The gene structure of TYR and its associated MethylRAD tags. The differentially DNA methylated sites are indicated by red tag names. b Comparison of DNA methylation levels between two groups for each tag in the tyr gene

Protein expression of three major pigment genes in skin tissues of Crucian carp

A western blotting analysis with commercially available antibodies was performed to determine whether the DEG findings were reflected at the protein level. The relative expression levels of MITFa, TYR, and DCT were higher in WCC skin tissues than those of RCC (Fig. 6), which was consistent with the results of the transcriptome mRNA differential expression analysis, and partially supports the reliability of mRNA differential expression analysis.

Fig. 6
figure 6

Verification of the expression of MITFa, TYR and DCT protein by western blot analysis. The “RCC-1”, “RCC-2,” and “RCC-3” represent the bands from the RCC; “WCC-1”, “WCC-2,” and “WCC-3” represent the bands from the WCC

Discussion

DNA methylation is essential for normal development and is associated with a number of key processes, including genomic imprinting, X-chromosome inactivation, repression of transposable elements, aging, and carcinogenesis [15, 16]. In this study, we constructed RNA-seq and DNA methylation libraries from skin tissues of RCC and WCC using transcriptome sequencing and the MethylRAD methods and discovered some DEGs and differentially methylated genes associated with pigmentation. There were 7117 DEGs found between RCC and WCC. A GO enrichment analysis of the DEGs revealed that the variation in pigmentation was related to cellular components and biological processes, such as pigmentation, biological regulation, and cellular process terms. Meanwhile, we found that some differentially methylated genes are enriched in some signaling pathways, including MAPK signaling pathway, cAMP signaling pathway, Wnt signaling pathway and Melanogenesis. This suggests that DNA methylation may be related to the formation of body color in crucian carp.

Melanocytes contain melanin, which is mainly synthesized from tyrosine. The three important enzymes were involved in melanin synthesis, including tyrosinases (TYR), tyrosinase-related protein 1 (TYRP1) and dopachrome tautomerase (DCT or TYRP2). Tyrosinase (TYR) carries out tyrosine hydroxylation to L-DOPA, which is rapidly oxidized to DOPAquinone. The DOPAquinone spontaneously undergoes cyclization to DOPAchrome. Under the action of the dopachrome tautomerase (DCT or TYRP2) and tyrosinase-related protein 1 (TYRP1), the DOPAchrome form DHICA-melanin is insoluble, which rapidly oxidizes and polymerizes to form dark or brown-black pigments. Our results show that the mRNA and protein expression levels of tyr and dct in RCC were significantly lower than those in WCC, and that the tyr and dct gene have modified DNA methylation sites, which were negatively correlated with their mRNA expression levels. This suggests that the tyrosine metabolism pathway was associated with different skin color in red crucian carp and white crucian carp. Besides, we found that the mitfa mRNA and protein level in skin tissue of RCC was significantly lower than that of WCC, and DNA methylation level was negatively correlated with gene expression level. The mitfa gene is a key factor regulating survival and development of melanocytes. Mitfa directly regulates the expression of multiple genes necessary for melanophore development, including dct, tyr, tyrp1, c-kit, and bcl2 [17]. Ectopic mis-expression of mitf is sufficient to confer a melanoblast phenotype [18, 19]. Some transcription factors can act on the promoter region of MITF, thereby promoting the expression of MITF, such as beta-catenin, Paired box 3 (PAX3), SRY-box containing gene10 (SOX10) [20]. In zebrafish, null mutations in the microphthalmia-associated transcription factor (MITF) ortholog mitfa lack all neural crest-derived melanocytes but retain the other two neural crest pigment cell types, xanthophores and iridophores [21]. This suggests that mitfa may be play important roles information of body color in crucian carp.

We also discovered some differentially DNA methylated genes involved in the MAPK and cAMP signaling pathways, which might be related to pigmentation in crucian carp. Extracellular signal regulated kinase 2 (ERK2)/MAPKs are reported to play an important role during post-translational modification of mitf [22, 23]. Activation of c-KIT in melanocytes results in phosphorylation of MITF at Ser73 by ERK2. Phosphorylation at Ser73 increases recruitment of transcriptional coactivator p300 [24], while simultaneously targeting MITF for ubiquitin dependent proteolysis [25]. Our results indicate that mRNA expression of the erk2 gene was lower in RCC than that in WCC, leading to the inability of mitf to be phosphorylated and move to the nucleus, which may have decreased mitf mRNA level. This could also affect the expression of some phosphorylation genes, such as cell cycle and apoptosis-related genes. In the cAMP signaling pathways, Alpha-melanocyte stimulating hormone (a-MSH) is one such peptide hormone that has been suggested to modulate the pigment response. Studies have shown a role of a-MSH in promoting pigment synthesis a-MSH binds to the melanocortin 1 receptor (MC1r), which is a member of the seven-transmembrane receptor family that activates adenylcyclase via G-protein signaling [26]. Although cAMP triggers numerous downstream effects, one important target is the MITF gene, which is transcriptionally up-regulated by cAMP signaling in a melanocyte-restricted fashion, thus linking extracellular signals to MITF expression and the transcriptional regulation of pigmentation [27]. Meanwhile, cAMP may be a second messenger in the dispersion of pigment induced by prolactin and that a novel protein receptor coupled with a Gs protein may be present in the membrane of erythrophores and xanthophores of teleost fish [28]. In this study, we found that the mRNA level of MC1R was down-regulate in the skin from RCC to WCC. In summary, DNA methylation may be involved in the formation and development of different body color patterns in crucian carp.

In addition, we found that 4-hydroxyphenylpyruvate dioxygenase (HPDA) in the tyrosine metabolic pathway is highly expressed in the skin of RCC, which is the key enzyme in tyrosine metabolism. HPDA catalyzes conversion of hydroxyphenyl pyruvate to ursolic acid (HGA), and HGA enters the TCA cycle through a series of reactions [29]. On the other hand, HPDA reduces blood tyrosine level, thereby reducing the source of melanin [30]. Furthermore, foxd3 mRNA levels in WCC were significantly higher thanthose in RCC. The forkhead transcription factor, Foxd3, represses expression of Mitf, and its negative regulation is also important for the cell fate specification [31,32,33,34]. The DNA methylation level of the Foxd3 gene was higher in WCC than that in RCC, and negatively correlated with its mRNA level. Foxd3 represses expression of mitfa in red and white koi, which not only inhibits the expression of melanocytes but also contributes to the formation of xanthophores and erythrophores [35]. In addition, foxd3 interacts with pax3a to activate the MITF promoter [36] and might respond to interleukin (IL) signaling via a poorly understood pathway [37]. Our results indicate that the pax3a mRNA level was reduced in RCC compared to that in WCC. The mRNA levels of IL3 and IL8 were higher in RCC than those in WCC. These results indicate that foxd3 reduces the formation of melanocytes and may increase the formation of xanthophores and erythrophores. Therefore, the high foxd3 mRNA level in RCC may be associated with the different coloration patterns established in crucian carp.

In addition, the BH4 de novo biosynthetic pathway has been associated with the synthesis of pteridine pigments in goldfish and zebrafish [38]. Three enzymes are involved in this pathway, such as GTP cyclohydrolase I (EC 3.5.4.16; GCH), 6-pyruvoyl tetrahydropterin synthase (EC4.6.1.10; PTPS), and sepiapterin reductase (EC 1.1.1.153). These enzymes catalyze the conversion of GTP to dihydroneopterin triphosphate, dihydroneopterin triphosphate to 6-pyruvoyl tetrahydropterin, and 6-pyruvoyl tetrahydropterin to BH4, respectively [39]. The gch1 and ptps mRNA expression levels were significantly upregulated in RCC compared to those in WCC. This includes introduction of pteridine signaling pathways in the xanthophores and erythrophores of RCC.

Conclusions

We evaluated the global gene expression patterns and DNA methylation patterns in the skin of crucian carp with different coloration patterns. Functional analyses of differentially expressed transcripts revealed that the molecular mechanism of body coloration in crucian carp is strongly related to disruptions in gene expression during pigmentation. Specifically, the expression of genes is associated with different coloration. Our results provide a foundation for further characterization of gene expression during pigmentation in different colored crucian carp. Further biochemical and physiological research on these candidate marker genes should be implemented in the future.

Methods

Material

RCC and WCC, 1 year old, were collected from the Engineering Center for Fish Breeding of the National Education Ministry, Hunan Normal University. Pieces of skin tissue were surgically excised from RCC and WCC, being anesthetized with 100 mg/L MS-222 (Sigma-Aldrich) before dissection. Each sample of skin tissue was collected from siblings, a mixture of six fishes, which were subsequently used for RNA and DNA isolation. All the sampling procedures were conducted according to the standards and ethical guidelines established by the Animal Ethical Review Committee of Hunan Normal University, Changsha, China.

cDNA library construction and sequencing

Total RNA was isolated from skin tissues using Trizol Reagent (Invitrogen) following the manufacturer’s protocol. After removing genomic DNA using DNase I (Fermentas, Vilnius, Lithuania), cDNA was constructed from 2 μg of total RNA per sample following the protocol for the Illumina HiSeqTM 2500 (Illumina Corp., San Diego, CA, USA). Total RNA was purified using RNA Trizol (Invitrogen) and quantified with the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). The cDNA libraries were synthesized using the mixed mRNA fragments as templates and the TruSeq RNA sample preparation kit V2 (Illumina), according to the manufacturer’s instructions. The libraries were sequenced using the Illumina HiSeq 2500 at the BGI Shenzhen, Guangdong, China.

De novo assembly and functional annotation of Transcriptome

Raw reads were filtered using FastQC software to obtain the paired-end clean reads, and all clean reads were used for assembly with Trinity [40, 41]. The assembled transcripts were annotated by BLASTx (National Center for Biotechnology Information [NCBI], Bethesda, MD, USA) against the protein databases of the NCBI, nonredundant (Nr), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG), and Gene Ontology (GO) with an e-value of 1e-5. GO annotation was performed according to the Nr annotation using Blast2GO software [42], and the GO functional classification was classified using WEGO software [43]. The relevant biological pathways were identified through gene enrichment analyses of the KEGG categories and annotation [44]. The number of genes included in each KEGG pathway was counted and the significance of gene enrichment for each KEGG pathway was calculated using the hypergeometric distribution test.

Levels of transcript expression and detection of differentially expressed transcripts

CD-HIT software was used to cluster transcripts and remove redundant sequences in transcriptomes of the two crucian carp strains and the referenced sequence data obtained [45]. Then, the numbers of fragments that mapped to each reference transcript were calculated and BLAT was used to align the paired-end reads to these reference sequences in each transcriptome [46]. The transcript expression level was calculated using the fragments per kilobase per million fragments mapped (FPKM) method [47]. The reference transcripts were filtered out when FPKM <0.07 in all samples of each group. Unigenes with FDR ≤ 0.001 and |log2Ratio| ≥ 2 were considered the DEGs for further analysis.

Quantitative real-time polymerase chain reaction (PCR) analysis

The real-time PCR analysis [48] was performed using the Prism 7500 Sequence Detection System (Applied Biosystems, Foster City, CA, USA) with a miScript SYBR Green PCR kit (Qiagen, Valencia, CA, USA). The amplification conditions were as follows: 50 °C for 5 min and 95 °C for 10 min, followed by 40 cycles at 95 °C for 15 s and 60°Cfor 45 s. The average threshold cycle (Ct) was calculated for each sample using the 2-∆∆Ct method, and normalized to β-actin. Finally, a melting curve analysis was completed to validate the specific generation of the expected product. Primers are available in Additional file 1: Table S3. For each sample, the qRT-PCR analysis was performed by three biological replicates.

DNA sample isolation and MethylRAD library preparation and sequencing

Genomic DNA of the two crucian carp strains was extracted from skin tissues using the conventional cetyltrimethyl ammonium bromide method following the manufacturer’s protocol. The MethylRAD library was prepared by digesting genomic DNA using FspEI (New England Biolabs, Ipswich, MA, USA) at 37 °C for 4 h. Then, the ligated products were amplified in 20 μl reactions. PCR products were purified using a QIAquick PCR purification kit (Qiagen) and were subjected to single end sequencing on an Illumina HiSeq2500 sequencer. The methyl-RAD methylated tag libraries were constructed using the MethylRAD technique and sequenced on the Hiseq X Ten platform. Raw reads were first subject to quality filtering and adaptor trimming. The FspEI sites (CmC/mCDS) extracted from the crucian carp genome were built as reference sites, and high-quality reads were mapped against these reference sites using the SOAP program. The number of reads supporting methylcytosines was calculated to quantifythe MethylRAD data and normalized as RPM (RPM = (read coverage per site/high-quality reads per library) × 1,000,000).

DNA methylation data analysis

Raw reads were trimmed to remove adaptor sequences as well as the terminal 2-bp length from each site to eliminate artifacts. The MethylRAD data were analyzed using reference-based and de novo approaches. The Methyl-RAD sequencing tag for each sample was compared to the reference genome (Reference genome download link: https://ftp.ncbi.nlm.nih.gov/Genomes/all/GCF./000/951/615/GCF_000951615.1_common_carp_genome/GCF_000951615.1_common_carp_genome_genomic.fna.gz) using the SOAP software (ver. 2.21; the primary alignment parameters included: -r 0, −M 4). DNA methylation sites with a sequence depth of no less than 3 were judged to be reliable. The function of the gene was described by a GO function enrichment analysis of the gene where the differential methylation site was located. The number of genes included in each GO entry was counted and the significance of gene enrichment for each GO entry was calculated using the hypergeometric distribution test [49].

Western blotting

We analyzed conservation of some proteins encoded by the DEGs and selected microphthalmia associated transcription factor a (MITFa), tyrosinase protein (TYR), and the dopachrome tautomerase protein (DCT), for which commercially available antibodies were available for western blotting. First, samples containing 50 μg of total lysed protein from RCC and WCC skin tissues was separated by15% sodium dodecyl sulfate-polyacrylamide gel electrophoresis for 1–2 h at 100 V, transferredto a polyvinylidene fluoride membrane (GE Healthcare Bio-Sciences, Uppsala, Sweden), and blocked in TBST containing 5% nonfat milk for 1 h at room temperature. The membrane was incubated with primary antibodies (rabbit polyclonal antibody to MITFa, TYR and DCT, 1:2000; GeneTex Inc., Irvine, CA, USA) overnight at 4 °C. The blots were then washed extensively and incubated with the corresponding peroxidase-conjugated anti-rabbit IgG (1:2000, Sigma St. Louis, MO, USA) for 1 h at room temperature. Proteins were detected using an enhanced chemiluminescence kit (Thermo Fisher Pierce, Rockford, IL, USA) and the ChemiDoc XRS+ imaging system (Bio-Rad, Hercules, CA, USA).

Abbreviations

COG:

Clusters of Orthologous Groups

DCT:

Dopachrome tautomerase

DEGs:

Differentially expressed genes qRT-PCR:quantitative polymerase chain reaction

GO:

Gene Ontology

HPDA:

4-hydroxyphenylpyruvate dioxygenase

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MAPK:

Mitogen activated protein kinase

MITFa:

Microphthalmia associated transcription factor a

NCBI:

National Center for Biotechnology Information

RCC:

Red crucian carp

TYR:

Tyrosinase

TYRP1:

Tyrosinase-related protein 1

VEGF:

vascular endothelial growth factor

WCC:

White crucian carp

References

  1. Komiyama T, Kobayashi H, Tateno Y, Inoko H, Gojobori T, Ikeo K. An evolutionary origin and selection process of goldfish. Gene. 2009;430:5–11.

    Article  CAS  PubMed  Google Scholar 

  2. Wang J, Liu S, Xiao J, Tao M, Zhang C, Luo K, Liu Y. Evidence for the evolutionary origin of goldfish derived from the distant crossing of red crucian carp × common carp. BMC Genet. 2014;15:33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Wang SY, Luo J, Murphy RW, SY W, Zhu CL, Gao Y, Zhang Y. Origin of Chinese goldfish and sequential loss of genetic diversity accompanies new breeds. PLoS One. 2013;8:e59571.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Fujii R. The regulation of motile activity in fish chromatophores. Pigment Cell Res. 2000;13:300–19.

    Article  CAS  PubMed  Google Scholar 

  5. Oshima N. Direct reception of light by chromatophores of lower vertebrates. Pigment Cell Res. 2001;14:312–9.

    Article  CAS  PubMed  Google Scholar 

  6. Kelsh RN. Genetics and evolution of pigment patterns in fish. Pigment Cell Res. 2004;17:326–36.

    Article  CAS  PubMed  Google Scholar 

  7. Liu JH, Zhang YQ, Gui SY, Liu WB, Xiao J, Xiao YM. Observation and regression models on body colour inheritance and development in crucian carp and carp. Aquacult Int. 2016;24:1191–9.

    Article  Google Scholar 

  8. Zhang Y, Liu J, Peng L, Li R, Zhang H, Zou L, Liu W, Xiao Y. Comparative transcriptome analysis of molecular mechanism underlying gray-to-red body color formation in red crucian carp (Carassius Auratus, red var.). Fish Physiol Biochem. 2017;43(5):1387–98.

  9. Koh KP, Rao ADNA. Methylation and methylcytosine oxidation incell fate decisions. Curr Opin Cell Biol. 2013;25:152–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Gifford CA, Ziller MJ, Gu H, Trapnell C, et al. Transcriptional and epigenetic dynamics during specification of human embryonic stem cells. Cell. 2013;153:1149–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Li XM, Song YN, Xiao GB, Zhu BH, Xu GC, Sun XW and Li JT. Gene Expression Variations of Red—White Skin Coloration in Common Carp (Cyprinus carpio). International Journal of Molecular Sciences. 2015;16:21310–29.

  12. Wang S, Lv J, Zhang L, Dou J. MethylRAD: a simple and scalable method for genome-wide DNA methylation profiling using methylation dependent restriction enzymes. Open Biol. 2015;150130:5.

    Google Scholar 

  13. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Statist Soc B. 1995;57:289–300.

    Google Scholar 

  15. Mhanni AA, Mcgowan RA. Variations in DNA(cytosine-5)-methyltrasferase-1 expression during oogenesis and early development of the zebrafish. Dev Genes Evol. 2002;212(11):530–3.

    Article  CAS  PubMed  Google Scholar 

  16. Varriale A, Bernardi G. DNA methylation and body temperature in fishes. Gene. 2006;385:111–21.

    Article  CAS  PubMed  Google Scholar 

  17. Steingrimsson E, Copeland NG, Jenkins NA. Melanocytes and the microphthalmia transcription factornetwork. Annu Rev Genet. 2004;38:365–411.

    Article  CAS  PubMed  Google Scholar 

  18. Planque N, Raposo G, Leconte L, Anezo O, Martin P, Saule S. Microphthalmia transcription factor induces both retinal pigmented epithelium and neural crest melanocytes from neuroretina cells. J Biol Chem. 2004;279:41911–7.

    Article  CAS  PubMed  Google Scholar 

  19. Tachibana M, Takeda K, Nobukuni Y, Urabe K, Long JE, Meyers KA, Aaronson SA, Miki T. Ectopic expression of MITF, a gene for Waardenburg syndrome type 2, converts fibroblasts to cells with melanocyte characteristics. Nat Genet. 1996;14:50–4.

    Article  CAS  PubMed  Google Scholar 

  20. Cheli Y, Ohanna M, Ballotti R, Bertolotto C. Fifteen-year quest for Microphthalmia associated transcription factor target genes. Pigment Cell Melanoma Res. 2010;23:27–40.

    Article  CAS  PubMed  Google Scholar 

  21. Lister JA, Robertson CP, Lepage T, Johnson SL, Raible DW. Nacre encodes a zebrafish microphthalmia-related protein that regulates neural-crest-derived pigment cell fate. Development. 1999;126:3757–67.

    CAS  PubMed  Google Scholar 

  22. Hemesath TJ, Price ER, Takemoto C, Badalian T, Fisher DE. MAP Kinase links the transcription factor Microphthalmia to c-kit signalling in melanocytes. Nature. 1998;391:298–301.

    Article  CAS  PubMed  Google Scholar 

  23. Wu M, Hemesath TJ, Takemoto CM, Horstmann MA, Wells AG, Price ER, Fisher DZ, Fisher DE. C-kit triggers dual phosphorylations, which couple activation and degradation of the essential melanocyte factor mi. Genes Dev. 2000;14:301–12.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. Mansky KC, Sankar U, Han J, Ostrowski MC. Microphthalmia transcription factor is atarget of the p38 MAPK pathway in response to receptor activator of NF-kB ligand signaling. J Biol Chem. 2002;277:11077–83.

    Article  CAS  PubMed  Google Scholar 

  25. Price ER, Ding HF, Badalian T, Bhattacharya S, Takemoto C, Yao TP, Hemesath TJ, Fisher DE. Lineage-specific signaling in melanocytes.C-kit stimulation recruits p300/CBP to microphthalmia. J. Biol.Chem. 1998;273:17983–6.

    Article  CAS  PubMed  Google Scholar 

  26. Voisey J, Box NF, van Daal AA. Polymorphism study of the human agouti gene and its association with MC1R. Pigment Cell Res. 2001;14:264–7.

    Article  CAS  PubMed  Google Scholar 

  27. Slominski A, Tobin DJ, Shibahara S, et al. Melanin pigmentation in mammalian skin and its hormonal regulation. Physiol Rev. 2004;84:1155–48.

    Article  CAS  PubMed  Google Scholar 

  28. oshima N, goto M. Prolactin signaling in Erythrophores and Xanthophores of teleost fish. Pigment Cell Res. 2000;13:35–40.

    Article  PubMed  Google Scholar 

  29. Bradle FC, Lindstedt S, Lipscomb JD, et al. 4-Hydroxyphe-nylpyruvate dioxygenase is an iron-tyrosinate protein. J Biol Chem. 1986;261(25):11693–6.

    Google Scholar 

  30. Crouch NP, Adlington RM, Baldwin JE, et al. A mechanistic rationalisation for the subst rate specificity of recombinant mammalian 4-hydroxyphenylpyruvate dioxygenase (4-HPPD). Tet rahedron. 1997;53(20):6993–7010.

    Article  CAS  Google Scholar 

  31. Curran K, Raible DW, Lister JA. Foxd3 controls melanophore specification in the zebrafish neural crest by regulation of Mitf. Dev Biol. 2009;332:408–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Yan B, Liu B, Zhu C-D, Li K-L, Gong X-L, Wang C-H, et al. microRNA regulation of skin pigmentation in fish. J Cell Sci. 2013;126:3401–8.

    Article  CAS  PubMed  Google Scholar 

  33. Ignatius MS, Moose HE, El-Hodiri HM, Henion PD. Colgate/hdac1 repression of foxd3 expression is required to permit mitfa-dependent melanogenesis. Dev Biol. 2008;313:568–83.

    Article  CAS  PubMed  Google Scholar 

  34. Thomas AJ, Erickson CA. FOXD3 regulates the lineage switch between neural crest-derived glial cells and pigment cells by repressing MITF through a non-canonical mechanism. Development. 2009;136:1849–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Jiang Y, Jiang Y, Zhang S, Xu J, Feng J, Mahboob S, Al-Ghanim KA, Sun X, Xu P. Comparative transcriptome analysis reveals the genetic basis of skin color variation in common carp. PLoS One. 2014;9:e108200.

  36. Watanabe A, Takeda K, Ploplis B, Tachibana M. Epistatic relationship between Waardenburg syndrome genes MITF and PAX3. Nat Genet. 1998;18:283–6.

    Article  CAS  PubMed  Google Scholar 

  37. Kamaraju AK, Bertolotto C, Chebath J, Revel M.  Pax3 down-regulation and shut-off of melanogenesis in melanoma B16/F10.9 by interleukin-6 receptor signaling. J Biol Chem. 2002;277:15132–41.

  38. Ziegler I. The pteridine pathway in zebrafish: regulation and specification during the determination of neural crest cell-fate. Pigment Cell Res. 2003;16:172–82.

    Article  CAS  PubMed  Google Scholar 

  39. Woo HJ, Hwang YK, Kim YJ, Kang JY, Choi YK, Kim C, Park YS. Escherichia Coli 6-pyruvoyltetrahydropterin synthase ortholog encoded by ygcM has a new catalytic activity for conversion of sepiapterin to 7, 8-dihydropterin. FEBS Lett. 2002;523:234–8.

    Article  CAS  PubMed  Google Scholar 

  40. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, MacManes MD, Ott M, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protocols. 2013;8:1494–512.

    Article  CAS  PubMed  Google Scholar 

  41. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotech. 2011;29:644–52.

    Article  CAS  Google Scholar 

  42. Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, Wang J, Li S, Li R, Bolund L, Wang J. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;34:W293–7.

  43. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:D480–4.

  44. Conesa A, Götz S, García-Gómez J M, et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

  45. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.

    Article  CAS  PubMed  Google Scholar 

  46. Kent WJ. BLAT—the BLAST-like alignment tool. Genome Res. 2002;12(4):656–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Meth. 2008;5:621–8.

    Article  CAS  Google Scholar 

  48. Xu K, Wen M, Duan W, Ren L, Hu F, Xiao J, Wang J, Tao M, Zhang C, Wang J, Zhou Y, Zhang Y, Liu Y, Liu S. Comparative analysis of testis Transcriptomes from triploid and fertile diploid cyprinid fish. Biol Reprod. 2015;92(4):1–12.

  49. Cohen-Karni D, Xu, D, Apone L, Fomenkov A, Sun Z, Davis PJ. et al. The MspJI family of modification-dependent restriction endonucleases for epigenetic studies. Proc Natl Acad Sci U S A. 2011;108:11040–5.

Download references

Acknowledgments

The authors would like to express their thanks to the anonymous referees for their constructive comments on the paper, which have greatly improved its presentation.

Funding

This work has been supported by the National Science Foundation of China (No. 31472272), the National Basic Research Program of China (No. 2012CB722305), Hunan Provincial Natural Science Foundation of China (No. 14JJ2058), the Scientific Research Fund of Hunan Provincial Education Department (No.14FJ4261), the Cooperative Innovation Center of Engineering and New Products for Developmental Biology of Hunan Province (No. 20134486).

Availability of data and materials

The complete clean reads have been uploaded to the SRA site (http://www.ncbi.nlm.nih.gov /sra/).

The differential DNA methylation sites were located (link: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/951/615/GCF_000951615.1_common_carp_genome/GCF_000951615.1_common_carp_genome_genomic.gff.gz).

Author information

Authors and Affiliations

Authors

Contributions

YMX and YQZ designed the experiments. YQZ, JHL, WF, WTX, HQZ, SJC and WBL carried out the experiments. YQZ, YMX, JHL, LYP, and WBL conducted the statistical analysis and discussion. YQZ, YMX and LYP organized and wrote the manuscript. All the authors have read the paper and agreed to list their names as coauthors.

Corresponding author

Correspondence to Yamei Xiao.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the Animal Ethical Review Committee (AERC) of Hunan Normal University, Changsha, China. The collection of fish samples was permitted by the Engineering Center of Fish Breeding of the Education Ministry, Hunan Normal University. All sampling procedures were conducted according to the standards and ethical guidelines established by the AERC. And the researchers involved in animal experiment were certified under a professional training course for laboratory animal practitioners held by the Institute of Experimental Animals, Hunan Province, China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1: Table S1.

Statistics of function annotation of Red and White skin crucian carp. Table S2. Sample sequencing data volume and comparison rate. Table S3. List of primers for the qRT-PCR validation of Unigenes identified. Figure S1. The GO function is classified by histogram of differences in methylation sites. (DOCX 107 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, Y., Liu, J., Fu, W. et al. Comparative Transcriptome and DNA methylation analyses of the molecular mechanisms underlying skin color variations in Crucian carp (Carassius carassius L.). BMC Genet 18, 95 (2017). https://doi.org/10.1186/s12863-017-0564-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12863-017-0564-9

Keywords