Genetic Characteristic and RNA-Seq Analysis in Transparent Mutant of Carp–Goldfish Nucleocytoplasmic Hybrid

In teleost, pigment in the skin and scales played important roles in various biological processes. Iridophores, one of the main pigment cells in teleost, could produce silver pigments to reflect light. However, the specific mechanism of the formation of silver pigments is still unclear. In our previous study, some transparent mutant individuals were found in the carp–goldfish nucleocytoplasmic hybrid (CyCa hybrid) population. In the present study, using transparent mutants (TM) and wild type (WT) of the CyCa hybrid as a model, firstly, microscopic observations showed that the silver pigments and melanin were both lost in the scales of transparent mutants compared to that in wild types. Secondly, genetic study demonstrated that the transparent trait in the CyCa hybrid was recessively inherent, and controlled by an allele in line with Mendelism. Thirdly, RNA-Seq analysis showed that differential expression genes (DEGs) between wild type and transparent mutants were mainly enriched in the metabolism of guanine, such as hydrolase, guanyl nucleotide binding, guanyl ribonucleotide binding, and GTPase activity. Among the DEGs, purine nucleoside phosphorylase 4a (pnp4a) and endothelin receptor B (ednrb) were more highly expressed in the wild type compared to the transparent mutant (p < 0.05). Finally, miRNA-Seq analysis showed that miRNA-146a and miR-153b were both more highly expressed in the transparent mutant compared to that in wild type (p < 0.05). Interaction analysis between miRNAs and mRNAs indicated that miRNA-146a was associated with six DEGs (MGAT5B, MFAP4, GP2, htt, Sema6b, Obscn) that might be involved in silver pigmentation.


Introduction
In teleost, the pigment chromophores in the skin and scales could be involved in several important biological processes, such as camouflage, thermoregulation, mate choice, mimicry, warning of toxicity, and immunity [1]. The pigments are produced by several pigment cells in fish, including melanophores, xanthophores, and iridophores [2]. Pigmentation is a very complicated process, which is regulated by a series of cellular, genetic, endocrine, environmental, and physiological factors [3]. Different pigment cells in fish display different colors due to their unique organelles. The melanin particles produced by melanocytes could absorb the incident light of a specific wavelength and make fish appear black or brown. The pigments produced by xanthophores are mainly pteridine, which can filter out certain wavelengths of light to make fish appear yellow or red [4]. In addition, iridophores could produce silver pigments, in which guanine could combine with water to form purine-containing reflecting platelets [5].
The transparent trait in fish is closely related to the development of pigment cells. Generally, most fish are transparent in their early lives, while the body in adult fish is no longer transparent due to the accumulation of pigment in skin and peritoneum. The cytoplasm of iridophores contains a large number of crystalline plates, which could reflect a certain wavelength of light and make them iridescent. The color presented by iridophores belongs to the configuration color, which is different from other kinds of pigment cells [5]. The absence of silver pigments produced by iridophores could allow light to penetrate body and make fish transparent. Previous studies have reported that zebrafish shady, rose, and barley transparent mutations are all associated with the lost of iridophores. Mutations in shady could result in a lack of iridophores throughout the whole life [6,7]. Mutations in rose, a gene encoding the endothelin receptor b1a (ednrb1a), could cause an adult-specific reduction of iridophores [8,9]. In addition, the transparent mutant not only shows the lost of iridophores, but also a significant decrease in the number of melanocytes.
To date, a series of genes have been reported to be involved in the determination of skin color, such as pro-opiomelanocortin (POMC), melanocyte stimulating hormone (MSH), microphthalmia-associated transcription factor (MITF), kit oncogene (KIT), tyrosinase (TYR), tyrosine related protein-1 (TYRP1), and tyrosine related protein-2 (TYRP2) [2,10]. In melanocytes, MITF could induce the biosynthesis of melanin through the up-regulation of melanogenic enzyme expression [11]. Further studies found that endothelin receptor B (ednrb) could induce the phosphorylation of MITF and its mutation could reduce the pigmentation in melanocytes [12,13]. In addition, previous studies reported that purine nucleoside phosphorylase 4a (pnp4a) was also involved in guanine synthesis and its mutation could cause the deficiency of iridophores [14,15]. However, most of the genetic resources involved in the color of fish remain to be discovered.
MicroRNAs (miRNAs) are a group of single-stranded, non-coding RNA molecules with an average size of approximately 22 nucleotides, which play crucial roles in numerous biological processes through promoting their target gene's degradation or inhibiting mRNA translation. It is important for the crosstalk between miRNAs and mRNAs to steady the signal transduction and the transcriptional activities as well as maintain homeostasis in many organs [16][17][18]. In mammals, it has been reported that the miR-137 could affect body color pattern in mice [19]. In Drosophila, the loss of miR-8 could reduce pigmentation in dorsal abdomen [20]. In addition, miR-429 silencing in common carp could increase the transcript level of Foxd3, and influence the melanin production [21]. These studies indicate that miRNAs could also play an important role in body color formation [15]. However, little is known about the functional role of miRNA in the regulation of silver pigment synthesis.
In our previous studies, some transparent mutant (TM) individuals were found in the wild type (WT) of the carp-goldfish nucleocytoplasmic hybrid (CyCa hybrid) population [22]. In the WT, the silver pigments and melanin in the scales could block the light to pass through fish body, so they were opaque. However, there were no reflective substances in the skin and peritoneum of TM, so the internal organs of TM could be observed directly. In the present study, using the WT and TM as the model, (i) the hybridization experiment was used to further explore the genetic characteristics of the transparent trait in common carp; (ii) the RNA-Seq was used to examine the differentially expressed genes between the TM and WT; (iii) the miRNA-Seq was also recruited to detect the functional role of miRNA in the formation of pigment in common carp. This study will help us to understand the formation mechanism of the transparent trait in teleost.

