Regulatory pathway analysis of coat color genes in Mongolian horses

Background Studies on the molecular genetics of horse skin pigmentation have typically focused on very few genes and proteins. In this study, we used Illumina sequencing to determine the global gene expression profiles in horses with white-colored coats and those with black-colored coats, with the goal of identifying novel genes that could regulate horse coat color. Results Genes encoding ribosomal-associated proteins were highly expressed in horse skin. We found a total of 231 unigenes that were differentially expressed between horses with white coats and horses with black coats; 119 were down-regulated, and 112 were up-regulated. Many of the up-regulated genes in black horses, such as genes related to tyrosine metabolism, may directly regulate dark coat color. Keratin genes, MIA family genes, fatty acid-related genes, and melanoma-associated genes were also differentially regulated, which suggests that they may play important roles in coat color formation. Conclusions These findings show that the transcription profiles from white and black horse skin provide useful information to understand the genetics underlying the control of skin melanin synthesis in horses, which may enhance our knowledge of human skin diseases, such as melanoma and albinism.


Background
Mammalian coat colors are determined by the quantities and distributions of melanin, which are dependent on the interaction between the genotype and the environment [1]. The most important pigments, including melanin and its derivatives, are synthesized in melanocytes by oxidation of tyrosine or tyrosine-related materials.
Mutation analyses have identified various genes that are involved in determining coat color in the horse [2]. Many of these genes regulate the expression and distribution of melanin, and their mutation can cause different phenotypes of the coat, skin, and eyes. Many of these gene mutations are situated within more than 60 loci that affect phenotype, and most are highly conserved in mammals, though the extent of their effect on pigment deposition varies [3]. Genes that commonly regulate skin and coat color in different mammalian species can be separated into 2 categories: one regulates the production, proliferation, or migration of different types of melanocytes, and the other affects pigment synthesis directly. Therefore, the formation of different skin and coat colors is determined by the regulation of genes that can change the progression/differentiation of melanocytes or the process of melanin synthesis [2].
At the cellular level, melanocytes put pigment granules into hair and skin cells; the presence and function of melanocytes therefore determine the amount, type, and character of pigmentation. Melanocytes originate along the neural crest, which also gives rise to the spinal cord and brain, and then migrate to the skin during embryogenesis [4]. The importance of this is that the pigmentary and nervous systems are closely allied in the embryo, and certain genes affect both.
Melanocytes can produce either pheomelanin or eumelanin as determined by a receptor on their surface, MC1R [5]. This receptor is activated by melanocyte-stimulating hormone, which is secreted by the pituitary gland. When activation occurs, the cells form eumelanin, while in the absence of activation, the cells form pheomelanin [6]. Whether the receptor becomes activated is determined by 2 factors: the presence or absence of melanocyte-stimulating hormone, and the presence or absence of cell surface receptors. Because melanocyte-stimulating hormone is consistently available to most cells in horses, regulation at the level of surface receptors is more important.
Previous studies to identify genes involved in skin pigmentation have focused on genetic polymorphisms. In the present study, we generated transcriptome profiles for horses with black or white skin utilizing highthroughput RNA deep-sequencing technology. Black and white skin from Mongolian horses was collected, and differentially expressed genes were identified by RNA-Seq, a high-throughput sequencing platform that allows for the detection and quantification of transcripts at low abundance, including novel transcript units [7]. The identification of genes for melanin production, distribution, and formation can provide a theoretical basis for the selection of skin traits during the selective breeding of horses. Additionally, increased understanding of the molecular mechanisms involved in skin pigmentation may have significance for other animals, including humans, in terms of skin-related diseases such as melanoma and albinism.

Animal selection
This investigation involved 2 Mongolian horses with black coat color (black1, black2) and 2 Mongolian horses with white coat color (white1, white2). All guidelines of the Institutional Animal Ethics Committee and the Animal Care Guidelines of the Inner Mongolia Agricultural University were followed when conducting experiments on the Mongolian horses.

Feeding and management of animals
The black and white Mongolian horses were raised under equivalent conditions. All animals had free access to the same natural pasture. The animals were provided the same feeding regime for 6 months prior to skin collection.

Sampling and excision biopsy
Skin biopsies were taken from the 4 Mongolian horses to analyze the enzymatic activity and metabolic status of melanocytes in the skin. The samples were collected in an exam room at the Inner Mongolia Agricultural University according to internal protocols and procedures. To achieve both sedation and analgesia, each animal received an intravenous bolus of 0.4 ml of detomidine hydrochloride per 100 kg of body weight. The inside of the hind legs underwent a wide cleaning and disinfection. After shaving to remove the hair, an excisional biopsy was carried out, and a small area of ≤1 cm 2 was removed and stored in a tube with a sample protector for tissues. The wound was then immediately sutured with reabsorbent closures.

