
==== Front
BMC Genomics
BMC Genomics
BMC Genomics
1471-2164
BioMed Central London

10772
10.1186/s12864-024-10772-7
Research
Genomic and transcriptomic landscape to decipher the genetic basis of hyperpigmentation in Lanping black-boned sheep (Ovis aries)
Chong Yuqing 1
Xiong Heli 12
Gao Zhendong 1
Lu Ying 1
Hong Jieyun 1
Wu Jiao 1
He Xiaoming 1
Xi Dongmei 1
Tu Xiaolong tuxiaolong@mail.kiz.ac.cn

34
Deng Weidong dengwd@ynau.edu.cn

1
1 https://ror.org/04dpa3g90 grid.410696.c 0000 0004 1761 2898 Yunnan Provincial Key Laboratory of Animal Nutrition and Feed, Faculty of Animal Science and Technology, Yunnan Agricultural University, Kunming, 650201 China
2 https://ror.org/010paq956 grid.464487.d Yunnan Academy of Animal Husbandry and Veterinary Sciences, Kunming, 650224 China
3 grid.9227.e 0000000119573309 State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming, 650204 China
4 https://ror.org/05qbk4x57 grid.410726.6 0000 0004 1797 8419 Kunming College of Life Science, University of the Chinese Academy of Sciences, Kunming, 650204 China
9 9 2024
9 9 2024
2024
25 84528 11 2023
4 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Background

Lanping black-boned sheep (LPB) represent a distinctive mammalian species characterized by hyperpigmentation, resulting in black bone and muscle features, in contrast to their conventional counterparts exhibiting red muscle and white bone. The genetic basis underlying LPB hyperpigmentation has remained enigmatic.

Methods

In this study, we conducted whole-genome sequencing of 100 LPB and 50 Lanping normal sheep (LPN), and integrated this data with 421 sequenced datasets from wild and domestic sheep, shedding light on the genetic backdrop and genomic variations associated with LPB. Furthermore, we performed comparative RNA-Seq analysis using liver sample to pinpoint genes implicated in the pigmentation process. We generated a comprehensive dataset comprising 97,944,357 SNPs from 571 sheep, facilitating an in-depth exploration of genetic factors.

Results

Population genetic structure analysis revealed that the LPB breed traces its origin back to LPN, having evolved into a distinct breed. The integration of positively selected genes with differentially expressed genes identified two candidates, ERBB4 and ROR1, potentially linked to LPB hyperpigmentation. Comparative analysis of ERBB4 and ROR1 mRNA relative expression levels in liver, spleen, and kidney tissues of LPB, in comparison to Diqing sheep, revealed significant upregulation, except for ERBB4 in the liver. Gene expression heatmaps further underscored marked allelic frequency disparities in different populations.

Conclusion

Our findings establish the evolutionary lineage of the LPB breed from LPN and underscore the involvement of ERBB4 and ROR1 genes in melanin synthesis. These results enhance our comprehension of the molecular basis of hyperpigmentation and contribute to a more comprehensive depiction of sheep diversity.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12864-024-10772-7.

Keywords

Lanping black-boned sheep (LPB)
Hyperpigmentation
Whole-genome sequence
RNA-Seq
Integrative analysis
http://dx.doi.org/10.13039/501100001809 National Natural Science Foundation of China 32302707 32160771 Yunnan Province Basic Research Program Project202301AU070117 http://dx.doi.org/10.13039/501100019413 Yunnan Provincial Department of Education Science Research Fund Project 2023J0392 Yunnan Provincial Agricultural Union Foundation202101BD070001-006 http://dx.doi.org/10.13039/501100018531 Major Science and Technology Projects in Yunnan Province 202202AE090005 Yunnan Provincial Academician and Expert Workstation202105AF150043 issue-copyright-statement© BioMed Central Ltd., part of Springer Nature 2024
==== Body
pmcIntroduction

Hyperpigmentation in animals has long piqued the curiosity of evolutionary and developmental biologists, offering valuable insights into melanocyte development and hyperpigmentation-related disease research [1–3]. Recently discovered in Lanping County, China, Lanping black-boned sheep (LPB) exhibit striking hyperpigmentation in their skin and internal tissues, resembling the world-renowned silky fowl [4]. They are unique among mammals, displaying extensive pigmentation throughout their bodies, including the eye conjunctiva, muscles, periosteum, kidney, heart, lung, liver, and trachea. This distinctive trait makes them highly prized for their superior meat quality, particularly appealing to oriental preferences. LPB are found emerging from Lanping normal sheep (LPN), which lacks hyperpigmentation in its skin and internal organs, like most sheep breeds. However, the exact genetic lineage of LPB, its relationship with LPN, and the underlying genetic mechanisms responsible for hyperpigmentation in LPB remain poorly understood.

Melanocytes and the melanin they produce play critical biological roles, including protection against UV radiation, thermal regulation, and camouflage [5]. Dysfunction of melanocytes can lead to various diseases such as albinism, melanoma, microphthalmia, and Waardenburg syndrome [6]. Melanocytes produce two types of melanin: eumelanin and pheomelanin. Melanin synthesis involves a series of complex enzymatic and biochemical reactions. The key regulatory factors in melanin production are tyrosinase (TYR), tyrosinase-related protein 1 (TYRP1), and tyrosinase-related protein 2 (TYRP2), with TYR playing a particularly crucial role. Melanin synthesis begins with the oxidation of L-tyrosine and/or L-dihydroxyphenylalanine (L-DOPA) to dopaquinone (DQ), which serves as the substrate to produce both eumelanin and pheomelanin. Dopaquinone undergoes polymerization to form colorless dopachrome. Dopachrome is highly unstable and is rapidly oxidized to dopachrome by another dopaquinone molecule. Under the action of TYRP2, dopachrome is hydroxylated to 5,6-dihydroxyindole-2-carboxylic acid or spontaneously decarboxylated to 5,6-dihydroxyindole. Subsequently, 5,6-dihydroxyindole is catalyzed to 5,6-indolequinone-2-carboxylic acid by TYRP1. 5,6-dihydroxyindole is catalyzed by TYR to form 5,6-indolequinone, which, along with 5,6-indolequinone-2-carboxylic acid and 5,6-indolequinone, ultimately contributes to eumelanin formation. In the presence of cysteine or other thiol-containing compounds, dopaquinone reacts with cysteine to form cysteine-dopachrome, which is further oxidized to pheomelanin. Among these enzymatic reactions, TYR acts as the rate-limiting enzyme in melanin synthesis, as all other reactions can occur under physiological pH conditions [7]. The process of melanin synthesis in melanocytes is regulated by autocrine and paracrine factors. During this process, the microphthalmia-associated transcription factor (MITF) serves as a central regulator by upregulating the two key enzymes, TYR and TYRP1, to promote melanin production. Studies have shown that paracrine factors secreted by keratinocytes, such as MC1R, EDN, and SCF, are major regulators of melanin synthesis. In the skin, in addition to keratinocytes, fibroblasts and immune cells in the dermis also regulate melanin synthesis in melanocytes through paracrine factors [8]. Research on melanocytes has been continuous over the past few decades. To date, studies have primarily focused on the formation, localization, and melanin synthesis of melanocytes, as well as the key cytokines and signaling pathways involved at different stages. These investigations have progressively elucidated the molecular mechanisms underlying melanocyte development and their biological functions, providing a theoretical basis for related disease research and valuable insights for studies on economically significant traits associated with pigmentation in animals.

To date, more than 688 genes associated with pigmentation have been identified [9] (last update: 8 December 2022), These genes encompass melanocyte development-related genes such as MITF, SOX10, PAX3, WNT, EDN3, EDNRB, KIT, MC1R, and melanogenic genes including TYR, TYRP1, DCT, and PEML [6, 10]. In several studies, gene duplication of Endothelin 3 (EDN3) has been identified as the cause of hyperpigmentation in silky fowl [1, 3] and other black-boned chicken breeds, like the Ayam Cemani in Indonesia [11]. However, a study by Darwish et al. [12] challenged this association, revealing that END3 duplication is not linked to hyperpigmentation in LPB. Notably, pigmentation in LPB and silky fowl followed a similar pattern of progressively darkening with age. While no pigmentation was observed at birth, it became visible at around one year of age and pronounced at approximately two years of age. Interestingly, melanin content in the liver of LPB was higher than in other tissues, including the trachea, lung, heart, kidney, and skin, but it was less pronounced and sometimes undetectable in the liver tissue of silky fowl [5, 13]. This highlights distinct mechanisms underlying hyperpigmentation in LPB compared to silky fowl. Although several pigmentation-related genes, such as TYR [14], TYRP1 [15], TYRP2 [16], and MC1R [17], have been explored for their potential involvement in LPB hyperpigmentation, but no conclusive links have been established.

