Characterization of Chicken Skin Yellowness and Exploration of Genes Involved in Skin Yellowness Deposition in Chicken

Skin color is an important economic trait in meat-type chickens. A uniform bright skin color can increase the sales value of chicken. Chickens with bright yellow skin are more popular in China, especially in the broiler market of South China. However, the skin color of chickens can vary because of differences in breeds, diet, health, and individual genetics. To obtain greater insight into the genetic factors associated with the process of skin pigmentation in chickens, we used a colorimeter and high-resolution skin photographs to measure and analyze the skin color of chickens. By analyzing 534 chickens of the same breed, age, and feed condition, we found that the yellowness values of the chickens varied within this population. A significant positive correlation was found between the cloacal skin yellowness values before and after slaughter, and the cloacal skin yellowness value of live chickens was positively correlated with the overall body skin yellowness value. Additionally, chicken skin yellowness exhibited low heritability, ranging from 0.07 to 0.27. Through RNA sequencing, 882 genes were found to be differentially expressed between the skin with the highest and lowest yellowness values. Some of these differentially expressed genes may play an important role in yellow pigment deposition in chicken skin, which included TLR2B, IYD, SMOC1, ALDH1A3, CYP11A1, FHL2, TECRL, ACACB, TYR, PMEL, and GPR143. In addition, we found that the expression and variations of the BCO2 gene, which is referred to as the yellow skin gene, cannot be used to estimate the skin yellowness value of chickens in this population. These data will help to further our understanding of chicken skin yellowness and might contribute to the selection of specific chicken strains with consistent skin coloration.


INTRODUCTION
Ma-Huang chicken is a high-quality broiler chicken with golden skin that is preferred by consumers and is one of the most popular chickens in the frozen chicken market of South China. Under the condition of an increasing frozen fresh chicken market, the first impression that a broiler presents to consumers is made by its skin. The yellowness and uniformity of skin color are two important factors influencing consumer choice. Many factors can regulate skin color, and skin yellowness value varies greatly between individuals. The market price of chickens with white or light yellow skin is lower than that of chickens with yellow skin. How to select for uniform skin color and generate chickens with higher yellowness is a critical problem limiting the development of commercial Ma-Huang chickens.
Coloration is a phenotypic trait related to a variety of adaptive functions, including temperature regulation, camouflage, mate choice, etc (Hill, 1991;Hamilton et al., 2013). Skin pigmentation is a complex biological process. Skin color can be influenced by environmental factors, nutrition, and heredity. Carotenoids are compounds responsible for skin color in animals. Most animals can only acquire carotenoid pigments from the environment. Carotenoids can be synthesized by plants, bacteria, and fungi and affect the red, orange, and yellow characteristics of the skin (Goodwin, 1954;Walsh et al., 2012;Toews et al., 2017). Chickens obtain carotenoids through eating feeds. Once ingested, these pigments are transported to the mucosal cells of the duodenum. Scavenger receptor class B, type 1 can recognize carotenoids incorporate with lipid micelles and facilitate carotenoids transportation (Eroglu and Harrison, 2013;Toomey et al., 2017). Once absorbed, carotenoids are transported by lipoproteins through the bloodstream to the target tissues. The abundance and nature of lipoproteins can limit the capacity to obtain carotenoids during this process (Galvan et al., 2016). Furthermore, many species can metabolize carotenoids into alternate, more oxidized forms that are deposited into the integument. The result is that several factors, such as amounts and species of carotenoids ingested, nutrient and health status of the host, genetic background, and food matrix composition, influence the bioavailability of carotenoids (West and Castenmiller, 1998). Therefore, there are large differences between species or individuals in their capacity to obtain these pigments from their diets, resulting in skin color variation. However, unlike endogenously synthesized pigments such as melanin, which have well-characterized genotype-phenotype connections, researchers were poorly understood the genes and pathways involved in exogenously derived carotenoid metabolism and deposition (Toews et al., 2017).
β-Carotene oxygenases are the best-characterized set of enzymes that are involved in carotenoid metabolism, including β-carotene 15, 15 oxygenase (BCO1) and β-carotene 9 ,10 oxygenase (BCO2) (Toews et al., 2017). BCO2 is a mitochondrial protein with a set of varied targets, including zeaxanthin, lutein, and lycopene. BCO2 is able to directly influence coloration by asymmetrically decomposing colored carotenoids into colorless defatted carotenoids, which hinder the deposition of colored pigments in fat and skin (Amengual et al., 2011). In domesticated chickens, lower expression of BCO2 is associated with carotenoid degradation into colorless derivatives (Eriksson et al., 2008). Eriksson et al. (2008) analyzed gene sequences associated with a specific phenotype found in domestic chickens across numerous wild junglefowl and domestic breeds and found that three SNPs in BCO2 were significantly related to yellow skin [located at chr24: 6,264,085 (G < A), chr24: 6,273,428 (A < G), and chr24: 6.287.900 (G < A)]. The presence of BCO2 initiates the degradation of carotenoids and protects mitochondria from oxidative damage (Harrison and Kopec, 2020). Although the above studies have indicated the importance of BCO2 in chicken yellow skin regulation, its function in the regulation of skin pigmentation and skin yellowness value still remains unclear.
In the present study, we systematically analyzed the variation in skin yellowness in yellow-skinned broiler chickens under intensive conditions. We found that the yellowness values of skin varied within the chicken population under intensive conditions. To identify the candidate genes that influence yellowness in chicken skin, we used RNA-seq to investigate differences in skin transcriptomes between the skin tissues of the Ma-Huang chicken showing high yellowness and low yellowness. We identified some candidate genes and metabolic pathways involved in yellow pigment deposition in chickens. All of these findings will help to establish a foundation for breeding chickens with high yellow skin coloration.

