Identification and characterization of skin color microRNAs in Koi carp (Cyprinus carpio L.) by Illumina sequencing

Background MicroRNAs (miRNAs) are endogenous, small (21–25 nucleotide), non-coding RNAs that play important roles in numerous biological processes. Koi carp exhibit diverse color patterns, making it an ideal subject for studying the genetics of pigmentation. However, the influence of miRNAs on skin color regulation and variation in Koi carp is poorly understood. Results Herein, we performed small RNA (sRNA) analysis of the three main skin colors in Koi carp by Illumina sequencing. The results revealed 330, 397, and 335 conserved miRNAs (belonging to 81 families) and 340, 353, and 351 candidate miRNAs in black, red, and white libraries, respectively. A total of 164 differentially expressed miRNAs (DEMs) and 14 overlapping DEMs were identified, including miR-196a, miR-125b, miR-202, miR-205-5p, miR-200b, and etc. Target prediction and functional analysis of color-related miRNAs such as miR-200b, miR-206, and miR-196a highlighted putative target genes, including Mitf, Mc1r, Foxd3, and Sox10 that are potentially related to pigmentation. Determination of reference miRNAs for relative quantification showed that let-7a was the most abundant single reference gene, and let-7a and miR-26b was the most abundant combination. Conclusions The findings provide novel insight into the molecular mechanisms determining skin color differentiation in Koi carp, and serve as a valuable reference for future studies on tissue-specific miRNA abundance in Koi carp. Electronic supplementary material The online version of this article (10.1186/s12864-018-5189-5) contains supplementary material, which is available to authorized users.


Background
Koi carp, a colorful variant of common carp (Cyprinus carpio L.), has been artificially selected for centuries, and is a popular pet and hobby species that exhibits a wide variety of colors and phenotypes [1]. Unlike the slow rate of phenotypic changes occurring in wild populations, the rapid establishment of diverse colors and coloration patterns under artificial selection makes Koi carp an ideal subject for studying the genetics of pigmentation [2]. Determination of skin color is a complicated process in fish, linked to various cellular, genetic, nutritional, physiological, and environmental factors [3]. Exactly how skin color patterns form is a long standing question among biologists. Extensive studies have identified a series of genes involved in the determination of skin pigmentation, including Melanocortin receptor 1 (Mc1r), Microphthalmia-associated transcription factor (Mitf ), Tyrosinase (Tyr), Tyrosine related protein-1 (Tyrp1), and Melanocyte-stimulating hormone (Msh) [4][5][6]. Additional genes related to the determination of skin color in fish remain to be elucidated.
MicroRNAs (miRNAs) are single-stranded, non-coding, highly-conserved 19-24 nucleotides (nt) molecules that regulate gene expression at the post-transcriptional level by directly targeting RNA-inducing silencing complex (RISC) to cognate messenger RNA targets [7,8]. Crosstalk between miRNAs and mRNAs is important for transcriptional and signal transduction events involved in multiple biological processes such as apoptosis, cell proliferation, cancer, embryo development, and skin pigmentation [9][10][11][12]. Therefore, skin-expressed miR-NAs might play a vital role in skin differentiation, color formation and skin diseases. Previous studies have investigated the expression and functions of some miRNAs in animal skin. For example, the skin-specific miR-203 was found to define a molecular boundary between proliferative basal progenitors and terminally differentiating suprabasal cells, ensuring proper identify of neighboring layers [13]. In Drosophila, dorsal abdomen pigmentation is decreased when miR-8 lost [14], and coat color in mice is affected when miR-137 is over-expressed [15]. In white alpaca (WA) and brown alpaca (BA) skin, deep sequencing identified 35 and 13 conserved differentially expressed miRNAs (DEMs), indicating potential functions in coat skin color regulation [16]. Meanwhile, Yan et al. identified 13 DEMs between red and white skin in common carp, miR-429 is a potential regulator because it's silencing resulted in an obvious change in skin pigmentation [12]. Thus, miRNAs could be involved in the regulation of skin color. However, the molecular and cellular mechanisms regulating skin color variation in fish, especially a variety of colors appeared on a single fish, such as Koi carp (Fig. 1), remain poorly understood.
Reference genes are employed as internal controls for normalizing real-time quantitative reverse transcription PCR (qRT-PCR) data during analysis of target gene expression, and this approach can capture all non-biological variation [17]. At present, β-actin, 18S rRNA, 5S rRNA and glycer-aldehyde-3-phosphate dehydrogenase (GAPDH) are most widely used as candidate reference genes for mRNA expression [18]. In typical miRNA expression studies, 18S rRNA and small nuclear RNAs such as U6 snRNA are commonly used as reference genes [19,20]. The selection of appropriate reference genes is not trivial, and studies indicate that a single gene is unlikely to behave consistently in all tissues and/or in all physiological and pathological situations [21,22]. Therefore, internal genes should be verified and selected prudently. However, to date, few reports have analyzed the expression stability of miRNAs used as reference in teleosts.
Herein, we investigated the potential functions of miRNAs in post-transcriptional regulation of pigmentation by constructing Koi carp miRNA libraries from white, red, and black fish skin, identified hundreds of conserved miRNAs and candidate novel miRNAs, and revealed evidence of differential expression of miR-NAs in different colors. We examined candidate reference miRNAs in 15 tissues (blood, brain, muscle, liver, kidney, spleen, intestine, heart, gill, fin, eye, ovary, and white/red/black skin) to minimize the risk of co-regulation artifacts. The results provide fundamental knowledge for unraveling the genetic mechanisms underlying color trait variation, and could assist the exploitation of gene regulation linked to skin color determination.