In this study, we aimed to provide a comprehensive understanding of the genetic phylogeny and genetic basis associated with hyperpigmentation in LPB. We achieved this by conducting whole-genome sequencing of 100 LPB and 50 LPN, followed by the integration of this data with a dataset encompassing 421 counterparts representing seven wild and 64 domestic sheep breeds from diverse geographic regions. Additionally, RNA-Seq analysis of liver tissues at two developmental stages was employed to investigate genes associated with melanocyte development and melanogenesis. Through a meticulous combination of genomic and transcriptomic analyses, we identified several potential gene implicated in hyperpigmentation.

Materials and methods

Samples collection

LPB and LPN are both local sheep breeds from Lanping County, Nujiang Prefecture, Yunnan Province, China. They share similar body size and appearance characteristics, with no apparent differences. The distinguishing feature between them lies in LPB having black trait such as black bones and black meat, while LPN lacks these traits. Both populations are raised under consistent conditions, and have no genetic relationship.

Among the 5 local sheep farms in Lanping County, a total of 100 LPB (half male and half female) with no genetic relationship, particularly obvious black-bone traits were selected from the 443 adult LPB, while 50 LPN (half male and half female) with no genetic relationship were selected from the 298 Lanping indigenous sheep. Each sheep was verified to meet the breed standards. Blood samples were drawn from the jugular vein and promptly frozen at -20 °C for further analysis.

The LPB is currently the only known mammal with the black-trait, which is absent in other sheep breeds. To improve the efficiency of screening functional genes related to the black trait, we additionally included 421 sheep of diverse backgrounds. These data encompassed 37 wild sheep individuals from seven different breeds, including Mouflon, Urial, Argali, Bighorn, Thin horn, and Barbary. Furthermore, there were 384 domestic sheep samples, comprising 247 individuals from 14 different breeds in Tibet, Northern and Eastern China, and 137 individuals from 50 diverse breeds across South Asia (SA), Europe, Africa, the Middle East (ME), and Oceania. These additional samples were obtained from 227 publicly available datasets from the NCBI database and 194 datasets shared by Dr. Yu Jiang from Northwest A&F University (Supplementary table S1).

For transcriptome analysis, a total of 23 sheep were slaughtered, including 6 LPB (3 male and 3 female) and 5 LPN (2 male and 3 female) at 2 months of age, and 6 LPB (3 male and 3 female) and 6 LPN (4 male and 2 female) at about 2 years of age. Liver tissue from each individual was sampled for RNA-Seq.

Whole-genome sequence and SNPs calling

Total genomic DNA from 150 blood samples were extracted using the Takara Blood Genome DNA Extraction Kit following the manufacturer’s guidelines. DNA quality and quantity was examined using a NanoDrop device and 1% agarose gel electrophoresis. Paired-end sequencing libraries were generated with an insert size of about 350 bp, according to the manufacturer’s instructions, and sequenced on Illumina Hiseq X Ten platform to an average raw read sequence coverage of ~ 10X. After removing the low-quality raw reads (N content ratio > 10%, low quality base ratio > 50%), high quality reads were mapped to the sheep reference genome assembly Oar_v4.0 using Burrows-Wheeler Aligner (BWA v0.7.15). Picard tools (v1.119) was utilized to identify and remove duplicate reads from the BAM file. After mapping, SNPs were called by GATK (v3.7.0) [18] from the BAM file and the output SNPs were further filtered using vcftools (v4.2) [19]. Raw SNPs were filtered using the following criteria, Quality Depth (QD) < 2.0, root mean square of Mapping Quality (MQ) < 40.0, Fisher Stand (FS) > 60.0, HaplotypeScore > 13.0, and MQRankSum ≤ 12.5, and biallelic SNPs with no more than 10% missing rate. Filtered SNPs were further annotated by ANNOVAR [20] based on the gene annotation of the sheep reference genome Oar_v4.0 and then classified as variations in exonic, intronic, intergenic, upstream, downstream, 5’UTR, 3’UTR, and splicing sites.

Population Structural analysis

Principle component analysis (PCA) was performed based on autosomal SNPs for 571 individuals using GCTA (v1.92.4 beta) [21]. PLINK (v1.9) [22] was used to prune SNPs in linkage disequilibrium in a sliding 100-SNPs window at 10-SNPs steps. After pruning, 302,088 representative SNPs of 571 sheep were extracted for NJ-tree construction and structure analysis. NJ-tree was constructed based on a matrix of pairwise genetic distances using PHYLIP (v3.698) [23] and visualized using the online software iTOL (https://itol.embl.de/). Admixture (v1.2.3) [24] was used to assess the admixture proportions with the default setting; the number of possible genetic clusters K ranged from 2 to 8.

Introgression analysis

Based on the Patterson’s D-statistic (ABBA-BABA test) [25], we used the Pythonscript ABBABABAwindows.py [26] to compute the introgression among two populations (LPB, and LPN). Windows with absolute D-statistic values close to zero (e.g., low 1% windows) were regarded as the genomic regions least affected by gene flow, whereas the highest D-statistic values (e.g., top 1% windows) were regarded as those most affected by gene flow.

Selective sweep analysis

We calculated the genome-wide distribution of FST and LSBL values between LPB and none black-bone domestic sheep population with a 20 kb sliding window using vcftools (v4.2) [19], and the Cross Population Extended Haplotype Homozygosity (XPEHH) values was calculated using Selscan (v2.0.0) [27] .The top 1% of values for the FST, LSBL, and XPEHH were considered as the candidate regions under strong selective sweeps and visualized in the R package. All the candidate regions were assigned to corresponding SNPs and genes using in-house scripts. We further performed KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway and GO (Gene Ontology) terms analysis using the online software g: Profiler [28], and p < 0.05 was used as the threshold for significantly enriched pathways and GO terms [29].

RNA-Seq analysis

RNA-Seq was performed on the liver tissues of 6 LPB and 5 LPN at the pigmentation-not-visible-stage, and in 6 LPB and 6 LPN at the pigmentation-evident stage. The total RNA of 23 liver samples was isolated using Trizol RNA Reagent (Takara), according to the manufacturer’s instructions. The integrity of total RNA was assessed by agarose gel (1%) electrophoresis, prepared in the same manner as for DNA detection. Electrophoresis was conducted at 100 V for 15 min. For concentration and purity analysis, 1 µL of total RNA solution was diluted with 99 µL of sterile water, using sterile water as a blank control. Measurements were performed using a UV-visible spectrophotometer. Then, mRNA libraries were constructed using a NEBNext Ultra Directional RNA Library Prep Kit. The mRNA libraries were then sequenced on an Illumina Hiseq X Ten platform with 100 bp paired-end reads. The raw reads were filtered by removing the low-quality reads and then mapped to the sheep reference genome Oar_v4.0 using TopHat (v2.0.9) [30]. The expression abundance of each assembled transcript was measured through Fragments per Kilobase of exon model per Million mapped reads (FPKM) values. The edgeR 4.0 was used to test the significance of inter group differences in RNA-Seq data processing. The genes that exhibited |log2(fold change) | ≥ 1 and adjusted p ≤ 0.05 in the comparisons between LPB and LPN individuals from two pigmentation stages were considered as differentially expressed genes (DEGs). KEGG pathway and GO term enrichment analysis of differentially expressed genes were performed using the online OmicShare tools.

qRT-PCR

Total RNA was isolated from liver, spleen, and kidney tissues of 6 LPB and 6 Diqing sheep individuals, utilizing TRIzol reagent (Thermo Scientific). cDNA synthesis was achieved with the PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara), and was synthesized by Sangon Biotech (Shanghai) Co., Ltd. Uniform cDNA concentration of 1ug/ul during reverse transcription. Quantitative RT-PCR (qRT-PCR) experiments were conducted with TB Green® Premix Ex Taq™ Takara, the cycle program of qRT-PCR consisted of pre-denaturation (95°C, 3 min), denaturation (95°C, 30 s), annealing (59°C, 30 s) and extension (72°C, 30 s), finally melt curve stage (Ramp rate: +0.5°C /5 s). Simultaneous amplification efficiency E in the range of 90%-110%. Raw data were processed with the ∆∆ Ct method (Kenneth & Thomas, 2001), where ∆Ct = Ct (gene of interest) – Ct (housekeeping gene); ∆∆Ct = ∆Ct (treated sample) – ∆Ct (untreated sample). Primer sequence (5’-3’) for the ERBB4 gene were forward: CGGAACAGTGTGATGGCAGA, reverse: TGTGGGCACTTCTTGACACA, and for ROR1 gene (forward: ACAGAGTTGTCAGTCACTGCT, reverse: CTAATCCGCAGTCGAGAGCC), and reference gene is GAPDH (forward: GCCCTCAACGACCACTTTGT, reverse: TCGGGAGATTCTCAGTGTGG).