Ethics Statement
Ma-Huang chickens were obtained from a chicken farm in Guangzhou, Guangdong Province, China. Animal experiments were performed in accordance with the regulations and guidelines established by the Animal Care Committee of South China Agricultural University (approval number: SCAU#0106; 25 November 2018).

Animals and Traits Measurement
Commercial Ma-Huang chicken pure line broiler chickens were used in the present study. Chickens at 120 days of age raised under the same feeding management conditions were chosen. Slaughter tests and traits measurements were carried out in a total of 534 Ma-Huang chickens (24 males and 510 females). These chickens were from one generation with complete pedigree. All chickens were slaughtered by bloodletting in the jugular vein, then scalded in tanks filled with water at a temperature of 64 • C for 220 s and plucked with a chicken defeatherer. The slaughter determination was analyzed using methods adapted from previous studies (Sirri et al., 2010). Growth traits and carcass traits were measured according to the NY/T 823-2004 standard (Chen et al., 2004).
The skin yellowness was measured using a 3nh-NH310 colorimeter (3NH, Guangzhou, China). The colorimeter can directly obtain a * value (which represents redness), b * value (which represents yellowness), and L * value (which represents luminance) from a single measurement. We measure the skin yellowness by b * value. To reduce measurement errors, skin color values were determined by the same person, and quantification was performed three times at each position. The skin b * values (yellowness) of the cloaca (for living chickens) was measured before bloodletting. This area was chosen because it contains a few feathers or blood vessels. After bloodletting, scalding and plucking, the skin b * values (yellowness) of eight different parts were measured. All color measurement areas were free from obvious defects (plucking damage, bruises, full blood vessels, hemorrhage, discoloration, and any other condition that might affect a uniform color reading).

Digital Skin Color Value Quantification
Digital photographs of 349 slaughtered chickens, placed on a white disk as the background, were successfully obtained with a Canon EOS 80D color camera (Canon, Tokyo, Japan). We developed an in-house program in MATLAB 2018a (The MathWorks, Inc., Natick, MA, United States) to automatically retrieve the color values of the chicken skin. RGB values were converted into HSV space according to standard formulas implemented by the rgb2hsv function of MATLAB (Supplementary Figure 1). Median HSV values, H values, S values, and V values were then reported for each image. In the HSV model, H represents hue (H represents the variation of the color type), S is saturation (S represents the purity or intensity of the color), and V is brightness or luminance. The V values were observed under environmental light, which could not be fully controlled, so we focused on the H and S color dimensions. A photograph is provided in Supplementary Figure 2 as an example. The MATLAB binary "get_skin_color" is provided in the Supplementary Material.

DNA and RNA Extraction
Genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen, CA, United States). Total RNA was extracted from tissues using RNAiso reagent (Takara, Otsu, Japan) and quantified using a NanoDrop spectrophotometer and an Agilent 2100 Bioanalyzer (Thermo Fisher Scientific, MA, United States).

