Effects of SLC45A2 and GPNMB on Melanin Deposition Based on Transcriptome Sequencing in Chicken Feather Follicles

Simple Summary Understanding the molecular mechanisms involved in feather color formation is important because coloration is one of the most recognizable characteristics in chickens. In this study, after transcrip-tome sequencing of the wing and neck feather follicle tissues of chickens with different plumage colors, we retrieved differentially expressed genes (DEGs) with the same trends in both the wing and neck and then identified DEGs that may be associated with melanin deposition through GO and KEGG annotation and PPI analysis. Finally, we verified that two genes in chicken melanocytes, SLC45A2 and GPNMB, promote melanocyte melanin deposition in chickens. Abstract As an essential genetic and economic trait, chicken feather color has long been an important research topic. To further understand the mechanism of melanin deposition associated with coloration in chicken feathers, we selected feather follicle tissues from the neck and wings of chickens with differently colored feathers (yellow, sub-Columbian, and silver) for transcriptome analysis. We focused on genes that were expressed in both the wings and neck and were expressed with the same trends in breeds with two different plumage colors, specifically, SLC45A2, GPNMB, MLPH, TYR, KIT, WNT11, and FZD1. GO and KEGG enrichment analyses showed the DEGs were enriched in melanin-related pathways, such as tyrosine metabolic pathway and melanogenesis, and PPI analysis highlighted the genes SLC45A2 and GPNMB as associated with melanin deposition. Verification experiments in chicken melanocytes demonstrated that these two genes promote melanocyte melanin deposition. These data enrich our knowledge of the mechanisms that regulate chicken feather color.


Introduction
There are a variety of color combinations and patterns in chicken plumage, making it one of the most colorful terrestrial vertebrates; therefore, these color variations have attracted interest from researchers. In addition to attracting individuals of the opposite sex, feather color helps keep chickens away from predators [1]. In addition, the color of the plumage affects the first impression consumers have when buying chickens, and so plumage color is considered an important economic trait in the industry.
The variation in plumage colors is attributed to the distribution and deposition of melanin, including eumelanin (gray or black) and pheomelanin (yellow or red), produced within melanosomes in melanocytes, and the presence of structural colors resulting from the reflection, refraction, and scattering of light within the feathers [2,3]. Studies have shown that the biosynthesis process of melanin is a series of biochemical reactions initiated by

