Comparative microRNA-seq Analysis Depicts Candidate miRNAs Involved in Skin Color Differentiation in Red Tilapia

Differentiation and variation in body color has been a growing limitation to the commercial value of red tilapia. Limited microRNA (miRNA) information is available on skin color differentiation and variation in fish so far. In this study, a high-throughput Illumina sequencing of sRNAs was conducted on three color varieties of red tilapia and 81,394,491 raw reads were generated. A total of 158 differentially expressed miRNAs (|log2(fold change)| ≥ 1 and q-value ≤ 0.001) were identified. Target prediction and functional analysis of color-related miRNAs showed that a variety of putative target genes—including slc7a11, mc1r and asip—played potential roles in pigmentation. Moreover; the miRNA-mRNA regulatory network was illustrated to elucidate the pigmentation differentiation, in which miR-138-5p and miR-722 were predicted to play important roles in regulating the pigmentation process. These results advance our understanding of the molecular mechanisms of skin pigmentation differentiation in red tilapia.


Introduction
Tilapia (genus Oreochromis), which originated in Africa, is one of the most important food fishes in the world. Due to its fast growth and high adaptability to complex cultural systems, tilapia is widely accepted and has emerged as the second most produced fish in aquaculture around the world [1,2]. There exist various kinds of tilapia species, including Nile tilapia (O. niloticus), blue tilapia (O. aureus), mutant reddish-orange Mozambique tilapia (O. mossambicus), red tilapia (Oreochromis sp. red tilapia) and others originating from a multi-crossbreeding strategy. Among those shown above, red tilapia, largely identifiable by its uniform red skin and bright pink peritoneum, is a very attractive commercial breed in certain markets.
Red tilapia usually display a variety of striking skin colors as a result of long-term natural selection, genetic regulation, feeding conditions, camouflage for threatening behaviors, mate-choice, and adaptation to low temperature, which play important roles in numerous biological processes [3][4][5][6][7]. The pigmentation differentiation in genetic breeding and skin color variation during the overwintering period are the main problems limiting the development of commercial red tilapia cultures. Coloration patterns including whole pink (WP), pink with scattered black spots (PB) and pink with scattered red spots (PR) have been found in the Malaysian red tilapia breeding population. The pigmentation differentiation is not reversible and skin color variation is reversible with the environmental temperature increasing [3][4][5][6][7]. Genetically, the formation of skin color has been proved to be a rather complicated bioprocess which mainly involves diverse pigments' biosynthetic pathways determined by chromatophores or pigment cells associating with a series of cellular, nutritional, physiological and environmental factors [5,8]. It was reported that dozens of genes have been involved in the color pattern formation of zebrafish [9][10][11][12][13]. Several transcriptome analyses of non-model fish with various skin colors have been conducted to explore coloration-related genes and the genetic basis effectively [5,12,13].
MicroRNAs (miRNAs) are non-coding RNA molecules ranging from 21-25 nt which play important roles in regulating gene expression in biological processes. miRNAs have also been suggested to have a crucial role in the coloration of body skin. For instance, miR-137 affected body color pattern in mice [14], and loss of miR-8 decreased pigmentation of the dorsal abdomen in the fruit fly [15]. It has been reported that several miRNAs were related to the regulation of skin color in the common carp [16], and miR-429 seemed to be a potential regulator of pigmentation by targeting the 3'untransliated region (3 UTR) region of foxd3 [17]. However, there are still few reports on the functions of miRNA regarding skin color differentiation.
In our previous study, an Illumina RNA-seq of transcriptome were conducted on three red tilapia varieties born with different colors (WP, PB and PR). We screened 148 differentially expressed genes (DEGs) and identified several key genes related to pigment synthesis, such as tyr, tyrp1, pmel/silv, sox10, slc24a5, cbs and slc7a11 [5]. In this study, a high-throughput sequencing strategy was conducted to identify skin color-related miRNAs among these three types of coloration, and several differentially expressed miRNAs were validated by quantitative real-time polymerase chain reaction (qRT-PCR). Furthermore, putative target genes of differentially expressed miRNAs were predicted and a subsequent enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was employed to gain insight into the function of miRNAs in the regulation of skin color in fish.