Genetic Parameter Estimates
The restricted maximum likelihood method implemented by the DMU package was used to obtain estimates of the phenotypic and genetic (co)variance and heritability (Madsen and Jensen, 2013), and the following linear model was used for the analysis of the data: where y was the phenotypic value of a trait, a was the vector of the animal additive genetic effect, b was the vector of the fixed effects, including gender (two levels) and pen (six levels), e was the vector of random residuals, and X and Z were the incidence matrices. A 534 individuals with complete pedigree and phenotypic trait records were included.

RNA-seq and Bioinformatics Analysis
Among the 534 Ma-Huang chickens, three chickens with high cloacal skin yellowness values were chosen and defined as the high yellowness (s_deep) group; another three chickens with low cloacal skin yellowness values were chosen and defined as the low yellowness (s_light) group. There was a significant difference in the average yellowness of the s_deep and s_light groups (12.53 ± 0.69 vs. 6.99 ± 0.66, p < 0.001, n = 3). A piece of skin (approximately 8 mm in diameter) was isolated from the cloaca, and the subcutaneous fat was scraped off and then snap-frozen in liquid nitrogen before RNA-seq was conducted. Oligo(dT)attached magnetic beads were used to purify mRNA. The purified mRNA was fragmented into small pieces with fragmentation buffer at the appropriate temperature. Then, first-strand cDNA synthesis was performed. Thereafter, A-Tailing Mix and RNA-Index Adapters were added via incubation for end repair. The cDNA fragments obtained from the previous step were amplified by PCR, and the products were homogenized with Ampure XP Beads and then dissolved in EB solution. The products were validated in an Agilent Technologies 2100 system for quality control. The double-stranded PCR products from the previous  step were heated, denatured and circularized via the splint oligo sequence to obtain the final library. The single-stranded circular DNA was formatted as the final library. The final library was amplified with phi29 to generate DNA nanoballs (DNBs), which included more than 300 copies of a single molecule. The DNBs were loaded into the patterned nanoarray, and single-end 50-base reads were generated on the BGI-500 platform (BGI-Shenzhen, China). Adaptors, reads with unknown base(N) > 5%, and lowquality sequences were removed using the FASTX tool, which is used to reduce data noise and to clean raw reads. Furthermore, SOPA2(version 2.2.1) was used to filter out potential rRNA reads.

Statistical Analysis
Descriptive statistical analysis was performed using all available yellowness values records by the SPSS 21.0 (SPSS Inc., Chicago, IL, United States). Pearson correlation coefficients (r) and linear regression were used to evaluate the correlation among traits by the SPSS 21.0. A p-value < 0.05 was considered statistically significant. For statistical analysis of two contrast, we using onesample t-test through SPSS. We considered p < 0.05 to be statistically significant. * p < 0.05; * * p < 0.01.