Results

Genome sequencing and variations

We conducted genome sequencing for a cohort of 100 LPB sheep and 50 LPN sheep (Supplementary table S1), resulting in a substantial dataset comprising 27.69 billion raw paired-end reads. This dataset yielded an average sequencing depth of 9.30X per individual, with an average genome coverage of 96.47% (Supplementary table S2). To ensure a comprehensive genomic comparison, we augmented our dataset by merging it with a larger collection of whole-genome sequences. This expanded dataset included individuals from seven wild sheep breeds (n = 37 sheep) and 64 domestic sheep breeds (n = 384 sheep), originating from diverse regions, including Tibet, Northern and Eastern China, South Asia, the Middle East, Africa, Europe, and Oceania [31, 32] (Supplementary table S1). The augmented dataset exhibited an average coverage depth of 9.70X per individual.

Following stringent variant filtering criteria, we identified a total of 97,944,357 high-quality single nucleotide polymorphisms (SNPs) across the entire cohort of 571 sheep (Supplementary table S3). These SNPs were categorized into various functional classes, with the majority falling into intergenic regions (65.15%) and intronic regions (32.71%) (Supplementary figure S1). Notably, LPB displayed a lower SNP density (11.81/1000 bp) and nucleotide diversity (π) (2.50 × 10− 3) in comparison to non-black-boned sheep, which encompassed LPN and 64 domestic sheep breeds (exhibiting 19.88 SNPs per 1000 bp and 2.94 × 10− 3 nucleotide diversity) (Supplementary table S4).

Genetic relationship of LPB with LPN and worldwide sheep breeds

We delved into the genetic relationships of LPB with other sheep breeds in a global context. To accomplish this, we employed Principal Component Analysis (PCA) (Fig. 1A), Neighbor-Joining (NJ) phylogeny (Fig. 1B), and ADMIXTURE analysis (Fig. 1C), all of which were based on the extensive dataset of whole-genome SNPs.

The PCA plot and NJ-tree unveiled a distinct separation of Lanping local sheep, comprising both LPB and LPN, from domestic and wild sheep breeds. However, they exhibited a closer genetic affinity with Chinese sheep populations, particularly Tibetan sheep, and those from Northern and Eastern China. This observation suggests that Lanping sheep can be regarded as a sub-population within the broader Chinese sheep genetic landscape. Tibetan sheep occupied an intermediate position between Lanping sheep and Northern and Eastern Chinese sheep, and have strong genetic affinity with the latter two. ADMIXTURE analysis further substantiated this relationship, revealing that Tibetan and Northern and Eastern Chinese sheep exhibited higher levels of shared ancestral components than they did with LPB and LPN when analyzed at K = 4. This aligns with a prior study, which reported a closer genetic connection between Tibetan and Northern and Eastern Chinese sheep in contrast to the relationship between Tibetan sheep and Yunnan-Kweichow sheep populations [31].

Fig. 1 Population genetic analysis of LPB and 9 subgroups of 571 domestic and wild sheep individuals from around the world.A: Principal components 1 and 2 for 571 individual sheep. B: Neighbor-joining phylogenetic tree of 571 individual sheep. C: Population structure of individual 571 with K = 2–8. The abbreviations and subgroups of the breeds and individuals are shown in Supplementary table S1

The Principal Component Analysis (PCA) plot did not exhibit a clear demarcation between LPB and LPN populations. However, the Neighbor-Joining (NJ) phylogenetic tree revealed that while LPB and LPN clustered together in one clade, they further separated into two distinct sub-clusters. Notably, within this sub-cluster delineation, ten LPB individuals displaying hyperpigmentation traits were found to be nested within the LPN branch (Fig. 1B). This intriguing observation raises the possibility that these ten individuals may represent hybrid progeny resulting from interbreeding between LPB and LPN.

ADMIXTURE analysis also supports the results of the PCA and physiological tree mentioned above (Fig. 1C), and the admixture levels and genetic components indicated a high ancestral coefficient and a strong genetic resemblance between the two populations. Nevertheless, a more nuanced perspective emerged when assessing genetic divergence through pairwise fixation index (FST) calculations. The FST value of 0.0210 between LPB and LPN surpassed the genetic differentiation observed between three separate breeds of Tibetan lineage, with FST(TIB/OLA) = 0.0150, FST(TIB/PRT) = 0.0184, and FST(PRT/OLA) = 0.0154 (Table 1). Therefore, while LPB and LPN displayed robust genetic similarities, the notable genetic divergence, as indicated by the FST value, suggests that LPB is indeed derived from LPN and has evolved into a distinct breed.

Table 1 Population divergence measured as FST

Breeds	LPB	LPN	TIB	VLT	OLA	PRT	HUS	CNM	ORI	AGL	
Lanping black-boned sheep (LPB)		0.0210	0.0565	0.0797	0.0555	0.0611	0.0718	0.0883	0.1512	0.3006	
Lanping normal sheep (LPN)	0.0191		0.0508	0.0739	0.0497	0.0552	0.0647	0.0788	0.1452	0.3038	
Tibetan sheep (TIB)	0.0389	0.0372		0.0349	0.0150	0.0184	0.0376	0.0647	0.1283	0.3089	
Valley Tibetan sheep (VLT)	0.0584	0.0567	0.0407		0.0404	0.0436	0.0591	0.0840	0.1557	0.3873	
Oula sheep (OLA)	0.0426	0.0409	0.0259	0.0491		0.0154	0.0343	0.0622	0.1312	0.3359	
Prairie Tibetan sheep (PRT)	0.0494	0.0475	0.0311	0.0552	0.0344		0.0404	0.0689	0.1393	0.3676	
Hu sheep (HUS)	0.0512	0.0492	0.0408	0.0610	0.0448	0.0509		0.0549	0.1298	0.3577	
Chinese merino sheep (CNM)	0.0488	0.0455	0.0445	0.0596	0.0475	0.0527	0.0485		0.1194	0.3165	
Mouflon wild sheep (ORI)	0.0622	0.0586	0.0575	0.0712	0.0594	0.0644	0.0623	0.0539		0.3283	
Argali wild sheep (AGL)	0.0845	0.0853	0.0825	0.1079	0.0910	0.0998	0.0988	0.0892	0.0904		
* FST is given above the diagonal and standard deviation (SD) for each combination is given below.

Genomic signatures associated with hyperpigmentation in LPB

To identify genomic signatures associated with hyperpigmentation in LPB, we conducted a comprehensive analysis of the genome. We assessed allele frequency differentiation by calculating FST between the LPB population and none-black-boned domestic sheep. Using a window size of 20 kb and a step size of 10 kb, we detected a total of 576 putatively positively selected genes (PSGs) based on empirical distributions (with a cutoff of the top 1st percentile) (Fig. 2, Supplementary table S5).

Gene Ontology (GO) enrichment analysis of these PSGs revealed an enrichment in processes related to development. Additionally, a significant Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, the MAPK signaling pathway, was identified, which is crucial for melanocyte development and melanogenesis. This pathway includes 16 genes: MAP3K1, ANGPT1, DUSP1, TGFB3, NFKB1, FGF7, ERBB4, KIT, NF1, MAPKAPK5, FGFR4, PTPN5, TP53, MAP3K12, FGFR1, MAP4K4 (Supplementary table S6 and S7).

Furthermore, we employed locus-specific branch length statistics (LSBL) to analyze the population of LPB, none black-boned sheep, and wild sheep using LPB as the target population. By applying a window size of 20 kb and a step size of 10 kb, we identified 504 genes as PSGs based on empirical distributions (with a cutoff of the top 1st percentile) (Fig. 2, Supplementary table S8). KEGG and GO terms analysis of these PSGs revealed 4 pathways and 39 terms that were over-represented (Supplementary table S9 and S10). The MAPK signaling pathway, which includes 14 genes (ANGPT1, IL1R1, TGFB3, RASGRF2, PAK1, MECOM, ERBB4, KIT, NF1, MAPKAPK5, FGFR4, PAK2, MAP3K7, MAP4K4), was one of the pathways identified.