Signature of Sequencing Data in Red Tilapia.
To obtain the miRNA expression signature related to red tilapia skin colors, we constructed small RNA sequencing libraries from 9 sample pools isolated from the skins of adult red tilapia with three different colors: WP, PB and PR (Supplementary Figure S1, triplicate for each color). An average 9,872,682 raw reads, corresponding to 9,043,832 clean reads after removing reads containing adaptor sequences and low-quality reads (Table 1), were obtained from each library. We summarized the length distribution of all clean reads and found that 8 samples, except one sample in the WP group, showed a peak of miRNA of typical length at 22 nt ( Figure 1A). Approximately 53.32% to 68.73% clean reads were matched to the Nile tilapia genome (Supplementary Table S1). All clean reads were used for functional annotation ( Figure 1B) and the remnant sRNA reads with rRNA removed were used for searching published mature miRNAs. The unannotated sRNA reads were processed to predict novel miRNAs using MIREAP software (http://sourceforge.net/projects/mireap). A total of 275 miRNAs were identified to be the same or similar to the mature miRNAs (known), and 205 miRNAs were predicted to be novel miRNAs with hairpin structures (novel, Supplementary Table S2).

MicroRNAs (miRNAs) Show Differential Expression in Tilapia with Different Skin Colors
With the criteria of |log 2 (Fold change)| ≥ 1 and q-value ≤ 0.001, we identified 103 significant differentially expressed miRNAs (DEMs) in PB (including 64 up-regulated and 39 down-regulated miRNAs) compared with WP ( Figure 2A); 79 significant DEMs in PR (containing 38 up-regulated and 41 down-regulated miRNAs) compared with WP; and 74 significant DEMs in PB (containing 44 up-regulated and 30 down-regulated miRNAs) compared with PR ( Figure 2B; Supplementary Table  S3). Among these miRNAs, 158 miRNAs including 7 overlapped DEMs were identified ( Figure 2C). Moreover, 46 miRNAs (including 16 known and 30 novel DEMs) were DEMs in the PB and PR groups compared with the WP group ( Figure 2C). The heatmap of these 16 known DEMs indicated that most of them, except down-regulated miR-92b-3p, miR-141-5p, miR-2188-3p, were up-regulated in PB and PR skins compared with WP skins (Figure 3). We then randomly selected 15 DEMs (including 3 down-regulated and 12 up-regulated miRNAs) and validated the expression patterns using qRT-PCR. The results showed that the qRT-PCR expression patterns were absolutely in agreement with those from RNA-seq analysis (Figure 4), indicating the RNA-seq data was reliable.
Previous observations have suggested that miR-107b was differentially expressed in the brain tissue of rainbow trout and found to target sox6 [20], which might act as a transcription factor indirectly regulating embryonic development and melanin biosynthesis by modulating sox10 in B16 melanoma cells [21]. The melanocortin system exerted its multiple functions via G protein-coupled receptors (GPCRs) [22]. In non-small cell lung cancer cells, GPR124, one of the GPCRs, was proved to be a direct target of miR-138-5p and its expression was suppressed by miR-138-5p [23]. Moreover, miR-141-3p was shown to indirectly interact with the GPCR signaling pathway via its targets [24]. MiR-27d, miR-16a and miR-722 had been identified as playing crucial roles in sex differentiation of tilapia, the cell cycle of whitefish liver cells, and immune system development and response, respectively [25][26][27][28]. However, their specific functions in pigmentation in red tilapia were poorly reported. In this present study, the fact that these four DEMs, miR-107b, miR-138-5p, miR-141-3p and miR-27d, were up-regulated in fish with PB and PR skins compared with fish with WP skins suggested that they played important roles in pigmentation in red tilapia. To predict the metabolism pathway of the target genes of these 16 interested DEMs between fishes with PB and PR skins, KEGG pathway analysis was conducted on their target genes. Enrichment analysis showed that their target genes were classified into several pigmentation-related pathways, including the Wnt signaling pathway, the MAPK (mitogen-activated protein kinase) signaling pathway, and melanogenesis ( Figure 6; Supplementary Table S5). Among the targets, mc1r, the most important regulator of melanogenesis, positively regulated cyclic AMP (cAMP) response-element binding protein (creb) and downstream mitf to increase melanin synthesis [29]; wnt4 was also reported to be maintained at higher levels in normal peripheral retinal pigment epithelium cells in dark agouti rats [30]. In normal melanocytes, the activation of Mc1r-cAMP-Mift pigmentation signaling cascade has been proved to induce the production of transforming growth factor (TGF)-β1 in transformed melanocytes in human melanoma cells and play crucial roles in ultraviolet (UV) light exposure-induced skin pigmentation [31,32]. Moreover, the inactivation of the Mc1r-cAMP-Mift signaling cascade was essential for sesamol-decreased pigment formation and melanogenesis in melanocyte cells and zebrafish [33].
Recently, a transcriptome analysis showed wnt4 related to melanogenesis in common buzzards with different color [34]. The gallic acid-induced inhibition of melanogenesis and hypopigmentation in B16F10 cells was associated with the activation of Akt, MEK/ERK (extracellular regulated protein kinases), and Wnt/β-catenin signaling [35]. Also, mift was reported to be required for β-catenin-induced melanoma growth [19]. These studies suggested that all these DEMs played crucial roles in pigmentation in red tilapia by modulating the related pathways via its target genes.

Sample Collection
Red tilapias born with WP, PB and PR skin colors (53.1 ± 1.30 g, Supplementary Figure S1), were obtained from the Qiting Pilot Research Station (Yixing, Jiangsu, China), affiliated to the Freshwater Fisheries Research Center (FFRC), Chinese Academy of Fishery Sciences. Fishes were maintained in conical fibreglass tanks (water depth: 50 cm, volume: 256 L) in a flow-through water system during the acclimation and experimental periods. Water temperature was maintained at 27 ± 1 • C, pH = 7-8, with dissolved oxygen (DO) > 6 mg·L −1 and NH 4 -N < 0.5 mg·L −1 . Aeration was supplied 24 h per day and the photoperiod was 12D:12L.
Skin tissues were collected from 18 WP, 18 PB (pink, scattered 2/3 or above with black spots on all skin) and 18 PR (pink, scattered 2/3 or above with red spots of all skin) red tilapia individuals, respectively. All samples were immediately frozen in liquid nitrogen and then stored at −80 • C before RNA isolation.

Small RNA Library Preparation and Sequencing
Total RNA was extracted from skin tissues using TRIZOL (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. Then genomic DNA was removed from RNA samples using DNase (New England Biolabs, Ipswich, MA, USA), and RNA purity was assessed using a Nanodrop2000 spectrophotometer (Thermo Fisher Scientific, BRIMS, Cambridge, MA, USA). An equal amount of total RNA from 6 individuals from each color group (WP, PB and PR red tilapia) was pooled, and 9 samples or RNA pools, including 54 individuals, were prepared accordingly. The purified RNA pools were ligated with Illumina 3 and 5 adapters (Illumina, San Diego, CA, USA) using T4 ligase (New England Biolabs, Ipswich, MA, USA). RNA was reverse-transcribed into the first strand cDNA using reverse transcriptase and amplified by PCR for 15 cycles using primers complementary to the adaptor sequences. Then, nucleotide fractions 140-150 bp in length were purified for Illumina sequencing library preparation. For sequencing, each library was loaded into one lane in the mode of 50 bp single-end Illunima Hiseq 2500 (Illumina, San Diego, CA, USA).

Sequencing Analysis for miRNA Profiling
Quality control of all sequencing reads was conducted by FastQC software (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/). We removed sequences in disorder, including poor quality reads, 3 adaptor null reads, reads with 5 adaptor contaminants, and reads shorter than 18 nt. The remaining sequences were used for mapping to the whole genome sequence of Nile tilapia [36], using the SOAP2 program (http://soap.genomics.org.cn) with a tolerance of 2 mismatches [37]. Then, rRNAs, tRNAs, snRNAs and snoRNAs in the matched sequences were filtered out by blasting against Rfam (http://rfam.sanger.ac.uk/) and National Center for Biotechnology Information (NCBI) GenBank (http://www.ncbi.nlm.nih.gov/genbank/) databases. After being classified into different categories, the remaining reads were aligned to the miRNA precursors in zebrafish in the miRBase version 21 database (http://www.mirbase.org/) to identify conserved miRNAs [38]. The unannotated sRNAs reads were conducted using MIREAP software (http://sourceforge.net/projects/mireap) to predict candidate novel miRNAs.

Data Normalization, Processing and Identification of DEMs
Sequencing data was normalized using the following formula: normalized expression = (actual miRNA sequencing readscount/total miRNAs reads count) × 1,000,000. Finally, DEGseq [39] was used to identify DEMs based on the miRNAs expression values by pairwise comparison. The cut-off criteria for DEMs were |log 2 (Fold change)| ≥ 1 and q-value ≤ 0.001.

Target Prediction for DEMs
The target genes of DEMs were identified using the RNAhybrid [40], miRanda [41] and Targetscan 7.0 [42] with default settings. The genes supported by either two algorithms were considered as the targets of DEMs in this study. The target reference sequences were whole genome sequences of Nile tilapia.

Enrichment Analysis for DEM Targets
The KEGG pathway database informs people of how molecules or genes work [43]. In order to investigate the related biological function of DEM targets in the pigmentation, the pathways of biochemical and signal transduction significantly associated with the DEM targets were determined through the KEGG pathway (available online: http://www.kegg.jp/) analysis using the KOBAS server [44].

MiRNA-mRNA Regulatory Network Analysis
According to the regulatory relationships between key DEMs and their targets, the regulatory network containing key miRNAs and target mRNAs was constructed and visualized by Cytoscape software, a standard tool for integrated analysis and visualization of biological networks [45].

QRT-PCR Analysis
Total RNA was extracted as described above. For the reverse transcription of miRNAs, the Prime Script RT Reagent Kit (Takara Bio, Dalian, China) were used. QRT-PCR was performed on the ABI PRISM 7500 Real-time PCR System (Applied Biosystems, Foster City, CA, USA) using the 2X SYBR Green Master Mix reagent (Takara Bio, Dalian, China) following the below steps: initial denaturation at 95 • C for 5 min, followed by 95 • C for 15 s and 60 • C for 45 s for 40 cycles, and one cycle of 95 • C to 65 • C. All reactions were conducted in triplicate. The relative expression levels of the DEMs were measured in terms of threshold cycle value (Ct) and were normalized to 5S rRNA using the equation 2 −∆∆Ct . The miRNA specific primers (Supplementary Table S6) were synthesized by Genscript Inc. (Nanjing, China).

Conclusions
In this study, we screened the DEMs related to red tilapia pigmentation and predicted their target genes which might be associated with pigmentation process and melanogenesis by being involved in signaling pathways, including the Mc1r-cAMP-Mift signaling cascade, via controlling target genes. Moreover, the Wnt and MAPK signaling pathways were predicted to be associated with the pigment formation in red tilapia. Our findings provided insights into the mechanisms underlying miRNA-mediated pigmentation in red tilapia. Further experiments are needed to validate the exact roles and molecular mechanisms of miRNAs and their targets in pigment formation in red tilapia.