microRNA regulation of skin pigmentation in golden-back mutant of crucian carp from a rice-fish integrated farming system
BMC Genomics volume 24, Article number: 70 (2023)
MicroRNAs (miRNAs) are endogenous small non-coding RNAs (21–25 nucleotides) that act as essential components of several biological processes. Golden-back crucian carp (GBCrC, Carassius auratus) is a naturally mutant species of carp that has two distinct body skin color types (golden and greenish-grey), making it an excellent model for research on the genetic basis of pigmentation. Here, we performed small RNA (sRNA) analysis on the two different skin colors via Illumina sequencing.
A total of 679 known miRNAs and 254 novel miRNAs were identified, of which 32 were detected as miRNAs with significant differential expression (DEMs). 23,577 genes were projected to be the targets of 32 DEMs, primarily those involved in melanogenesis, adrenergic signaling in cardiomyocytes, MAPK signaling pathway and wnt signaling pathway by functional enrichment. Furthermore, we built an interaction module of mRNAs, proteins and miRNAs based on 10 up-regulated and 13 down-regulated miRNAs in golden skin. In addition to transcriptional destabilization and translational suppression, we discovered that miRNAs and their target genes were expressed in the same trend at both the transcriptional and translational levels. Finally, we discovered that miR-196d could be indirectly implicated in regulating melanocyte synthesis and motility in the skin by targeting to myh7 (myosin-7) gene through the luciferase reporter assay, antagomir silencing in vivo and qRT-PCR techniques.
Our study gives a systematic examination of the miRNA profiles expressed in the skin of GBCrC, assisting in the comprehension of the intricate molecular regulation of body color polymorphism and providing insights for C. auratus breeding research.
MicroRNAs (miRNAs) are small, single-stranded, endogenous non-coding RNA molecules containing 21 to 23 nucleotides . Mature miRNAs are produced from longer primary transcripts that are sheared and processed by a series of nucleases and subsequently assembled into an RNA-induced silencing complex (RISC) that recognizes the 3′-untranslated regions (UTRs) of the target mRNA, thereby inducing cleavage or blocking translation of the gene . Recent studies have demonstrated that miRNA-mRNA interactions are significantly related to the transcriptional and signal transduction events involved in a wide range of biological processes, including apoptosis, cell proliferation, embryo development and skin pigmentation [3,4,5]. Via Illumina sequencing, Wang et al.  found that miR-107b and miR-141-3p were up-regulated in pink/black and pink/red skin of red tilapia (Oreochromis spp.) compared to whole pink skin. Luo et al.  discovered 164 differentially expressed miRNAs (DEMs) between the three skin colors (black, red and white) of Koi carp (Cyprinus carpio L.) such as miR-196a, miR-202 and miR-206. Hao et al.  identified 60 DEMs between the red-colored skin and black-colored of leopard coral grouper (Plectropomus leopardus). Meanwhile, Yin et al.  investigated that miR-196a regulates the skin pigmentation of koi carp by targeting transcription factor mitfa (melanocyte inducing transcription factor a). Other research also found that miR-196 family played a critical role in eye formation , fin bud initiation and regeneration  during the early development stages. These many pieces of information point to a potential role for miRNAs in the regulation of fish skin color.
Fish skin and scales frequently exhibit fascinating color patterns, and these patterns play crucial roles in a variety of biological processes, such as mating selection, camouflage, thermoregulation, photoprotection, mimicry and warning off predators . Fish skin color is determined by a complicated process involving a series of cellular, genetic, environmental, nutritional and physiological factors . As of now, many genes associated with the melanogenesis, the wnt signaling pathway and the MAPK signaling pathway, such as oca2 (oculocutaneous albinism II), mitf, map2k4 (Mitogen-Activated Protein Kinase Kinase 4), α-msh (α-melanocyte stimulating hormone) and foxd3 (forkhead box d3), have been identified and are thought to be involved in the regulation of skin color [5, 14, 15]. But undoubtedly, the majority of the gene resources involved in determining fish skin color patterns are still under study.
Crucian carp (C. auratus), one of the most well-liked freshwater species in China with highly economic values, has undergone extensive artificial and natural selections throughout its evolution, and was able to adapt to a variety of environments after being bred into numerous strains for commercial production [16, 17]. The golden-back crucian carp (GBCrC) are distributed in the Dong’s Rice Fish Duck System (DRFDS), one of the Globally Important Agricultural Heritage Systems (GIAHS) in the world, which is located in Qiandongnan Miao and Dong Autonomous Prefecture, Guizhou Province, China. Due to its golden scales and skin on the dorsum, this fish was easier to observe in the paddy field. It has tender meat, soft bones and spines, and its appetizing soup is devoid of muck and fishy flavor. All of these make it worthy of scientific and reasonable breeding and exploitation. In our previous study, we investigated and examined its external morphological traits, distribution, paddy field selection, reproductive mode, genetic diversity and origin. The results showed that GBCrC has a gynogenesis reproductive mode and that the golden coloration of scales and skin is a mutation of native C. auratus that live under special conditions of altitude and annual temperature accumulation, and that this golden color does not disappear under other environmental conditions or through aquaculture methods . Meanwhile, we also identified mRNAs and proteins that are closely associated with two different skin colors (golden and greenish grey) (data not yet published). However, there is still a gap in our knowledge of the genetics in the two skin differentiation types as well as the molecular regulatory mechanisms, which calls for further research employing advanced techniques.
Herein, we further identified skin color related miRNAs using the RNA-Seq approach and screened the DEMs between golden and greenish grey skin, and initially investigated how miR-196d contributed to the skin color differentiation and pigmentation in GBCrC. The findings will provide fundamental knowledge for understanding the genetic mechanisms underlying variation in color traits and may enable the exploitation of gene regulation in skin color determination.
Overview miRNA expression signature in crucian carp skin
To identify miRNAs expressed in GBCrC skin, GO and GR skin small RNA libraries were analyzed by deep sequencing. A total of 11,555,738 and 10,941,036 raw reads were yielded from the two libraries, respectively (Table s1). After filtering out adaptor sequences, low quality reads, insert-null reads, and reads < 18 nt, the reads were further matched with databases for Repeat-associated RNA, GenBank, and Rfam (e.g., rRNA, tRNA, snRNA and snoRNA). Finally, 5,133,894 and 4,002,749 valid reads were retrieved for analysis (Table s1). Of these, 91.28% and 90.29% of the small RNA sequences obtained were 21–24 nt in size, which is the typical size range for Dicer-derived products (Fig. 1).
To further identify the conserved miRNAs and predicted novel miRNAs in the skin of GBCrC, blast search with an E-value cutoff of 10 was used to search for the remnant small RNAs with predicted hairpin structures against the miRbase 22.0. Following alignment and subsequent sequence analysis, 679 miRNAs were discovered to be identical to the mature miRNAs (known), 505 conservative miRNAs were found to match known precursor miRNAs of particular species and 254 miRNAs were projected to be novel miRNAs with hairpin structure (Table s2). Additionally, we counted 15 miRNAs (five highly expressed in both the golden skin and the greenish grey skin, five highly expressed in the golden skin only and five highly expressed in the greenish grey skin only) to better understand the miRNAs that are highly expressed in the skin of GBCrC (Table 1).
Differential expression of miRNAs in different skin
To gain insight into the function of miRNAs in pigmentation processes, it is essential to have precise information on their expression patterns. With the criteria of p-value ≤ 0.05, we identified 32 significant DEMs in GO (including 12 up-regulated and 20 down-regulated miRNAs) compared with GR samples (Table s3). By drawing a volcano plot (Fig. 2A), we were able to better understand the general distribution of DEMs and discovered that certain miRNAs with point overlap due to their relatively close p values (e.g., gmo-miR-126-3p and dre-miR-126a-3p, gmo-let-7a-3-3p and hsa-let-7a-3p, etc.). Six novel predicted miRNAs were discovered among these miRNAs, three of which were strongly up-regulated in golden skin (PC-3p-28150, dre-mir-196d-p3 and ssa-miR-196a), and three of which were highly expressed in greenish grey skin (PC-3p-33311, PC-3p-29800, and PC-3p-24422) (Fig. 2B). Following that, we used qRT-PCR to confirm the expression patterns in 12 DEMs (including 6 down-regulated miRNAs and 6 up-regulated miRNAs). Our qRT-PCR expression patterns matched the results of RNA-seq analysis, indicating that the validity of the RNA-seq data (Fig. 3 and Table s4). For example, based on the small RNA sequencing results, the expression level of dre-miR-196d in GO was almost 1.894 times higher than in GR, compared with a 3.168-fold difference in the qRT-PCR results.
Target genes prediction and enrichment analysis
Target prediction analysis was performed using the miRanda and TargetScan algorithms to further explore the potential functions of DEMs in various colored skin. The outcomes demonstrated that 46,011 transcripts and 23,577 genes could be targeted by 32 miRNAs with considerably differential expression. The miRNAs might target one or more transcripts, and some target transcripts might be regulated by multiple miRNAs (Table s5). For instance, ssa-miR-196a, ccr-miR-139, PC-5p-18590 and gmo-miR-489 can target 4402, 5443, 8420 and 5593 transcripts, respectively. Among these, genes like styk1 (serine/threonine/tyrosine kinase 1), fgfr1 (fibroblast growth factor receptor 1), celf2 (CUGBP elav-like family member 2) and slco4a1 (solute carrier organic anion transporter family member 4A1) can be regulated by numerous miRNAs simultaneously, indicating a complicated regulatory network of these DEMs in different skin differentiation and pigmentation (Fig. 4).
Based on the functional annotation of the target transcripts, enrichment analysis (GO and KEGG) was carried out, and a number of biological processes and pathways were noticeably identified (Table s6). Intriguingly, we discovered that miRNAs can be markedly enriched to biological processes such as protein phosphorylation (GO:0,006,468), regulation of transcription, DNA-templated (GO:0,006,355) and ion transport (GO:0,006,811), molecular function such as transferase activity (GO:0,016,740), protein kinase activity (GO:0,004,672) and metal ion binding (GO:0,046,872), and cellular component such as nucleus (GO:0,005,634) and cytoplasm (GO:0,005,737) (Fig. 5A). Additionally, some KEGG pathways, such as Wnt signaling network, MAPK signaling pathway, melanogenesis and adrenergic signaling in cardiomyocytes, appear to be involved in the color differentiation and pigmentation of the C.auratus skin (Fig. 5B).
Comprehensive regulation of miRNAs on their target expression
In order to search for the targets of DEMs and to understand the impacts of post-transcriptional regulation of miRNAs on gene expression in different skin color pigmentations, we further combined transcriptomic and proteomic data to analysis the molecular networks of miRNA targets. For miRNA-mRNA-protein association analysis, data that met the criteria were screened. Then, we determined 52 groups of reciprocal interaction networks matching the criteria and separated them into two main groups based on the 23 up-regulation and 29 down-regulation of miRNAs (Table s7). Meanwhile, we also discovered some features of miRNA-regulated expression, such as: 1) there were six expression patterns of miRNAs-mRNAs-proteins (down-down-down, down-down-up, down-up-down, down-up-up, up-down-down and up-up-up), of which down-up-up and up-down-down regulation ways were relatively high; 2) one miRNA could target and regulate multiple mRNAs simultaneously, such as ssa-miR-196a could target tnnt2 (troponin T, cardiac muscle isoforms-like), rtn2 (reticulon-2-like), col6a6 (collagen alpha-6 (VI) chain-like) and myoz1 (myozenin-1-like) (Fig. 6); 3) one mRNA or protein could be regulated by numerous miRNAs, such as the scarb2 (lysosome membrane protein 2-like) gene and LIMP-2 (lysosomal integral membrane protein type 2) protein are simultaneously regulated by PC-3p-29800, PC-3p-28510, gmo-let-7a, has-let-7a and ssa-let-7e (Fig. 6). As a result, miRNAs exhibited a more complicated regulation of their target expression, indicated by both repression and activation at the transcriptional and translational levels.
Then, we conducted GO annotation and KEGG enrichment analysis on these highly differentially expressed miRNA-mRNA-protein interactions (Table s7). Several significant biology processes (e.g., sarcomere organization, ion transport, and calcium ion transport), cellular component (e.g., membrane, integral component of membrane, and cytoplasm), molecular function (e.g., actin binding, actin filament binding, and protein binding) (Fig. s1) and pathways (e.g., MAPK signaling pathway, focal adhesion, and adrenergic signaling in cardiomyocytes) (Fig. s2) were enriched.
The role of miR-196d involved in the skin pigmentation
In our previous studies, we conducted in vivo investigations on miRNA function using the miRNA antagomir approach. Here, we also discovered that miR-196d expression in fish could be efficiently suppressed by its antagomir at a level of 40 mg/kg through the tail vein injection, but not by the mismatched antagomir (Fig. s3). Then, we used miR-196d antagomir (positive control), negative antagomir (negative control), and PBS (blank control) to treat wild-type GBCrC. We discovered that the GBCrC skin had a higher melanin content after miR-196d antagomir injection relative to its negative controls (Fig. 7A). In addition, we also examined the tissue-specific expression profiles of miR-196d, and discovered that it was strongly expressed in the gonad and brain, followed by golden skin, heart and fin, but not in the liver, intestine or gill (Fig. s4). These findings suggested that miR-196d expression levels might have an impact on fish skin pigmentation.
To elucidate the role of miR-196d in the skin pigmentation process, we further explored it in conjunction with data from sRNA sequencing and identified the myh7 as a potential target gene in the adrenergic signaling pathway. The alignment of miR-196d with myh7 3′-UTR is illustrated in Fig. 7B based on the sequence complementary. Additionally, we used two luciferase reporters, which were the wild-type 3′-UTR and the mutant 3′-UTR of myh7. These luciferase reporters were co-transfected into HEK293T cells together with a miR-196d mimic. To regulate the non-specific effects of expression, a scrambled miRNA mimic that has no similarity to the C.auratus genome was utilized. The findings showed that the miR-196d mimic dramatically reduced the luciferase activity of the wild-type myh7 3′-UTR but had no effect on the luciferase activity of the mutant 3′-UTR, indicating that miR-196d directly inhibited myh7 expression by binding to its 3′-UTR regions (Fig. 7C).
Meanwhile, we further detected the expression levels of myh7 and tpm1 (tropomyosin alpha-1 chain isoform X16) genes, which are related to the adrenergic signaling pathway. The findings indicated that the expression levels in the miR-196d antagomir injection group were all significantly higher than those in the controls within 1 day, 5 days, and 10 days (p < 0.05), suggesting that miR-196d silencing can affect the expression of genes involved in the adrenergic signaling pathway and regulate the production of melanin (Fig. 8).
In fish, the formation of body color is caused by the light absorption of biological pigments or by the interference of chromatophores with light-reflecting substances, whereas color patterns are produced by the combination and distribution of different pigment cell types [19, 20]. In a previous study, we identified the distinct phenotype of GBCrC and suggested that two golden regions (one is from the 3rd squama on the dorsum to the origin of dorsal fin, and another is from the base of 6th-10th fin ray of the dorsal fin to the base of caudal fin) might result from special altitudes and annual accumulated temperatures . In order to better understand the molecular roles underlying the differentiation and pigmentation of different skin colors, the current study sorted the DEMs that were closely related to them and conducted a preliminary study on miR-196d function, with the aim of laying the theoretical groundwork for the subsequent deeper investigation of the regulatory mechanisms.
MiRNAs play pivotal roles in various biological processes by promoting mRNA degradation or inhibiting mRNA translation, and increasing evidences suggest that their dysregulations have an impact on the skin differentiation and pigmentation [21, 22]. Herein, our sequencing results revealed a similar read length distribution of 21–24 nt, which is comparable in size to typical Dicer-derived products . Most of the identified miRNAs including let-7a, miR-199a, miR-214, miR-1, etc., were also shown to be abundantly expressed in other species, such as koi carp , rainbow trout (Oncorhynchus mykiss)  and Japanese flounder (Paralichthys olivaceus) . Similarly, miR-138-5p and let-7a were also predicted to be highly expressed in the miRNA sequencing investigation of the diverse body colors of red tilapia  and Manila clam (Ruditapes philippinarum) . Studies have shown that miR-33b, which specifically targets hypoxia inducible factor (hif)-1α, can reduce malignant melanoma cell proliferation and glycolysis . Meanwhile, it has been demonstrated that miRNA-139 could bind to its target gene igf1r (insulin-like growth factor receptor type 1) to regulate the PI3K/AKT signaling pathway, which in turn controls the proliferation and spread of malignant melanoma cells . As a result, it is likely that our subsequent research will concentrate on how these miRNAs affect the skin-color regulation of C. auratus.
MiRNAs paly a variety of biological functions, including targeting transcripts and suppressing post-transcriptional gene expression . By target prediction and examination of the correlation between their expression levels, we further investigated the potential relationship between miRNAs and genes. The results showed the complexity of the miRNA-mRNA network and the potential for a single miRNA to target multiple genes and for numerous miRNAs to target the same gene. Similar findings in other research have suggested that particular miRNAs might be involved in the gene regulation of numerous physiologically active processes, and this observation supports the idea of functional redundancy among miRNAs [25, 29]. MiRNAs in various cell types form a comprehensive, multi-layered network system through interactions with signaling pathways and regulatory elements .
Additionally, we performed enrichment analysis on the targeted genes to better understand the molecular function. Among these, GO function analysis is frequently used to explain the molecular function, cellular component and associated biological process. We discovered that several GO terms, including protein phosphorylation (GO:0,006,468), transcription regulation and DNA-templated (GO:0,006,355), nucleus (GO: 0,005,634) and cytoplasm (GO: 0,005,737), can be significantly enriched in biological processes (Fig. 5A), suggesting that the various biological functions underlying the color variations in C.auratus skin. Furthermore, the transcription factor binding (GO: 0,008,134), which is crucial to controlling the transcription of proteins in cells, primarily reflects the enrichment difference between two skin colors. The MAPK, PI3K/Akt, Wnt/β-catenin signaling pathway and melanogenesis pathway are the important regulatory pathways that promote or inhibit the melanin synthesis . KEGG pathway analysis showed that several different pathways for DEMs to target genes were classified, including the adrenergic signaling in cardiomyocyte (ko04261), MAPK signaling pathway (ko04010), wnt signaling pathway (ko04310) and melanogenesis (ko04916), which might be involved in pigment formation and regulation (Table s6). The MAPK signaling pathway, which involved 630 genes, was the main one among them. Previous research has shown that members of the MAPK family proteins, including p38, ERK and JNK, are essential for melanogenesis and p38 also activates MITF, which can increase the expression of melanogenic enzymes . While through suppressing MITF, the ERK and/or JNK/SAPK pathways cause the down-regulation of melanin synthesis . Additionally, the wnt signaling pathway has been shown to be involved in regulating pigment cell differentiation in zebrafish (Danio rerio) , and the relationship between Wnt/β-catenin signaling and MITF has been reported as a key feature of melanocyte development and subsequent pigmentation . It is therefore obvious that different genes in each of these pathways may be involved in regulating skin production in various color phenotypes, but the precise mechanism of action requires further research.
In general, the type of complementarity between miRNAs and their target genes dictates the regulatory interaction, which typically involves two post-transcriptional regulatory mechanisms, mRNA cleavage or translational repression, both of which negatively influence the expression of target genes . However, the mRNA expression profile and protein level of the same genes can be different due to various factors in regulatory mechanisms (e.g., RNA secondary structure, the effects of regulatory proteins, and codon bias) . In our study, we found there were six expression patterns of miRNAs-mRNAs-proteins (down-down-down, down-down-up, down-up-down, down-up-up, up-down-down, and up-up-up). What is more, the expression direction of some miRNAs, their target genes and proteins (e.g., miR-196d-myh7-MYH7, let-7e-scarb2-LIMP-2 and PC-3p-33311-podnl1 (podocan-like 1)-PODNL1) was the same. This might be under two potential mechanisms: 1) counter-regulation of different post-transcriptional mechanisms and 2) the down-regulation of protein production being offset by a feedback mechanism . On the other hand, a small number of miRNAs can regulate protein expression by pairing with proteins to the opposite-expression direction (e.g., miR-142-mlpha (melanophilin a)-MLPH, miR-458-5p-mlpha-MLPH and PC-3p-24422-flnc (filamin C)-FLNC). This is in line with previous study that most animal miRNAs do not precisely complement their targets, leading to their primary mode of action being translational suppression rather than mRNA breakage and destruction . According to our research, it is critical to take into account both aspects of miRNA regulation, including translational repression that results in activation during skin color differentiation and pigmentation.
Pigment cells typically have small muscle fibers and nerve endings distributed on the membrane. The primary cause of the variation in body color in fish is due to these cells changing their morphology when muscle fibers contract, causing the pigment particles to shift under the influence of kinesin, and thus scattering or accumulating in different areas of the body surface . It has been demonstrated that melanocytes had β-adrenergic receptors, and that the adrenaline (AD) was responsible for the diffusion of melanin granules to darken the body color of animals, whereas noradrenaline (NE) was in charge of the aggregation of pigments within melanophores to lighten animals’ bodies, both belonged to the category of neuromodulation of body color changes . Adrenergic receptor (AR) belongs to G protein-coupled receptors, which are receptor proteins for the action of catecholamines and are divided into four subtypes, including α1, α2, β1 and β2 . Among them, α2-AR agonist regulates the production of pigment in melanocytes and erythrophores, includes the aggregation of melanin granules through mediating NE and affects the melanocyte apoptosis; β2-AR stimulates adenylate cyclase (cAMP), a key step comparable to the regulation of melanogenesis pathway, and will mediate AD to inhibit melanocyte aggregation . In this study, we discovered that some miRNAs with differential expression (such as miR-196d, miR-139 and PC-3p-24422) were significantly enriched to key members of the adrenergic signaling pathway, including tpm1, myoz1, myh7 and act3 (actin related gene 3), suggesting that these miRNAs may be involved in regulating the development of various pigment cells. Therefore, we further explored how miR-196d involved in regulating the differentiation and pigmentation process of skin color by targeting the myh7 gene.
Antagomirs are single-stranded RNA molecules (21–23 nt) conjugated to cholesterol that are complementary to a particular miRNA target and contain either a base change or a mispairing at the Ago2 cleavage site to prevent Ago2 cleavage . In our previous studies, we employed a miRNA antagomir technique to carry out in vivo miRNA loss-of-function assays [8, 45]. Herein, we also found that miR-196d expression could be effectively inhibited in vivo by its equivalent antagomir but not the mismatched miRNA. The α-MSH peptide participates in the melanogenesis process by binding to mc1r, raising intracellular cAMP levels, promoting adenylyl cyclase, and activating protein kinase A, which in turn activates tyrosinase to create melanin (Fig. 9A) . Levodopa, an intermediate product of melanin synthesis, is another metabolic pathway that forms catecholamines, which are converted to dopamine, NE and AD by the action of the corresponding transferase enzymes, which can phosphorylate tyrosine and thus regulate tyrosinase activity through the AC/cAMP kinase system and affect melanin content . It has been demonstrated that in the adrenergic signaling pathway, AD binded to β1-AR or β2-AR receptors to activate Gs proteins, which influenced the concentration of the second messenger cAMP by activating adenylate cyclase (AC), and ultimately affected the expression of genes like act3, tpm1 and myh7 through activating the protein kinase A system (PKA) and Ca2+ signaling pathway . In the present study, we found a binding site of miR-196d in the 3′-UTR regions of myh7, and characterized their effects on myh7 using a reporter assay. Meanwhile, we also discovered that the miR-196d antagonist injection group had higher levels of myh7 and tpm1 gene expression than control groups, implying that miRNA-196d, which targeted the myh7 gene in the β2AR-GS-AC-cAMP-PKA signal transduction pathway, may indirectly regulate the synthesis of melanin in crucian carp (Fig. 9B).
In summary, the differentiation and pigmentation of C.auratus skin is a complex event involving several genes and pathways that interact cross-talk. We provided evidence that miRNAs emerged as crucial roles in these processes, and comprehending the regulatory mechanism at the post-transcription level will offer fresh perspectives on the regulation of skin color development and formation.
Our study provides basic information on miRNAs expression in GBCrC, including knowledge of known and potential novel miRNAs, as well as understanding of DEMs related to different skin colors and their associated target genes. We built an interaction module of mRNAs, proteins and miRNAs based on 10 up-regulated and 13 down-regulated miRNAs in golden skin, and found that miR-196d, which targeted the myh7 gene may be indirectly implicated in regulating skin melanocyte synthesis and motility. The findings would provide to a better understanding of the complex molecular regulation of C.auratus body color polymorphism.
One hundred eight-month-old female GBCrC (average weight: 40 ± 3 g) were collected from paddy fields in Pingzheng village, which is a part of the DRFDS and is located at 25°38′N-108°39′E in Congjiang County, Qiandongnan Miao and Dong Autonomous Prefecture, Guizhou, China. They were kept in a water circulation system in 200-L tanks under 12-h light/dark photoperiod at 24 ± 1 °C. Then, three fish were tranquilized in 15 ~ 20 mg/L MS-222 buffered to pH 7.0 ~ 7.5, and golden skin (GO) and greenish grey skin (GR) were collected for small RNA library construction. Meanwhile, the other three fish samples (GO, GR, gonad, gill, liver, muscle, blood, heart, caudal fin, spleen, intestine, eye and brain) were collected for functional analysis. All samples were immediately snap-frozen in liquid nitrogen and stored at -80℃ until use.
Small RNA library construction and sequencing analysis
RNA samples were obtained from different colored skin of GBCrC using TRIzol (Invitrogen, USA) according to the manufacturer’s protocol. The purity and quantity were analyzed by RNA Nano LabChip (Agilent, USA) and Bioanalyzer 2100 with RIN value > 7.0, then kept in -80℃ refrigerator. Small RNA libraries were constructed following the protocol of Illumina’s TruSeq Small RNA Sample Preparation Kits (San Diego, CA, USA). After that, single-end sequencing (50 bp) was run followed the manufacturer’s instructions on an Illumina Hiseq 2500 at the LC-BIO (Hangzhou, China).
Raw reads were subjected to ACGT101-miR program (LC Sciences, Houston, USA)  to remove adapter dimers, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats. Subsequently, unique sequences (18 ~ 26 nucleotide) were mapped to the C.auratus genome (https://www.ncbi.nlm.nih.gov/genome/?term=goldfish) by blast with a tolerance of one mismatch. The mapped sequences were separated into two categories: known miRNAs, which are found at the hairpin arms, and novel miRNAs, which are found at the arm that is opposite the annotated arm. After that, the remaining sequences were mapped to other fish species precursors in miRBase 22.0 by BLAST search, and the mapped pre-miRNAs were further blasted against the C.auratus genome to determine their genomic locations. To identify the novel predicted miRNAs, the unmapped sequences were blasted against the C.auratus genome database, and the hairpin RNA structures comprising sequences were predicated by RNAfold software (http://rna.tbi.univie.ac. at/cgi-bin/RNAfold.cgi).
Differentially expressed miRNAs and target genes prediction
The methods used for data normalization were those outlined in a prior study . miRNAs were regarded as differentially expressed based on normalized deep-sequencing levels in greenish grey and golden skin. The p-value was calculated by using Chi-square (X2) test and Fisher exact test. The significance level for each test was set at 0.05. Those miRNAs with p-value ≤ 0.05 were considered as DEMs.
Two computational target prediction algorithms, TargetScan  and Miranda , were utilized to locate miRNA binding sites in order to anticipate the genes that the most prevalent miRNAs will target. The overlaps were then determined using the combined data predicted by the two methods. Gene Ontology (GO; http://www.geneontology.org) biological process categories were used to identify functions strongly associated with the predicted target genes of miRNAs  and KEGG enrichment analysis was carried out using the pathway database (http://www.genome.jp/kegg/pathway.html) in a statistical test on each DEMs [53, 54], and the results were plotted using R program v3.0.1 (https://cran-archive.r-project.org/bin/windows/base/old/3.0.1/).
Based on the results of the previous transcriptome and proteome sequencing of golden back crucian carp with GO and GR samples (data not yet published), we performed a triple association analysis for differential mRNAs satisfying FC > 2 or FC < 0.5 with p < 0.05, differential proteins satisfying FC > 1.2 or FC < 0.833 with p < 0.05, and differential miRNAs satisfying p < 0.05. For analysis, Cytoscape v2.8.1 (https://cytoscape.org/download.html) plotted and networked the outcome files.
Silencing miR-196d in vivo with the antagomir
The antagomirs employed are single-stranded RNA that consists of 21–23 nucleotides by chemically modified as follows: the miR-196d antagomir is CsCsCCACAACACGAAACUGUUAAsCsCsAs-Chol-3′. Additionally, we set up miRNA negative antagomir (CsCsACACACACUUCCUUACAUUsCsCsAs-Chol-3′) and phosphate-buffered saline (PBS) groups for the experiment’s control, respectively. The subscript “s” represents for phosphorothioate linkage and “Chol” stands for cholesterol connected via hydroxyproline in chemical modifications of miRNA antagomir, and all nucleotides are 2’-OMe-modified. Following that, 60 GBCrCs that were being temporarily raised in 200 l tanks were evenly divided into three groups that subjected to the tail vein injection of miR-196d antagomir, negative antagomir and PBS, respectively. Golden skin samples were taken at 0-, 1-, 5-, 10-, and 14-day after injection, then instantly frozen in liquid nitrogen and kept at -80 °C.
Luciferase reporter assay
Myh7’s 3′-UTR was synthesized based on the crucian carp sequence (XM_026241629.1) and was then separately cloned into pSE3575 vector (Sunbio Medical Biotech Co., Ltd., China) by directional cloning. By altering the seed region of the predicted miR-196d site, the mutant myh7 3′-UTR reporters were produced. HEK293 cells, which do not express miR-196d, were seeded into 48-well plates at a density of 1 × 105 cells per well the day before transfection. Then, using lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA), 25 ng of the luciferase reporter vector carrying 50 nM miR-196d mimic, miR-196d wild-type (wt), or the 3′-UTR mutant was co-transfected with 5 ng of the Renilla luciferase control vector (pRL-TK, Promega, Madison, WI, USA) in 24-well plates. Using the Dual Luciferase Reporter Assay System (Promega, Madison, WI, USA), luciferase assays were performed in accordance with the manufacturer’s instructions at 48 h after transfection. A liquid scintillation counter was used to detect luciferase activity, the firefly luciferase activity was normalized to Renilla luciferase activity.
Quantitative Real-Time PCR (qRT-PCR) analysis
Total RNA extraction procedures followed the same procedures as those above mentioned for small RNA libraries construction. The step-by-step process for qRT-PCR and reverse transcription is consistent with our earlier study (Dong et al., 2020). All the primers (Table s8) were designed using Primer Premier v5.0 and synthesized by Sangon Biotech., China. U6 snRNA and β-actin were employed as the internal controls of the DEM and target genes, respectively, in three biological replications of each reaction. The relative expression values were determined based on the 2-△△Ct method .
The ELISA kit (Zhenke, Shanghai, China) was used to quantify the activities of melanin (ME), and the precise steps were the same as those described in the article by Wang et al. . All results were presented as means ± standard deviation, and SPSS v21.0. (SPSS Inc., IL, USA) was used to analyze the data. By using a one-way ANOVA and post hoc Duncan’s multiple range tests, the sample points of various treatments were examined. P < 0.05 was deemed to be significant.
Availability of data and materials
The original contributions presented in the study are publicly available. All raw sRNA sequencing data have been deposited in the NCBI Gene Expression Omnibus (GEO) with the accession numbers GSE212497.
National Center for Biotechnology Information
- Mitf :
Microphthalmia-associated transcription factor
- Mc1r :
Melanocortin receptor 1
- Msh :
RNA-inducing silencing complex
Real-time quantitative reverse transcription polymerase chain reaction
Freshwater Fisheries Research Center
Differentially expressed miRNAs
Kyoto Encyclopedia of Genes and Genomes
- Wnt/β-catenin signaling:
Wingless and INT-1 signaling transduction
Mitogen-activated protein kinase
Cyclic adenosine monophosphate
Cai Y, Yu X, Hu S, Yu J. A brief review on the mechanisms of miRNA regulation. Genom Proteom Bioinf. 2009;7:147–54.
Fabian MR, Sonenberg N, Filipowicz W. Regulation of mRNA translation and stability by microRNAs. Annu Rev Biochem. 2010;79:351–79.
Di Leva G, Croce CM. miRNA profiling of cancer. Curr Opin Genet Dev. 2013;23:3–11.
Wang F, Jia Y, Wang P, Yang Q, Du Q, Chang Z. Identification and profiling of Cyprinus carpio microRNAs during ovary differentiation by deep sequencing. BMC Genom. 2017;18:333.
Yan B, Liu B, Zhu CD, Li KL, Yue LJ, Zhao JL, et al. microRNA regulation of skin pigmentation in fish. J Cell Sci. 2013;126:3401–8.
Wang L, Zhu W, Dong Z, Song F, Dong J, Fu J. Comparative microRNA-seq analysis depicts candidate miRNAs involved in skin color differentiation in red tilapia. Int J Mol Sci. 2018;19:1209.
Luo M, Wang L, Zhu W, Fu J, Song F, Fang M, et al. Identification and characterization of skin color microRNAs in koi carp (Cyprinus carpio L.) by Illumina sequencing. BMC Genom. 2018;19:779.
Hao R, Zhu X, Tian C, Jiang M, Huang Y, Zhu C. Integrated analysis of the role of miRNA-mRNA in determining different body colors of leopard coral grouper (Plectropomus leopardus). Aquaculture. 2022;548: 737575.
Yin H, Luo M, Luo W, Wang L, Zhu W, Fu J, et al. miR-196a regulates the skin pigmentation of koi carp (Cyprinus carpio L.) by targeting transcription factor mitfa. Aquac Res. 2021;52:229–36.
Qiu R, Liu Y, Wu JY, Liu K, Mo W, He R. Misexpression of miR-196a induces eye anomaly in Xenopus laevis. Brain Res Bull. 2009;79:26–31.
He X, Yan YL, Eberhart JK, Herpin A, Wagner TU, Schartl M, et al. MiR-196 regulates axial patterning and pectoral appendage initiation. Dev Biol. 2011;357:463–77.
Protas ME, Patel NH. Evolution of coloration patterns. Annu Rev Cell Dev Biol. 2008;24:425–46.
Nüsslein-Volhard C, Singh AP. How fish color their skin: a paradigm for development and evolution of adult patterns: multipotency, plasticity, and cell competition regulate proliferation and spreading of pigment cells in zebrafish coloration. BioEssays. 2017;39:1600231.
Klaassen H, Wang Y, Adamski K, Rohner N, Kowalko JE. CRISPR mutagenesis confirms the role of oca2 in melanin pigmentation in Astyanax mexicanus. Dev Biol. 2018;441:313–8.
Liu JH, Wen S, Luo C, Zhang YQ, Tao M, Wang DW, et al. Involvement of the mitfa gene in the development of pigment cell in Japanese ornamental Koi carp (Cyprinus carpio L). Genet Mol Res. 2015;14:2775–84.
Zhou L, Gui J. Natural and artificial polyploids in aquaculture. Aquac Fish. 2017;2:103–11.
Gui JF, Zhou L, Li XY. Rethinking fish biology and biotechnologies in the challenge era for burgeoning genome resources and strengthening food security. Water Biology and Security. 2021. https://doi.org/10.1016/j.watbs.2021.11.001.
Zhang X, Fu J, Luo M, Cai-rang Z, Min Q, Hu J, et al. Genetic and reproductive mode analyses of a golden-back mutant of crucian carp from a rice-fish integrated farming system. Aquac Rep. 2022;24: 101146.
Parichy DM. Pigment patterns: fish in stripes and spots. Curr Biol. 2003;13:R947–50.
Olsson M, Stuart-Fox D, Ballen C. Genetics and evolution of colour patterns in reptiles. Seminar Cell Dev Biol. 2013;24:529–41.
Yi R, Fuchs E. MicroRNA-mediated control in the skin. Cell Death Differ. 2010;17:229–35.
Li X, Ponandai-Srinivasan S, Nandakumar KS, Fabre S, Xu Landén N, Mavon A, et al. Targeting microRNA for improved skin health. Health Sci Rep. 2021;4:e374.
Starega-Roslan J, Krol J, Koscianska E, Kozlowski P, Szlachcic WJ, Sobczak K, et al. Structural basis of microRNA length variety. Nucleic Acids Res. 2011;39:257–68.
Juanchich A, Bardou P, Rué O, Gabillard JC, Gaspin C, Bobe J, et al. Characterization of an extensive rainbow trout miRNA transcriptome by next generation sequencing. BMC Genom. 2016;17:164.
Wang N, Wang R, Wang R, Tian Y, Shao C, Jia X, et al. The integrated analysis of RNA-seq and microRNA-seq depicts miRNA-mRNA networks involved in Japanese flounder (Paralichthys olivaceus) albinism. PLoS ONE. 2017;12: e0181761.
Ding J, Wen Q, Huo Z, Nie H, Qin Y, Yan X. Identification of shell-color-related microRNAs in the Manila clam Ruditapes philippinarum using high-throughput sequencing of small RNA transcriptomes. Sci Rep. 2021;11:8044.
Zhao Y, Wu C, Li L. MicroRNA-33b inhibits cell proliferation and glycolysis by targeting hypoxia-inducible factor-1α in malignant melanoma. Exp Ther Med. 2021;22:1924.
Yang C, Xia Z, Zhu L, Li Y, Zheng Z, Liang J, et al. MicroRNA-139-5p modulates the growth and metastasis of malignant melanoma cells via the PI3K/AKT signaling pathway by binding to IGF1R. Cell Cycle. 2019;18:3513–24.
Alberti C, Cochella L. A framework for understanding the roles of miRNAs in animal development. Development. 2017;144:2548–59.
Horsburgh S, Fullard N, Roger M, Degnan A, Todryk S, Przyborski S, et al. MicroRNAs in the skin: role in development, homoeostasis and regeneration. Clin Sci. 2017;131:1923–40.
Luo M, Lu G, Yin H, Atuganile M, Dong Z. Fish pigmentation and coloration: Molecular mechanisms and aquaculture perspectives. Rev Aquac. 2021;13:2395–412.
Bellei B, Maresca V, Flori E, Pitisci A, Larue L, Picardo M. P38 regulates pigmentation via proteasomal degradation of tyrosinae. J Biol Chem. 2010;285:7288–99.
Kim DS, Jeong YM, Park IK, Hahn HG, Lee HK, Kwon SB, et al. A new 2-imino-1,3-thiazoline derivative, KHG22394, inhibits melanin synthesis in mouse B16 melanoma cells. Biol Pharm Bull. 2007;30:180–3.
Dorsky RI, Raible DW, Moon RT. Direct regulation of narce, a zebrafish MITF homolog required for pigment cell formation, by the Wnt pathway. Genes Dev. 2000;14:158–62.
Gajos-Michniewicz A, Czyz M. WNT signaling in melanoma. Int J Mol Sci. 2000;21:4852.
Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136:215–33.
Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, et al. Global quantification of mammalian gene expression control. Nature. 2011;473:337–42.
Bartel DP, Chen CZ. Micromanagers of gene expression: the potentially widespread influence of metazoan microRNAs. Nat Rev Genet. 2004;5:396–400.
Cloonan N. Re-thinking miRNA-mRNA interactions: Intertwining issues confound target discovery. BioEssays. 2015;37:379–88.
Nascimento AA, Roland JT, Gelfand VI. Pigment cells: a model for the study of organelle transport. Annu Rev Cell Dev Biol. 2003;19:469–91.
Sivamani RK, Pullar CE, Manabat-Hidalgo CG, Rocke DM, Carlsen RC, Greenhalgh DG, et al. Stress-mediated increases in systemic and local epinephrine impair skin wound healing: potential new indication for beta blockers. PLoS Med. 2009;6: e12.
Li F, De Godoy M, Rattan S. Role of adenylate and guanylate cyclases in beta1-, beta2-, and beta3-adrenoceptor-mediated relaxation of internal anal sphincter smooth muscle. J Pharmacol Exp Ther. 2004;308:1111–20.
Fujii R. The regulation of motile activity in fish chromatophores. Pigment Cell Res. 2000;13:300–19.
Krützfeldt J, Rajewsky N, Braich R, Rajeev KG, Tuschl T, Manoharan M, et al. Silencing of microRNAs in vivo with ‘antagomirs.’ Nature. 2005;438:685–9.
Dong Z, Luo M, Wang L, Yin H, Zhu W, Fu J. MicroRNA-206 regulation of skin pigmentation in koi carp (Cyprinus carpio L.). Front Genet. 2020;11:47.
Yamanome T, Chiba H, Takahashi A. Melanocyte-stimulating hormone facilitates hypermelanosis on the non-eyed side of the barfin flouder, a pleuronectiform fish. Aquaculture. 2007;270:505–11.
Zhang X, Chien EY, Chalmers MJ, Pascal BD, Gatchalian J, Stevens RC, et al. Dynamics of the β2-adrenergic G-protein coupled receptor revealed by hydrogen-deuterium exchange. Anal Chem. 2010;82:1100–8.
Han X, Yin H, Song X, Zhang Y, Liu M, Sang J, et al. Integration of small RNAs, degradome and transcriptome sequencing in hyperaccumulator Sedum alfredii uncovers a complex regulatory network and provides insights into cadmium phytoremediation. Plant Biotechnol J. 2016;14:1470–83.
Li X, Shahid MQ, Wu J, Wang L, Liu X, Lu Y. Comparative small RNA analysis of pollen development in autotetraploid and diploid rice. Int J Mol Sci. 2016;17:499.
Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. Elife. 2015;4: e05005.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5:R1. https://doi.org/10.1186/gb-2003-5-1-r1.
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.
Kanehisa M. Torward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28:1947–51.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51:D587–92.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-DDCT method. Meth. 2001;25:402–8.
We would like to thank Hangzhou LC-Bio Technology CO., Ltd., (Hangzhou, China) for providing transcriptomics and proteomics services. Also, we want to express our gratitude to Dr. Yalun Dong of the University of the Sunshine Coast for his help with the manuscript and adjustments to the grammar. Additionally, the relevant pathway maps in Figure 9 are referenced from the KEGG database and are hereby acknowledged.
This work was supported by grant  05 from the Germplasm Resources Project of Guizhou Academy of Agriculture Sciences; grants CARS-45 and CARS-46 from the Earmarked Fund for China Agriculture Research System.
Ethics approval and consent to participate
The fish samples used in this study were obtained from Guizhou Fisheries Research Institute in China with permission from the owner. All the experimental procedures were reviewed and approved by the Bioethical Committee of the Freshwater Fisheries Research Center of the Chinese Academy of Fishery Sciences.
The animals used in this study were handled in accordance with relevant guidelines and regulations of the Ministry of Science and Technology, Beijing, China (No. 398, 2006). Also, the study was conducted in accordance with the ARRIVE guidelines (https://arriveguidelines.org) for the reporting of animal experiments.
Consent for publication
The authors have declared that no competing interest exists.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Gene ontology (GO) enrichment analysis of mRNAs targeted by miRNAs that were significantly differentially expressed in the GO and GR groups.
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of the targets of the total differentially expressed miRNAs (DEMs). Gene number, number of target genes in each pathway.
Effect of antagomir treatment on miR-196d expression. Fish was injected with miR-206 antagomir, miR-196d antagomir or left untreated, respectively. Asterisk (*) indicates a significant difference compared with the control group. Each sample was analyzed in triplicate.
Expression pattern of miR-196d in different tissues. Go, gonad; Br, brain; He, heart; Fi, fin; Gos, golden skin; Grs, greenish grey skin; Sp, spleen; Mu, muscle; In, intestine; Li, liver; Gi, gill; Ey, eye; Bl, blood. * p <0.05, ** p <0.01, *** p <0.001.
Distribution of sequenced reads from raw data to cleaned sequences. Table s2. Summary of known and predicted miRNAs in this study. Table s3. Differentially expressed miRNAs (DEMs) in different skin color of crucian carp. Table s4. The results of qRT-PCR and small RNA sequencing. Table s5. Results of target gene analysis of various miRNAs. Table s6. GO function and KEGG pathway enrichment analysis of differentially expressed miRNAs. Table s7. Information on the correspondence between differentially miRNAs, mRNAs, and proteins. Table s8. Primers for skin DEMs and reference miRNAs in crucian carp.
About this article
Cite this article
Zhang, X., Luo, M., Jiang, B. et al. microRNA regulation of skin pigmentation in golden-back mutant of crucian carp from a rice-fish integrated farming system. BMC Genomics 24, 70 (2023). https://doi.org/10.1186/s12864-023-09168-w
- Skin color
- Golden-back crucian carp