Finally, we employed the XPEHH test to analyze the LPB and non-black-bone domestic sheep populations. Based on the empirical distribution (previously 1%), we identified 1193 genes as PSGs (Fig. 2, Supplementary table S11). According to the statistical results, a total of 151 genes have been strongly subjected to selection pressure across the Fst, LSBL, and XPEHH methods (Supplementary figure S2). KEGG and GO term analysis of these PSGs revealed that 41 pathways and 2894 terms were significantly overrepresented (Supplementary tables S12 and S13). The MAPK signaling pathway included 28 genes (TNIK, NOTCH2, SPAG9, PDE6G, SLC30A10, PTPRC, NRP1, PLCB1, PLCG2, FBXW7, PTPN11, LIF, CSPG4, EGFR, ERBB4, PTK2B, ADRB3, FGFR1, INSR, FLT4, FGFR4, ADRA1B, ADRA2C, KIT, FERMT2, BMP4, TPD52L1, SASH1). A total of 12 genes (ABCA12, CSNK1A1, KIT, ADAM10, RB1, CEP131, CDH11, PEPD, PTPN11, IFT122, TBX10, and STAR) were both detected by FST, LSBL, and XPEHH analysis. These results provide valuable insights into the genetic mechanisms associated with hyperpigmentation in LPB and shed light on the involvement of specific genes and pathways in this unique trait.

Fig. 2 Genomic regions with strong selective signals (FST, LSBL, andXPEHH) in LPB identified by comparing genomes of refined groups of LPB and none black-bone domestic sheep population based on population structure analysis. The population genetic differentiation FST values, and the LSBL values are were plotted in 20 kb genomic bins with a 10 kb step. The significance threshold of the selection signature was arbitrarily assigned to the top 1% percentile outliers for each individual test and is indicated with red horizontal solid line

DEGs in two pigmentation developmental stages between LPB and LPN

Our previous studies demonstrated that hyperpigmentation in the LPB is the result of ectopic distribution of melanocytes in skin and internal tissues and subsequent melanin synthesized in melanocytes (unpublished data). To understand the differential gene expression associated with pigmentation development in LPB and LPN, RNA-Seq was performed on liver tissues at two different developmental stages: when pigmentation is not yet visible (2-month-old, stage 1) and when pigmentation becomes evident (2-years-old, stage 2). (Supplementary table S14). Liver tissue was chosen for this analysis due to its high melanin content in adult LPB. The total RNA extracted from sheep liver tissue showed an OD260/OD280 ratio between 1.8 and 2.1, and a total RNA concentration between 0.8 and 3 µg/µL (Supplementary table S15), indicating good purity and concentration of the extracted RNA. Concurrently, the gel electrophoresis results showed that three distinct bands corresponding to 28 S, 18 S, and 5 S rRNA were visible without DNA contamination or RNA degradation bands (Supplementary figure S3), indicating that the extracted RNA had good integrity. Approximately 1.35 billion raw reads were generated for the analysis (Supplementary table S16).

At stage 1, there were 13,959 and 13,402 expressed genes in the liver tissues of LPB and LPN, respectively. A total of 619 DEGs were identified through comparisons between LPB and LPN (Fig. 3A, Supplementary table S17). Notably, six of these DEGs were enriched in the melanogenesis pathway, with five of them being up-regulated (PRKACB, KIT, EDNRB, CAM, and CAMK) and one being down-regulated (Wnt). (Fig. 3B, Supplementary table S18). Additionally, GO analysis revealed significant enrichment in the molecular function of lipid binding (Supplementary table S19).

At stage 2, there were 13,917 and 13,518 expressed genes in the liver tissues of LPB and LPN, respectively, and 224 DEGs were identified. Among these, 127 DEGs were up-regulated (including 41 exclusively expressed in LPB), and 97 genes were down-regulated in LPB (Fig. 3C, Supplementary table S20). These DEGs were enriched in 18 significantly over-represented KEGG pathways (Fig. 3D, Supplementary table S21), but no pigmentation-related pathways were identified. Over-represented GO clusters were related to the regulation of cellular processes, oxidoreductase activity, and the extracellular region (Supplementary table S22).

Finally, we conducted a further analysis of the interaction between breed and age factors. A total of 35 DEGs were identified, including 7 significantly upregulated genes and 28 significantly downregulated genes (p < 0.05) (Fig. 3E, Supplementary table S23). Due to the small number, these DEGs were not significantly enriched in any KEGG pathways. We then performed GO enrichment analysis, which revealed 153 enriched pathways (p < 0.05). These pathways were primarily related to immune responses, including innate immune response, defense response to other organisms, response to biotic stimulus, immune response, negative regulation of viral process, and biological processes involved in interspecies interactions between organisms, and no pathways related to melanin were enriched (Supplementary table S24).

Fig. 3 The volcano plots and KEGG pathways of DEGs between LPB and LPN at stages 1 and 2. A: Volcano plot of DEGs identified by comparing the expression levels in liver tissues of LPB and LPN at stage 1. B: The significant 21 KEGG pathways of DEGs at stage (1) C: Volcano plot of DEGs identified by comparing the expression levels in liver tissues of LPB and LPN at stage (2) D: The significant 18 KEGG pathways of DEGs at stage 2. E: Volcano plot of DEGs identified by comparing the expression levels in liver tissues of breed vs. age. F: The significant 24 GO pathways of DEGs for breed vs. age

Screening of candidate genes related to melanin pathway

To further assess the expression abundance of PSGs identified by FST and LSBL statistics, the overlapping genes between PSGs and DEGs in two pigmentation developmental stages were investigated. One gene, KIT, was found to be shared between PSGs and DEGs. KIT is a well-known key gene involved in melanocyte development and melanogenesis. Interestingly, KIT exhibited significantly higher FST and LSBL values, compared to adjacent genomic regions. This suggests a strong selective sweep acting on this gene, indicating its importance in hyperpigmentation. However, the study did not find significant differences in allele frequencies between LPB and LPN populations (Supplementary figure S4).

To identify candidate genes related to melanin synthesis in LPB, we screened all positive selected genes among four populations (LPB, LPN, none black bone Ovis aries, and wild sheep). Finally, two candidate genes (ERBB4 and ROR1) that may be related to melanin synthesis were identified. Both ERBB4 and ROR1 exhibited significantly higher FST, LSBL, and XPEHH values (Fig. 4A), indicating strong selection signals. Among all the positive selection genes screened, only ERBB4 and ROR1 showed significant gene frequency differentiation in different populations. Gene frequency analysis revealed that the frequency of mutant alleles was highest in LPB, followed by LPN, and was significantly lower in other populations (none black bone Ovis aries, and wild sheep). Notably, the frequencies of mutant alleles in LPB and LPN were very similar, suggesting a close genetic relationship (Fig. 4B). Even more noteworthy is that these two genes, ERBB4 and ROR1, represent novel findings not previously documented within the list of 688 known melanoma-associated genes [5]. To account for this genetic similarity, the sliding window gene flow analysis was performed, and the results showed that the ERBB4 and ROR1 genes were within the window D value range of TOP5%, indicating a relatively large gene flow in the LPB and LPN populations (Dataset S1). This genetic closeness may explain why ERBB4 and ROR1 genes were not identified in the transcriptome sequencing results of these two populations.

Fig. 4 The relationship betweenERBB4andROR1genes and hyperpigmentation in sheep. A: FST, LSBL, and XPEHH values around the ERBB4 and ROR1 genes. The parts within the two black dashed lines correspond to gene positions. B: Genotype frequency heatmap of ERBB4 and ROR1 genes. C: Gene relative expression levels of ERBB4 and ROR1 genes in the liver, spleen, and kidneys of Lanping black-boned sheep (LPB) and Diqing sheep. Red represents LPB, blue represents Diqing Sheep

To eliminate the impact of genetic relatedness on transcriptional results, we further examined two genetically distant populations, LPB and Diqing sheep, using qRT-PCR to detect ERBB4 and ROR1 gene relative expression in various tissues. The results demonstrated significant differences in the relative expression levels of ERBB4 in the spleen and kidney between LPB and Diqing sheep (p < 0.0001), with ROR1 showing similar results in the liver, spleen, and kidney (Fig. 4C), indicating that these genes may play a role in hyperpigmentation in LPB.