Sample collection
The Koi carp used in this study were obtained from the aquaculture base of Freshwater Fisheries Research Center (FFRC) affiliated with the Chinese Academy of Fishery Sciences. Animals (average weight: 1000 ± 80 g) were raised at 24 ± 1°C in 256 L tanks in a circulation water system for 1 week before experiments and fed twice daily (in the morning and late afternoon). Aeration was supplied constantly and a 12 h light / 12 h dark photoperiod was employed.
All of the surviving fish were tranquilized in 20-30 mg/L MS-222 buffered to pH 7.0-7.5, and three sets of skin tissues were collected from nine Koi carp (black, white and red spots in one fish synchronous); three black skin (BS), white skin (WS), and red skin (RS) samples per group, and another 12 tissue samples (blood, brain, muscle, liver, kidney, spleen, intestine, heart, gill, fin, eye, and ovary) were also collected. All fresh tissue samples were immediately snap-frozen in liquid nitrogen and stored at − 80°C until use.

RNA isolation, library preparation and sequencing
Total RNA was obtained from Koi carp samples using TRIzol (Invitrogen, UK) according to the manufacturer's protocol, and genomic DNA was removed using DNase I (New England Biolabs). RNA purity was assessed using the Nanodrop-2000 instrument (Thermo Scientific, USA), and all RNA samples had an A260:A280 ratio > 1.9. An equal amount of total RNA from three fish individuals from each color group (black, white, and red spots) was pooled, and a total of nine RNA pools were prepared and used for library construction. A pair of Illunima adapters was ligated to the RNA 5′ and 3′ termini of RNAs, and the resultant products were reverse-transcribed and  was mapped to only one annotation, we followed the rule rRNAetc (Genbank> Rfam) > known miRNA >repeat >exon >intron3, as previously described [23]. Subsequently, the remaining reads were identified as conserved miRNAs by BLAST searching against miRBase 21.0, allowing no more than two mismatches [24]. Sequences that did not match existing or conserved miRNAs were used to identify potentially novel miRNA candidates [25,26]. The enrichment degree of each miRNA was identified by counting the number of reads in each sample, and miRNAs were family clustered in each library through miRBase 21.0 database.

