Identifying and profiling the microRNAs associated with skin colour in the Muchuan black-bone chicken

Abstract MicroRNAs (miRNAs) are an abundant class of small non-protein-coding RNAs that regulate gene expression in plants and animals. Skin colour is an important economic consideration in chicken production, and chickens with black skin have high market value. Nevertheless, little research has been conducted on miRNA regulation of melanogenesis in chicken skin. In this study, we sampled the dorsal skins of chickens with black (BS) and white (WS) skin to construct six small RNA libraries. High-throughput sequencing technology was then used to identify which miRNAs were expressed differentially between the BS and WS phenotypes. A total of 645 known and 64 novel miRNAs were identified from the six sequencing libraries. Additionally, the expression of 18 miRNAs was significantly different between the two phenotypes, including 9 miRNAs that were up-regulated, and 9 that were down-regulated. We identified 2 miRNAs, i.e. miR-204 and miR-6631-5p, that may be important for melanogenesis in chicken skin. Our data contribute to the understanding of the molecular mechanism, through miRNA regulation, of melanin formation in chicken skin. Highlights Differentially expressed miRNAs were identified in black skin (BS) and white skin (WS) of black-bone chickens Eighteen DEMs were detected in BS and WS groups. miR-204 and miR-6631-5p were potential candidates influencing skin melanogenesis in chicken.