RNA extraction
Total RNA was extracted using Trizol (TaKaRa) according to the manufacturer's instructions, and its purity and integrity were assessed by 1% agarose gel electrophoresis. Subsequently, genomic DNA was removed by treating the RNA sample with RNase-free DNase I for 30 min at 37°C.

cDNA library construction and Illumina sequencing
To obtain poly (A) mRNA from the total RNA, oligo (dT) magnetic beads were used (Illumina). The RNA was broken down into short fragments by the addition of fragmentation buffer. These short fragments were then used as templates for first-strand cDNA synthesis with random hexamers and reverse transcriptase (Illumina). To synthesize the second strand of the cDNA, a solution of RNase H (Illumina), DNA polymerase I (Illumina), dNTPs, and buffer was used. The ends of the resulting double-stranded cDNA fragments were repaired, and adapters were ligated. The final version of the cDNA library was prepared from these products by purification and subsequent amplification by PCR. Using the Illumina Hiseq 2000 platform, the four prepared cDNA libraries were sequenced, resulting in 100-bp paired-end reads.

Sequence preprocessing and functional annotation
The raw sequence data was processed using a Perl script developed in our lab; it removed the adapter sequences and filtered the low-quality reads and genes with N ≥ 10%. The resultant clean sequence data were mapped to the horse genome (Equus caballus) using TopHat2 without the discordant or mixed options. The reads that were uniquely mapped to the horse genome were analyzed to determine the approximate gene abundance, and gene expression levels were calculated by the reads per kilobase per million mapped (FPKM). Additionally, differentially expressed genes between the 2 Mongolian horses with black coat color and 2 Mongolian horses with white coat color were estimated by edgR using FPKM based on multiple significance tests. If |log2 (fold change)| > 1.4 and FDR < 0.05, then these genes were considered as differently expressed.