Differential expression of miRNAs
To compare the expression levels of miRNAs among the black, white, and red skin libraries, the frequency of miRNA counts was normalized as transcripts per million (TPM). Calculated as follows: normalized expression (TPM) = (actual miRNA count/number of total clean reads) × 1,000,000. Only miRNAs with altered by more than two-fold relative to controls were considered significantly differentially expressed (p-value < 0.05) [27] and classified as DEMs. A positive value indicated up-regulation, while a negative value indicated down-regulation. Fold-change was calculated as follows: Fold-change = log2 (treatment/control).
The p-value was calculated as follows: Prediction and analysis of miRNA target genes Target genes of DEMs were identified using the miRanda [28], RNAhybrid [29] and TargetScan [30], and overlapping DEMs from all three programs provided the final results. The rules used for target prediction were based on those suggested previously by Allen et al... [31] and Schwab et al [32], as follows: (1) No more than four mismatches between sRNA and target (G-U bases count as 0. Functions strongly associated with the predicted target genes of miRNAs were determined using Gene Ontology (GO; http://www.geneontology.org) [33] biological process categories, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was conducted using pathway database (http://www.genome.jp/kegg/pathway.html) [34] to statistically test DEMs.

QRT-PCR and analysis of reference miRNAs
Total RNA was extracted as described above. For reverse-transcription of miRNAs, the PrimeScript RT Reagent Kit (Takara Bio, Dalian, China) were used. qRT-PCR was performed on a CFX-96 Real-time PCR System (Bio-Rad, CA, USA) in 25 μL reactions containing 12.5 μL SYBR Advantage Premix (2×) reagent (Takara Bio), 0.5 μL miRNA-specific forward primer (10 μM), 0.5 μL miScript universal primer (10 μM), and 2 μL PCR template (cDNA). Amplification was performed with an initial denaturation at 95°C for 10s, followed by 40 cycles at 95°C for 5 s and 60°C for 20s, and a final cycle from 95°C to 65°C. All reactions were conducted in triplicate, which six biological replicates. The relative expression levels of DEMs were measured in terms of threshold cycle value (Ct) and were normalized to U6 snRNA using the eq. 2-ΔΔCt method [35]. Candidate reference miRNAs used for evaluation included let-7a, miR-140-3p, miR-21, miR-15c, miR-26b, miR-92-5p, and miR-145a-5p due to their low expression differences according to the Illumina deep sequencing results (data not shown). Additionally, 5 s rRNA, 18 s rRNA and U6 snRNA were also selected for evaluation as reference genes due to their broad use in the literature. QRT-PCR methods were the same as described above for DEMs, and the stability of reference miRNAs was tested using geNorm [36], NormFinder [37] and Bestkeeper [38] programs. All primers used in qPCR experiments (Additional file 1: Table S1) were synthesized by Sangon Biotech. (Shanghai, China). 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).

Identification of Koi carp miRNAs via sRNA sequencing
To identify miRNAs expressed in the skin of Koi carp, BS, RS, and WS sRNA libraries were determined and analyzed by Illumina deep sequencing. A Total of 20,502,959, 16,916,198, and 19,839,224 reads were obtained from the BS, RS, and WS libraries, respectively. After filtering out low-quality reads, insert-null, reads < 18 nt, and removing adaptor sequences, 19,649,441 (96.01%), 16,193,916 (95.91%), and 18,981,507 (95.87%) clean reads were retrieved for further analysis (Additional file 2: Table S2). In order to annotate the different classes of sRNAs, clean reads were compared with GenBank, Rfam, miRbase, exon and intron, and repeat-associated RNA databases. Representations of different types of RNA sequences, including rRNAs, tRNAs, snRNAs, snoRNAs, and scRNAs are shown in Table 1.
The length distribution of sRNA reads is shown in Fig. 2. The results were similar among the three BS, RS, and WS groups, with 22-or 23-nt sRNAs accounting for 73.82, 74.20, and 73.10% of total sequences, respectively (Fig. 2).

Characterization of miRNAs in different skin tissues
To further identify conserved miRNAs and predict novel miRNAs in the three skin color tissues, all sRNA sequences were mapped to known miRNAs in miRBase 21.0. After alignment and additional sequence analysis, 330, 397, and 335 conserved miRNAs were identified from BS, RS, and WS libraries, respectively, that are highly conserved in other species (Additional file 3). These miRNAs exhibited a broad range of expression levels, ranging from 32,428,443 counts for the most abundant, to a single count. The results including known carp miRNAs in the miRBase database for common carp, including 137, 144, and 138 miRNAs for BS, RS, and WS groups, respectively, of which miR-199-5p, miR-199-3p, let-7a, miR-21, miR-26a, miR-99, miR-100, miR-125b, miR-22a, miR-126-3p, and miR-146a were the most abundant (Additional file 4). To obtain a clearer perspective of the most abundant known carp miRNAs, we compared those with the 10 highest read numbers ( Table 2) and~81 miRNA families were identified in each library. Several miRNA families were expressed at levels > 100,000 reads, such as let-7, miR-199, miR-126, miR-27, miR-21, etc. The degree of sequence conservation of miRNAs was assessed to investigate phylogenetic evolution, and 20 conserved families were clustered into three groups based on phylogenetic distribution (Fig. 3). Five miRNA families (let-7, miR-124, miR-184, miR-375, and miR-31) are shared in protostomes and deuterostomes, 13 families are presented only in vertebrate, and the other two families (miR-135 and miR-727) have only been identified in fish.
We also found 1563 miRNAs not previously identified in Koi carp, comprising 340, 353, and 351 in BS, RS, and WS groups, respectively (Additional file 5). These were considered specific candidate miRNAs for further investigation. The seven most abundant miRNAs (novel-miRn0378, novel-miRn0924, novel-miRn1323, novel-miRn0321, novel-miRn0902, novel-miRn0547, and novel-miRn1422) all had > 300 reads (Table 3). We also investigated the most abundant novel miRNAs in each skin tissue, and novel-miRn0022 and novel-miRn0032 were dominantly expressed in white skin rather than black or red skin, with 1165 and 1108 reads, respectively (Table 4).

Expression profiling of miRNAs in different skin patterns
In the high-throughput sequencing analysis, many miR-NAs exhibited dissimilar expression levels among the three skin color samples. Using the criteria |log 2 (Foldchange)| ≥ 1 and q-value ≤0.001, volcano plots of the three pairwise comparisons (BS vs. WS, RS vs. BS, and RS vs. WS) revealed the expression trends for each pair (Fig. 4a). We also constructed a histogram of DEMs in the three skin tissues (Fig. 4b). Compared with black skin, 53 miRNAs were up-regulated in white skin, including miR-200b and miR-135c, while 36 miRNAs were down-regulated. A total of 50 miRNAs were down-regulated in white skin compared with red skin, including miR-217, while 45 miRNAs were up-regulated. There were 47 miRNAs displaying greater abundance in black skin compared with red skin, including miR-196a, while 34 miRNAs were down-regulated (Additional file 6). The proportion of DEMs reflects their specific functions and related biological mechanisms in different skin tissues. A Venn diagram (Fig. 4c) was generated to visually compare the expression of miRNAs in Koi carp. Among the significant DEMs, 164 miRNAs (including 14 overlapping DEMs) were identified, and 30 miRNAs (including 17 conserve and 13 novel DEMs) were significant DEMs in WS and RS group compared with the BS group, suggesting that they likely play a key role in the color variation process (Fig. 4c). The heatmap of these 17 known DEMs (Fig. 5) indicates that eight miRNAs including miR-203b-3p, miR-26-3p, miR-205-5p, etc. were up-regulated in WS and RS groups compared with the BS group, while six miRNAs, including aca-miR-210-3p, miR-125c, miR-206, etc. were down-regulated in WS and RS groups compared with BS group. The remaining miRNAs (ccr-miR-196a, miR-202, and ipu-miR-196a) exhibited different expression trends (Fig. 5).
In order to validate the reliability of the Illumina sequencing results, expression in different skin tissues was examined by qRT-PCR for 18 miRNAs, including  conserved and specific miRNAs, selected randomly from the dataset. The results showed that qRT-PCR expression patterns of all selected miRNAs were absolutely in agreement with the results of RNA-seq analysis, indicating that the RNA-seq data was reliable (Fig. 6). For example, based on the deep sequencing result, the expression level of let-7f-5p in RS was almost 10 times higher than in WS, compared with an 8-fold difference in the qRT-PCR results.

Target genes prediction and signaling pathways analysis
Target genes were predicted based on the common carp (C. carpio) genome sequence (http://www.fishbrowser.org/database/Commoncarp_genome/) using a combination of TargetScan, RNAhybrid, and miRandato refine and improve the prediction results. A total of 2495 putative target genes were predicted, and some unique miR-NAs, such as miR-206, miR-125c, miR-217, miR-205-5p, and novel-miRn1185, appear to target hundreds of genes. Moreover, our results indicate that some genes may be regulated by more than one miRNAs, such as cytoglobin-1-like (No: DN51764) that is targeted by nine miRNAs. Target genes were subjected to GO analysis for functional classification (Additional file 7). The result showed that many miRNAs were assigned to pigmentation-related terms such as melanocyte migration (GO: 0097324), retinal pigment epithelium development (GO: 0010669), and melanosome transport (GO: 0032402). These genes enriched in pigmentation-related processes are informative and worthy a further study. KEGG pathway analysis was also performed to classify putative miRNA targets, and 275 annotated signaling pathways were found to be involved in a wide range of biological functions. For instance, at least 10 pathways involved in pigmentation-related biology were identified, including Wnt/β-catenin signaling, MAPK signaling, melanogenesis, oxidative phosphorylation, cell cycle, TGF-beta signaling, melanoma, cAMP signaling, epidermal growth factor receptor (EGFR) tyrosine kinase inhibitor resistance, and Notch signaling (Fig. 7). Herein, we mainly focused on melanogenesis, Wnt signaling, and MAPK signaling pathways, and putative gene    results showed that miR-133c, miR-26a-3p, miR-196a, and miR-181a-5p all target the Mitf gene. Foxd3 is one of the earliest molecular makers in the Guanine pathway, and altering the expression level of this gene can control lineage choice between neural or glial and pigment cells, as demonstrated by repressing Mitf during the early phase of neural crest migration [42]. Herein, Foxd3 was found to be targeted by miR-429b, miR-26a-3p, miR-200b, miR-10b-5p, and miR-181a-5p.

Identification of optimal reference genes
As shown in Fig. 9 Analysis of miRNA expression stability (M) values using the geNorm algorithm was applied to each reference gene based on average pairwise variation between all candidate genes (M value is negatively correlated with miRNA stability; the lower the value, the higher the stability, and vice versa). The results showed that comparison of let-7a and miR-26b gave an M value of 0.661, which indicates high stability. The stability of the 10 reference genes in different tissues was ranked let-7a and miR-26b > U6 > miR-145a-5p > miR-21 > miR-92a-5p > miR-15c > miR-140-3p > 5 s rRNA >  Table 5). The NormFinder algorithm is similar to that used by geNorm, and the lower the expression of internal reference genes, the higher the stability. As shown in Fig. 10, the NormFinder reference gene stability in different Koi carp tissues was ordered let-7a > U6 > miR-26b > miR-15c > miR-145a-5p > miR-92a-5p > miR-21 > 5 s rRNA > miR-140-3p > 18 s rRNA, hence let-7a was the most stable gene. In BestKeeper analysis, based on average Ct values, the stability of internal reference genes should have a low standard deviation (SD) and coefficient of variation (CV), as well as a high Pearson correlation coefficient (R). The BestKeeper results showed that the best-performing internal reference gene was let-7a, followed by miR-26b, while the most unstable gene was miR-140-3p, followed by 18 s rRNA (Table 6). Taking the results of these three analyses together, let-7a was the most stable single reference gene, let-7a and miR-26b was the most stable reference combination in different tissues, while 18 s rRNA was the most unstable gene.

Discussion
Body color and coloration patterns are prominent features related to intraspecific communication, interspecific interactions, photoprotection, photoreception, thermoregulation, and market value in many animals, including Koi carp. Thus, in the present work we explored the molecular mechanisms of pigmentation regulation and evolution to facilitate improving coloration and ornamental value [43]. Increasing evidence suggests that miRNAs play a vital role in skin pigmentation by controlling the expression of a vast number of genes [44,45]. Herein, we generated miRNA transcriptome data from black, white, and red skin samples of Koi carp using a deep sequencing approach, and the findings provide a foundation for future skin pigmentation studies.
The BS, WS, RS miRNA libraries constructed from different colored skin tissue displayed a similar read length distribution of 22-23 nt miRNAs, equivalent in size to typical Dicer-derived products. The most abundant miRNAs were let-7a, miR-21, miR-26a, miR-15c, etc. with > 10,000 read each. These miRNAs are also highly expressed in other species; let-7 family was first characterized in Caenorhabditis elegans, and plays a key role in regulating late developmental events by down-regulating lin-41 and possibly other genes [46]; miR-21 is an abundantly expressed miRNA in mammalian cells, and is evolutionarily conserved across a wide range of vertebrate species. Additionally, it is significantly up-regulated in gastric cancer, and it targets 15-PGDH, suggesting the PGE2/PI3K/Akt/Wnt/β-catenin axis might be a novel pathway for gastric cancer treatment [47]. Thus, these miRNAs are considered to be housekeeping molecules involved in the maintenance of basic cellular activities.
To gain insight into the possible significance of skin miRNAs in Koi Carp, we determined expression profile of miRNAs in different colored (red, black, and white) skin tissues. We found that 30 miRNAs (17 conserved and 13 novel) were significant DEMs in WS and RS groups compared with the BS group, and 14 overlapping DEMs were shared by all three skin tissues, including miR-196a, miR-125b, miR-202, miR-125c, miR-205-5p, miR-200b, and etc. implying that these DEMs may therefore be involved in skin pigmentation. Knowledge of tissue-and cell-specific expression patterns of miRNAs could directly inform functional studies [48]. A previous study has revealed miR-200b and miR-196a as candidates that may be important for hair follicle (HF) development because their expression was reduced in the Dkk1 transgenic skin [49]. Meanwhile, miR-125b serves as a rheostat that controls stem cell proliferation, fate commitment, and differentiation [50]. Serum levels of miR-205-5p were determined by real-time PCR in 11 patients with metastatic melanoma and 16 patients without metastasis, and differences indicated that the miRNA network may be involved in the pathogenesis of melanoma metastasis [51]. Also, a significant correlation was observed between miR-9-5p and miR-205-5p, and between miR-203a-3p and miR-205-5p, suggesting that expression levels of these miRNA were not independently regulated in melanoma patients [51]. Herein, we found that miR-196a, miR-202, miR-125b andmiR-125c, etc. were highly expressed in white and red skin compared with black skin, indicating specific functions in difference tissues. Thus, specific DEMs may regulate skin pigmentation in fish. MiRNA target identification is important for predicting their functions. Although computational approaches have been widely used to predict targets, most methods suffer from a high false-positive rate [52]. Therefore, we used three software packages for analysis of target genes to refine and improve the results. We found that some miRNAs such as miR-206, miR-125c, miR-217, etc. targeted hundreds of genes, and some genes were regulated by more than one miRNA, indicating a complicated regulatory network between miRNAs and their targets gene. Skin pigmentation in fish is a complex process which involves numerous physiological, cellular and genetic factors [53]. Previous studies found that signaling pathway such as the Wnt, MAPK, cAMP, melanogenesis, Notch, and Fgf pathways strongly influence skin pigmentation process [41]. Herein, we constructed the putative gene pathways with some miRNAs involved in red, black, and white color skins pigmentation process. The molecular mechanism of melanin biosynthesis was extensively studied due to its biomedical significance, in mammals and birds, the ratio of eumelanin and pheomelanin largely determines an animal's overall color: darker (black to brown) phenotypes result from the increased deposition of eumelanin, while lighter (red to yellow) phenotypes result from the increased deposition of pheomelanin [54]. Xu (2014) also elucidated the pigmentation mechanisms of Hebao and Songpu carp, in the melanogenesis pathway, tyrosine is oxidized to form dopaquinone, and then intracellularly catalyzed to become eumelanin; However, cystine and dopaquinone can switch off the eumelanin synthesis pathway and promote   [54,55,57]. (1) represents the pheomelanin biosynthesis pathway; (2) represents the eumelanin pathway. (3) represents the guanine pathway the synthesis of pheomelanin; And slc7a11 encodes the transmembrane cystine/glutamate exchanger (xCT), which transports cystine into melanocytes to synthesize pheomelanin [55]. Both miR-200 and miR-205 are highly expressed in normal hair follicles compared epithelia tumors in human [56]. The Mc1r gene is a key factor involved in the production of pheomelanin and eumelanin in melanocytes, and we found that it is targeted by miR-200a, miR-200b, and miR-206, and miR-206 presented up-regulated in black skins compared with red and white skins suggested that it played important roles in pigmentation in Koi carp. Wnt and MAPK are implicated in numerous development and physiological process. Sun reported that both Wnt and MAPK signaling pathways are highly likely to be involved in melanin biosynthesis in common carp [57]. Mitf is a member of the Myc-related family of basic helix-loop-helix leucine zipper (bHLH-Zip) transcription factors, and Mitf was found in cells related to the retinal pigment epithelium (RPE), as well as cells behind the optic cup that were probably derived from neural crest tissue and could develop into iridophore pigment cells [58]. Mutations of the Mitf gene causes a variety of phenotypes, most notably in pigmented cells [59]. Foxd3 can affect the lineage between neural or glial and pigment cells by repressing Mitf during the early phase of neural crest migration Fig. 10 The stability of candidate reference genes analyzed by NormFinder. Values with a different superscript letter denote significant differences between genes (P < 0.05)  [42]. Kevin et al... reported a model which melanophores and iridophores descend from a common precursor whose fate is regulated by the interplay transcription factors regulating Mitf and Foxd3, and loss of Mitfa resulted in supernumerary ectopic iridophores, while loss of Foxd3, a Mitfa repressor, resulted in fewer iridophores [60]. Double mutants displayed restoration of iridophores, suggesting that one of Foxd3 is to suppress Mitfa to promote iridophore development [60]. In our present study, we also observed miR-429b, miR-26a-3p, and miR-181a-5p can target the Foxd3 gene, and miR-133c, miR-196a act on the Mitf gene. Therefore, we have carried out a partial analysis of the regulatory network of the three color skins, while the specific mechanisms and functional need to be deeply excavated. A number of reference genes for mRNA qRT-PCR analysis have been established, but miRNA qRT-PCR studies are fewer in number. Thus, we attempted to select suitable reference genes for normalization to eliminate or decrease non-biological variation between test samples. The two most commonly used reference genes for normalization in miRNA analyses are U6 and 18S rRNA, but there has been some debate as to the validity since they are expressed at levels greater than target genes [61,62]. . found that miR-101a was the most stable miRNA when all tissue types were considered separately [63]. In Siniperca chuatsi, all six miR-NAs (miR-101a, miR-146a, miR-22a, miR-23a, miR-26a and let-7a) exhibited better expression consistency than U6 in most of the conditions examined, and the best combination of reference genes was miR-22a and miR-23a [64]. Herein, we systematically evaluated the reliability of seven miRNA genes across 15 tissues. Our results suggest that let-7a was the most stable single reference gene, and the best combination of reference genes was let-7a and miR-26b, followed by U6, while 18 s rRNA was the most unstable gene.

Conclusions
Our findings provide fundamental information on the expression of conserved and potential novel miRNAs in Koi carp skin, and identified DEMs and their respective target genes that may be related to skin color. We identified 137, 144, and 138 miRNAs for BS, RS, and WS groups which already present in the C.carpio miRBase database, and miR-203b-3p, miR-125c, ccr-miR-196a, miR-202, and etc. were presented differences in expression in the three different skin color tissue samples. The best combination of reference genes was let-7a and miR-26b. Validation of the functional role of DEMs and reference genes will provide a better understanding of miRNA-mediated regulatory mechanisms influencing skin color in fish.