Discussion

This study conducted genomic and transcriptomic analyses to explore the genetic basis of hyperpigmentation in LPB. Through genetic structure analysis, we identified both genetic similarities and differences between LPB and LPN, confirming that LPB is a distinct breed derived from LPN. Additionally, the genomic analysis identified several functional genes associated with pigmentation in LPB.

In this study, we performed RNA-Seq analysis to examine the DEGs in liver tissues of LPB and LPN at two distinct developmental stages to elucidate the genetic basis underlying pigmentation development. At the first stage (2-month-old, when pigmentation is not yet visible), five genes (PRKACB, KIT, EDNRB, CAM, and CAMK) were upregulated, suggesting their potential roles in the early stages of pigmentation. The downregulation of the Wnt gene may indicate its inhibitory role in melanogenesis. The significant enrichment in lipid binding functions further implies that these genes may could influence melanin biosynthesis through the regulation of lipid metabolism, offering new insights into the mechanisms underlying early pigmentation in LPB. In contrast, at the second stage (2-years-old, when pigmentation becomes evident), no gene were significantly enriched in pigmentation-related KEGG pathways. This suggests that by this stage, the genetic regulation of pigmentation might have already occurred, or could be governed by factors specific to other tissues. The overrepresented GO terms related to oxidoreductase activity and the regulation of cellular processes imply that the gene expression changes at this stage might be more related to maintaining cellular homeostasis and responding to environmental stressors, rather than directly contributing to melanin production. In summary, our findings highlight the crucial role of gene expression changes, particularly during early developmental stages, in the pigmentation process of LPB.

By combining the selective signature with DEGs, this study identified the KIT gene as potentially involved in hyperpigmentation in LPB. The biosynthesis of melanin is meticulously governed by a multifaceted interplay of intracellular signaling cascades, among which the most commonly implicated ones are MC1R/α-MSH [33], cAMP/PKA [34], MAPK/ERK [35], PI3K/Akt [36], and Wnt/β-catenin [37] signaling pathway. The KIT gene we detected is the member of MAPK signaling pathway. In terms of pigmentation, KIT gene mainly regulates the migration of melanocyte precursors to the dermis, epidermis, inner ear, and eye choroid along a specific path [38]. Some studies show that fox white coat was associated with deletion of KIT gene Exon 12 [39], white spots in Arabian camels were associated with KIT gene mutations [40], the translocation of KIT gene led to varying degrees of white spots in cattle [41]. Although the role of KIT gene in pigmentation is well-established in other species, no significant frequency differences in alleles were observed between the LPB and LPN populations for this gene. This finding suggests that additional factors may contribute to hyperpigmentation in LPB.

The ERBB4 and ROR1 genes were further identified, based on their strong selection signals and significant differences in allele frequencies among four populations: LPB, LPN, none black bone Ovis aries, and wild sheep. The ERBB4 gene encodes human epidermal growth factor 4 (HER-4) and plays a crucial role in cell signaling through autophosphorylation, which, in turn, regulates cell growth and division. Studies have reported that the expression rate of ERBB4 protein in melanoma is 66% [42], and the somatic mutation rate of ERBB4 in metastatic malignant melanoma is 18.99% [43], suggesting that ERBB4 may play an important role in melanoma. Receptor tyrosine kinase like orphan receptor 1 (ROR1) is a member of the RORs family that plays a role in the development process, including the growth and development of bones and neurons, cell motility, and cell polarity [44, 45]. ROR1 is highly expressed in melanoma cells and regulates the phosphorylation of Dvl-2, an upstream component of the Wnt pathway, thereby activating both classical and non-classical Wnt signaling pathways. ROR1 can positively regulate the expression of Wnt5a, while indirectly activating STAT3 to inhibit the transcription of melanocyte antigens (MART1, GP100, and promoters PAX3 and MITF), leading to decreased pigment proliferation and antigenicity of melanocytes [46]. Furthermore, in melanoma, the ROR1 gene sustainably undergoes phosphorylation, thus bypassing the Wnt signaling pathway initiated by Wnt5a to promote the occurrence and development of melanoma cells [47], and during the progression of melanoma, ROR1 can promote tumor cell invasion and metastasis by increasing the expression of ROR2 [48, 49].

In particular, we found similar allele frequency distributions in both LPB and LPN populations. These two breeds of sheep both live in the same region (Lanping County, Nujiang Prefecture, Yunnan Province, China), and their body shape and appearance are very similar, except for their black traits (such as black eyelids and gums). According to the survey, there is hybridization between the two populations, and our calculated gene flow results for both populations also confirm it, and this provides an explanation for the similarity in gene frequency distribution between those two populations. To mitigate the impact of gene flow, we selected two distant populations (LPB and Diqing sheep), which distributed in Nujiang Prefecture and Diqing Prefecture of Yunnan Province, respectively. Further qRT-PCR detection results showed that in the LPB population, the relative expression levels of ERBB4 and ROR1 genes in the liver, spleen, and kidney were significantly higher than those in the Diqing sheep population (except the liver for ERBB4 gene). The strong selection signals and allele frequency differences among populations suggest that ERBB4 and ROR1 could be key regulators of melanin synthesis and contribute to the black traits of LPB.

Conclusion

In conclusion, Lanping black-boned sheep (LPB) represent a valuable model for investigating the genetic mechanisms underlying hyperpigmentation in mammals. Through high-throughput whole-genome sequencing, we generated a vast dataset of single nucleotide polymorphisms (SNPs) within the LPB genome. The analysis of these autosomal SNPs provides evidence that LPB have evolved as a distinct breed from their none-black-boned counterparts (LPN). Furthermore, we pinpointed two genes, ERBB4 and ROR1, which are potentially associated with melanin synthesis and hyperpigmentation. These findings contribute to a more comprehensive understanding of non-cutaneous melanocyte development and the genetic basis of hyperpigmentation in LPB. The genetic insights obtained from this study not only shed light on the intriguing phenomenon of hyperpigmentation in LPB but also hold broader implications for our understanding of the genetic underpinnings of pigmentation in mammals.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Acknowledgements

We sincerely appreciated that Dr. Yu Jiang from Northwest A&F University and Dr. Dongdong Wu from Kunming Institute of Zoology, Chinese Academy of Sciences freely provided unpublished 421 sheep data and altruistically technical help.

Author contributions

Y.Q.C. the acquisition, analysis, interpretation of data, have drafted the work; H.L.X. the analysis, interpretation of data, have drafted the work; Z.D.G. the analysis; Y.L. the analysis; J.Y.H. have revised the work; J.W. have revised the work; X.M.H. have revised the work; D.M.X. have revised the work; X.L.T. design of the work, analysis, interpretation of data, and have revised the work; W.D.D. the acquisition, interpretation of data, and have revised the work.

Funding

This work was supported by grants from the National Nature Science Foundation of China (32160771, 32302707), Yunnan Provincial Agricultural Union Foundation (202101BD070001-006), Major Science and Technology Projects in Yunnan Province (202202AE090005), Yunnan Province Basic Research Program Project (202301AU070117), Yunnan Provincial Department of Education Science Research Fund Project (2023J0392), “Xingdian Talent” Industry Innovation Talent Program in Yunnan Province (XDYC-CYCX-2022-0029), and Yunnan Provincial Academician and Expert Workstation (202105AF150043).

Data availability