Enrichment analysis
The functional annotation and pathway enrichment of genes that were differentially expressed between blackand white-skinned horses were performed using the online analysis tool DAVID (DAVID 6.7: https:// david.ncifcrf.gov/tools.jsp), which is a program that manages the enriched Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that characterize genes. The differentially expressed genes (DEGs) were mapped to the GO database, and the hypergeometric test was utilized to determine which GO terms were significantly enriched among the DEGs against the background of the horse genome; GO terms were considered to be significantly enriched in the DEGs when they had corrected p-values <0.05. To determine which KEGG pathways, and therefore which complex biological behaviors, were enriched in the DEG data, a similar method was used (threshold: corrected p-value <0.05).

Quantitative real time PCR (qRT-PCR) validation
Immediately after its removal, a piece of skin (approximately 0.5 cm 2 ) was frozen in liquid nitrogen until needed for later qRT-PCR analysis. Its total RNA was obtained using Trizol (TaKaRa) according to the manufacturer's protocol. After its extraction, the total RNA was dissolved in nuclease-free water. Approximately 0.5 μg of total RNA was used as a template for first-strand cDNA preparation using the PrimerScript RT reagent kit (TaKaRa) according to the manufacturer's instructions. The synthesized cDNA was diluted to 0.1 μg/μl for analysis by qRT-PCR (Bio-Rad) using the SYBR Green Realtime PCR Master Mix (TaKaRa). The housekeeping gene GAPDH was selected as the control. To determine the relative levels of gene expression, the 2 −ΔCt (ΔCt = Ct target gene − Ct GAPDH ) method was used. ANOVA (using SAS software 9.0) was used to determine which genes were differentially expressed between black-and white-colored horse skins.

Results
Assembly and functional classification of unigenes from horse skin A total of 24,301,563, 22,691,201, 19,423,074 and 20,465,943 pair-end reads were generated in the black1, black2, white1 and white2 libraries, respectively, after the raw reads were filtered from the skins of white and black horses. Of the 231 known genes, 202 were marked using GO analysis. This collection of genes was sorted into 63 classes using their putative functions as a basis for categorization. The largest collection of genes was sorted entirely by general function. Additionally, we annotated the known genes utilizing GO classification analysis and sorted them into 3 groups (74.6% were "biological processes," 11.1% were "cellular components," and 14.3% were "molecular function") based on their presumed functions (Fig. 1).

Highly expressed genes in horse skins
The keratin-family gene latherin (LATH) was the most highly expressed gene in horse skin (Table 1). Twentyfive of the 30 most highly expressed genes encoded ribosomal proteins; Equus caballus eukaryotic translation elongation factor 1 alpha 1 (EEF1A1) and three unknown genes were also highly expressed in the horse skin.
GO enrichment analysis of in genes that are differentially expressed in black-and white-colored skins We assessed the 202 unigenes that were differentially expressed in black and white coat color horses by GO enrichment analysis. Among 97 known genes that were classified as "biological processes," 15 genes were related to tyrosine metabolism; additionally, there were 58 known genes related to skin development. Among the 54 unigenes that were classified as "cellular components," 17 genes were related to the endoplasmic reticulum; and among 41 known "molecular function" genes, 13 genes were related to acyltransferase activity and carboxylesterase activity ( Table 2).

Analysis of differentially expressed genes in black and white coat color skin
When gene expression in white-and black-colored coats was compared, we found 231 genes that were differentially expressed; 119 DEGs were down-regulated in white coat color, and 112 DEGs were up-regulated (Fig. 2). The 20 genes with the most significantly different expressions are listed in Table 4; this list shows 10 up-regulated genes and 10 down-regulated genes.

Discussion
The mechanisms of melanogenesis are intricate, and recent publications have furthered our understanding of melanin production in the skin [8,9]. The interactions of several basic factors, most of which are from single genetic loci, are responsible for the various colors of most horses. With the development of advanced sequencing technology, skin and coat color genetics have been well studied [10]; however, most studies have focused on gene polymorphisms, whereas the impact of gene expression on coat color has not been well characterized.
The GO and KEGG pathway analyses indicated that the vast majority of the DEGs were related to "biological processes". The pathways related to tyrosine metabolic and catabolic processes were of particular interest in our dataset. Tyrosine is a non-essential amino acid responsible for melanin production via its oxidation and polymerization [11].
To confirm the Illumina results, we performed qRT-PCR for specific genes. The DCT, HPD, and TAT genes were verified to be differentially expressed in white and black horses. These genes each participate in tyrosine metabolism. Dopachrome tautomerase (DCT) catalyzes the rearrangement of dopachrome to the carboxylated derivative DHICA [12]. The results of DCT on UV DNA damage and survival pathways in human melanocytic cells were determined by knockdown tests in melanoma cells, neonatal foreskin melanoblasts in monoculture, and co-cultures with human keratinocytes [13,14]. DCT plays a major role in the coat color of cattle [15,16] and rhesus macaques [17]. At the DCT locus, 3 mutations are known that affect pigmentation phenotypes. Guyonneau (2004) generated a knockout of the Dct gene in mice with effects restricted to pigment production and coat color [18]. As shown in Fig. 4, DCT is involved in regulating eumelanin and pheomelanin levels [19].
The Fe(II)-containing non-heme oxygenase 4-hydroxy phenylpyruvate dioxygenase (HPD) catalyzes the conversion of 4-hydroxyphenylpyruvate into homogentisate in the catabolism of tyrosine. HPD is linked to 1 of the oldest known inherited metabolic disorders, alkaptonuria, which is caused by low levels of homogentisate in the bloodstream. HPD is also directly linked to type III tyrosinemia [20]. When its concentration is low in the human body, high levels of tyrosine occur in the blood, which can cause mild mental retardation at birth and subsequent degradation in vision [21]. The liver enzyme tyrosine aminotransferase (TAT) catalyzes the transformation of tyrosine into 4-hydroxyphenylpyruvate [22]. It is the rate-limiting enzyme in the tyrosine catabolic pathway that also includes HPD [23]. Because the accumulation of tyrosine in blood causes toxic effects to tissues and organs [24], the breakdown of tyrosine by TAT is very significant for human health. Type II tyrosinemia is caused by a deficiency in TAT [25]. Thus, each of these differentially regulated genes has functions in important processes affected by tyrosine metabolism. These 3 genes were up-regulated in the tyrosine metabolism pathway in horses of black coat color, indicating that the DCT, HPD, and TAT genes affect the formation of dark coat color.
In our study, 4 keratin genes were indicated to be differentially expressed in black/white horse skin: keratin 35 was up-regulated in white coat color skin, and keratin 34, keratin 26, and keratin 39 were up-regulated in black coat color skin. Keratinocytes generate a large number of paracrine factors that affect the growth of pigment cells and their proliferation and behavior; changes in the expression of these factors affect melanocytosis via receptor-mediated signaling pathways [26]. This result suggests that the keratin family also affects dark coat color formation.
We also identified 2 differentially expressed genes from the melanoma inhibitory activity (MIA) family: MIA and otoraplin (OTOR). MIA has been shown to have growthinhibitory activity on malignant melanoma cells in vitro [27,28]. MIA is expressed and secreted by melanoma cells, but not melanocytes, and its mRNA levels parallel progressive malignancy of melanocytic tumors [29,30]. In fact, increased MIA serum levels are considered to be a reliable tumor marker in detecting and monitoring metastatic disease and responses to therapy [31,32]. OTOR is secreted via the Golgi apparatus and may function in cartilage development and maintenance. A frequent polymorphism in the translation start codon of this gene, potentially associated with alternative polyA sites, is associated with forms of deafness [33]. Gray horses are at an increased risk for melanoma, with 70-80% over the age of 15 presenting with melanomas [34]. It is possible that melanoma in light coat color horses is associated with the expression of MIA family genes. Other genes identified as differentially expressed in our study include PTN, RGS13, and SPARC. The secreted heparin-binding protein pleiotrophin (PTN)  has been shown to be involved in cell growth and differentiation [35]. Not much is known about the effects of PTN on skin pigmentation and melanocyte function. Transfection studies, however, have shown that PTN decreases melanogenesis through MITF degradation via Erk1/2 activation [36]. In vitro, the chemotaxis of B cells is controlled by the regulator of G-protein signaling 13 (RGS13); this control is thought to be carried out by increasing the GTPase activity of G α proteins that are coupled to chemokine receptors [37]. RGS13 expression also reduces cAMP production [38], which plays an important role in melanoma even though genetic alterations in components of this pathway are not commonly found in melanomas [39,40]. The incorporation of collagen into the skin is controlled by secreted protein acidic and rich in cysteine (SPARC) (or osteonectin or BM-40) [41]. SPARC has various roles that it carries out in cooperation with many extracellular matrix units: it behaves as a de-adhesive molecule, regulates cytokine and growth factor activities, and inhibits the cell cycle [42]. SPARC has also been shown to be strongly expressed in advanced primary and metastatic melanomas [43]. Although the molecular events responsible for their activity remain to be defined, the upregulation of PTN, RGS13, and SPARC in black color skin    [44,45]. KIT encodes stem cell factor (steel factor) that is involved in stem cell differentiation and subsequent melanocytic migration [46]. Failure of melanocyte migration results in white markings on horses, including all-white coats. We did not find any differences in gene expression of KIT in the present study, though this could be due to the small sample size. However, similar global gene expression profiling of black and white sheep skin using Illumina sequencing [47] and transcriptome profiling of black and white rabbits [48] also did not find differential expression of KIT. Similar to our results, these studies did find up-regulated genes related to tyrosine metabolism and melanogenesis, including DCT in sheep, as well as keratin family genes. Zhang and coworkers [49] cataloged global gene expression profiles in Lueyang chickens with white versus black skin by Illumina2000 sequencing, and found differential expression of KIT. They also discovered upregulation of tyrosine metabolism genes in black skinned chickens, but not changes in the expression of keratin family members. Overall, these results and those of the present study indicate that while there are many common differences in gene expression between white and black skinned vertebrates, they also vary considerably by species. More work is needed to determine the specific factors involved and their mechanisms.

Conclusion
In this study, differentially expressed genes in horses with white and black coat colors were screened using high-throughput sequencing. The genes identified in this study may affect the formation of horse coat color directly or indirectly. Although the genes controlling horse coat color formation are not completely known, the transcription analysis presented in this study provides valuable information. For the first time, we present a collection of genes that are differentially expressed in horse skins of different colors. Some genes have not been described previously, and others are known, but they are likely to be involved in skin pigmentation and other physiological functions. The description of these differentially expressed genes will allow further study of the molecular regulation of horse coat color.

Availability of data and materials
We have provided detailed information about materials and method in our manuscript, so we will not provide data and supporting materials in this section.
Authors' contributions BL and XH conceived the research. QZ and DB collected the skin of Mongolian horses. BL, YZ, and XH performed RNA extraction, cDNA library construction, and Illumina sequencing. WS carried out sequence preprocessing and functional annotation, BL carried out qRT-PCR validation. BL and DM wrote the paper. All authors read and approved the final manuscript.

Ethics approval
All animal experiments were conducted in accordance with the Institutional Animal Ethics Committee and Animal Care Guidelines of the Inner Mongolia Agricultural University, which governed the use of experimental animals.

Consent for publication
Not applicable.