Ethics Statement
In accordance with the protocol approved by the Institutional Animal Care and Use Committee (IACUC) of Henan Agricultural University, China (Permit Number: 11-0085;

Ethics Statement
In accordance with the protocol approved by the Institutional Animal Care and Use Committee (IACUC) of Henan Agricultural University, China (Permit Number: 11-0085; date: 13 June 2011), all sample collection and treatment procedures were carried out with strict adherence to ethical guidelines to ensure animal welfare and minimize suffering.

Laboratory Animal Samples and Collection
Hens with different plumage colors were selected as experimental subjects (H line chicken: sub-Columbian feathers; Gushi chicken: yellow feathers; Hy-Line chicken: silver feathers). Chickens were housed and raised in our laboratory using standard methods. The tissues were categorized into six groups: HN (sub-Columbian plumage nuchal follicle tissue), HW (sub-Columbian plumage wing follicle tissue), YN (yellow plumage nuchal follicle tissue), YW (yellow plumage wing follicle tissue), SN (white plumage nuchal follicle tissue), and SW (white plumage wing follicle tissue). Each group consisted of three biological replicates, and each replicate contained three follicle tissue samples.
Nine 5-week-old hens were selected for each group, and nuchal follicle tissue and wing follicle tissue were collected immediately. First, the surface of the feather was cleaned with 75% alcohol cotton, the feather root was clamped with tweezers, and the feather follicles were quickly removed; then, another small forceps was used to separate the feather follicle tissue at the root of the feather. The feather was cleaned with PBS, quickly placed in a cryopreservation tube, and finally frozen in a liquid nitrogen tank.

RNA Extraction, Library Construction, and High-Throughput Sequencing
RNA was extracted from the follicles was performed using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA). The RNA quantity and quality were assessed through 1% agarose gel electrophoresis. Eighteen RNA sequence libraries were constructed-HW1, HW2, HW3, HN1, HN2, HN3, YW1, YW2, YW3, YN1, YN2, YN3, SW1, SW2, SW3, SN1, SN2, and SN3. The sequencing libraries for the NEBNext ® UltraTM Directional RNA Library Prep KIT (New England Biolabs, Ipswich, MA, USA) were generated on the Illumina platform. Initially, the raw data in fastq format were processed using an in-house Perl script. During this step, the raw data were filtered to obtain clean reads by excluding those containing adapters, poly-N sequences, and low-quality reads. The clean reads were mapped to the Gallus_gallus-6.0 genome (http://ftp.ensembl.org/pub/release-96 /fasta/gallus_gallus/dna/) URL (accessed on 11 March 2019) using HISAT2 [24], and the expression levels of each gene were calculated using fragments per kilobase of transcript per million mapped fragments (FPKM) [25]. The formula is shown as follows: FPKM= cDNAFragment MappedFragment(Millions) * TranscriptLength(kb) [25]. Differential expression analysis of the 6 groups was performed using DESeq2 [4]. To control the false discovery rate (FDR), only genes with p adjusted < 0.05 were considered differentially expressed.
To screen for alternative splicing events of DEGs associated with melanin deposition, equal amounts of RNA were collected from the nuchal and wing of sub-Colombian and yellow feathers into four mixing pools: the HN group, HW group, YN group, and YW group. The four libraries were then analyzed for full-length transcriptomes using the PacBio Iso-seq platform.

GO and KEGG Enrichment Analysis
Gene Ontology (GO) enrichment analysis of the DEGs was conducted with the GOseq R package (https://www.rdocumentation.org/packages/goseq/versions/1.24.0/topics/ goseq, accessed on 23 June 2023). This analysis employed Wallenius noncentral hypergeometric distribution and adjusted for length bias in the DEGs [26]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource that provides insights into the advanced functions and utilities of biological systems according to the pathways to which

Protein-Protein Interaction Analysis
To examine the interactive connections between DEGs, we obtained protein-protein interaction (PPI) data for the DEGs from the STRING database (http://string-db.org, Organism: Chicken). Subsequently, the DEG PPI network was constructed and visualized using Cytoscape version 3.6.1 (http://www.cytoscape.org/).

Plasmid Construction
To construct an overexpression plasmid for the SLC45A2 and GPNMB genes, the coding DNA sequences (CDSs) of the SLC45A2 and GPNMB genes were amplified from H-line chicken feather follicle cDNA by PCR using gene-specific homologous recombinant clone primers (Supplementary Materials Table S1).

Cell Culture, Identification, and Transfection
Melanocytes were isolated from the peritoneum of 20-day-old Xichuan black-bone chickens. The peritoneum tissue was washed multiple times with phosphate-buffered saline (PBS) containing penicillin and streptomycin. Next, the tissue was cut into pieces and digested with dispase II (Roche, China) and trypsin-EDTA (Gibco, Waltham, MA, USA) at 37 • C for 1 h with gentle shaking every 15 min. To stop the digestion, Medium 254 (Gibco,Waltham, MA, USA) supplemented with 1% Human Melanocyte Growth Supplement (HMGS; Gibco; Waltham, MA, USA) was added at a 1:3 ratio. The resulting cell suspension was then filtered through screens of varying mesh sizes (100, 200, and 500). Next, the filtered suspension was centrifuged at 1000 rpm for 5 min. The obtained cell pellet was resuspended to obtain melanocytes. These melanocytes were cultured at 37 • C and 5% CO 2 [27].
Immunofluorescence analysis was used to confirm that isolated primary cells were melanocytes. Primary cells were fixed with 4% paraformaldehyde for 10 min and then permeabilized with 0.1% Triton X-100 for 1 h, followed by 3 washes with PBS. Then, 3% bovine serum albumin (BSA) was added for 1 h to block nonspecific binding, and melanocyte inducing transcription factor-melanocyte inducing transcription factor (MITF; 1:200, ab122982) (Abcam, Cambridge, UK) was added and incubated at 4 • C for 16 h, followed by 3 washes with PBS. The secondary antibodies used for staining were goat anti-mouse IgG conjugated to Alexa ® Fluor 488 (ab150078; Abcam) and goat anti-rabbit IgG conjugated to Alexa Fluor 555 (ab150117; Abcam). The samples were then incubated with the secondary antibody for 1 h in the dark and washed 3 times with PBS. The cells were incubated with DAPI (Abcam) for 10 min at room temperature in the dark and quickly washed 3 times with PBS to stain the nuclei. Finally, the cells were observed under an inverted fluorescence microscope.
Transient transfections were carried out using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions, and the cells were collected after 48 h of incubation.

Tyrosinase and Melanin Content Measurements
After transfection for 48 h, the melanocytes were washed twice with ice-cold PBS, digested with trypsin and centrifuged for 10 min to obtain the pellet. Then, intracellular tyrosinase was measured by enzyme-linked immunosorbent assay (ELISA) using a Tyrosinase ELISA Assay KIT (Nanjing Jiancheng Bioengineering Institute) according to the manufacturer's instructions. Melanin content was quantitated using the Amplite™ Fluorimetric Melanin Assay Kit (ATT Bioquest, CA, USA).

Statistical Analyses
Statistical analysis was conducted using SPSS Statistics 27 to assess differences among groups. One-way ANOVA was performed, and significance levels were indicated by * (p < 0.05) and ** (p < 0.01). The data are presented as the mean ± SEM (standard error of the mean).

Transcriptome Sequencing and Quality Control
Eighteen libraries were sequenced on an Illumina HiSeq platform. As shown in Supplementary Materials Table S2, clean reads, mapped reads, uniquely mapped reads, and Q20, Q30, and GC contents for each library were identified. A total of 476.78 Mb of clean data was obtained, with an average of 26.49 Mb of clean data for each sample. Of these, the uniquely mapped reads to the chicken reference genome (Gallus_gallus-6.0) ranged from 77.06% to 84.86%. The sequencing data obtained from the experiment exhibited high quality, as indicated by Q20 (>95%) and Q30 (>90%) values and the similar GC content. To ensure accessibility and transparency, the raw sequence data have been deposited in the NCBI Sequence Read Archive under accession number PRJNA859098.

Differential mRNA Expression in Feather Follicles
DESeq2 was used to analyze the DEGs between the groups. Fold change ≥ 1.5 and FDR < 0.05 were used as screening criteria for DEG detection. A total of 12104 genes were detected in this study; 276 DEGs were obtained from the YW group and the HW group, of which 149 genes were upregulated and 127 genes were downregulated; and 212 DEGs were obtained from the YN group and the HN group, of which 98 genes were upregulated and 114 genes were downregulated (Figure 2A,B). Compared to the SW group, the HW group had 4198 DEGs identified, with 2591 genes upregulated and 1607 genes downregulated. A total of 7418 DEGs were obtained in the SN group and the HN group, including 4334 upregulated and 3084 downregulated genes ( Figure 2C,D). To identify the DEGs responsible for the formation of plumage color, we targeted DEGs that had the same trend in the YW vs. HW groups and in the YN vs. HN groups or that had the same trend in the SW vs. HW groups and in the SN vs. HN groups (p < 0.05) ( Figure 2E). If a DEG had the same trend in the YN vs. HN group and the YW vs. HW group, we denoted it as the Y vs. H group and similarly for the S vs. H group.
To some extent, alternative splicing allows the same gene to produce multiple isoforms. Four AS forms were detected, namely, the alternative 3 splice site (A3SS), alternative 5 splice site (A5SS), retained intron (RI), and skipped exon (SE) [28]. We identified AS forms of genes expressed in feather follicle tissues of the HN vs. YN group and HW vs. YW group, which showed similar proportions for the two groups. Unfortunately, we did not find different alternative splicing events in the DEGs of the HN vs. YN groups or in the HW vs. YW groups.
isoforms. Four AS forms were detected, namely, the alternative 3′ splice site (A3SS), alte native 5′ splice site (A5SS), retained intron (RI), and skipped exon (SE) [28]. We identifi AS forms of genes expressed in feather follicle tissues of the HN vs. YN group and H vs. YW group, which showed similar proportions for the two groups. Unfortunately, w did not find different alternative splicing events in the DEGs of the HN vs. YN groups in the HW vs. YW groups.

Data Validation
We randomly selected nine DEGs for validation using qRT-PCR and found that the expression levels of all genes were consistent with the results of RNA-seq, indicating that the detection and expression of our RNA-seq are accurate ( Figure 3). Venn diagram showing the overlap of the DEGs in four comparisons.

Data Validation
We randomly selected nine DEGs for validation using qRT-PCR and found that the expression levels of all genes were consistent with the results of RNA-seq, indicating that the detection and expression of our RNA-seq are accurate ( Figure 3). The left Y-axis displays the FPKM derived from the RNA-seq, while the data from RT-qPCR are shown on the Y-axis on the right. The data are represented as the means ± SEDs; red represents the FPKM value of RNA-seq, blue represents the RT-qPCR with GAPDH as the internal reference gene; * indicates p < 0.05; ** indicates p < 0.01.

GO Analysis and KEGG Pathway Enrichment Analysis
Gene enrichment analysis was conducted to investigate the distribution of GO items and reveal differential enrichment in the biological process, cellular component, and molecular function categories. The GO classification of the DEGs can be observed in Figure 4 and Supplementary Materials Table S3. Within the biological process category, the predominant subcategories identified in all four groups included cellular process, single-organism process, biological regulation, and metabolic process. The left Y-axis displays the FPKM derived from the RNA-seq, while the data from RT-qPCR are shown on the Y-axis on the right. The data are represented as the means ± SEDs; red represents the FPKM value of RNA-seq, blue represents the RT-qPCR with GAPDH as the internal reference gene; * indicates p < 0.05; ** indicates p < 0.01.

GO Analysis and KEGG Pathway Enrichment Analysis
Gene enrichment analysis was conducted to investigate the distribution of GO items and reveal differential enrichment in the biological process, cellular component, and molecular function categories. The GO classification of the DEGs can be observed in Figure 4 and Supplementary Materials Table S3. Within the biological process category, the predominant subcategories identified in all four groups included cellular process, single-organism process, biological regulation, and metabolic process.
The cellular process was the most affected of the subcategories, including upregulated genes AQP1 in the Y vs. H groups, ras homolog family member C (RHOC), frizzled class receptor 1 (FZD1), frizzled class receptor 6 (FOXO6), dihydroorotate dehydrogenase    Table S4), listing the top 20 pathways with the lowest Q values. Among the pathways that are significantly associated with melanin are tyrosine metabolism, the MAPK signaling pathway, adrenergic signaling in cardiomyocytes, and the calcium signaling pathway. Among these pathways, the tyrosine metabolism pathway was the most significantly enriched pathway in the Y vs. H groups (p < 0.05). By analyzing the top 20 pathways with the lowest Q values in the four comparison groups, we found the distinctly expressed gene TYR in the tyrosine metabolic pathway. the lowest Q values. Among the pathways that are significantly associated with melanin are tyrosine metabolism, the MAPK signaling pathway, adrenergic signaling in cardiomyocytes, and the calcium signaling pathway. Among these pathways, the tyrosine metabolism pathway was the most significantly enriched pathway in the Y vs. H groups (p < 0.05). By analyzing the top 20 pathways with the lowest Q values in the four comparison groups, we found the distinctly expressed gene TYR in the tyrosine metabolic pathway. HN groups. The x-axis represents the rich factor, which is the ratio of the DEGs annotated with the pathway term to the total number of genes annotated with the pathway term. The greater the rich factor is, the greater the degree of enrichment. The y-axis shows each KEGG pathway name. Each round point represents a specific KEGG pathway. The circle size indicates the number of DEGs associated with each significantly enriched pathway. The circle color indicates the significance level (q-value). A q-value < 0.05 was considered to indicate significant enrichment. Light purple indicates the least significant, and orange represents the most significant. HN groups. The x-axis represents the rich factor, which is the ratio of the DEGs annotated with the pathway term to the total number of genes annotated with the pathway term. The greater the rich factor is, the greater the degree of enrichment. The y-axis shows each KEGG pathway name. Each round point represents a specific KEGG pathway. The circle size indicates the number of DEGs associated with each significantly enriched pathway. The circle color indicates the significance level (q-value). A q-value < 0.05 was considered to indicate significant enrichment. Light purple indicates the least significant, and orange represents the most significant.
In addition, the KEGG classification table shows that there were other melanin-related pathways, such as melanogenesis, Wnt signaling, Notch signaling, calcium signaling, adherens junction, and mTOR signaling in the four comparison groups (Supplementary Materials  Table S5). Interestingly, we found that the melanogenesis pathway was common in all four groups, so we focused on that pathway ( Figure 6). We found that nine DEGs were enriched in the melanogenesis pathway-TYR, KIT, WNT11, FZD1, frizzled class receptor 2 (FZD2), frizzled class receptor 7 (FZD7), frizzled class receptor 9 (FZD9), adenylate cyclase 3 (ADCY3), and adenylate cyclase 3 (ADCY8), which had the same trend in the YW vs. HW groups and the YN vs. HN groups, or had the same trend in the SW vs. HW groups and SN vs. HN group (Supplementary Materials Table S6). common in all four groups, so we focused on that pathway ( Figure 6). We found that nine DEGs were enriched in the melanogenesis pathway-TYR, KIT, WNT11, FZD1, frizzled class receptor 2 (FZD2), frizzled class receptor 7 (FZD7), frizzled class receptor 9 (FZD9), adenylate cyclase 3 (ADCY3), and adenylate cyclase 3 (ADCY8), which had the same trend in the YW vs. HW groups and the YN vs. HN groups, or had the same trend in the SW vs. HW groups and SN vs. HN group (Supplementary Materials Table S6).

PPI analysis of DEGs
In the Y vs. H groups, the DEG PPI network was composed of nine proteins and nine pairs of PPIs. Moreover, TYR interacted with SLC45A2, KIT, GPNMB, and MFSD12. Melanophilin (MLPH) interacted with SLC45A2, TYR, and MFSD12 ( Figure 7A). In addition, TYR was mainly enriched in tyrosine metabolism and melanogenesis. Furthermore, the PPI network in the S vs. H groups consisted of 161 proteins and 107 pairs of PPIs (Supplementary Materials Table S7). The key core node FZD1 exceeded seven interactions, including WNT11, transcription factor 7-like 2 (TCF7L2), and RHOC. Among them, both heat shock protein family A (Hsp70) member 8 (HSPA8) and heat shock protein family B (small) member 1 (HSPB1) were enriched in the MAPK signaling pathway, while WNT11 and TCF7L2 were enriched in the melanogenesis pathway and the Wnt signaling pathway ( Figure 7B).

PPI Analysis of DEGs
In the Y vs. H groups, the DEG PPI network was composed of nine proteins and nine pairs of PPIs. Moreover, TYR interacted with SLC45A2, KIT, GPNMB, and MFSD12. Melanophilin (MLPH) interacted with SLC45A2, TYR, and MFSD12 ( Figure 7A). In addition, TYR was mainly enriched in tyrosine metabolism and melanogenesis. Furthermore, the PPI network in the S vs. H groups consisted of 161 proteins and 107 pairs of PPIs (Supplementary Materials Table S7). The key core node FZD1 exceeded seven interactions, including WNT11, transcription factor 7-like 2 (TCF7L2), and RHOC. Among them, both heat shock protein family A (Hsp70) member 8 (HSPA8) and heat shock protein family B (small) member 1 (HSPB1) were enriched in the MAPK signaling pathway, while WNT11 and TCF7L2 were enriched in the melanogenesis pathway and the Wnt signaling pathway ( Figure 7B).

Effects of SLC45A2 and GPNMB Overexpression on Melanin Deposition
Notably, we found that the expression of SLC45A2 and GPNMB in the H line was significantly lower than that in the Gushi line using transcriptome sequencing and qRT-PCR. Moreover, in the transcriptome sequencing, the eumelanin-related genes MLPH and

Effects of SLC45A2 and GPNMB Overexpression on Melanin Deposition
Notably, we found that the expression of SLC45A2 and GPNMB in the H line was significantly lower than that in the Gushi line using transcriptome sequencing and qRT-PCR. Moreover, in the transcriptome sequencing, the eumelanin-related genes MLPH and TYR were significantly reduced in Group H compared to Group Y (p < 0.05). Given that SLC45A2 has the same expression trend as MLPH and TYR, high expression of SLC45A2 and GPNMB is also speculated to affect the deposition of melanin.
To explore the biological role of SLC45A2 and GPNMB in eumelanin in chickens, we overexpressed the two genes separately in chicken primary melanocytes. Primary melanocytes were identified by immunofluorescence staining for MITF ( Figure 8A). The mRNA expression levels of the SLC45A2 and GPNMB genes in melanocytes were increased by approximately 18000-fold and 3-fold after transfection with pcDNA3.1-SLC45A2-EGFP and pcDNA3.1-GPNMB-EGFP, respectively, compared to the control group that was transfected with the pcDNA3.1-EGFP plasmid (Figure 8 B,C).  After overexpression of SLC45A2 or GPNMB, the mRNA expression levels of the eumelanin marker genes TYR, MLPH, and MITF were significantly increased (p < 0.01; Figure 8D,E). Consistently, western blotting also confirmed that the protein levels of MITF were upregulated following melanocytes' transfection with pcDNA3.1-SLC45A2-EGFP and pcDNA3.1-GPNMB-EGFP, respectively, compared to the control group ( Figure 8F,G). Full WB images can be found in Supplementary Materials Figure S1. Furthermore, the tyrosinase and melanin contents in cells transfected with the recombinant plasmids pcDNA3.1-SLC45A2-EGFP or pcDNA3.1-GPNMB-EGFP were significantly higher than those in the control group (p < 0.01) ( Figure 8H-K).

Discussion
In chicken feathers, melanin distribution and deposition are responsible for plumage color variations, and feather follicles are the only place where plumage melanin is produced [29]. During feather follicle melanogenesis, mature melanocytes produce melanosomes; newly synthesized melanin then migrates to the feather keratinocytes, giving the feathers their particular color [30]. Pigmentation is a complex feature, and many studies have shown that melanin pigmentation is strongly genetically controlled [31]. In recent years, several key genes involved in melanin deposition have been identified. These include MLPH, TYR, KIT, PMEL, MITF, agouti signaling (ASIP), melanocortin 1 receptor (MC1R), and endothelin 3 (EDN3) [32]. However, there has been little research on the association of pigmentation on chicken feathers. In particular, because most of the feather color is caused by eumelanin and pheomelanin, and most of the marker genes related to melanin deposition, for example TYR, affect both eumelanin and pheomelanin, it is particularly important to further study the colorful plume color due to total melanin deposition [33,34]. In our transcriptome data, there were significant differences in the expression of MLPH, TYR, and KIT (p < 0.05). MLPH is an important structural protein that assists in the transport of melanosomes, ensuring that melanosome bodies accumulate at the dendritic borders of melanocytes and release them into surrounding tissues for deposition, regulating the color of the animal's skin and coat. MLPH gene mutations result in chickens with the Japanese quail feather lavender diluted phenotype [35,36]. The expression of the MLPH gene in the goose has tissue specificity, with the highest expression in black eyes and the lowest expression in back skin, abdominal skin, and fins, but it was not detected in the heart, liver, and other organs [37]. This expression pattern is consistent with the melanin deposition pattern of geese, showing that the MLPH gene is related to melanin deposition [37]. TYR is a key gene in the biosynthesis of melanin in animals, and the tyrosinase it encodes has the catalytic activity of tyrosine hydroxylase and dopa oxidase and is a rate-limiting enzyme in the process of melanin biosynthesis. The lack of tyrosinase catalytic function hinders melanin synthesis, resulting in corresponding albinism traits in the animal's coat and skin [38]. At present, many research results show that functional mutations in the TYR gene and differences in gene expression levels in animals such as cattle, sheep, dogs, cats, martens, rabbits, rats, and poultry will lead to albino traits in their coats [14,39,40]. The insertion of an intact avian retroviral sequence in intron 4 in the chicken TYR gene causes deletion of transcript exon 5, eventually leading to chicken recessive white feather mutation [8]. The mast stem cell factor receptor encoded by the KIT gene is a member of the tyrosine kinase receptor family, which acts on blood cells and melanin, affecting the growth and development of blood cells and melanocytes [41]. When KIT is expressed in a melanocyte precursor, its ligand SCF activates receptors on the cell membrane through the MAPK signaling pathway, ultimately initiating melanin synthesis [42]. Together with the previous results, we speculate that the expression of MLPH, TYR, and KIT plays an important role in the deposition of melanin in chicken feathers.
To understand the potential functions of the identified mRNAs, the DEGs were subjected to functional enrichment analyses. KEGG pathway analysis showed that the differentially expressed genes TYR, KIT, WNT11, FZD1, FZD2, FZD7, FZD9, ADCY3, and ADCY8 were significantly enriched in tyrosine metabolic pathway or melanogenesis. This result suggested that the DEGs identified in the feather follicles play an important role in melanin biosynthesis and pigmentation. In previous studies, we collected dorsal feather follicles (sub-Columbian feathers) and ventral feather follicles (white feathers) of the nuchal area of the H line for RNA-seq and found that the expression of WNT11 and SLC45A2 significantly differed in these two sites, while WNT11 was enriched in the melanogenesis pathway [22]. WNT11 also plays a crucial role in the Wnt signaling pathway. Several studies have demonstrated that the Wnt pathway is involved in melanin synthesis [43,44]. In addition, we found that the expression of WNT11 in the feather follicles of the H line was significantly increased after feeding 1.0% tyrosine daily [22]. However, there are very few reports on the synthesis of WNT11 in melanin. In this study, the expression level of WNT11 in sub-Columbian feathers was significantly higher than that in white feathers (p < 0.05). Using KEGG enrichment analysis, we found that WNT11 was enriched in the melanogenesis pathway and Wnt signaling pathway. In the PPI network, WNT11 interacted with FZD1, which is also involved in the melanogenesis pathway. We speculate that FZD1 and WNT11 are responsible for the formation of dark spots (melanin) at specific sites in sub-Columbian feathers, but the exact genetic mechanism needs to be further investigated.
In our transcriptome data, we focused on DEGs (p < 0.05) between different plumage colors and selected DEGs with the same trend in both wing and neck. Moreover, through PPI analysis, genes that are consistent with the expression trend of MLPH, TYR, and KIT, such as SLC45A2 and GPNMB, are likely to exhibit similar functions to melanin marker genes.
SLC45A2 is often found to be overexpressed in black tissues compared with white tissues [45]. Moreover, after transfection of SLC45A2 into mouse melanocytes, MITF, TYR, TYRP1, DCT, and tyrosine contents were significantly increased (p < 0.05) [46]. This is similar to the result of significantly elevated MITF and TYR after we overexpressed SLC45A2 on chicken primary melanocytes (p < 0.05). According to the results of previous experiments in our laboratory, the expression of SLC45A2 on the dorsal side of the nuchal (sub-Columbian plumage) area of the H line was significantly higher than that on the ventral side of the nuchal (white plumage) area, and the addition of different concentrations of tyrosine to melanocytes significantly affected the expression of SLC45A2 [22]. In our experimental results, after overexpressing SLC45A2, the tyrosinase content and melanin content in cells were significantly increased, so we hypothesized that the high expression of SLC45A2 could promote the deposition of melanin. Notably, SLC45A2 is the known pathogenic gene of an autosomal recessive hypopigmentary disorder, oculocutaneous albinism type 4 (OCA4) [45]. Moreover, the missense mutation of SLC45A2 results in a dominant silver locus that dilutes pheomelanin [47,48]. Similarly, no effect on eumelanin was observed in SLC45A2 mutant horses, with some effect on pheomelanin dilution [49]. Cysteine is a key amino acid in the synthesis of pheomelanin, so we speculate that this mutation prevents cysteine transport to melanosomes. In addition, SLC45A2 transports sugar in yeast cells in an acid-dependent manner [50]. The anti-melanogenic effects of sugar and sugar derivatives are probably due to three different mechanisms: by increasing melanosomal pH, by interfering with melanosome maturation, or by inhibiting TYR maturation by blocking N-glycosylation [51]. It has been previously demonstrated that the proton/glucose exporter SLC45A2 mediates melanin synthesis and melanoma metastasis primarily by modulating melanosomal glucose metabolism. In human melanocytes and melanoma cell lines, the SLC45A2 gene is highly enriched, and MATP, the gene's protein, is found in melanosomes [52]. By knocking down MATP using siRNAs, melanin content and tyrosinase activity were reduced. It helps maintain an appropriate melanosome pH, so that copper ions can be incorporated into the tyrosinase active site correctly [52]. Therefore, we speculate that in the melanocytes of chicken feather follicles, SLC45A2 may regulate the pH of melanosomes by regulating the glucose metabolism of melanosomes, leading to the deposition of melanin.
At the cyclin-dependent kinase inhibitor 2A (CDKN2A) tumor suppressor locus, there are two non-coding and two missense mutations associated with sex-linked barring [53]. Schwochow Thalmann found that one or both of the non-coding changes cause an upregulation of CDKN2A expression in feather follicles during feather growth [10]. Melanocyte progenitor cells may stop dividing prematurely when CDKN2A expression is high, distinguishing into pigment-producing cells instead [10]. In the KASP typing results published earlier in our laboratory, the silver feather only had the mutations of SLC45A2 and no mutation of cyclin-dependent kinase inhibitor 2A (CDKN2A), while sub-Columbian plumage have both SLC45A2 mutations and CDKN2A mutations. We speculate that CDKN2A is also a key gene influencing melanin deposition in chicken feathers. But CDKN2A did not appear in the differentially expressed genes in our transcriptome data. However, we verified that the gene was significantly different in the sub-Columbian plumage vs. yellow plumage by q-PCR (p < 0.05) (Supplementary Materials Figure S2). The results were consistent with the previous studies.
Melanin is synthesized in the melanosomes of melanocytes, and its formation goes through four stages [54]. GPNMB is a type I transmembrane glycoprotein, which is present in all stages of melanosome formation (I-IV), and especially enriched in mature (stage III and IV) melanosomes [55,56]. The expression level of GPNMB in melanocytes was found to be inversely correlated with the metastatic capacity of human melanomas and was shown to be associated with the development of the retinal pigment epithelium and the iris [57,58]. In mice with pigmentary glaucoma, a premature stop codon mutation in the GPNMB gene results in iris pigment dispersal [59]. Combined with the results of this experiment, in chicken melanocytes, overexpression of GPNMB can increase intracellular tyrosinase content and melanin content, we speculated that the high expression of GPNMB could promote the deposition of melanin. During development, GPNMB is expressed in a pattern similar to that of MITF, DCT, and PMEL [60]. GPNMB and PMEL are highly homologous proteins that play an important role in melanosome biogenesis [61,62]. In addition, GPNMB adheres to PAM212 keratinized cells in an RGD-dependent manner and is possibly involved in melanocyte development, melanin synthesis, and melanoma production [63]. MITF has been regarded as an essential transcription factor for the regulation of more than 25 pigmentation genes, such as TYR, TYRP1, DCT, PMEL17, and so on [64]. In addition, MITF is involved in the Wnt/β-catenin signaling pathway, MC1R/α-MSH signaling pathway, and SCF/C-KIT signaling pathway [65][66][67]. The first two pathways can increase MITF expression, and these pathways affect melanin production by regulating the transcription of TYR, TYRP1, and DCT through MITF [68]. Based on this experiment, after overexpression of GPNMB in chicken melanocytes, the expression of TYR and MITF increased. From this, we speculate that GPNMB affects the deposition of melanin by upregulating the transcription of TYR, TYRP1, and DCT by upregulating MITF.
It is worth noting that by screening for differentially expressed genes with the same expression trend as marker genes such as TYR in the two groups of this experiment, combined with PPI correlation analysis, we finally speculated that SLC45A2 and GPNMB genes may affect the deposition of total melanin on chicken feather follicles. However, since the experimental cells were obtained from the peritoneal tissue of black chickens, the cells can only secrete eumelanin, so we can only prove that SLC45A2 and GPNMB have an effect on the deposition of eumelanin, and whether they have an effect on pheomelanin needs further experimental verification.

Conclusions
In this study, transcriptome sequencing and subsequent GO, KEGG, and PPI analyses were used to screen for genes with the same trend as melanin deposition marker genes and DEGs (SLC45A2, GPNMB, MLPH, TYR, KIT, WNT11, and FZD1) in the neck and wing feather follicles of chickens. Moreover, we found that SLC45A2 and GPNMB can promote the deposition of melanin in chicken feathers. Overall, the candidate genes detected in these experiments can help provide novel insights into the mechanisms regulating feather color development in chickens.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ani13162608/s1, Table S1. List of overexpression vector construction and qPCR primers; Table S2. Statistics of quality control and reads mapping of sequencing date; Table S3. GO classification of the differentially expressed genes; Table S4. The top 20 pathways in KEGG enrichment in the four comparisons; Table S5. KEGG classification of the differentially expressed genes; Table S6. DEGs in the Melanogenesis pathway; Table S7. PPI network of DEGs; Figure S1: Full Western Blots representing images used in Figure 8F,G; Figure S2: Differences in the relative expression of CDKN2A in the two comparison groups.