Genetic Analysis of the Transparent Trait
Experimental CyCa hybrids were bred in the Yaowan Testing Ground of the Yangtze River Fisheries Research Institute (Hubei Province). The wild types (WTs) in CyCa hybrids (three males and three females) were used as broodstocks for crossing to produce the F1 generation ( Figure 1A). After the F1 generation matured, the transparent mutants (TMs) in the F1 generation (six males and six females) were used as parents for self-crossing to produce the F2 generation ( Figure 1B). In addition, three male TMs in F1 generation and three female WTs in the F0 generation were used as parents for back-crossing to produce the F2 generation, respectively ( Figure 1C). Similarly, three female TMs in the F1 generation and three male WTs in the F0 generation were also recruited to do the back-crossing experiment ( Figure 1D). When the phenotypes were easily distinguished, the number of WTs and TMs in the F1 and F2 generation were counted, respectively. In the present study, the common carps were used according to the protocol approved by the committee for animal use at Huazhong Agricultural University (ethical approval no. HBAC20091138; date: 15 November 2009).

Genetic Analysis of the Transparent Trait
Experimental CyCa hybrids were bred in the Yaowan Testing Ground of the Yangtze River Fisheries Research Institute (Hubei Province). The wild types (WTs) in CyCa hybrids (three males and three females) were used as broodstocks for crossing to produce the F1 generation ( Figure 1A). After the F1 generation matured, the transparent mutants (TMs) in the F1 generation (six males and six females) were used as parents for self-crossing to produce the F2 generation ( Figure 1B). In addition, three male TMs in F1 generation and three female WTs in the F0 generation were used as parents for back-crossing to produce the F2 generation, respectively ( Figure 1C). Similarly, three female TMs in the F1 generation and three male WTs in the F0 generation were also recruited to do the back-crossing experiment ( Figure 1D). When the phenotypes were easily distinguished, the number of WTs and TMs in the F1 and F2 generation were counted, respectively. In the present study, the common carps were used according to the protocol approved by the committee for animal use at Huazhong Agricultural University (ethical approval no. HBAC20091138; date: 15 November 2009). The wild types (WTs) in the carp-goldfish nucleocytoplasmic hybrids (CyCa hybrids) were used as broodstocks for crossing to produce F1 generation. (B) After the F1 generation matured, the transparent mutants (TMs) in the F1 generation were used as parents for self-crossing to produce the F2 generation. (C) TMs in the F1 generation and WTs in the F0 generation were used as the male parent and female parent for back-crossing to produce the F2 generation, respectively. (D) TMs in the F1 generation and WTs in the F0 generation were used as the female parent and male parent for back-crossing to produce the F2 generation, respectively. When the phenotypes were easily distinguished, the number of WTs and TMs in the F1 and F2 generation were counted, respectively.