Introduction
Skin colour is an important economic consideration in chicken production. The Muchuan black-bone chicken has a deep black skin and is deemed to produce high quality meat (Yu et al. 2017). As such, black-bone chickens sell for a higher price than chickens with white skin. However, birds with white or lighter skin can be produced during the breeding process for black-bone chickens, which leads to economic losses. Research on chicken genetics and breeding has largely focussed on identifying genes associated with skin colour with the aim of increasing the amount of melanin in the skin using genetic techniques (Zhang et al. 2015a(Zhang et al. , 2015b. To date, more than 150 genes and loci associated with skin colour that influence pigmentation have been discovered (Cieslak et al. 2011). Still, the key regulators and regulatory mechanisms that govern skin colour remain unknown.
MicoRNAs (miRNAs) are approximately 19-24 nucleotides (nt) in length and are a type of endogenous non-coding RNA that regulate gene expression by targeting mRNAs for cleavage or translational repression (Bartel 2004). Endogenous miRNAs help to regulate various metabolic processes, including cell cycle progression and cell proliferation, differentiation, development, apoptosis, and migration (Brennecke et al. 2003;Shenoy and Blelloch 2014). Previous research has indicated that miRNAs can target genes or pathways involved in melanocyte biology (Mione and Bosserhoff 2015). In other vertebrates, such as ducks (Apopo et al. 2015), alpacas (Xue et al. 2012), goats (Wu et al. 2014), mice (Dynoodt et al. 2013), and fish (Yan et al. 2013), miRNAs are known to be responsible for pigmentation. However, no studies to date have reported on the role of miRNAs in determining skin colour in chickens. The black-bone chicken has intense pigmentation in its dermal skin layer across its entire body. Because skin colour in the black-bone chicken can vary (including white skin), the black-bone chicken is an ideal model for investigating the molecular mechanism of melanin pigmentation during melanogenesis in chicken skin.
In this study, RNA sequencing (RNA-Seq) was used to investigate miRNA expression profiles in chickens with black and white skin. We were able to identify miRNAs associated with melanogenesis that were differentially expressed during skin colour formation, extending the knowledge of miRNAs involved with melanin formation in chicken skin.

Ethics statement
All experimental procedures were approved by the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and the Animal Care and Use Committee of the Leshan Normal University.

Animal selection and phenotype measurement
A total of 133 female Muchuan black-bone chickens (40 weeks old) were purchased from the Muchuan County Black Phoenix Black-Bone Chicken Industry Co., Ltd. (Leshan, China) for use in this study. All chickens had free access to feed and water and were exposed to natural light (12 h/day) and temperatures of 17-25 C during the study. The dorsal skin colour of each chicken was measured using a portable colorimeter (3NH NR10QC; Shenzhen 3NH Technology Co., Ltd., Shenzhen, China). Colours were recorded as L Ã , a lightness value ranging from 0 (completely black) to 100 (completely white); a Ã , a colour scale ranging from green (negative) to red (positive); and b Ã , a colour scale ranging from blue (negative) to yellow (positive) ( Skrlep and Candek-Potokar 2007). While measuring, the colorimeter was in full contact with the skin to avoid losing any light from the instrument. Of the 133 chickens, we selected three individuals with black skin (group BS) and three with white skin (group WS) for RNA samples (Supplementary Table  S1). Phenotypic differences between the two groups were confirmed by quantifying the expression of the tyrosinase (TYR) gene.

Tissue collection and RNA isolation
The chickens in the BS and WS groups were killed by exsanguination to obtain skin samples. The dorsal skins were immediately frozen in liquid nitrogen and then stored at À80 C. Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA) following the manufacturer's protocol. RNA degradation and contamination were monitored using 1% agarose gels. RNA purity and concentrations were checked using the NanoPhotometer V R spectrophotometer (Implen, Westlake Village, CA) and the Qubit V R RNA Assay Kit of the Qubit V R 2.0 Fluorometer (Life Technologies, Carlsbad, CA). RNA integrity was also assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA).

Small RNA library construction and sequencing
Six small RNA libraries were constructed using the NEBNext V R Multiplex Small RNA Library Prep Set for Illumina V R (New England Biolabs, Ipswich, MA) following the manufacturer's instructions. Index codes were added so that sequences could be attributed to the proper sample. From each sample, 3 lg total RNA was used for library construction. The small RNAs were first ligated with Illumina 3t and 5a adapters. Then, first strand cDNAs were synthesised using M-MuLV Reverse Transcriptase (reverse ribonuclease H). Polymerase chain reaction (PCR) amplification was performed using the LongAmp Taq 2X Master Mix. PCR products were purified using 8% polyacrylamide gel (100 V, 80 min). DNA fragments of length 140-160 bp were recovered and dissolved in 8 lL elution buffer. Finally, library quality was assessed using the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. All small RNA libraries were sequenced with 50bp paired end reads using the Illumina HiSeq2500 platform. The raw sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) with accession number PRJNA515587.

Data processing
High-quality clean reads were obtained from raw data after removing reads with poly-N tails, 5i adapter contaminants, poly-A, -T, -G, or -C sequences, reads without 3a adapters or insert tags, and other low-quality reads. Clean reads were mapped to the Gallus_gallus-5.0 assembly of the chicken genome using Bowtie v. 1.2.2 (Langmead et al. 2009). After mapping, the small RNA reads were matched against known miRNAs in the miRBase 21.0 database (Ana and Sam 2014). Small RNA reads originating from protein-coding genes, repeat sequences, rRNA, tRNA, small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA) were also identified by alignment with the publicly available Rfam (http://rfam.sanger.ac.uk) and RepBase (http:// www.girinst.org/repbase/) databases. Potential novel miRNAs were predicted from the remaining unannotated clean reads using miREvo (Wen et al. 2012) and mirdeep2 (Friedl€ ander et al. 2012) software. Small RNA tags were aligned with mature miRNA with an allowance of one positional mismatch to detect potential miRNA base edits. Additionally, the miFam.dat data file from miRBase (http://www.mirbase.org/ftp.shtml) was used to look for miRNA families, whereas searches for the families of novel miRNA precursors were performed using the Rfam database (http://rfam.sanger. ac.uk/search/).

Identification of differentially expressed miRNAs
Relative miRNA expression was estimated by transcript per million clean tags (TPM) as follows: normalised expression ¼ (mapped readcount/total reads)Â1,000,000 (Zhou et al. 2010). The expression levels of the BS and WS groups were compared using the DESeq2 package (v. 1.8.3) in R (Anders and Huber, 2010). Significance thresholds for differences in expression were set at p < .05 and jlog 2 (fold change)j >1.

Predicting target genes and functional annotation
To understand the functional role of miRNAs that were differentially expressed in skin melanogenesis for the two chicken groups, the target genes of these miRNAs were predicted using the RNAhybrid (Rehmsmeier et al. 2004), PITA (Kertesz et al. 2007), and miRanda (Betel et al. 2008) analysis tools. Gene ontology (GO) analysis of the target genes was then conducted using the GOseq method (Young et al. 2010). In this analysis, p < .001 represents significantly enriched GO terms. Pathway enrichment analysis based on the Kyoto Encyclopaedia of Genes and Genomes (KEGG) Pathway database (http://www.genome.jp/kegg/pathway.html) was carried out using KOBAS v. 2.0 (Mao et al. 2005).
Validation of candidate miRNA and their target genes expression using quantitative real-time PCR Quantitative real-time PCR (qRT-PCR) was used to detect miRNAs that were differentially expressed in the two chicken groups. The total RNA from the same samples used for RNA sequencing was reverse-transcribed to cDNA using the TransScript V R miRNA First-Strand Synthesis Kit (Transgen, Beijing, China). Subsequently, miRNA and gene quantification were conducted using the Lightcycler Roche 480 system (Roche, Basel, Switzerland) with the LightCycler 480 SYBR Green I Master kit. Forward primers were  Table S2), and reverse primers were produced using a first-strand cDNA synthesis kit (Transgen, Beijing, China). The qRT-PCR primers of miRNA-target genes were listed in Supplementary  Table S2. All reactions were performed in triplicate for each sample, and the chicken U6 small-nuclear RNA (U6 snRNA) and b-actin were used to normalise the qRT-PCR results for each miRNA and gene, respectively (Sun et al. 2016). Relative miRNA expression was calculated using the 2 ÀDDCT method (Livak and Schmittgen 2001).

Phenotypic confirmation
TYR codes for tyrosinase, a multifunctional enzyme that is essential for melanin biosynthesis in melanocytes (Yu et al. 2019). TYR expression was used to confirm that the differences in skin colour had a genetic basis (Figure 1(a)). Our results revealed that TYR expression was significantly different (p < .01) between chickens with black and those with white skins (Figure  1(b)). Therefore, the chickens selected for the BS and WS groups were suitable candidates for screening the miRNAs involved in skin colouration.

Overview of miRNA sequencing
A total of 57.36 million and 59.76 million raw reads were obtained from the skins of BS and WS chickens, respectively. After filtering, we obtained more than 18 million clean reads for each library, with a mean clean read ratio of 99.09% (Supplementary Table S3). Most of the reads (>79%) were 21-23 nt in length, and the 22-nt small RNA was the most abundant (>51%). Almost all the small RNAs (>89%) mapped to the reference genome. The Q20 and Q30 values, which represent sequencing accuracy, were greater than 99% and 98% for each sample, respectively (Supplementary  Table S3). These results indicate that small RNA sequences were reliably produced in our study, and the mapped reads were suitable miRNA candidates in the subsequent analyses.

Annotation of small RNA sequences
The mapped small RNA reads were annotated according to references in the miRBase21.0 database. Most small RNAs were categorised as miRNAs (mean 53.14%) or were unannotated (mean 54.98%; Supplementary Table S4). The percentages of small RNAs annotated as rRNA, tRNA, snRNA, snoRNA, repeat sequences, exons, and introns were relatively low, ranging from 3.24 to 8.47%.

Identification of known chicken miRNAs and candidates for novel miRNAs
To identify the miRNAs in chicken skin, the sequences of mapped reads were compared with those of known mature miRNAs and their precursors in miRBase v21.0. A total of 645 known miRNAs and 64 novel miRNAs were identified from the six chicken sequence libraries (Supplementary Tables S5 and S6). In terms of miRNA expression, the read numbers of the 10 most prolific miRNAs comprised more than 65% of the total reads for each sequence library. Our analysis of miRNA expression revealed that most miRNAs were expressed by a small proportion of the miRNA genes, consistent with other chicken studies (Wu et al. 2017a(Wu et al. , 2017b. The miRNA gga-miR-143-3p had the highest expression level in both groups. Each miRNA was expressed at different levels among the six libraries (Supplementary Table S5). Most of the novel miRNAs were expressed at low levels compared to the known miRNAs. For example, the highest TPM of novel-353, a novel miRNA, were 25 and 492 in the BS and WS groups, respectively. Similar expression patterns have been reported in studies of chickens (Wu et al. 2017a) and other species (Zhiqiang et al. 2013;Wei et al. 2018). Additionally, results from the base bias analysis of known miRNAs in each library revealed that uracil (U) was the most common base at the first-nucleotide position in miRNAs with lengths of 20, 22, and 24 nt (Supplementary Figure S1). U was also most common  Figure S2).

Identification of miRNAs differentially expressed between black and white chickens
To identify candidate miRNAs associated with skin melanogenesis, relative expression was first calculated and normalised with reference to TPM. The miRNAs that were differentially expressed between the BS and WS groups were then identified using DESeq2. In total, 18 miRNAs had significantly different expression levels between the two groups, including 8 up-regulated and 8 down-regulated miRNAs (Table 1).

Target gene prediction and gene functional annotation
Target gene prediction helps to determine the potential downstream biological functions of miRNAs. We predicted that the miRNAs that were differentially expressed between the two chicken groups target a total of 755 genes. Some of these predicted target genes and their related families, such as the solute carrier family, the Ras oncogene family, and the tumour necrosis factor receptor superfamily, are known to be involved in melanogenesis. Furthermore, we conducted GO annotation and KEGG pathway analysis to identify the potential functional roles of these differentially expressed miRNAs. From these analyses, 29 significantly enriched GO terms and 5 KEGG pathways were identified. In GO terms, the target genes were strongly involved with 'cellular component' and 'biological process' (Supplementary Table  S7). Pathway analysis results revealed that the target genes were mainly involved with purine and pyrimidine metabolism, base excision repair, and glycosaminoglycan biosynthesis of chondroitin sulphate and RNA polymerase (Supplementary Table S7).

Validation of miRNAs differentially expressed and their target genes between black and white chickens
We used qRT-PCR to validate the expression levels of 6 miRNAs that were differentially expressed between the two groups. We selected 3 up-regulated miRNAs.
i.e. gga-miR-194, gga-miR-204, and gga-miR-215-5p, and 3 down-regulated miRNAs, i.e. gga-miR-1731-5p, gga-miR-1769-3p, and gga-miR-12264-3p. The results from qRT-PCR matched those from our RNA-Seq analysis (Figure 2), confirming that similar expression patterns were observed using both methods. Five candidate genes: SLC35F6, SLC18A3, SLC25A39, RAB24 and RABL2B, were predicted to be the target genes of the miRNAs gga-miR-1454, gga-miR-204, gga-miR-12264-3p, gga-miR-6631-5p and gga-miR-1731-5p, respectively. Expression levels of these miRNA targets in the skin tissues of BS and WS were assessed by qRT-PCR. The results showed that the expression levels of SLC18A3 and RAB24 had significantly different between the BS and WS groups (Figure 3). Figure 2. Validation of the microRNAs that were differentially expressed between chickens with black and white skin using quantitative real-time polymerase chain reaction (qPCR). Transcript per million clean tags values were used to calculate gene expression in RNA sequencing (RNA-seq) and to normalise gene expression in the WS group to the value of '1'. Ã p < .05, ÃÃ p < .01. All data are presented as mean ± standard error. BS: chickens with black skin; WS: chickens with white skin.

Discussion
We investigated miRNA expression profiles in the skins of black and white phenotypes of the black-bone chicken to elucidate the regulation and roles of miRNAs in skin pigmentation, particularly melanogenesis. In total, 18 miRNAs were found to be differentially expressed between the skins of black and white chickens. These include gga-miR-204, which is a miRNA known to be associated with melanogenesis. The expression of gga-miR-204 is significantly different between black and white feather bulbs in ducks (Apopo et al. 2015), and is also associated with skin colour in fish (Yan et al. 2013). In our study, gga-miR-204 was up-regulated in chickens with black skin, similar to results from previous studies (Adijanto et al. 2012b;Apopo et al. 2015). ; In humans, the microphthalmia-associated transcription factor (MITF) gene and gga-miR-204 are necessary for the maintenance of retinal pigment epithelial (RPE) cells. MITF promotes the differentiation of RPE cells by regulating gga-miR-204 expression (Adijanto et al. 2012a(Adijanto et al. , 2012bKutty et al. 2016), using the transcription factor PAX6 to upregulate gga-miR-204 (Shaham et al. 2013). At the start of melanogenesis in the RPE, PAX6 also synergizes with MITF to activate the expression of the genes (TYR, TYRP1, and dopachrome tautomerase) involved in pigment biogenesis (Shaul et al. 2014). MITF has been reported to regulate PAX6 expression in chicken RPE cells as well (Mochii et al. 1998). Additionally, our previous research has revealed that MITF expression is associated with skin colour in the black-bone chicken . These findings suggest that gga-miR-204 may play an important role in melanogenesis in chicken skin. Further studies are required to elucidate the molecular mechanisms that link gga-miR-204 to skin colour in chickens.
MiRNAs are key regulators in the post-transcriptional modifications of several biological processes (Krol et al. 2010). The majority of miRNAs have several to several thousand gene targets (Shiwei et al. 2014). Thus, it is challenging to identify specific targets in the study of miRNA function. To elucidate the potential regulatory network that affects skin colour in chickens, we predicted the target genes of the miRNAs that were differentially expressed between the black and white chickens. SLC18A3 was predicted to be the target gene of the miRNA gga-miR-204. Expression levels of SLC18A3 in the skin tissues of BS and WS were assessed by qRT-PCR. The results showed that SLC18A3 had significantly lower expression levels in BS than in WS. In humans and other vertebrates (fish, mice, birds, horses, sheep, and Xenopus laevis), members of the SLC family such as SLC24A5, SLC24A4, SLC7A11, SLC2A11B, SLC36A1, and SLC45A2 have also been reported to be associated with melanin synthesis or pigmentation (Mariat et al. 2003;Gunnarsson et al. 2007;Ginger et al. 2008;Sturm Figure 3. Comparison of expression levels for 5 target mRNAs using qRT-PCR. Ã p < .05, ÃÃ p < .01. All data are presented as mean ± standard error. BS: chickens with black skin; WS: chickens with white skin. 2009; Liu et al. 2011;He et al. 2012;Yuan et al. 2015). Results from our previous study revealed that five SLC members, SLC6A9, SLC38A4, SLC22A5, SLC35F3, and SLC16A3, were also associated with melanogenesis in the muscles of black-bone chickens ). The SLC superfamily comprises a major group of membrane transport proteins that are present in various cells (Amber et al. 2009). However, only a few SLC members are known to be associated with pigmentation and melanogenesis. Our results indicate that gga-miR-204 and its SLC target gene may be important for pigmentation in chicken skin. Further research will help to clarify the relationship of gga-miR-204 and SLC18A3 with melanogenesis in chicken skin.
Ras-related proteins not only are critical regulators of cellular membrane trafficking (Pfeffer 2001) but also are involved in a variety of processes such as skin pigmentation (Wasmeier et al. 2006). To date, several Rab proteins such as RAB27A (Strom et al. 2002;Hume et al. 2007;Yoshida-Amano et al. 2012), RABRP1 (Fujikawa et al. 2002), RAB8 (Chabrillat et al. 2005), RAB11B (Tarafder et al. 2014), RAB29 ), RAB7 (Jordens et al. 2006), RAB38 (Brooks et al. 2007), RAB9A (Mahanty et al. 2016), and RAB32 (Wasmeier et al. 2006) have been reported to have crucial roles in the pigmentation process, such as melanosome biogenesis, degradation, and transport. In this study, we predicted that gga-miR-6631-5p target RAB24, that belong to the Ras oncogene family. We have also found that RAB24 was expressed significantly higher in BS than that in WS. RAB27A can regulate the peripheral distribution of melanosomes in melanocytes, resulting in tissue colouration (Hume et al. 2001). Mutations in RAB27A can also disrupt melanosome transport, thereby affecting melanosome transfer to the surrounding tissues such as skin, feathers, or hair (Cieslak et al. 2011). Additionally, RAB32 and RAB38 are involved in the transport of melanogenic enzymes, including TYR and TYRP1, to melanosomes in melanocytes (Wasmeier et al. 2006;Marubashi et al. 2016). The expression of RAB29 mRNA was reported to be associated with muscle pigmentation in chickens as well . Therefore, we postulated that gga-miR-6631-5p may affect skin colour in the blackbone chicken by regulating the expression of its respective target gene (RAB24).

Conclusions
We used RNA-Seq to identify melanin-related miRNAs and their target genes in the black-bone chicken. In total, 18 differentially expressed miRNAs were identified between the skins of black and white chickens. MiR-204 and miR-6631-5p, were potential candidates influencing skin colour of chicken. Although additional studies are needed to characterise the functional roles of these miRNAs and their target genes in the regulation of skin melanogenesis, our results contribute an initial understanding of the molecular mechanisms governing the production of black or white skin in chickens.