The raw sequences of 150 sheep were deposited NCBI with SRA accession nos. SRR17839521 - SRR18086562. The 23 liver transcriptomic data reported in this article have been deposited in the China National GeneBank DataBase with the project ID PRJCA026513 (BioSample: SAMC3759895-SAMC3759918, https://download.cncb.ac.cn/gsa2/CRA016699/).

Declarations

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Ethical approval

All the experimental procedures were approved by the Animal Care and Use Committee of Yunnan Agricultural University (Approval Code: 2022718, Approval Data: 8, July, 2022). The euthanasia process fully complies with local animal welfare laws, guidelines, and policies, and the process is as follows: (1) Anesthesia and Sedation: To minimize pain and distress, all sheep were anesthetized prior to euthanasia. The anesthetic used was xylazine, a common sedative in veterinary practice, and a dosage of 0.5 mg/kg body weight was administered intramuscularly. The dosage was calculated based on the animal’s weight to ensure effective sedation. (2) Monitoring Anesthetic Depth: After administering xylazine, the sheep were monitored to ensure they were fully sedated. We checked for signs of deep anesthesia, such as the absence of response to painful stimuli (e.g., pinch reflex), lack of corneal reflex, and muscle relaxation. (3) Euthanasia Procedure: Once deep anesthesia was confirmed, the sheep were euthanized by exsanguination.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Yuqing Chong and Heli Xiong contributed equally to this work.
==== Refs
References

1. Dorshorst B Molin AM Rubin CJ Johansson AM Stromstedt L Pham MH Chen CF Hallbook F Ashwell C Andersson L A complex genomic rearrangement involving the endothelin 3 locus causes dermal hyperpigmentation in the chicken PLoS Genet 2011 7 12 e1002412 10.1371/journal.pgen.1002412 22216010
Dorshorst B, Molin AM, Rubin CJ, Johansson AM, Stromstedt L, Pham MH, Chen CF, Hallbook F, Ashwell C, Andersson L. A complex genomic rearrangement involving the endothelin 3 locus causes dermal hyperpigmentation in the chicken. PLoS Genet. 2011;7(12):e1002412.22216010 10.1371/journal.pgen.1002412
2. Li Y Zhu X Yang L Li J Lian Z Li N Deng X Expression and network analysis of genes related to melanocyte development in the Silky Fowl and White Leghorn embryos Mol Biol Rep 2011 38 2 1433 41 10.1007/s11033-010-0248-2 20848220
Li Y, Zhu X, Yang L, Li J, Lian Z, Li N, Deng X. Expression and network analysis of genes related to melanocyte development in the Silky Fowl and White Leghorn embryos. Mol Biol Rep. 2011;38(2):1433–41.20848220 10.1007/s11033-010-0248-2
3. Shinomiya A Kayashima Y Kinoshita K Mizutani M Namikawa T Matsuda Y Akiyama T Gene duplication of endothelin 3 is closely correlated with the hyperpigmentation of the internal organs (fibromelanosis) in silky chickens Genetics 2012 190 2 627 38 10.1534/genetics.111.136705 22135351
Shinomiya A, Kayashima Y, Kinoshita K, Mizutani M, Namikawa T, Matsuda Y, Akiyama T. Gene duplication of endothelin 3 is closely correlated with the hyperpigmentation of the internal organs (fibromelanosis) in silky chickens. Genetics. 2012;190(2):627–38.22135351 10.1534/genetics.111.136705
4. Deng WD Yang SL Huo YQ Gou X Shi XW Mao HM Physiological and genetic characteristics of black-boned sheep (Ovis aries) Anim Genet 2006 37 6 586 8 10.1111/j.1365-2052.2006.01530.x 17121605
Deng WD, Yang SL, Huo YQ, Gou X, Shi XW, Mao HM. Physiological and genetic characteristics of black-boned sheep (Ovis aries). Anim Genet. 2006;37(6):586–8.17121605 10.1111/j.1365-2052.2006.01530.x
5. Faraco CD Vaz SA Pastor MV Erickson CA Hyperpigmentation in the silkie fowl correlates with abnormal migration of fate-restricted melanoblasts and loss of environmental barrier molecules Dev Dyn 2001 220 3 212 25 10.1002/1097-0177(20010301)220:3<212::AID-DVDY1105>3.0.CO;2-9 11241830
Faraco CD, Vaz SA, Pastor MV, Erickson CA. Hyperpigmentation in the silkie fowl correlates with abnormal migration of fate-restricted melanoblasts and loss of environmental barrier molecules. Dev Dyn. 2001;220(3):212–25.11241830 10.1002/1097-0177(20010301)220:3<212::AID-DVDY1105>3.0.CO;2-9
6. Saleem MD Biology of human melanocyte development, Piebaldism, and Waardenburg syndrome Pediatr Dermatol 2019 36 1 72 84 10.1111/pde.13713 30561083
Saleem MD. Biology of human melanocyte development, Piebaldism, and Waardenburg syndrome. Pediatr Dermatol. 2019;36(1):72–84.30561083 10.1111/pde.13713
7. Pillaiyar T Namasivayam V Manickam M Jung SH Inhibitors of Melanogenesis: an updated review J Med Chem 2018 61 17 7395 418 10.1021/acs.jmedchem.7b00967 29763564
Pillaiyar T, Namasivayam V, Manickam M, Jung SH. Inhibitors of Melanogenesis: an updated review. J Med Chem. 2018;61(17):7395–418.29763564 10.1021/acs.jmedchem.7b00967
8. Yuan XH Jin ZH Paracrine regulation of melanogenesis Br J Dermatol 2018 178 3 632 9 10.1111/bjd.15651 28494100
Yuan XH, Jin ZH. Paracrine regulation of melanogenesis. Br J Dermatol. 2018;178(3):632–9.28494100 10.1111/bjd.15651
9. Baxter LL Watkins-Chow DE Pavan WJ Loftus SK A curated gene list for expanding the horizons of pigmentation biology Pigment Cell Melanoma Res 2019 32 3 348 58 10.1111/pcmr.12743 30339321
Baxter LL, Watkins-Chow DE, Pavan WJ, Loftus SK. A curated gene list for expanding the horizons of pigmentation biology. Pigment Cell Melanoma Res. 2019;32(3):348–58.30339321 10.1111/pcmr.12743
10. Pavan WJ Sturm RA The Genetics of Human skin and hair pigmentation Annu Rev Genomics Hum Genet 2019 20 41 72 10.1146/annurev-genom-083118-015230 31100995
Pavan WJ, Sturm RA. The Genetics of Human skin and hair pigmentation. Annu Rev Genomics Hum Genet. 2019;20:41–72.31100995 10.1146/annurev-genom-083118-015230
11. Dharmayanthi AB Terai Y Sulandari S Zein MS Akiyama T Satta Y The origin and evolution of fibromelanosis in domesticated chickens: genomic comparison of Indonesian Cemani and Chinese silkie breeds PLoS ONE 2017 12 4 e0173147 10.1371/journal.pone.0173147 28379963
Dharmayanthi AB, Terai Y, Sulandari S, Zein MS, Akiyama T, Satta Y. The origin and evolution of fibromelanosis in domesticated chickens: genomic comparison of Indonesian Cemani and Chinese silkie breeds. PLoS ONE. 2017;12(4):e0173147.28379963 10.1371/journal.pone.0173147
12. Darwish HYA Zhang Y Cui K Yang Z Han D Dong X Mao H Deng W Deng X Molecular cloning and characterization of the endothelin 3 gene in black bone sheep J Anim Sci Biotechnol 2018 9 57 10.1186/s40104-018-0272-y 29988351
Darwish HYA, Zhang Y, Cui K, Yang Z, Han D, Dong X, Mao H, Deng W, Deng X. Molecular cloning and characterization of the endothelin 3 gene in black bone sheep. J Anim Sci Biotechnol. 2018;9:57.29988351 10.1186/s40104-018-0272-y
13. Han D Wang S Hu Y Zhang Y Dong X Yang Z Wang J Li J Deng X Hyperpigmentation results in aberrant Immune Development in Silky Fowl (Gallus gallus Domesticus Brisson) PLoS ONE 2015 10 6 e0125686 10.1371/journal.pone.0125686 26047316
Han D, Wang S, Hu Y, Zhang Y, Dong X, Yang Z, Wang J, Li J, Deng X. Hyperpigmentation results in aberrant Immune Development in Silky Fowl (Gallus gallus Domesticus Brisson). PLoS ONE. 2015;10(6):e0125686.26047316 10.1371/journal.pone.0125686
14. Deng WD Xi DM Gou X Yang SL Shi XW Mao HM Pigmentation in Black-boned sheep (Ovis aries): association with polymorphism of the tyrosinase gene Mol Biol Rep 2008 35 3 379 85 10.1007/s11033-007-9097-z 17520341
Deng WD, Xi DM, Gou X, Yang SL, Shi XW, Mao HM. Pigmentation in Black-boned sheep (Ovis aries): association with polymorphism of the tyrosinase gene. Mol Biol Rep. 2008;35(3):379–85.17520341 10.1007/s11033-007-9097-z
15. Li G Xiong H Xi D Memon S Deng W An examination of melanogenic traits and TYRP1 polymorphism in Nanping and Romney Marsh sheep breeds Archives Anim Breed 2018 61 1 131 41 10.5194/aab-61-131-2018
Li G, Xiong H, Xi D, Memon S, Deng W. An examination of melanogenic traits and TYRP1 polymorphism in Nanping and Romney Marsh sheep breeds. Archives Anim Breed. 2018;61(1):131–41.10.5194/aab-61-131-2018
16. Deng WD Tan Y Wang X Xi D He Y Yang S Mao H Gao S Molecular cloning, sequence characteristics, and polymorphism analyses of the tyrosinase-related protein 2 / DOPAchrome tautomerase gene of black-boned sheep (Ovis aries) Genome 2009 52 12 1001 11 10.1139/G09-078 19953128
Deng WD, Tan Y, Wang X, Xi D, He Y, Yang S, Mao H, Gao S. Molecular cloning, sequence characteristics, and polymorphism analyses of the tyrosinase-related protein 2 / DOPAchrome tautomerase gene of black-boned sheep (Ovis aries). Genome. 2009;52(12):1001–11.19953128 10.1139/G09-078
17. Deng WD Shu W Yang SL Shi XW Mao HM Pigmentation in Black-boned sheep (Ovis aries): association with polymorphism of the MC1R gene Mol Biol Rep 2009 36 3 431 6 10.1007/s11033-007-9197-9 18075782
Deng WD, Shu W, Yang SL, Shi XW, Mao HM. Pigmentation in Black-boned sheep (Ovis aries): association with polymorphism of the MC1R gene. Mol Biol Rep. 2009;36(3):431–6.18075782 10.1007/s11033-007-9197-9
18. McKenna A Hanna M Banks E Sivachenko A Cibulskis K Kernytsky A Garimella K Altshuler D Gabriel S Daly M The genome analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data Genome Res 2010 20 9 1297 303 10.1101/gr.107524.110 20644199
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, et al. The genome analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303.20644199 10.1101/gr.107524.110
19. Danecek P Auton A Abecasis G Albers CA Banks E DePristo MA Handsaker RE Lunter G Marth GT Sherry ST The variant call format and VCFtools Bioinformatics 2011 27 15 2156 8 10.1093/bioinformatics/btr330 21653522
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8.21653522 10.1093/bioinformatics/btr330
20. Wang K Li M Hakonarson H ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data Nucleic Acids Res 2010 38 16 e164 10.1093/nar/gkq603 20601685
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16):e164.20601685 10.1093/nar/gkq603
21. Yang J Lee SH Goddard ME Visscher PM GCTA: a tool for genome-wide complex trait analysis Am J Hum Genet 2011 88 1 76 82 10.1016/j.ajhg.2010.11.011 21167468
Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88(1):76–82.21167468 10.1016/j.ajhg.2010.11.011
22. Purcell S Neale B Todd-Brown K Thomas L Ferreira MAR Bender D Maller J Sklar P Bakker PIWD Daly MJ PLINK: a tool set for whole-genome association and population-based linkage analyses Am J Hum Genet 2007 81 3 559 75 10.1086/519795 17701901
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, Maller J, Sklar P, Bakker PIWD, Daly MJ. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75.17701901 10.1086/519795
23. Retief JD Phylogenetic analysis using PHYLIP Methods Mol Biol 2000 132 243 58 10547839
Retief JD. Phylogenetic analysis using PHYLIP. Methods Mol Biol. 2000;132:243–58.10547839
24. Alexander DH Novembre J Lange K Fast model-based estimation of ancestry in unrelated individuals Genome Res 2009 9 19
Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;9:19.
25. Patterson N Moorjani P Luo Y Mallick S Rohland N Zhan Y Genschoreck T Webster T Reich D Ancient admixture in human history Genetics 2012 192 3 1065 93 10.1534/genetics.112.145037 22960212
Patterson N, Moorjani P, Luo Y, Mallick S, Rohland N, Zhan Y, Genschoreck T, Webster T, Reich D. Ancient admixture in human history. Genetics. 2012;192(3):1065–93.22960212 10.1534/genetics.112.145037
26. Martin SH Davey JW Jiggins CD Evaluating the use of ABBA-BABA statistics to locate introgressed loci Mol Biol Evol 2015 32 1 244 57 10.1093/molbev/msu269 25246699
Martin SH, Davey JW, Jiggins CD. Evaluating the use of ABBA-BABA statistics to locate introgressed loci. Mol Biol Evol. 2015;32(1):244–57.25246699 10.1093/molbev/msu269
27. Szpiech ZA Hernandez RD Selscan: an efficient multithreaded program to perform EHH-based scans for positive selection Mol Biol Evol 2014 31 10 2824 7 10.1093/molbev/msu211 25015648
Szpiech ZA, Hernandez RD. Selscan: an efficient multithreaded program to perform EHH-based scans for positive selection. Mol Biol Evol. 2014;31(10):2824–7.25015648 10.1093/molbev/msu211
28. Raudvere U Kolberg L Kuzmin I Arak T Adler P Peterson H Vilo J G:profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update) Nucleic Acids Res 2019 47 W1 W191 8 10.1093/nar/gkz369 31066453
Raudvere U, Kolberg L, Kuzmin I, Arak T, Adler P, Peterson H, Vilo J. G:profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019;47(W1):W191–8.31066453 10.1093/nar/gkz369
29. Zhang DY Zhang XX Li FD Yuan LF Li XL Zhang YK Zhao Y Zhao LM Wang JH Xu D Whole-genome resequencing reveals molecular imprints of anthropogenic and natural selection in wild and domesticated sheep Zool Res 2022 43 5 695 705 10.24272/j.issn.2095-8137.2022.124 35843722
Zhang DY, Zhang XX, Li FD, Yuan LF, Li XL, Zhang YK, Zhao Y, Zhao LM, Wang JH, Xu D, et al. Whole-genome resequencing reveals molecular imprints of anthropogenic and natural selection in wild and domesticated sheep. Zool Res. 2022;43(5):695–705.35843722 10.24272/j.issn.2095-8137.2022.124
30. Kim D Pertea G Trapnell C Pimentel H Kelley R Salzberg SL TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions Genome Biol 2013 14 4 R36 10.1186/gb-2013-14-4-r36 23618408
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.23618408 10.1186/gb-2013-14-4-r36
31. Hu XJ Yang J Xie XL Lv FH Cao YH Li WR Liu MJ Wang YT Li JQ Liu YG The Genome Landscape of Tibetan Sheep reveals adaptive introgression from Argali and the history of early human settlements on the Qinghai-Tibetan Plateau Mol Biol Evol 2019 36 2 283 303 10.1093/molbev/msy208 30445533
Hu XJ, Yang J, Xie XL, Lv FH, Cao YH, Li WR, Liu MJ, Wang YT, Li JQ, Liu YG, et al. The Genome Landscape of Tibetan Sheep reveals adaptive introgression from Argali and the history of early human settlements on the Qinghai-Tibetan Plateau. Mol Biol Evol. 2019;36(2):283–303.30445533 10.1093/molbev/msy208
32. Yang J Li WR Lv FH He SG Tian SL Peng WF Sun YW Zhao YX Tu XL Zhang M Whole-genome sequencing of native Sheep provides insights into Rapid adaptations to Extreme environments Mol Biol Evol 2016 33 10 2576 92 10.1093/molbev/msw129 27401233
Yang J, Li WR, Lv FH, He SG, Tian SL, Peng WF, Sun YW, Zhao YX, Tu XL, Zhang M, et al. Whole-genome sequencing of native Sheep provides insights into Rapid adaptations to Extreme environments. Mol Biol Evol. 2016;33(10):2576–92.27401233 10.1093/molbev/msw129
33. Hsiao JJ Fisher DE The roles of microphthalmia-associated transcription factor and pigmentation in melanoma Archives Biochem Biophys 2014 563 28 34 10.1016/j.abb.2014.07.019
Hsiao JJ, Fisher DE. The roles of microphthalmia-associated transcription factor and pigmentation in melanoma. Archives Biochem Biophys. 2014;563:28–34.10.1016/j.abb.2014.07.019
34. Bellei B, Maresca V, Flori E, Pitisci A, Larue L, Picardo M. p38 regulates pigmentation via Proteasomal Degradation of Tyrosinase. J Biol Chem 2010, 285.
35. Yu-Mi, Kim E-C, Lim LH-M, Young-Kwon. Seo: Rice Bran Ash Mineral Extract increases pigmentation through the p-ERK pathway in zebrafish (Danio rerio). Int J Mol Sci 2019.
36. Moon H-R, Jo SY, Kim HT, Lee WJ, Won CH, Lee MW, Choi JH, Chang SE. Loratadine, an H 1 antihistamine, inhibits melanogenesis in human melanocytes. BioMed Research International 2019, 2019.
37. Su T-R Lin J-J Tsai C-C Huang T-K Yang Z-Y Wu M-O Zheng Y-Q Su C-C Wu Y-J Inhibition of melanogenesis by gallic acid: possible involvement of the PI3K/Akt, MEK/ERK and Wnt/β-catenin signaling pathways in B16F10 cells Int J Mol Sci 2013 14 10 20443 58 10.3390/ijms141020443 24129178
Su T-R, Lin J-J, Tsai C-C, Huang T-K, Yang Z-Y, Wu M-O, Zheng Y-Q, Su C-C, Wu Y-J. Inhibition of melanogenesis by gallic acid: possible involvement of the PI3K/Akt, MEK/ERK and Wnt/β-catenin signaling pathways in B16F10 cells. Int J Mol Sci. 2013;14(10):20443–58.24129178 10.3390/ijms141020443
38. Vandamme N Berx G From neural crest cells to melanocytes: cellular plasticity during development and beyond Cell Mol Life Sci 2019 76 10 1919 34 10.1007/s00018-019-03049-w 30830237
Vandamme N, Berx G. From neural crest cells to melanocytes: cellular plasticity during development and beyond. Cell Mol Life Sci. 2019;76(10):1919–34.30830237 10.1007/s00018-019-03049-w
39. Yan SQ Hou JN Bai CY Jiang Y Zhang XJ Ren HL Sun BX Zhao ZH Sun JH A base substitution in the donor site of intron 12 of KIT gene is responsible for the dominant white coat colour of blue fox (Alopex lagopus) Anim Genet 2014 45 2 293 6 10.1111/age.12105 24308634
Yan SQ, Hou JN, Bai CY, Jiang Y, Zhang XJ, Ren HL, Sun BX, Zhao ZH, Sun JH. A base substitution in the donor site of intron 12 of KIT gene is responsible for the dominant white coat colour of blue fox (Alopex lagopus). Anim Genet. 2014;45(2):293–6.24308634 10.1111/age.12105
40. Holl H, Isaza R, Mohamoud Y, Ahmed A, Almathen F, Youcef C, Gaouar S, Antczak DF, Brooks S. A Frameshift Mutation in KIT is Associated with White spotting in the Arabian Camel. Genes (Basel) 2017, 8(3).
41. Hoque MR, Lee SH, Lim DJ, Cho IC, Lee JH. Investigation of KIT gene polymorphisms in Korean cattle. J Anim Sci Technol 2012, 54(6).
42. Feinmesser M Veltman V Morgenstern S Tobar A Gutman H Kaganovsky E Tzabari C Sulkes J Okon E Different patterns of expression of the erbB family of receptor tyrosine kinases in common nevi, dysplastic nevi, and primary malignant melanomas: an immunohistochemical study Am J Dermatopathol 2010 32 7 665 75 10.1097/DAD.0b013e3181d1e6f0 20559114
Feinmesser M, Veltman V, Morgenstern S, Tobar A, Gutman H, Kaganovsky E, Tzabari C, Sulkes J, Okon E. Different patterns of expression of the erbB family of receptor tyrosine kinases in common nevi, dysplastic nevi, and primary malignant melanomas: an immunohistochemical study. Am J Dermatopathol. 2010;32(7):665–75.20559114 10.1097/DAD.0b013e3181d1e6f0
43. Prickett TD Agrawal NS Wei X Yates KE Lin JC Wunderlich JR Cronin JC Cruz P Rosenberg SA Samuels Y Analysis of the tyrosine kinome in melanoma reveals recurrent mutations in ERBB4 Nat Genet 2009 41 10 1127 32 10.1038/ng.438 19718025
Prickett TD, Agrawal NS, Wei X, Yates KE, Lin JC, Wunderlich JR, Cronin JC, Cruz P, Rosenberg SA, Samuels Y. Analysis of the tyrosine kinome in melanoma reveals recurrent mutations in ERBB4. Nat Genet. 2009;41(10):1127–32.19718025 10.1038/ng.438
44. Green JL Kuntz SG Sternberg PW Ror receptor tyrosine kinases: orphans no more Trends Cell Biol 2008 18 11 536 44 10.1016/j.tcb.2008.08.006 18848778
Green JL, Kuntz SG, Sternberg PW. Ror receptor tyrosine kinases: orphans no more. Trends Cell Biol. 2008;18(11):536–44.18848778 10.1016/j.tcb.2008.08.006
45. Valenzuela DM Stitt TN DiStefano PS Rojas E Mattsson K Compton DL Nuñez L Park JS Stark JL Gies DR Receptor tyrosine kinase specific for the skeletal muscle lineage: expression in embryonic muscle, at the neuromuscular junction, and after injury Neuron 1995 15 3 573 84 10.1016/0896-6273(95)90146-9 7546737
Valenzuela DM, Stitt TN, DiStefano PS, Rojas E, Mattsson K, Compton DL, Nuñez L, Park JS, Stark JL, Gies DR, et al. Receptor tyrosine kinase specific for the skeletal muscle lineage: expression in embryonic muscle, at the neuromuscular junction, and after injury. Neuron. 1995;15(3):573–84.7546737 10.1016/0896-6273(95)90146-9
46. Dissanayake SK Olkhanud PB O’Connell MP Carter A French AD Camilli TC Emeche CD Hewitt KJ Rosenthal DT Leotlela PD Wnt5A regulates expression of tumor-associated antigens in melanoma via changes in signal transducers and activators of transcription 3 phosphorylation Cancer Res 2008 68 24 10205 14 10.1158/0008-5472.CAN-08-2149 19074888
Dissanayake SK, Olkhanud PB, O’Connell MP, Carter A, French AD, Camilli TC, Emeche CD, Hewitt KJ, Rosenthal DT, Leotlela PD, et al. Wnt5A regulates expression of tumor-associated antigens in melanoma via changes in signal transducers and activators of transcription 3 phosphorylation. Cancer Res. 2008;68(24):10205–14.19074888 10.1158/0008-5472.CAN-08-2149
47. Hojjat-Farsangi M Ghaemimanesh F Daneshmanesh AH Bayat AA Mahmoudian J Jeddi-Tehrani M Rabbani H Mellstedt H Inhibition of the receptor tyrosine kinase ROR1 by anti-ROR1 monoclonal antibodies and siRNA induced apoptosis of melanoma cells PLoS ONE 2013 8 4 e61167 10.1371/journal.pone.0061167 23593420
Hojjat-Farsangi M, Ghaemimanesh F, Daneshmanesh AH, Bayat AA, Mahmoudian J, Jeddi-Tehrani M, Rabbani H, Mellstedt H. Inhibition of the receptor tyrosine kinase ROR1 by anti-ROR1 monoclonal antibodies and siRNA induced apoptosis of melanoma cells. PLoS ONE. 2013;8(4):e61167.23593420 10.1371/journal.pone.0061167
48. O’Connell MP Marchbank K Webster MR Valiga AA Kaur A Vultur A Li L Herlyn M Villanueva J Liu Q Hypoxia induces phenotypic plasticity and therapy resistance in melanoma via the tyrosine kinase receptors ROR1 and ROR2 Cancer Discov 2013 3 12 1378 93 10.1158/2159-8290.CD-13-0005 24104062
O’Connell MP, Marchbank K, Webster MR, Valiga AA, Kaur A, Vultur A, Li L, Herlyn M, Villanueva J, Liu Q, et al. Hypoxia induces phenotypic plasticity and therapy resistance in melanoma via the tyrosine kinase receptors ROR1 and ROR2. Cancer Discov. 2013;3(12):1378–93.24104062 10.1158/2159-8290.CD-13-0005
49. Matsuda T Suzuki H Oishi I Kani S Kuroda Y Komori T Sasaki A Watanabe K Minami Y The receptor tyrosine kinase Ror2 associates with the melanoma-associated antigen (MAGE) family protein Dlxin-1 and regulates its intracellular distribution J Biol Chem 2003 278 31 29057 64 10.1074/jbc.M302199200 12754255
Matsuda T, Suzuki H, Oishi I, Kani S, Kuroda Y, Komori T, Sasaki A, Watanabe K, Minami Y. The receptor tyrosine kinase Ror2 associates with the melanoma-associated antigen (MAGE) family protein Dlxin-1 and regulates its intracellular distribution. J Biol Chem. 2003;278(31):29057–64.12754255 10.1074/jbc.M302199200