RNA-Seq Library Construction and Illumina Sequencing
Total RNA was extracted from the whole carp fry of WTs and TMs which were two-weeks old by using a TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. To be specific, five fry were collected for one replicated group, and three replicated groups for TMs (total 15 carps) and WTs (total 15 carps) were all used for RNA extraction, respectively. A Nanodrop 2000 spectrophotometer was used to assess sample purity and RNA The wild types (WTs) in the carp-goldfish nucleocytoplasmic hybrids (CyCa hybrids) were used as broodstocks for crossing to produce F1 generation. (B) After the F1 generation matured, the transparent mutants (TMs) in the F1 generation were used as parents for self-crossing to produce the F2 generation. (C) TMs in the F1 generation and WTs in the F0 generation were used as the male parent and female parent for back-crossing to produce the F2 generation, respectively. (D) TMs in the F1 generation and WTs in the F0 generation were used as the female parent and male parent for back-crossing to produce the F2 generation, respectively. When the phenotypes were easily distinguished, the number of WTs and TMs in the F1 and F2 generation were counted, respectively.

RNA-Seq Library Construction and Illumina Sequencing
Total RNA was extracted from the whole carp fry of WTs and TMs which were two-weeks old by using a TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. To be specific, five fry were collected for one replicated group, and three replicated groups for TMs (total 15 carps) and WTs (total 15 carps) were all used for RNA extraction, respectively. A Nanodrop 2000 spectrophotometer was used to assess sample purity and RNA concentration, and the quality of RNA was analyzed by using an Agilent 2100 bioanalyzer with a RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA). Only high-quality RNA samples (RIN ≥ 8) were used to construct the sequencing library. RNA-Seq transcriptome libraries were prepared by using 5 µg of total RNA following TruSeq TM RNA sample preparation Kit (Illumina, San Diego, CA, USA). After quantification by TBS380, paired-end RNA-Seq sequencing libraries were sequenced with the Illumina HiSeq™ 4000 (Illumnia, San Diego, CA, USA). In this study, the cDNA libraries for TMs and WTs were constructed and sequenced in Nextomics Biosciences (Wuhan, China), respectively. The raw reads were filtered using NGSQCToolkit (v2.3.3). By removing reads containing adapters, reads containing poly-N, or low quality reads from raw data, the clean data were obtained. Using FastQC software to do quality control (QC) of clean data, if QC was qualified, clean reads were compared to the reference sequence (http://www.carpbase.org/CarpBase/download). Then, a series of subsequent analysis such as gene expression, variable shear, prediction of new transcripts, SNP detection, and gene structure optimization were conducted.

Small RNA Library Construction and Illumina Sequencing
A total amount of 10 µg of total RNA from WTs and TMs were used as input material to construct small RNA libraries. The small RNA libraries were constructed by using a Truseq TM Small RNA sample prep Kit (San Diego, CA, USA), following the Illumina protocol. Subsequently, SE50 sequencing of qualified sequencing libraries was performed on an Illumina HiSeq™ system [23]. In this study, three small RNA libraries for WTs and TMs were constructed and sequenced by Nextomics Biosciences (Wuhan, China), respectively. Then, a series of subsequent analysis such as sequencing data quality control, small RNA annotation analysis, and miRNA deep analysis would be carried out.

Differentially Expressed Genes (DEGs) Analysis and Functional Enrichment
To identify DEGs between TMs and WTs, the expression level of each transcript was calculated according to the fragments per kilobase of exon per million mapped reads (FRKM) method. The significance of DEGs between TMs and WTs was determined by using edgeR [24]. In this study, the threshold for DEGs was defined with p < 0.05 and |log2(Fold change)| ≥ 1. Functional annotation of gene ontology (GO) terms were analyzed by using Blast2GO software [25], and GO functional classification of unigenes was analyzed by using WEGO software [26]. Functional enrichment analysis including GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed by using Goatools (or KOBAS) software.

Differentially Expressed miRNAs (DEMs) Analysis and Functional Enrichment
Clean reads were filtered from raw reads by removing low quality reads and adaptor sequences, either shorter than 18 nt or longer than 32 nt. Then, the cleaned sRNA reads were mapped to the Rfam databases to discard rRNA-, snRNA-, scRNA-, snoRNA-, tRNA-, and ribozyme associated reads [27]. Then, the qualified clean reads were compared with the reference sequence by the genome comparison analysis software Bowtie (version 1.1.2) [28]. DEMs analysis was performed by using DEGseq software. The screening conditions for DEMs were p-adjust < 0.05 and |log2(Fold change)| ≥ 1. The potential known and novel miRNA targets were predicted by miRanda software [29], and expression levels of target genes were taken from DEGs identified in our transcriptome data. GO and KEGG analysis were performed for the target genes of DEMs by using Blast2GO and KASS.

Interaction Analysis of DEGs and DEMs
In order to define all the possible miRNA-mRNA interactions, we constructed the DE genes pool (p < 0.05 and |log2(Fold change)| ≥ 1) and DE miRNAs target genes pool (p-adjust < 0.05 and |log2(Fold change)| ≥ 1), respectively. Then the interaction analysis was carried out in the two gene pools. Briefly, the target genes of the DE miRNAs were firstly predicted basing on the miRanda (3.3a). Then, overlapped with the DE genes obtained by RNA-Seq, the target DE genes of DE miRNAs were further detected. According to the expression trend, the relationship between DE miRNAs and their corresponding target DE genes was determined to be up-up, up-down, down-down and down-up, respectively. The miRNA-mRNA regulatory networks were visualized by using Cytoscape software [30]. Given that miRNAs could negatively regulate the expression of their target mRNAs in most cases, only the data with a negative relationship (up-down and down-up) between miRNA and mRNA expression were selected for the miRNA-mRNA network mapping in the present study.

RT-PCR Analysis of DEGs
Total RNA was extracted from the same period carp fry as in the RNAseq experiment by using Trizol reagent (Invitrogen). Five fry were collected for one replicated group, and three replicated groups for TMs (total 15 carps) and WTs (total 15 carps) were all used for RNA extraction, respectively. The cDNA was synthesized by using 5 µg total RNA with the Hifair TM II 1st Strand cDNA Synthesis Kit following the manufacturer's protocol (Yesen, Shanghai). The primers of candidate genes were designed following the general principles with Premier 6.0 software. The housekeeping gene β-actin and glyceraldehyde-3-phosphate dehydrogenase (Gapdh) were used as the reference. The details of the primers used in RT-qPCR are listed in Table S1. RT-PCR was performed by using the Unique Aptamer TM qPCR SYBR ® Green Master Mix (Novogene, Tianjin, China) on an ABI PRISM 7500 Real-time Detection System (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's protocol. The relative expression levels of target genes were calculated by the 2 −∆∆CT method. ∆CT target gene = CT target gene − CT housekeeping gene , ∆∆CT = ∆CT TM − ∆CT WT , each target gene mRNA expression difference multiples with 2 −∆∆CT .

Genetic Analysis of Transparent Trait
In the CyCa hybrid carp inbred line, some individuals were found to present the transparent trait without reflective substances in their eyelids and scales, which were named as transparent mutants (TM). Since the TMs accounted for about 25% in F1 generation ( Figure 1A), it was presumed that the reflective trait in their broodstock was an "Aa" heterozygote, and the transparent trait was recessive over the reflective trait. In addition, all of the self-crossing F2 generation showed the transparent trait ( Figure 1B), while the ratio of TMs and WTs of back-crossing F2 generation was nearly 50% ( Figure 1C,D). These results further confirmed that the transparent trait in TM was "aa" and the heritance characteristic was consistent with Mendel's genetic law. Microscopic examination showed that there were a large number of melanin, yellow pigments, and silver pigments on the scales of WT ( Figure 2A). However, only a large number of yellow pigments, very few melanin and silver pigments, were found on the scales of TM ( Figure 2B). Due to the absence of silver pigments, the scales in TMs were transparent.

Transcriptome Analysis
After stringent quality assessment and data filtering, approximately 63,272,006 and 75,019,864 clean paired-end sequence reads with a length of 145 bp and a Q30 rate (those with a base quality greater than 20) over 96% were generated from TMs and WTs, respectively (Table S2). Compared with the reference genome of common carp, approximately 41,800,444 and 50,183,848 clean reads were mapped back to the transcriptome assembly in TM and WT libraries, respectively (Table S2). Compared to WT, a total of 331 transcript-derived unigenes were up-regulated while 792 unigenes were down-regulated in TM ( Figure S1A and Table S5). Here, 57 up-regulated DEGs were selected on the basis of FC TM/WT > 3 and 60 down-regulated DEGs were selected according to FC WT/TM > 15. Among them, the up-regulated DEGs in TM included epididymis-specific α-mannosidase, zinc-binding protein A33, FAM83H, interferon-induced very large GTPase 1, galactose-specific lectin nattectin, guanine nucleotide-binding protein subunit α-14, and others. (Table 1). For the down-regulated DEGs, they contained hemoglobin subunit α, 40S ribosomal protein S26, epididymis-specific α-mannosidase, purine nucleoside phosphorylase (pnp4a), purine nucleoside phosphorylase (pnp5a), and others (Table 2).

Transcriptome Analysis
After stringent quality assessment and data filtering, approximately 63,272,006 and 75,019,864 clean paired-end sequence reads with a length of 145 bp and a Q30 rate (those with a base quality greater than 20) over 96% were generated from TMs and WTs, respectively (Table S2). Compared with the reference genome of common carp, approximately 41,800,444 and 50,183,848 clean reads were mapped back to the transcriptome assembly in TM and WT libraries, respectively (Table S2). Compared to WT, a total of 331 transcript-derived unigenes were up-regulated while 792 unigenes were down-regulated in TM ( Figure S1A and Table S5). Here, 57 up-regulated DEGs were selected on the basis of FCTM/WT > 3 and 60 down-regulated DEGs were selected according to FCWT/TM > 15. Among them, the up-regulated DEGs in TM included epididymis-specific α-mannosidase, zinc-binding protein A33, FAM83H, interferon-induced very large GTPase 1, galactose-specific lectin nattectin, guanine nucleotide-binding protein subunit α-14, and others. (Table 1). For the down-regulated DEGs, they contained hemoglobin subunit α, 40S ribosomal protein S26, epididymis-specific α-mannosidase, purine nucleoside phosphorylase (pnp4a), purine nucleoside phosphorylase (pnp5a), and others (Table 2).   To identify the functional clusters and biochemical pathways of these DEGs, GO (gene ontology) enrichment analysis was performed by hypergeometric distribution testing. Based on the Nr annotation, 29,609 unigenes were assigned to 1739 GO terms. Using WEGO, 15,270 unigenes for the biological process category, 10,947 unigenes for the cellular component category, and 26,107 unigenes for the molecular function category were identified, respectively ( Figure S2). Based on the analysis of DEGs and GO annotation results, GO enrichment analysis was conducted under the condition of p < 0.05. Here, we listed 21 GO terms for up-regulated DEGs and 24 GO terms for down-regulated DEGs, respectively. For the up-regulated DEGs, the most abundant GO terms were annotated to peptidase activity, heme binding, iron ion binding, and other terms ( Figure 3A). For the down-regulated DEGs, the most abundant GO terms were annotated to molecular function, such as hydrolase activity, guanyl nucleotide binding, GTPase activity, nucleoside triphosphatase act, pyrophosphatase activity and GTP binding ( Figure 3B). interactions of these DEGs, the unigenes were mapped to the reference pathways recorded in the KEGG database. In total, 12,540 unigenes were annotated in KEGG and located to 5 known KEGG metabolic pathways, including cellular processes, environmental informational processing, genetic information processing, metabolism, organismal systems (Table 3). KEGG analysis of DEGs showed 4 different pathways, including other glycan degradation, antigen processing and presentation, microRNAs in cancer, and regulation of mitophagy-yeast (Table 4).   KEGG (Kyoto Encyclopedia of Genes and Genomes) is the main public database for systematic analysis of the metabolic pathways in cells. To further elucidate the biological functions and interactions of these DEGs, the unigenes were mapped to the reference pathways recorded in the KEGG database. In total, 12,540 unigenes were annotated in KEGG and located to 5 known KEGG metabolic pathways, including cellular processes, environmental informational processing, genetic information processing, metabolism, organismal systems (Table 3). KEGG analysis of DEGs showed 4 different pathways, including other glycan degradation, antigen processing and presentation, microRNAs in cancer, and regulation of mitophagy-yeast (Table 4).

miRNA Analysis
For the miRNA sequences, 10,979,661 and 12,296,112 total raw reads were obtained from TM and WT, respectively. After quality control, 9,876,199 and 11,368,562 total clean reads were left in TM and WT for further analysis, respectively. Among them, 7,898,353 in TM and 9,385,680 in WT clean reads were perfectly matched to the reference genome. (Table S3). To identify miRNAs in TM and WT, all of the mapped small RNA reads in the transcriptome were used as query against known and novel miRNAs in the database. In total, 54 known miRNAs and 396 novel miRNAs were obtained. Screening conditions for DEMs were p-adjust < 0.05 and |log2(Fold change)| ≥ 1. Firstly, we identified the expression profile of known miRNAs between TM and WT. The results indicated that only 2 differential expression known miRNAs, miR-146a and miR-153b, were highly up-regulated expression in TM compared to WT. Secondly, we identified the differential expression profile of novel miRNAs between TM and WT. The results indicated that 70 differential expression novel miRNAs were up-regulated and 125 differential expression novel miRNAs were down-regulated in TM vs. WT. Here, 2 up-regulated known DEMs, and 30 up-regulated and 30 down-regulated novel DEGs were listed according to the FC (fold change) value (Table S4). The overall distribution of the DEMs can be inferred from the volcanic map ( Figure S1B,C).
GO enrichment analysis was performed for the target genes of the DEMs, including two known and 195 novel DEMs. For the target genes (total 490) of two known DEMs, 15 GO terms from the three ontologies of molecular function, cellular component, and biological process were selected to describe their functions. As a result, some GO terms were identified in the guanyl-nucleotide exchange factor activity, GTPase activator activity, purine nucleotide metabolic process, GTP metabolic process, and other functions (Figure 4). Additionally, similar to known DEMs, the most abundant GO terms in the 64,546 target genes of 195 novel DEMs were annotated to metal ion binding, nucleic acid binding, purine nucleotide binding, and cellular process ( Figure 5).
Genes 2019, 10, x FOR PEER REVIEW 11 of 19 and WT, all of the mapped small RNA reads in the transcriptome were used as query against known and novel miRNAs in the database. In total, 54 known miRNAs and 396 novel miRNAs were obtained. Screening conditions for DEMs were p-adjust < 0.05 and |log2(Fold change)| ≥ 1. Firstly, we identified the expression profile of known miRNAs between TM and WT. The results indicated that only 2 differential expression known miRNAs, miR-146a and miR-153b, were highly up-regulated expression in TM compared to WT. Secondly, we identified the differential expression profile of novel miRNAs between TM and WT. The results indicated that 70 differential expression novel miRNAs were up-regulated and 125 differential expression novel miRNAs were down-regulated in TM vs. WT. Here, 2 up-regulated known DEMs, and 30 up-regulated and 30 down-regulated novel DEGs were listed according to the FC (fold change) value (Table S4). The overall distribution of the DEMs can be inferred from the volcanic map ( Figure S1B,C). GO enrichment analysis was performed for the target genes of the DEMs, including two known and 195 novel DEMs. For the target genes (total 490) of two known DEMs, 15 GO terms from the three ontologies of molecular function, cellular component, and biological process were selected to describe their functions. As a result, some GO terms were identified in the guanyl-nucleotide exchange factor activity, GTPase activator activity, purine nucleotide metabolic process, GTP metabolic process, and other functions (Figure 4). Additionally, similar to known DEMs, the most abundant GO terms in the 64,546 target genes of 195 novel DEMs were annotated to metal ion binding, nucleic acid binding, purine nucleotide binding, and cellular process ( Figure  5).

Guanyl-nucleoti de exchange factor activity
GTPase activator activity Ras guanyl-nucleotide exchange factor activity Sulfur compound binding Rho guanyl-nucleoti de exchange factor activity Protein tyrosine kinase activity Transcription coactivator activity Heparin binding Glycosaminoglycan binding Lipid transporter activity Transmembrane receptor protein tyrosine kinase activity Wide pore channel activity Rab GTPase binding Anion channel activity Alcohol binding Protein metabolic process Small molecule metabolic process Organic substance catabolic process Purine-containi ng compound metabolic process Purine nucleotide metabolic process Regulation of nucleotide metabolic process Regulation of purine nucleotide metabolic process Nucleoside metabolic process Ribonucleoside metabolic process Purine nucleoside metabolic process Purine ribonucleoside metabolic process Guanosine-cont aini ng compound metabolic process Purine nucleoside triphosphate metabolic process GTP metabolic process Regulation of GTPase   KEGG enrichment analysis was performed for the target genes of DEMs, including known and novel DEMs. For the known DEMs, 138 target genes were assigned to 344 KEGG pathways. For the novel DEMs, 9824 target genes were assigned to 344 KEGG pathways. For the known DEMs, target genes were mostly concentrated in metabolic pathways, Rap1 signaling pathway, pathways in cancer, biosynthesis of secondary metabolites, biosynthesis of antibiotics, and other pathways ( Figure 6A). For novel DEMs, target genes were mostly enriched in the mitogen-activated protein kinase (MAPK) signaling pathway, Human T-lymphotropicvirus 1 (HTLV-I) infection, oxytocin signaling pathway, axon guidance, microRNAs in cancer, and other pathways ( Figure 6B). KEGG enrichment analysis was performed for the target genes of DEMs, including known and novel DEMs. For the known DEMs, 138 target genes were assigned to 344 KEGG pathways. For the novel DEMs, 9824 target genes were assigned to 344 KEGG pathways. For the known DEMs, target genes were mostly concentrated in metabolic pathways, Rap1 signaling pathway, pathways in cancer, biosynthesis of secondary metabolites, biosynthesis of antibiotics, and other pathways ( Figure 6A). For novel DEMs, target genes were mostly enriched in the mitogen-activated protein kinase (MAPK) signaling pathway, Human T-lymphotropicvirus 1 (HTLV-I) infection, oxytocin signaling pathway, axon guidance, microRNAs in cancer, and other pathways ( Figure 6B).

Interaction Analysis Between miRNAs and mRNAs
In order to visualize the relationship between DEGs and DEMs, the RNA-Seq data and miRNA data were combined together to make the miRNA-mRNA regulatory network. For the known DEM, the up-regulated miR-146a was associated with 6 down-regulated DEGs including MGAT5B, MFAP4, GP2, htt, Sema6b, and Obscn (Table 5). In addition, miR-153b was also one of the known

Interaction Analysis Between miRNAs and mRNAs
In order to visualize the relationship between DEGs and DEMs, the RNA-Seq data and miRNA data were combined together to make the miRNA-mRNA regulatory network. For the known DEM, the up-regulated miR-146a was associated with 6 down-regulated DEGs including MGAT5B, MFAP4, GP2, htt, Sema6b, and Obscn (Table 5). In addition, miR-153b was also one of the known DEMs, but its target genes were not detected in the DEGs. For the novel miRNAs, correlation analysis was performed for 70 up-regulated miRNAs and their down-regulated target DEGs ( Figure 7A). The network diagram showed that some novel miRNAs, such as Novel_294, Novel_384, and Novel_272, were associated with multiple differentially expressed genes. Conversely, some target genes, such as TTN, Nlrc3, and NLRP, were regulated by multiple novel miRNAs ( Figure 7A). In addition, correlation analysis was also performed for 125 down-regulated miRNAs (such as Novel_212, Novel_377, and Novel_399) and their up-regulated target DEGs (such as Nrxn2, GALNT17, and MLRP3) ( Figure 7B). DEMs, but its target genes were not detected in the DEGs. For the novel miRNAs, correlation analysis was performed for 70 up-regulated miRNAs and their down-regulated target DEGs ( Figure 7A). The network diagram showed that some novel miRNAs, such as Novel_294, Novel_384, and Novel_272, were associated with multiple differentially expressed genes. Conversely, some target genes, such as TTN, Nlrc3, and NLRP, were regulated by multiple novel miRNAs ( Figure 7A). In addition, correlation analysis was also performed for 125 down-regulated miRNAs (such as Novel_212, Novel_377, and Novel_399) and their up-regulated target DEGs (such as Nrxn2, GALNT17, and MLRP3) ( Figure 7B).

Discussions
There are three main types of pigment cells in fish, namely melanophores, xanthophores, and iridophores [2]. Among them, iridophores could produce silver pigments to reflect light in fish scales, so the loss of silver pigments could make the fish transparent. In goldfish, the transparent mutants have been reported in China, which were named as crystal goldfish. Further genetic studies on crystal goldfish indicated that the transparent mutants were homozygotes of guanophore (-/-) [31]. Similarly, our present study demonstrated that the transparent trait in common carp was recessive over the reflective trait, and this phenomenon could completely conform to the Mendelian segregation rule of alleles. These results suggested that the transparent common carp could be a novel model to investigate the genetic basis of silver pigment synthesis.
In the present study, microscopic observations showed that the silver pigments and melanin were both lost in the scales of transparent mutants compared to that in wild types. Furthermore, GO and KEGG analysis showed that the DEGs between wild type and transparent mutants were mainly enriched in the metabolism of guanine, such as hydrolase, guanyl nucleotide binding, guanyl ribonucleotide binding, and GTPase activity. Consistent with these, previous studies have shown that the guanine was the major component of silver pigments, and it could combine with water to form a crystal plate to make scales reflect light [32]. The reduction in guanine deposition in iridophores could induce the establishment of the transparent medaka strains [33,34] (Wakamatsu et al., 2001;Ohshima et al., 2013). In addition, the differentiation of melanophores is coupled closely with the onset of pteridine synthesis which starts from GTP and is regulated through the control of GTP cyclohydrolase I activity [14].
Transcriptome analysis showed that several DEGs were detected between TMs and WTs, including purine nucleoside phosphorylase 4a (pnp4a) and endothelin receptor B (ednrb). As a member of pnp gene family in vertebrates, pnp4a was reported to have an important role in guanine synthesis in iridophores [14]. The previously reported zebrafish shady, rosy, and barley transparent mutations were all associated with the deficiency of iridophores caused by the mutant of pnp4a [7,32,35]. The CRISPR/CAS9 system has been used to further confirm that pnp4a was a causal gene on the guanineless locus of zebrafish [14,15]. In the present study, the slivery pigments in TM were completely lost, and pnp4a was less expressed in TM compared to WT. These results, taken together, suggested that pnp4a could also be involved in the formation of silver pigments in common carp. Similar to pnp4a, our results found that the pnp5a was also differentially expressed between WT and TM, which indicated that the pnp5a might be a guanine synthetic gene in teleost iridophores. Endothelin receptor b (ednrb) and its ligand, endothelin (edn), were vital in melanocyte development [36]. The signaling pathway of ednrb could promote the proliferation and differentiation of melanocyte stem cells and the regeneration of epidermal melanocytes, and ednrb mutations could give rise to pigmentation defects [12,13]. Further studies found that the activation of ednrb could induce the phosphorylation of microphthalmia-associated transcription factor (MITF), which was the key transcription factor for catalyzing tyrosine synthesis into melanin in melanocytes [15]. In our present study, the melanin in TM was completely lost and the expression level of ednrb in TM was significantly lower than that in WT, which suggested that ednrb might play an important role in the synthesis of melanin in common carp. Guanine nucleotide-binding protein was the subunit of the G protein-couple signaling pathway, and a previous study reported that it might be involved in the formation of coat color in mice [37]. Interestingly, our present study also found that guanine nucleotide-binding protein subunit α-14 (GNA14) was significantly differentially expressed between WT and TM, which suggested that GNA14 might be involved in the formation of pigment in teleost.
In mammals, miRNAs played pivotal roles in a variety of developmental processes, and their dysregulations were linked to several skin pigment diseases, such as melisma development [38] and melanoma [39]. In teleost, previous studies have reported that miR-429 was highly expressed in red skin of common carp, and inhibition of miR-429 could cause a substantial decrease in skin pigmentation [21], which suggested that the miRNAs also played an important role in the pigmentation of teleost. In our present study, DEMs analysis found that miR-146a and miR-153b were both more highly detected in TM compared to those in WT. In mammals, recent studies have reported that miR-146a played an important role in melanoma [39][40][41] and retinal pigment epithelial cells [42]. Similarly, miR-153 could suppress cell proliferation and invasion by targeting snail family transcriptional repressor 1 (SNAI1) in melanoma [43]. In the present study, further GO and KEGG analysis showed that the target genes of miR-146a and miR-153b were mainly enriched in the guanyl-nucleotide exchange factor activity, GTPase activator activity, purine nucleotide metabolic process, and GTP metabolic process. These results suggested that miR-146a and miR-153b might also play a vital role in the regulation of pigment synthesis in common carp. In addition, we also found several novel DEMs between TMs and WTs, such as Novel_294, Novel_384, Novel_272, Novel_377, Novel_212, and Novel_399. Interaction analysis showed that these DEMs were associated with several DEGs between TMs and WTs. Further GO and KEGG analysis showed that these target DEGs were mainly enriched in purine nucleotide binding and nucleic acid binding. These results suggested that these novel DEMs might also play an important role in the regulation of pigmentation.
In the present study, the transparent mutants were detected in the CyCa hybrid carp. Microscopic examination showed that the scales of transparent mutants were completely devoid of silver pigments and melanin. Genetic analysis indicated that the transparent trait was recessive over the opaque trait in common carp. To further examine the molecular mechanism of the transparent trait, RNA-Seq and miRNA-Seq were both used to detect the DEGs and DEMs between TMs and WTs. RNA-Seq analysis showed that several DEGs were detected between TMs and WTs, such as pnp4a and ednrb. Interestingly, previous studies have reported that the pnp4a and ednrb played an important role in the synthesis of silver pigments and melanin, respectively. In addition, our GO and KEGG analysis showed that the DEGs between TM and WT were mainly enriched in the metabolism of guanine, which was consistent with the previous studies that the silver pigment was the metabolic product of guanine. By using miRNA-Seq analysis, several known and novel DEMs were also detected between TM and WT. In addition, interaction analysis showed that the DEMs were associated with several DEGs. These results suggested that the miRNAs were also involved in carp pigmentation. In general, the present study could provide some new ideas for investigating the genetic mechanism of the transparent trait and the formation mechanism of sliver pigment in common carp.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/9/704/s1: Figure S1: Volcano plot of DEGs and DEMs in TM and WT: (A) Volcano plot of DEGs in TM and WT; (B) Volcano plot of DEMs (known) in TM and WT; (C) Volcano plot of DEMs (novel) in TM and WT. Splashes represent different genes or miRNAs, and the black splashes mean genes or miRNAs without significant differential expression. The red splashes mean significantly up-regulated genes or miRNAs, and the green splashes mean significantly down-regulated genes or RNAs. Figure S2: GO classification annotation chart: The abscissa represents the three ontologies of GO, which are molecular function, cellular component, and biological process. The ordinate represents the percent of genes and the number of genes. Table S1: Primer sequences. Table S2: Statistic of mRNA data. Table S3: Statistic of miRNA data.