Identification of Differentially Expressed Genes between White and Black Skin Tissues by RNA-Seq in the Tibetan Sheep (Ovis aries)

Zhen-Yang Wu1, Li Li1, Yu-Hua Fu2, Sheng Wang3, Qing-Ming An1, Xiao-Hui Tang4, Xiao-Yong Du3,5,* and Fei Zhou6,*

1Tongren University, Tongren 554300, Guizhou, PR China

2School of Computer Science and Technology, Wuhan University of Technology, Wuhan, 430070, Hubei, PR China

3Key Laboratory of Swine Genetics and Breeding of Ministry of Agriculture and Key Laboratory of Agriculture Animal Genetics, Breeding and Reproduction of Ministry of Education, College of Animal Science, Huazhong Agricultural University, Wuhan, 430070, Hubei, PR China

4Tibet Agricultural and Animal Husbandry University, Linzhi 860000, Tibet, PR China

5Hubei Key Laboratory of Agricultural Bioinformatics, College of Informatics, Huazhong Agricultural University, Wuhan, 430070, Hubei, PR China

6School of Civil Engineering and Architecture, Hubei University of Arts and Science, Xiangyang 441053, Hubei, PR China

ABSTRACT

Tibetan sheep is one of China’s three coarse-wool sheep, the wool is mainly used for the production of Tibetan carpets and its coat color as one of important economic traits which can affect wool price directly. In this study, to investigate the functional roles of genes in Tibetan sheep skin with different coat colors, we sequenced genes from six skin samples using Solexa sequencing. The RNA-Seq analysis generated 63,283,784 and 63,644,062 clean reads in black and white skin, respectively. A total of 60 differentially expressed genes (DEGs) were identified, providing evidence that the different coat color skin changed considerably. Among these differentially expressed genes, there were 18 up-regulated genes and 42 down-regulated genes in black skin tissue. Further analysis showed that these genes involved a series of biological processes such as “pigment biosynthetic process” and “pigment metabolic process”. KEGG pathways were analyzed for the differentially expressed genes and show that Melanogenesis signaling pathway may affect the process of coat color formation. Quantitative PCR confirmed differential expression of thirteen genes, TYR, TYRP1, TYRP2, MITF, β-catenin, Wnt3a and so on. These results will expand our understanding of the complex molecular mechanisms of coat color in Tibetan sheep and provide a foundation for future studies.

Article Information

Received 31 January 2020

Revised 01 March 2020

Accepted 09 March 2020

Available online 26 February 2021

Authors’ Contribution

XYD and FZ conceived and designed the experiment. LL and XHT collected samples. SW and QMA performed some experiments. YHF helped in data analysis and preparation of some figures. ZYW analyzed the data, prepared the manuscript and provided funding.

Key words

Tibetan sheep, Coat color, RNA-sequencing, Differentially expressed genes, Q-PCR.

DOI: https://dx.doi.org/10.17582/journal.pjz/20200122010145

* Corresponding authors: duxiaoyong@mail.hzau.edu.cn; zhoufei-0521@qq.com

0030-9923/2021/0002-0721 $ 9.00/0

Copyright 2021 Zoological Society of Pakistan

Introduction

Tibetan sheep are mainly distributed in the Qinghai-Tibet Plateau and adjacent to Gansu, Qinghai and other places, due to temperature, light, precipitation, altitude and other factors, Tibetan sheep formed many different ecological types (or subtypes). The wool as one of the important economic products, different quality wool has different uses and economic value. In addition to the quality and yield of wool, coat color is one of its important economic traits.

The study of the genetic mechanism of sheep coat color began at the beginning of the 20th century, and a large number of studies have shown that the coat color of sheep is controlled by multiple alleles at multiple loci and reveals certain inheritance laws (Terrill, 1947; Singh and Singh, 1971; Vsevolodov et al., 1987). Through hybridization experiments, a part of the gene loci controlling sheep coat color was revealed. So far, 11 gene loci have been used to control the color of sheep (Renieri et al., 2008). This work constructs the basic framework for controlling the coat color genes, which has a guiding role in the control of the genetic and breeding work of the coat color traits in sheep production.

With the continuous development of molecular biology and biological technology, coat color research change gradually from the phenotypic genotype to the molecular level of the study, such as Candidate gene strategy, gene linkage analysis and Genome-wide association study. Candidate gene strategy plays an important role in the regulation of important economic traits in sheep, such as sheep growth traits (An et al., 2017; Sahu et al., 2017), reproductive traits (Barakat et al., 2017; Hao et al., 2018), disease resistance (Sayre et al., 2012; Brown et al., 2013; Periasamy et al., 2014), milk production trait (Moioli et al., 2013) and hair follicle growth and development (Guo et al., 2009). In recent years, the use of candidate gene method for the study of coat color is increasing (Davila et al., 2014; Fontanesi et al., 2014; Zhang et al., 2017; Mahmoud et al., 2017). Candidate gene method is widely studied with mammalian coat color-related candidate genes such as MC1R, ASIP, KIT, TYR, TYRP1 and so on.

A large of genetic variants associated with complex diseases have been identified by the Genome-wide association analysis (GWAS). GWAS not only contributes significantly to disease research, but also plays an important role in locating and screening important economic traits of animals. Such as fleece weight (Fatemeh et al., 2017), sheep supernumerary nipple phenotype (Peng et al., 2017), disease resistance (Atlija et al., 2016; Ćalić et al., 2017), reproduction traits (Martinez-Royo et al., 2017; Abdoli et al., 2018) and sheep hair follicle growth (Goher et al., 2010). This method was also used in animal coat color analyzed. Fan et al. (2014) analyzed the two Chinese Holstein dairy cows using bovine 50K SNP chip, suggesting that KIT, IGFBP7, PDGFRA, MITF, ING3 and WNT16 genes may be related to coat color. On the sheep, Kijas et al. (2013) obtained 49034 SNP loci by GWAS and found that most of the loci were associated with the MC1R gene. Similarly, the research of coat color is still very little, especially in the research of sheep coat color, few reports.