Characteristic of Skin Yellowness in Ma-Huang Chicken
To understand the characteristics of skin pigmentation in chickens, we collected skin yellowness data from 534 Ma-Huang chickens. A total of eight different body skin positions were evaluated with a colorimeter (Figure 1A), and the discrete distribution of the yellowness value of each position is shown in Table 2. The results showed that the skin yellowness values varied among individuals ( Figure 1B). Thigh skin exhibited the lowest mean yellowness values among the measured positions, whereas abdominal fat exhibited the highest mean yellowness values. Next, we used hue and saturation values that were collected from photographs and analyzed with MATLAB software to assess the overall skin color of individual chickens. To better classify the skin yellowness of the chickens, the depth of the skin color of Ma-Huang chickens was graded into three levels according to the HSV values as follows: "nearly white, " "pale yellow, " and "yellow" (Figure 1C) (according to the HSV ranking, the individual's visual perception is divided into three parts. The specific skin color classification data are in the Supplementary  Table 1). Within a subsample, there were 120 chickens with "nearly white" skin, 120 chickens with "pale yellow" skin, and 120 chickens with "yellow" skin. The overall hue value of the skin ranged from 0.06 to 0.10, and the overall saturation value of the skin ranged from 0.15 (very light) to 0.58 (yellow).
To better understand the association between Hue-Saturation values and skin color, three representative chicken images were mapped to the 3D scatterplot space (Figure 1D). The darkness variation in saturation was easy to perceive, but there was little yellow-red variation in hue ( Figure 1D). These data indicate that under intensive farming conditions, the saturation of the skin color of broilers presented a greater difference than the hue.

The Correlation Analysis of Skin Yellowness Between Different Parts of Chicken Skin
To better understand the correlation of yellowness values among different positions and the correlation between skin yellowness and hue-saturation values, we conducted correlation analysis for each value. In general, there was a significant positive correlation among the yellowness values at each measured position, except for the V value and abdominal fat values (Figure 2A). Notably, the yellowness value of the cloaca before slaughter was positively correlated with the yellowness value of cloacal skin after slaughter as well as the yellowness values at other skin positions ( Figure 2B). These results suggested that the yellowness value of the cloacal skin in a live chicken can serve as a reference for the overall skin yellowness value of the carcass. Then, we analyzed the correlation between skin color values and other economic traits of chickens. The yellowness value of the cloaca before slaughter was positively correlated with many traits, especially growth traits such as live weight and tibia length ( Figure 2C). The HSV values were also positively correlated with these economic traits. However, we noted that many traits related to fat were uncorrelated with skin color traits, and the yellowness value of abdominal fat was significantly negatively correlated with abdominal fat weight. In summary, these results indicate that the skin yellowness values at different positions on the chicken are correlated with each other and that the abdominal fat weight is negatively correlated with skin color.

Estimation of Genetic Parameters of Skin Yellowness Value
Considering that skin yellowness has a crucial effect on chicken sales, we evaluated the genetic parameters associated with the skin yellowness value. The heritability of the skin yellowness value was low, ranging from 0.07 to 0.27 (Table 3). Shank skin yellowness presented the lowest heritability, whereas abdominal fat yellowness presented the highest heritability (Table 3). Among all carcass traits, the carcass weight exhibited the highest heritability, and the heritability of the eviscerated weight (EW), eviscerated weight with giblets (EWG), and abdominal fat weight (AFW) ranged from 0.29 to 0.39 (Table 3).

Differentiated Expressed DEG Between Low Yellowness Skin and High Yellowness Skin
To explore the genetic mechanism underlying the regulation of pigmentation, six RNA-seq libraries were obtained from the skin of chickens with high yellowness (s_deep) and low yellowness (s_light) (Figure 3A). There was a significant difference in the average yellowness of the s_deep and s_light groups (12.53 ± 0.69 vs. 6.99 ± 0.66, p < 0.001). During filtering, linker contamination, unknown base N content greater than 10%, and low-quality reads are sequentially removed (we define reads with a quality value of fewer than 15 bases that account for more than 50% of the total bases in the reads as low-quality reads). A total of 23.03 and 23.01 M valid reads were obtained from s_light and s_deep, respectively (Supplementary Table 2). According to a fold change ≥2 and an adjusted p ≤ 0.001, 882 unigenes were identified as differentially expressed between s_deep and s_light (Wang et al., 2010; Figure 3B). The 20 most significantly differentially expressed genes, including 10 upregulated and 10 downregulated, are listed in Table 4. Compared with the s_deep group, 612 genes were upregulated, and 270 genes were downregulated in the s_light group ( Figure 3B) (fold change ≥2, adjusted p ≤ 0.001). Through integrative analysis of these DEGs between the s_deep and s_light group, we found that 794 DEGs were commonly expressed in both samples, whereas 42 DEGs and 46 DEGs were specifically expressed in the s_deep chickens and the s_light chickens, respectively ( Figure 3C and Supplementary  Table 3). These group-specific DEGs may be important for resulting in the different performance of pigmentation between the two group chickens. According to the results of the analysis of differentially expressed genes, hierarchical clustering analysis (HCA) was performed on the combined DEGs using the heatmap function in R software. HCA demonstrated diverse expression profiles between s_deep skin tissues and s_light skin tissues ( Figure 3D).
To define the biological functions of the 882 DEGs, GO and KEGG enrichment analysis was performed with the phyper function of R software. Top 20 most significant pathways and GO terms were listed (Figures 3E,F). The KEGG pathway enrichment analysis showed that the DEGs were significantly involved in Focal adhesion, ECM-receptor interaction, olfactory transduction, drug metabolism-cytochrome P450, and cytokinecytokine receptor interaction ( Figure 3E). Notably, there are two enriched pathways related to cytochrome P450 in this result ( Figure 3E). Previous studies have shown that the gene encoding the cytochrome P450 enzyme was found to associate with carotenoid-based avian coloration phenotypes (Lopes et al., 2016;Mundy et al., 2016). Therefore, these cytochrome P450 related pathways may play an important role in the regulation of the chicken yellow skin trait. GO analysis showed that intermediate filament, structural constituent The X-axis represents the difference multiplied after log2 conversion, and the Y -axis represents the significant values after log 10 conversion. Red represents upregulated DEGs, blue represents downregulated DEGs, and gray represents non-DEGs. (C) venn diagram representing DEGs overlapping between the two groups. (D) A heat map clusters the FPKM values of the differential genes in each comparison group. The horizontal axis represents the log2 (FPKM + 1) of the sample, and the vertical axis represents genes. The redder the color of the color block, the higher the expression, and the bluer the color, the lower the expression. (E) KEGG enrichment analysis of the DEGs between s_deep and s_light. The y axis represents KEGG pathways and the x axis represents | log(p value)|. The number of genes enriched in the indicated pathway was marked. (F) GO enrichment analysis of DEGs between s_deep and s_light. The y axis represents GO terms and the x axis represents | log(p value)|. The number of genes enriched in the indicated GO terms was marked.
of cytoskeleton, intermediate filament cytoskeleton, muscle tissue development, and multicellular organismal process were enriched in the functions of DEGs ( Figure 3F). Furthermore, tyrosinase (TYR), pre-melanosomal protein (PMEL), and G protein-coupled receptor 143 (GPR143), which are involved in the regulation of pigmentation, were significantly differentially expressed between the two groups (Supplementary Table 3).

Gene Interaction Networks Involved in the Regulation of Skin Yellowness
Studies have shown that immunity and lipid metabolism are related to skin pigmentation. To validate our RNA-seq results and find some candidate genes involved in skin color regulation, we selected 11 DEGs implicated in immunity, lipid metabolism, and pigmentation for qPCR confirmation. We observed consistency between the qPCR results and RNAseq results (Figures 4A,B). Next, by using GeneMANIA, we generated a protein interaction network to show how these eleven DEGs communicated or interacted with each other (Figure 4C). This network contributes to our understanding of the regulatory mechanism linking pigment deposition and immunity in chickens. Additionally, we also used GeneMANIA software to predict the potential target genes of TYR, PMEL, GPR143, and BCO2, which were involved in pigmentation. Then we generated an interaction network consisting of these four genes and their potential targets (Figure 4D). This network may play roles on chicken skin carotenoid deposition. Within this network, TYR, PMEL, and GPR143 were differentially expressed between s_light and s_deep according to the RNA-sequencing results (Supplementary Table 3). PMEL and GPR143 are mainly involved in the regulation of melanin deposition, and their roles in carotenoid deposition have rarely been reported. Therefore, based on our interaction network, it will be an important research content to study how PMEL and GPR143 interact with TYR and BCO2, and ultimately affect the carotenoid deposition in chicken skin. The differential expression of these genes may lead to an imbalance in the network and ultimately affect yellow pigment deposition.

THE EXPRESSION AND THE THREE TYPICAL SNPs OF BCO2 CANNOT USE TO ESTIMATE THE CHICKEN SKIN YELLOWNESS IN THIS CHICKEN POPULATION
BCO2 plays an important role in determining the yellow skin of chickens. Mutations in the region of the BCO2 gene can result in a phenotypic change. Eriksson et al. (2008) reported that three single-nucleotide polymorphisms (SNPs: A, B, and C) around the BCO2 gene were associated with the yellow skin phenotype and that low BCO2 expression allows the deposition of yellow carotenoids in the skin. However, our RNA-seq and qPCR results showed that the BCO2 yellow skin gene was not significantly differentially expressed between the s_deep samples and the s_light samples (Figures 5A,B). To further test whether the three BCO2 SNPs are correlated with the skin yellowness value in this experimental group, we selected 40 two-tailed samples of skin yellowness values for verification (20 chickens with highest yellowness values and 20 chickens with lowest yellowness values). The genotyping results showed that SNP-A and SNP-B were the same in all 40 chickens, and only SNP-C showed polymorphism in these 40 chickens (Supplementary Table 4). Next, we compared the skin yellowness value between the chickens with GA genotype and the chickens with AA genotype. The result showed that there is no significant difference in skin yellowness value between chickens with different type of SNP-C ( Figure 5C). Correlation analysis results also showed that the SNP-C was not correlated with skin yellowness value in any area (Supplementary Table 5). Therefore, these results suggest that the expression and the three typical yellow skin SNPs of BCO2 may not be used to estimate chicken skin yellowness.

DISCUSSION
Skin color makes phenotypic selection difficult in chickens because the skin of chickens is covered with feathers. Several positions covered with few feathers can be used to measure the skin color directly, such as the skin around the cloaca and the skin of the underwing. However, there is a vein under the skin of the underwing, which can significantly affect the color value detection. Therefore, the skin around the cloaca may be the only position at which skin color can be directly measured in a live chicken. Importantly, there are no blood vessels under the skin around the cloaca, and our correlation analysis results showed that the skin yellowness value at this position was correlated with overall chicken skin color after slaughter. These results suggest   showed that the expression of BCO2 is no significant difference between s_deep and s_light samples. p value = 0.156815, one-sample t test was used to assess the difference between the s_deep and s_light samples. (B) qPCR results of the BCO2 gene (p = 0.11) n = 3. One-sample t test was used to assess the difference between the s_deep and s_light samples. (C) Skin yellowness value (B* value) in chickens with GA genotype and AA genotype. The skin yellowness has no significant difference between chickens with GA genotype and AA genotype.
that the cloacal skin yellowness of a live chicken may be a useful value for assessing the overall skin yellowness of a slaughtered chicken, and may imply for application in deep yellow skin chicken selection. At present, research on the skin color of poultry mainly depends on the use of a colorimeter and color plate to assess skin color after slaughter. For example, Sirri et al. (2010) used a CR-300 chromameter to determine the skin color of the breast, thighs, and shanks. Mueller et al. (2020) used a CR-300 chromameter to detect the skin color of the left breast to evaluate the skin color differences among different broiler breeds. The skin color detection position employed in the experiment of Rajput et al. was the same as employed by Sirri et al. (Rajput et al., 2013). However, a colorimeter can only determine the color within a small piece of skin, and the color plate measurement is influenced by subjective factors. To overcome these defects, we developed a method combining highresolution photos and MATLAB software valuation to quantify the global skin color of chickens. A similar method has been employed in human skin color research (Jacobs et al., 2013). On the basis of this method, we can use a single number to represent the skin color of a chicken. In addition, we found that many skin yellowness values determined by the colorimeter were positively correlated with the global skin color value analyzed by MATLAB, suggesting that the skin yellowness value of a single position is representative to some extent. Additionally, our results showed that the skin yellowness values detected at different positions were positively correlated with each other. Díaz-Gómez et al. (2017) also found that there was a correlation between breast skin yellowness and thigh skin yellowness in chickens. Therefore, Pigments may evenly deposit in different areas of the chicken body skin.
Broilers are prone to excessive fat deposition, but most of this fat is not physiologically necessary. Yellow fat is occasionally observed in cattle, sheep, pigs, and rabbits, but it occurs widely in poultry (Strychalski et al., 2016). Broilers deposit a portion of the adsorbed dietary pigments into the fat. Female broilers exhibit higher subcutaneous fat levels than males (Fletcher, 1992). Our experiment showed that the yellowness value of abdominal fat presents a very significant negative correlation with abdominal fat weight (Figure 2), which means that more abdominal fat deposition will lead to lighter skin color. The potential underlying mechanism may be the dilution of the pigment deposited in the fat by lipids, resulting in a decrease in the yellowness value.
The BCO2 gene encodes beta-carotene dioxygenase 2, which cleaves carotenoids (colorful substance) to apocarotenoids (colorless substance) via an asymmetrical cleavage reaction (Kiefer et al., 2001), is detected in many tissues and can affect color trait (Lindqvist et al., 2005). Kong et al. (2017) reported that the frequency of the reported SNPs (SNP A, SNP B, SNP C) of BCO2 was significantly different between yellow-and whiteskinned. Using CRISPR/Cas9 to knock out the EcBCO2 gene can lead to a color change in the hepatopancreas of prawns (Sun et al., 2020). However, there are also studies showing that BCO2 is not related to yellow skin traits. For example, Yang et al. (2012) found that one SNP of BCO2 (chr24: 6273428) was the same among three yellow skin chicken breeds. Xu et al. (2015) reported that the three classic SNPs of BCO2 and the yellow skin traits of Gushi chickens were not statistically correlated. In our experiment, both mRNA expression and association analysis showed that BCO2 and skin yellowness values were not correlated. Kiefer et al. (2001) proposed that when the expression of BCO2 is low in skin tissue, carotenoids remain intact and will be incorporated into keratinocytes. We speculate that BCO2 may determine whether the skin is yellow or white but cannot determine the yellowness value of skin. The skin yellowness of chickens may be a quantitative trait regulated by multiple genes.
Here, we found that the yellowness value of cloaca skin in live chicken is significantly correlated with skin yellowness of chicken after slaughter. However, this result still need further verified.
We are planning to use the yellowness value of cloaca skin in live chicken to select chickens and divide them into high and low yellowness groups. The chickens will then be slaughtered and detected for the yellowness value of their skin. If the high yellowness group has significantly higher yellowness value than those in low yellowness group after slaughtered, this method would be implied for application in yellow skin chicken selection. The lighter the skin color, the cheaper the chicken is sold in South China. Chicken breeding companies can weed out chickens with low yellowness value, therefore, resulting in more uniform skin color in the flock.
In addition, considering the low heritability of skin yellowness that our results showed, it would take much time and require much more records to improve this trait in chicken breeding. Phenotypes such as skin yellowness value, carotenoid levels in the blood, health condition, and composition of pigment additives in feed can be recorded and use to guide the breeding selection. Besides, molecular markers can also be used in chicken yellow skin breeding. Even though the three SNPs of BCO2 were not correlated with chicken skin yellowness in this study, other SNPs may be useful. DEGs found in our results, especially those genes included in the network, can be potential candidates to scan SNP and do the association analysis with chicken skin yellowness.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/genbank/, GSE144982.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Care Committee of South China Agricultural University (approval number: SCAU#0106; 25 November 2018).