In recent years, the rapid development of high-throughput sequencing technology makes it possible to quickly explore potential genes that control traits. Sheep and goat skin, hair follicle tissue research is also more and more. Kang et al. (2013) used RNA-Seq technique to analyze the transcriptome of the skin tissue of 1 month and 48 months of age, respectively. It was found that there were 37 differentially expressed genes in these two different periods. Probably related pathways are Notch, VEGF, Wnt, MAPK and so on. Fan et al. (2013) performed a transcriptome sequencing of 6 neonat sheep (3 in white and black) and found that 2235 known differential expression genes were associated with coat color. But the coat color research is still inadequate. Yue et al. (2015) performed de novo transcriptome sequencing of skin and annotated 22,164 unigenes. No studies have been attempted to identify mRNAs affecting coat color in the Tibetan sheep.

In this study, we sequenced mRNAs from black and white skin collected from 1-year old crossbreed white and black coat color sheep. These results provide new information on mRNA expression profiles in the sheep and identify possible mRNA regulated pathways related to pigmentation in hair follicles.

Materials and methods

Experimental animals and sample collection

In this study, we collected the samples of Tibetan sheep in Bayi Town, Linzhi area, Tibet Autonomous Region. To eliminate the genetic background, samples were collected from three 1-year old crossbreed white and black coat colored Tibetan sheep. In the living state, cut the sheep hair with the shears and make skin tissue completely naked. Cleaned and disinfected the naked skin surface with 75% alcohol. In the use of scalpels quickly cut 1cm × 1cm size of the skin tissue after slaughter. Divided skin into smaller tissue and stored quickly in the tissue preservation solution (TaKaRa). The preservation of the organization into the ice box for temporary storage, transported to the laboratory after 4°C refrigerator storage for 12h and then placed at -80°C permanent preservation.

RNA extraction and library construction

Removed the skin from the tissue preservation solution after thaw the frozen tissue in the ice box, cut the skin tissue with scissors into powder, add 1mL Trizol reagent, the use of tissue dispersion machine for homogenization. The supernatant was transferred to a new RNAse Free 1.5 mL centrifuge tube after centrifuged at 12,000 rpm. Total RNA was isolated from skin of Tibetan sheeps according to the manufacturer’s protocol.

Quality and quantity of RNA was examined using a NanoDrop 2000/2000C (Thermo Fisher Scientific Inc., Waltham, MA, USA) and integrity was detected using agarose gel electrophoresis. Six samples (three from white skin tissue and three from black skin tissue) from three Tibetan sheep each were sent to Genergy Biotechnology Co., Ltd. (Shanghai, China) for RNA library construction and sequencing. Sequencing was performed using an Illumina HiSeq 2000 Genome Analyzer (Illumina Inc.,Santiago, CA, USA). Libraries were constructed using a TruSeq Small RNA Samples Preparation kit (Illumina Inc., Santiago, CA, USA).

Data preprocessing

Sequencing-received raw image data were transformed by base calling into sequence data (raw data). The quality of the original raw sequencing data obtained by Solexa deep sequencing was assessed using FastQC software (http://www.bioinformatics.babraha m.ac.uk/projects/fastqc/). The clean reads were obtained by trimming the low-quality reads and eliminating reads with contaminants. Prior to mapping reads to the reference database, we filtered all sequences as follows: (i) remove adaptor sequences and low-quality sequences (the percentage of low quality bases with a quality value ≤5 was >50% in a read); (ii) remove the N >10% in a read (N means that cannot determine the base information); (iii) remove the 5’ adaptor contaminated reads; (iv) remove reads that don’t have 3’ adaptor sequences and inserts; (v) remove the 3’ adaptor sequences and (vi) remove the polyA /T/G/C reads.

Mapping reads to the reference genome

In this study, Tophat 2 (Kim et al., 2013) was applied to align the clean reads to the sheep genome (the sheep genomic DNA is from the Ensemble genome database, http://www.ensembl.org/) using the algorithm was a spliced mapping. The algorithm could segment the samples without full-length matching. In addition, Tophat 2 can align the whole sequence to the genomic exon. Matching allows two mismatches, each reads allow multi hits <= 1. Two databases were employed for sequence analysis, the ensemble sheep (Ovis aries) reference gene and reference genomic DNA.

Normalized expression levels of genes and screening of DEseq

The results of the Tophat 2 alignment were quantitatively analyzed using cufflink (Trapnell et al., 2013) and the FPKM values for each genes were calculated. The mapped read counts for each gene were normalized for RNA length and for the total read number in the lane according to reads per kilobase of exon model per million mapped reads (RPKM), which facilitates the comparison of transcript levels between samples (Mortazavi et al., 2008). The Cuffdiff under Cufflinks will be used to normalize the resulting FPKM values for the analysis of the differences between the samples.

Gene ontology (GO), pathway enrichment analysis and validation of differentially expressed genes

GO functional enrichment analysis and KEGG pathway analyses were performed using DAVID software (http://www.cytoscape.org/). Differentially expressed genes were confirmed with qPCR (Chen et al., 2005). Six samples (three white skin samples and four black skin samples) from four sheep were used in qPCR analysis. The primers are shown in Table I. One microgram of total RNA from each sample was reverse-transcribed into cDNA using the Thermo Scientific Revert Aid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA). GAPDH was used as the internal control. qPCR was performed using standard protocols on the Roche LightCycler 480 Real-Time PCR Detection System (Hoffmann-La Roche Ltd, Basel, Switzerland). The 2−ΔΔCt method was used to analyze the expression levels (Livak et al., 2001).

 

Table I.- The information of gene primer.

Gene Name	Primer sequence
OVIS-TYR-F	AGCCTGTGCCTCCTCCAAGA
OVIS-TYR-R	CAGTGGTGCCGTGGACAGAA
OVIS-TYRP1-F	GATGGCAGAGATGACCGAGAGG
OVIS-TYRP1-R	GCGGCAGGTTCCACAGTTGT
OVIS-TYRP2-F	GGGCAGCAAGGCAAGACGAT
OVIS-TYRP2-R	CAAGGTGACTCGGCGGTTGTAG
OVIS-β-catenin-F	TGAGGACCTACACTTACGAGAA
OVIS-β-catenin-R	ACAGCACCTTCAGCACTCT
OVIS-KIT-F	CAACCAAGGCGGCAGACAAGA
OVIS-KIT-R	TCAGCAAGTCCTCCAGGTCCAG
OVIS-EDNRB-F	CGTCCGCAACACCGCAAGTT
OVIS-EDNRB-R	CCAGCACGAACACTAGGCAGGA
OVIS-KITL-F	ACCCAGGCAGTTTACTCC
OVIS-KITL-R	AAGCCACTATCTCCAAGTCC
OVIS-MC1R-F	CATCGCCAAGAACCGCAACCT
OVIS-MC1R-R	AGCCTCCAGCAGCAGCATGA
OVIS-SOX10-F	CGGGAAGCCTCACATCGACT
OVIS-SOX10-R	TCCAACTCGGCCACGTC
OVIS-WNT3a-F	TCGCAGCCTGACTTCCG
OVIS-WNT3a-R	CACCATCTCTGAGGCGCTGT
OVIS-MITF-F	ACAGACGGCACCATCACCTTCA
OVIS-MITF-R	TGCTCCTTCGGCTGCTGGTT
OVIS-LEF1-F	ATCCTCTGGTCCACACTCAAGG
OVIS-LEF1-R	ACATCATCAGTCACACGGTTCC
OVIS-PKA-F	CGCAACCTCCTACAAGTG
OVIS-PKA-R	TTCCTCCTCATAGTCGTCAA
OVIS-GAPDH-F	ACGGCACAGTCAAGGCAGAGA
OVIS-GAPDH-R	CGCTCCTGGAAGATGGTGATGG

 

Results

Overview of sequencing data

To identify differentially expressed genes in the two types of skin (one is white and other is black), six mRNA libraries were constructed for Solexa sequencing. A total of 140,746,212 raw reads and 126,927,846 clean reads were obtained after eliminating the low quality reads, containing N sequences and adaptor sequences. The basic information of sequencing data is also listed in Table II

.

Table II.- The information of raw data filtering.

Sample Name	Total Reads	Clean Reads	Clean Ratio	Q20	Q30
1B	22,470,265	20,266,893	90.19%	98.18%	93.79%
2B	23,752,698	21,415,951	90.16%	98.05%	93.51%
3B	23,966,791	21,600,940	90.13%	98.10%	93.64%
1W	20,643,079	18,624,824	90.22%	97.96%	93.10%
2W	26,010,571	23,463,755	90.21%	97.97%	93.15%
3W	23,902,808	21,555,483	90.18%	97.97%	93.14%

 

Table III.- The most abundantly expressed gene in Tibetan sheep skin tissues (top 10).

Gene Names	Black skin	White skin	Sum
MT-ATP8	65192.9	130385.8	126843.1
KRTAP8-1	56621.7	113243.4	102791.2
KAP3.2	45360.2	90720.4	83371.3
COX1	30833.7	61667.46	60610.68
CST6	23010.4	46020.72	48337.67
KAP7	24826.9	49653.78	45179.5
MT-ATP6	19275	38550.06	37869.16
ZNFX1-AS1_2	18662.9	37325.78	37545.15
THYMB4X	18393.5	36787	36782.3
PRR9	18988.1	37976.12	36247.96

 

Gene expression levels analysis

Genetic quantification of the results of tophat alignment was performed using cufflink (VERSION: 2.1.1) and the FPKM values for each pair of genes were calculated. A total of 16,709 annotated genes were obtained. We list the top 10 genes with high expression levels in Table III. Found that these genes are ATP synthesis process, keratin and keratin associated protein, membrane-associated glycoprotein and other related.

Cuffdiff was used to analyze the normalized FPKM values between the samples. A total of 60 differentially expressed genes were obtained. Table IV lists only the top 10 up-regulated and the first 9 down-regulated genes.

Clustering of DEGs between the samples

To further analyze the differentially expressed genes, we do the heatmap of the differentially expressed genes (Fig. 1). It is convenient to see that there is no significant change between the two skin tissues (white skin and black skin). The clustering of the sample found that a black sheep skin sample and two other white skin tissue samples were clustered together to further indicate that there may be differences between individuals.

 

Table IV.- The information of differentially expressed gene between white skin and black skin in Tibetan sheep.

Gene names	White skin	Black skin	log2 (fold change)	p value
SNORD49	295.4	0	inf	0.00005
RETN	3.18	0	inf	0.00005
ZBP1	8.6	1.09	2.97	0.00005
MX2	41.91	5.36	2.97	0.00005
BST-2A	50.22	6.48	2.95	0.00005
HOXB6	10.16	1.53	2.74	0.00005
BST-2B	17.49	2.81	2.64	0.00005
ISG17	174.81	29.07	2.59	0.00005
RSAD2	25.75	4.44	2.53	0.00005
HOXB7	8.4	1.61	2.38	0.00005
PMEL	1.66	151.65	-6.51	0.00005
TYRP1	0.6	38.8	-6.02	0.00005
DCT	2.82	19.82	-2.81	0.00005
TYR	1.59	8.05	-2.34	0.00005
STMN2	2.37	11.52	-2.28	0.0001
VNN1	0.74	2.97	-2	0.00005
GPR143	3.4	13.08	-1.94	0.00005
OCA2	1.64	5.74	-1.81	0.00005
IL37	7.27	21.56	-1.57	0.00005

 

 

 

Analysis of differential gene function

Genotyping of the differential genes obtained by cuffdiff was performed using the DAVID software for gene function analysis. We listed biological process (Fig. 2; Supplementary Table I), cellular component (Fig. 2; Supplementary Table II) and molecular function (Fig. 2; Supplementary Table III), respectively. Although there are few differences in the expression of genes, but the results show that these genes in addition to “immune response”, “defense response” and “response to virus” and other biological processes related, there are many genes related to pigment deposition, such as “regulation of cell proliferation”, “pigment biosynthetic process”, “pigment metabolic process” and “pigment metabolic process”. In “cellular component” and “molecular function”, there are also genes associated with pigment deposition, such as “pigment granule” and “melanosome”.

The differentially expressed genes were analyzed by DAVID software to analyze the KEGG pathways of these differential genes. Only three KEGG pathways were obtained, namely “RIG-I-like receptor signaling pathway”, “Tyrosine metabolism” and “Melanogenesis”. The Melanogenesis pathway only covers three genes, TYR, TYRP1 and TYRP2.

Confirmation of differential gene expression by qPCR

We selected the most important three genes, TYR, TYRP1 and TYRP2 genes, to validate the expression profiles obtained by RNA-Seq. The qPCR results verified that these genes were differentially expressed at different color skins, consistent with the RNA-Seq findings (Fig. 3). Thus, RNA-Seq can provide reliable data for mRNA differential expression analysis from one side as with other similar findings.

 

 

We also selected candidate genes associated with melanocyte activity and melanin synthesis pathway for qPCR analysis. The genes included MC1R, KIT, KITL, MITF, WNT3a, SOX10, PKA, PAX3, EDNRB, β-catenin, and LEF1. Relatively high expression levels include PKA, MITF, EDNRB, Kit, LEF1 andβ-catenin (Fig. 4). Among these genes, KIT and EDNRB showed extremely significant differences (P<0.01), and β-catenin showed significant differences (P<0.05), while other differences were not significant. KIT is considered to be an important factor in the migration of melanocytes’ precursors, which can affect the activity of TYR in the later stage. MAPK and other signaling pathways can also be activated. However, in this study there was no difference in MITF gene, and a large number of studies have shown that MITF is an important transcription factor in the process of melanin synthesis, which can promote the transcription of tyrosine gene family. Therefore, it was speculated that the upstream genes with inhibitory effect on TYR might influence its expression in other ways.

The genes with relatively low expression in skin tissues are MC1R, SOX10, WNT3a and KITL. SOX10 (P< 0.01) was the most significant difference, while WNT3a and MC1R (P<0.05) were the most significant differences (Fig. 5). But based on differences in expression levels and multiples, these genes may not play a substantial role. Therefore, more genes need to be discovered to explain the formation mechanism of this complex trait.

Discussion

In this study, we used the mRNA-Seq sequencing technique to sequence the skin tissues of three crossed black and white sheep, and obtain 70,189,754 and 70,556,458 reads in white and black skin, respectively. These samples were compared to the sheep genome and annotated, and then differentially expressed genes were analyzed by cuffdiff, and 60 differentially expressed genes were obtained. Among these genes, there were 18 up-regulated genes and 56 down-regulated genes in black skin tissues. Further analysis showed that these genes involved a series of biological processes such as “pigment biosynthetic process”, “pigment metabolic process” and “pigment metabolic process”. Fan et al. (2013) used a high-throughput sequencing technique to analyze the skin tissues of 6 Sunit ewes (3 black and 3 white), and received a total of 2,235 differentially expressed genes, of which 479 genes were in black Skin uptake expression, among these genes, there are 49 known control of the coat color genes. The expression of differentially expressed genes in this study was not high, but there was a similar result to Fan’s work. We found that the distribution of data between the three groups may be different; it is possible to average the level of gene expression between the samples and impact the experimental results. Although the individual differences are large, but the results are still pointing to the sheep control the coat color by regulating the number of melanin. Melanin is produced by mature melanocytes, so there are two possible ways to reduce the production of melanin, one is through the regulation of melanoma stem cells into mature melanoma cells to achieve the reduction of the number of melanin, and other is the direct control of tyrosinase activity to control the amount of melanin (Yang et al., 2016). The mature melanocytes of the skin and hair follicles come from the melanocytic stem cells of the hair follicle, and the melanocytic cells in the hair follicles also change periodically as the hair follicles grow periodically, thus requiring the differentiation of melanocytes to supplement. The mature melanocytes of the skin and hair follicles come from the melanocytic stem cells of the hair follicle, and the melanocytic cells in the hair follicles also change periodically as the hair follicles grow periodically, thus requiring the differentiation of melanocytes to supplement. Melanocytes in the skin tissue are differentiated from the embryonic period of the neural crest cells, when the neural crest cell differentiation into melanocytes to the epidermal basal layer, the differentiation into mature melanoma cells, these melanocytes are generally not affected by hair follicle cycle the effect of growth changes has been stable in the basal layer of the skin (Kelsh et al., 2009; Mayo et al., 2014). But not all species of skin melanoma cells are so, so the epidermal melanoma cells and hair follicle melanoma cell regulation mechanism is a certain difference. We hypothesize that the reduction in the number of melanin in white hair follicles is likely to result in the lack of mature melanocytes in the hair follicles by controlling the differentiation of melanocyte stem cells. And skin tissue, it may be through the control of tyrosinase activity or the two controls caused by the combined effect. So it is necessary to carry out more in-depth study.

The expression levels of TYR, TYRP1 and TYRP2 were significant in the different skin tissues, and the expression level of TYR, TYRP1 and TYRP2 was very high in black skin tissue. The expression of TYR gene and TYRP1 gene in black skin was 46 times and 450 times higher than that of white skin, while that of TYRP2 was only 11 times. These three genes are melanoma cell-specific protein family members, respectively, catalyze the different stages of melanin biosynthesis. A large number of studies have confirmed that mutations in TYR can lead to whitening of many species such as sheep, mice, rabbit and so on (Song et al., 2017; Xing et al., 2015; Ming et al., 2016). The TYRP1 gene is not well known at present and is located in the brown locus of mice and is also known as the Brown gene. Through the study of the function of TYRP1 gene, it was proved that the protein had DHICA oxidase activity, which played an important role in the downstream pathway of melanin synthesis. In addition, studies have shown that TYRP1 in addition to control the synthesis of pigment process, but also on the stability of tyrosinase and to maintain the ultrastructure of melanoma cells and melanoma cell proliferation and differentiation plays an important role (Hrckova et al., 2017; Li and Li, 2015). The TYRP2 function has been identified as enabling the rapid incorporation of carboxylic acid precursors DHICA into biosynthetic melanin, which is important for reducing the cytotoxicity of these intermediates (Jiang et al., 2016; Andrzej et al., 2016; Nie et al., 2016). This study shows that these three genes are significantly different, especially the expression of TYR and TYRP1 genes. It is speculated that the difference between black and white coat color, whether in the melanin synthesis process, or in the process of proliferation and differentiation of melanoma cells are subject to varying degrees of block.

In black and white skin tissues of Tibetan sheep, β-catenin and Wnt3a were highly expressed in black skin, in which β-catenin in black skin was significantly higher than that in white skin (P<0.01), while Wnt3a in black skin was significantly higher than that in white skin (P<0.05), suggesting that the formation of different hair colors of Tibetan sheep may be related to the Wnt/β-catenin signaling pathway. The study of Yu et al. (2010) in brown and white alpacas also showed that both mRNA expression and protein expression of β-catenin were significantly higher in brown skin than in white skin. Song et al. (2012) also used qPCR and Western blot analysis in brown and white alpaca skin and found significant differences in Wnt3a in skin tissues of different colors. Wnt3a is the ligand of the classical WNT signaling pathway. Both Wnt3a and the activated WNT/β-catenin signaling pathway protein are expressed in the hair follicle melanocytes during the growth period, so they play an important role in the proliferation, differentiation and pigmentation of the hair follicle melanocytes. Studies in mouse melanin stem cells showed that Wnt3a activated the Wnt/β-catenin signaling pathway, thus promoting the differentiation of McSCs in this process, and Wnt3a induced hair follicle regeneration (Guo et al., 2016). A large number of studies have confirmed the important role of the Wnt/β-catenin signaling pathway in the differentiation, proliferation and melanin production of melanocytes (Zhao et al., 2014; Chung-Hsing et al., 2014; Guo et al., 2012). Among these genes, β-catenin is an important adhesion molecule in epithelial tissue and a widely distributed multifunctional protein involved in various developmental processes including neural crest. β-catenin also regulates the expression of MITF and can work with it to regulate the transcription of downstream genes (Cheong-Yong et al., 2015). MITF gene was first discovered in mice in 1942, located on chromosome 6 and able to bind to relevant downstream targets to regulate the expression of downstream proteins, thus affecting the proliferation, survival and pigmentation of melanocytes (Liu et al., 2017). However, the results showed that there was no significant difference in the expression of MITF gene in the skin tissues of Tibetan sheep with different colors (P>0.05), indicating that MITF gene might not be involved in the hair color regulation of Tibetan sheep. Han et al. (2013) came to a similar conclusion, pointing out that the mRNA expression level of MITF gene was significantly higher in the black skin tissues around the orbit of Tibetan sheep than in the brown-black skin tissues of the neck and the white skin tissues of the body side, while the mRNA expression level of MITF gene in the brown-black skin tissues of the neck and the white skin tissues of the body side was not significantly different. MITF promoter region contains binding regions of transcription regulatory factors such as SOX10 and LEF1, so it is also regulated by a variety of signal molecules. In the experiment, the expression of LEF1 gene in different hair colors was also not significantly different, while SOX10 gene was differentially expressed in skin tissues of different hair colors, but the expression level was very low. Therefore, it can be speculated that the hair color regulation of Tibetan sheep may be related to Wnt/β-catenin signaling pathway, but not to MITF gene and its related transcription factors. MC1R can activate intracellular adenylate cyclase system, and then activate the tyrosine kinase, tyrosine (TYR) involved in the synthesis of melanin. Zhang et al. (2019) transfected MC1R plasmid of alpaca into sheep melanocytes, and not only found that mRNA and protein expression levels of relevant downstream genes such as TYR and TYRP1 in sheep melanocytes were significantly up-regulated, but melanin synthesis was also significantly increased. Tang et al. (2018) also found that MC1R gene expression in dark skin tissues was higher than that in white skin tissues in different color yaks. The expression of MC1R gene in black skin tissues of Tibetan sheep is significantly higher than that in white skin, indicating that MC1R gene can increase the content of melanin by up-regulating downstream genes such as TYR in the process of melanin synthesis, thus playing a regulatory role in hair color.

The formation of hair color is not only related to the production of melanin, but also to the migration, distribution and proliferation of melanocytes. The KIT gene encodes mast/stem cell growth factor receptors, which together with their ligand (KITL) play a critical role in the survival and maintenance of hematopoietic stem cells and mast cells. It also affects the proliferation, differentiation and migration of melanocytes, so it plays a crucial role in the formation of melanin. Among them, abnormal expression of KIT gene can directly affect the number of melanocytes, which is also one of the reasons for the formation of white fur (Anello et al., 2019). The expression level of KIT gene in black skin of Tibetan sheep was significantly higher than that in white skin, while the expression level of KITL was not significantly different between the two skin tissues. The results showed that KIT gene plays an important role in the proliferation of melanocytes in Tibetan sheep skin. In addition to KIT gene, EDNRB gene was also highly expressed in the black skin tissue of Tibetan sheep and was highly expressed in the white skin. This is consistent with the research of Geng et al. (2010). The translation products of EDNRB gene and functional signals of EDN3 are necessary for the development of pigment cells derived from neural crest, which play a role in signal transmission during the development of neural crest derived melanocytes and are important candidate genes related to hair color (Liu et al., 2018).

Conclusion

In conclusion, the development of melanocytes is a complex event involving numerous genes and pathways that interact and show cross-talk. These studies will provide new insights into the regulation of coat color formation for understanding the regulatory mechanism at the post-transcription level.

Acknowledgments

This study was supported by the National Nature Scientific Foundation of China (Grant No. 31702097), the National Nature Scientific Foundation of China (Grant No. 31872978), Major innovation group project of education department of Guizhou province (No. (2016) 052), Project of science and technology department of Guizhou province (No. (2016) 7311) and National Key R&D Program of China (No. 2018YFD0502000).

Supplementary material

There is supplementary material associated with this article. Access the material online at: https://dx.doi.org/10.17582/journal.pjz/20200122010145

Statement of conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this article.

Reference

Abdoli, R., Mirhoseini, S.Z., Hossein-Zadeh, N.G. and Zamani, P., 2018. Screening for causative mutations of major prolificacy genes in Iranian fat-tailed sheep. Int. J. Fertil. Steril., 12: 51-55.

An, Q., Zhou, H., Hu, J., Luo, Y. and Jgh, H., 2017. Haplotypes of the ovine adiponectin gene and their association with growth and carcass traits in New Zealand romney lambs. Genes, 8: 160. https://doi.org/10.3390/genes8060160

Andrzej, J., Magdalena, G., Beata, H., Marek, K., Kornel, K., Katarzyna, G. and Grazyna, J.W., 2016. Single-nucleotide polymorphism of MC1R, ASIP, and TYRP2 genes in wild and farmed foxes (Vulpes vulpes). Can. J. Anim. Sci., 96: 172-179. https://doi.org/10.1139/cjas-2015-0066

Anello, M., Daverio, M.S., Silbestro, M.B., Vidal-Rioja, L. and Di, R.F., 2019. Characterization and expression analysis of KIT and MITF-M genes in llamas and their relation to white coat color. Anim. Genet., 50: 143-149. https://doi.org/10.1111/age.12769

Atlija, M., Arranz, J.J., Martinezvalladares, M. and Gutiérrezgil, B., 2016. Detection and replication of QTL underlying resistance to gastrointestinal nematodes in adult sheep using the ovine 50K SNP array. Genet. Sel. Evol., 48: 1-16. https://doi.org/10.1186/s12711-016-0182-4

Barakat, I.A.H., Salem, L.M., Daoud, N.M., Khalil, W.K.B. and Mahrous, K.F., 2017. Genetic polymorphism of candidate genes for fecundity traits in Egyptian sheep breeds. Biomed. Res., 28: 851-857.

Brown, E.A., Pilkington, J.G., Nussey, D.H., Watt, K.A., Hayward, A.D., Tucker, R., Graham, A.L., Paterson, S., Beraldi, D., Pemberton, J.M. and Slate, J., 2013. Detecting genes for variation in parasite burden and immunological traits in a wild population: testing the candidate gene approach. Mol. Ecol., 22: 757-773. https://doi.org/10.1111/j.1365-294X.2012.05757.x

Ćalić, I., Koch, J., Carey, D., Addo-Quaye, C. and Carlson, J.E., 2017. Genome-wide association study identifies a major gene for beech bark disease resistance in American beech (Fagus grandifolia Ehrh.). BMC Genom., 18: 547. https://doi.org/10.1186/s12864-017-3931-z

Chen, C., Ridzon, D.A., Broomer, A.J., Zhou, Z., Lee, D.H., Nguyen, J.T., Barbisin, M., Xu, N.L., Mahuvakar, V.R. and Andersen, M.R., 2005. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucl. Acids Res., 33: e179. https://doi.org/10.1093/nar/gni178

Cheong-Yong, Y., Soon-Tae, Y., Jin-Hwa, K., Jin, H.C., Sang, B.H., Eun, Y.S. and Eung, K., 2015. p21-Activated Kinase 4 critically regulates melanogenesis via activation of the CREB/MITF and β-Catenin/MITF pathways. J. Incest. Dermatol., 135: 1385-1394. https://doi.org/10.1038/jid.2014.548

Chung-Hsing, C., Rong-Kung, T., Ming-Hsien, T., Yi-Hsiung, L. and Tomohisa, H., 2014. The roles of Frizzled-3 and Wnt3a on melanocyte development: In vitro studies on neural crest cells and melanocyte precursor cell lines. J. Dermatol. Sci., 75: 100-108. https://doi.org/10.1016/j.jdermsci.2014.04.012

Davila, S.G., Gil, M.G., Resino-Talavan, P. and Campo, J.L., 2014. Association between polymorphism in the melanocortin 1 receptor gene and E locus plumage color phenotype. Poult. Sci., 93: 1089-1096. https://doi.org/10.3382/ps.2013-03611

Fan, R.W., Xie, J.S., Bai, J.M., Wang, H.D., Tian, X., Bai, R., Jia, X.Y., Yang, L., Song, Y.F., Herrid, M., Gao, W.J., He, X.Y., Yao, J.B., Smith, G.W. and Dong, C.S., 2013. Skin transcriptome profiles associated with coat color in sheep. BMC Genom., 14: 389. https://doi.org/10.1186/1471-2164-14-389

Fan, Y., Wang, P., Fu, W., Dong, T. and Qi, C., 2014. Genome-wide association study for pigmentation traits in Chinese Holstein population. Anim. Genet., 45: 740. https://doi.org/10.1111/age.12189

Fatemeh, E., Mohsen, G., Ghodrat, R.M. and Ayoub, F., 2017. Detection of QTL for greasy fleece weight in sheep using a 50 K single nucleotide polymorphism chip. Trop. Anim. Hlth. Prod., 49: 1657-1662. https://doi.org/10.1007/s11250-017-1373-x

Fontanesi, L., Vargiolu, M., Scotti, E., Latorre, R., Pellegrini, M.S.F., Mazzoni, M., Asti, M., Chiocchetti, R., Romeo, G., Clavenzani, P. and De, G.R., 2014. The KIT gene is associated with the English spotting coat color locus and congenital megacolon in checkered giant rabbits (Oryctolagus cuniculus). PLoS One, 9: e93750. https://doi.org/10.1371/journal.pone.0093750

Geng, J.J., Sun, L.T., Mu, X.L., Zhang, J., Jiang, J.B., Zhang, Y.L., H.Q. and Dong, C.S., 2010. Immunolocalization and expression analysis of endothel in receptor in sheep skin of different hair color. Sci. Agric. Sin., 43: 5147-5154.

Goher, M., Rauw, W., Thin, D. and Gomez-Raya, L., 2010. Selective genotyping using genome-wide association studies (GWAS) that are associated with fiber diameter in Merino sheep. J. Dairy Sci., 93: 169-169.

Guo, C.H., Jia, C.L., Zhang, W., Zhu, X.P. and Jia, Z.H., 2009. Polymorphisms in the 5’ flanking region of IGF-I gene are associated with cashmere fibre traits in liaoning cashmere goats. Progr. Biochem. Biophys., 36: 566-573. https://doi.org/10.3724/SP.J.1206.2008.00627

Guo, H.Y., Yang, K.Y., Deng, F., Ye, J.X., Xing, Y.Z., Li, Y.H., Lian, X.H. and Yang, T., 2012. Wnt3a promotes melanin synthesis of mouse hair follicle melanocytes. Biochem. Biophy. Res. Co., 420: 799-804. https://doi.org/10.1016/j.bbrc.2012.03.077

Guo, H.Y., Xing, Y.Z., Liu, Y.X., Luo, Y., Deng, F., Yang, T. and Yang, K., 2016. Wnt/β-catenin signaling pathway activates melanocyte stem cells in vitro and in vivo. J. Dermatol. Sci., 83: 45-51. https://doi.org/10.1016/j.jdermsci.2016.04.005

Han, J.L., Yue, Y.J., Pei, J., Guo, T.T., Liu, J.B., Guo, J., Sun, X.P., Feng, R.L. and Yang, B.H., 2013. The polymorphism of agouti signaling protein gene (Agouti) and the expression of agouti and microphthalmia-associated transcription factor gene MITF in different coat color skin tissues of Tibetan sheep (Ovis aries). J. agric. Biotechnol., 21: 338-347.

Hao, Z.Y., Wang, J.Q., Hu, J., Liu, X., Li, S.B., Zhang, Y.Q., Wang, G. and Luo, Y.Z., 2018. The tissue expression of the sheep (Ovis aries) STAT5a gene and the analysis of nucleotide sequence variation. J. agric. Biotechnol., 26: 96-103.

Hrckova, T.E., Majchrakova, Z., Bielikova, M., Soltys, K., Turna, J. and Dudas, A., 2017. A novel mutation in the TYRP1 gene associated with brown coat colour in the Australian shepherd dog breed. Anim. Genet., 48: 626. https://doi.org/10.1111/age.12563

Jiang, S., Yu, X. and Dong, C., 2016. MiR-137 affects melanin synthesis in mouse melanocyte by repressing the expression of c-Kit and Tyrp2 in SCF/c-Kit signaling pathway. Biosci. Biotechnol. Biochem., 21: 1-7. https://doi.org/10.1080/09168451.2016.1200455

Kang, X.L., Liu, G., Liu, Y.F., Xu, Q.Q., Zhang, M. and Fang, M.Y., 2013. Transcriptome profile at different physiological stages reveals potential mode for curly fleece in Chinese Tan sheep. PLoS One, 8: 111. https://doi.org/10.1371/journal.pone.0071763

Kelsh, R.N., Harris, M.L., Colanesi, S. and Erickson, C.A., 2009. Stripes and belly-spots-A review of pigment cell morphogenesis in vertebrates. Semin Cell Dev. Biol., 20: 90-104. https://doi.org/10.1016/j.semcdb.2008.10.001

Kijas, J.W., Serrano, M., McCulloch, R., Li, Y., Ortiz, J.S., Calvo, J.H. and Perez-Guzman, M.D., 2013. Genomewide association for a dominant pigmentation gene in sheep. J. Anim. Breed. Genet., 130: 468-475. https://doi.org/10.1111/jbg.12048

Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R. and Salzberg, S.L., 2013. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genom. Biol., 14: R36. https://doi.org/10.1186/gb-2013-14-4-r36

Li, L. and Li, X., 2015. Research progress of animal coat color candidate gene in tyrosinase related protein 1(TYRP1). J. Hebei Norm. Univ. Sci. Technol., 4: 2.

Liu, W., Li, Y., Cai Y.Q., Zhang, L.H., Li, H.F., Li, H. and Zhu, Z.W., 2017. Research progress of the regulators involved in MITF-related signaling pathways in melanocytes. Anim. Husband. Vet. Med., 49: 125-128.

Liu, Y., Xue, L.L., Chang, X.Y., Bai, Y., Li, L.X., Li, J.W., Yu, X.J., He, X.Y. and Wang, H.D., 2018. Effects of EDN3 on the expression of melanin-related genes in mouse melanocytes. China Anim. Husband. Vet. Med., 45: 581-589.

Livak, K.J. and Schmittgen, T.D., 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCt method. Methods, 25: 402-408. https://doi.org/10.1006/meth.2001.1262

Mahdavi, M., Nanekarani, S. and Hosseini, S.D., 2014. Mutation in BMPR-IB gene is associated with litter size in Iranian Kalehkoohi sheep. Anim. Reprod. Sci., 147: 93-98. https://doi.org/10.1016/j.anireprosci.2014.04.003

Mahmoud, A.H., Mashaly, A.M., Rady, A.M., Alanazi, K.M. and Saleh, A.A., 2017. Allelic variation of Melanocortin-1 receptor (MC1R) locus in Saudi indigenous sheep exhibiting different color coats. Asian-Austral. J. Anim., 30: 154. https://doi.org/10.5713/ajas.16.0138

Martinez-Royo, A., Alabart, J.L., Sarto, P., Serrano, M. and Lahoz, B., 2017. Genome-wide association studies for reproductive seasonality traits in Rasa Aragonesa sheep breed. Theriogenology, 99: 21-29. https://doi.org/10.1016/j.theriogenology.2017.05.011

Mayo, V., Sawatari, Y., Huang, C.Y.C. and Garcia-Godoy, F., 2014. Neural crest-derived dental stem cells-Where we are and where we are going? J. Dent., 42: 1043-1051. https://doi.org/10.1016/j.jdent.2014.04.007

Ming, L., Yi, L., Sa, R., Ji, R. and Ha, S., 2016. Polymorphisms of the tyrosinase (TYR) gene in Bactrian camel (Camelus bactrianus) with different coat colour. J. Camel Pract. Res., 23: 47-51. https://doi.org/10.5958/2277-8934.2016.00007.2

Moioli, B., Scata, M.C., De, M.G., Annicchiarico, G., Catillo, G. and Napolitano, F., 2013. The ACACA gene is a potential candidate gene for fat content in sheep milk. Anim. Genet., 44: 601-603. https://doi.org/10.1111/age.12036

Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L. and Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods, 5: 621-628. https://doi.org/10.1038/nmeth.1226

Nie, R.Q., Yang, Y.J., Xie, J.S., Fan, R.W., Xu, D.M., Yu, X.J. and Duan, Z.C., 2016. Influences of Pax6 PAI subdomain on MITF, TYR, TYRP1 and TYRP2 in melanocytes. Sci. Agric. Sin., 49: 3433-3442.

Peng, W.F., Xu, S.S., Ren, X., Lv, F.H. and Xie, X.L., 2017. A genome-wide association study reveals candidate genes for the supernumerary nipple phenotype in sheep (Ovis aries). Anim. Genet., 48: 570-579. https://doi.org/10.1111/age.12575

Periasamy, K., Pichler, R., Poli, M., Cristel, S., Cetra, B., Medus, D., Basar, M., Thiruvenkadan, A.K., Ramasamy, S., Ellahi, M.B., Mohammed, F., Teneva, A., Shamsuddin, M., Podesta, M.G. and Diallo, A., 2014. Candidate gene approach for parasite resistance in sheep - variation in immune pathway genes and association with fecal egg count. PLoS One, 9: e88337. https://doi.org/10.1371/journal.pone.0088337

Renieri, C., Valbonesi, A., La, M.V., Antonini, M. and Lauvergne, J.J., 2008. Inheritance of coat color in Merino sheep. Small Rumin. Res., 74: 23-29. https://doi.org/10.1016/j.smallrumres.2007.03.001

Sahu, A.R., Jeichitra, V., Rajendran, R. and Raja, A., 2017. Polymorphism in exon 3 of myostatin (MSTN) gene and its association with growth traits in Indian sheep breeds. Small Rumin. Res., 149: 81-84. https://doi.org/10.1016/j.smallrumres.2017.01.009

Sayre, B.L. and Harris, G.C., 2012. Systems genetics approach reveals candidate genes for parasite resistance from quantitative trait loci studies in agricultural species. Anim. Genet., 43: 190-198. https://doi.org/10.1111/j.1365-2052.2011.02231.x

Singh, L.B. and Singh, O.N., 1971. Inheritance of colour in Chokla and crossbred sheep. Indian Vet. J., 48: 260-266.

Song, Y., Xu, Y., Deng, J., Chen, M., Lu, Y., Wang, Y., Yao, H., Zhou, L., Liu, Z., Lai, L. and Li, Z., 2017. CRISPR/Cas9-mediated mutation of tyrosinase (Tyr) 3′UTR induce graying in rabbit. Scient. Rep., 7: 1569. https://doi.org/10.1038/s41598-017-01727-y

Song, Y.F., Tian, X., Lu, X.X., Guo, Q.Y., Zhang, D.L., Xie, X.J., Yu, X.J., He, J.P., Wang, H.D. and Dong, C.S., 2012. Expression and localization of Wnt3a in different colors of Alpaca skin. Acta Vet. Zootech. Sin., 43: 785-790.

Tang, P., Ling, X.X., Gao, Z.C., Jia, C.J., Liang, C.N., Wu, X.Y., Chu, M. and Yan, P., 2018. Skin histological observation of Yak with different coat colors and MC1R functional verification. Acta Vet. Zootech. Sin., 49: 1377-1386.

Terrill, C.E., 1947. Color on the legs of sheep; its inheritance in the Columbia and Targhee breeds. J. Heredity, 38: 89-92. https://doi.org/10.1093/oxfordjournals.jhered.a105697

Trapnell, C., Hendrickson, D.G., Sauvageau, M., Goff, L., Rinn, J.L. and Pachter, L., 2013. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat. Biotechnol., 31: 46. https://doi.org/10.1038/nbt.2450

Vsevolodov, E.B., Adalsteinsson, S. and Ryder, M.L., 1987. Electron spin resonance spectrometrical study of the melanins in the wool of some North European sheep in relation to their color inheritance. J. Heredity, 78: 120-122. https://doi.org/10.1093/oxfordjournals.jhered.a110330

Xing, S.Y., Yang, T.A., Guo, Q. and Yang, F.H., 2015. Correlation analysis of TYR gene mutation and coat color of mustela vison. J. Jilin Agric. Univ., 37: 232-236.

Yang, S., Zhang, J., Ji, K., Jiao, D. and Fan, R., 2016. Characterization and expression of soluble guanylate cyclase in skins and melanocytes of sheep. Acta Histochem., 118: 219. https://doi.org/10.1016/j.acthis.2016.01.002

Yu, X.J., Dong, C.S., Fan, K.H., He, J.P., Jiang, J.B., Zhu, Z.W. and Dong, Y.J., 2010. Study on the expression and localization of β-catenin in Alpaca skin of different colors. Acta Vet. Zootech. Sin., 41: 335-340.

Yue, Y.J., Liu, J.B., Yang, M., Han, J.L. and Guo T.T., 2015. De novo assembly and characterization of skin transcriptome using RNAseq in sheep (Ovis aries). Genet. Mol. Res., 14: 1371-1384. https://doi.org/10.4238/2015.February.13.16

Zhang, X., Li, W., Liu, C., Peng, X. and Lin, J., 2017. Alteration of sheep coat color pattern by disruption of ASIP gene via CRISPR Cas9. Scient. Rep., 7: 8149. https://doi.org/10.1038/s41598-017-08636-0

Zhang, Y.N., Zhang, Q.Y., Ren, Y.H. and Fan, R.W., 2019. Expression of Alp-MC1R in sheep melanocyte stimulates melanin production. Chinese J. Biochem. mol. Biol., 35: 901-906.

Zhao, J.R., Bai, X.F., Yang, K., Luo, Y. and Guo, H.Y., 2014. Expression of Wnt3a in hair follicle and melanocyte lineage. Progr. Modern Biomed., 14: 5814-5817.