The animal experiments were performed with the approval of the Animal Care Committee of the Institute of Animal Science, Guangdong Academy of Agricultural Sciences (Guangzhou, China) and Use Committee (GAAS-IAS-2009-73).

The experimental population was a three-generation cross-breeding population construction from two different meat-type chicken lines. The first founder line was Huiyang Bearded Chicken (HB), which is an indigenous breed in China. The second founder line was the commercial broiler breed “High Quality Chicken A Line” (HQLA). The ratio of HVP of Huiyang Bearded Chicken is much higher than that of HQLA. A total of 395 individuals from the F₂ population were used for GWAS, including 212 cocks and 183 hens from eight half-sib families. At 91 days of age, venous blood was collected and slaughtered. As shown in Figure 1, HVP is evaluated by peeling: chickens with visceral peritoneum pigmentation are HVP chickens, and chickens with no peritoneal pigmentation are normal chickens. Among the 395 birds, 124 individuals had HVP and 271 individuals had absent HVP. Heritability estimation was carried out by GAPIT (Wang and Zhang, 2021) software. We used a linear mixed model to estimate genetic variance based on genomic information.

Genomic DNA was extracted from venous blood using the phenol/chloroform method as described before (Jia et al., 2016). All the birds were genotyped using the 10X SLAF-seq (Specific-Locus Amplified Fragment Sequencing), which was performed by Biomarker Technologies (Beijing, China) Co., Ltd. For enzyme cutting site prediction, the genome of Gallus gallus 5.0 (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/315/GCF_000002315.4_Gallus_gallus-5.0) was used as the reference genome. The qualifying DNA of the sampled individuals was digested using HaeIII restriction enzymes. The length of paired-end reads was 125 bp, while the sequencing insert size was 364–444 bp. The consistent part of the variation obtained by GATK (DePristo et al., 2011) and samtools (Li et al., 2009) was used as the final variation site for subsequent analysis. After quality control, 115,706 SNP markers remained with call rate >70% and MAF> 0.01. BEAGLE version 5.2 was used to impute missing SNPs (Browning and Browning, 2013). The 30X SLAF-seq data were obtained from 22 F₀, with individuals as a reference panel. r ² of LD was calculated using PLINK 2.0 (Chang et al., 2015); then, regrouping was done, and the mean of LD and the LD decay plot were obtained using the ggplot2 package. The SNP distribution is shown in Figure 2A. The LD decay line quickly decreased first; after 50 kb, the distance tended to be slow (Figure 2B). These values thus indicate that the F₂ population has relatively short LD distances and quick LD decay.

GWAS was performed using the logistic mixed model (LMM); the statistical model is l o g i t ( π i ) = X i α + G i β + b i where π i = P ( y i = 1 ) is the probability of a HVP for animal I; X_i is the vector of covariates for animal i; α is the vector of fixed covariate effects including the intercept, gender, batch, and two principal components (PC1 and PC2); Gi is the genotype of animal i; and β is the genotype effect. It was assumed that b   ∼   N ( 0 ,   σ g 2 K ) , where σ g 2 is the additive genetic variance and K is the additive genetic relationship matrix based on genomic information. GWAS was implemented in the GMMAT (Chen et al., 2016) package, which is based on the R language.

The Bonferroni correction (Benjamini and Hochberg, 1995) was used to control the false positive rate in the GWAS analyses, namely, p = 0.05/N, where N is the number of SNP loci. In this study, we used a genome-wide significance threshold p-value of 4.32*10⁻⁷ (i.e., p = 0.05/115,706). Quantile–quantile (QQ) plots and Manhattan plots were plotted using the R package CMplot (Yin et al., 2021). Use of a web-based platform ensemble (URL:http://archive.ensembl.org/Gallus_gallus/Info/Index) translated the significant SNPs from Gallus gallus 5.0 to Gallus gallus 6.0.

The population structure was analyzed using Admixture software (Alexander et al., 2009), and admixture proportions were visualized with the R package pophelper (Francis, 2017). The genomic inflation factor (λ), defined as the ratio of the genome-wide chi-squared statistics mean (or median) to the expected distribution of this statistic under the null hypothesis, was used to assess the impact of population stratification on association analyses. Generally, the more the value of genomic inflation factor deviates from 1, the higher the rates of false positive errors in GWAS analyses (Devlin et al., 2001). The formula for the genomic inflation factor is as follows: λ = M e d i a n ( T i 2 ) 0.4549 where M e d i a n ( T i 2 ) is the median of the statistic from the real data and 0.4549 is the expected median under the null hypothesis.

Linkage disequilibrium (LD) of SNP loci significantly associated with HVP was performed with the software haploview (Barrett et al., 2005). After the haplotypes significantly associated with HVP are identified, candidate genes in the haplotype blocks are detected. The statistical analysis of haplotypes with HVP was performed with the R package haplo.stats (Schaid et al., 2002).

The nearest genes present within 50 kb upstream or 50 kb downstream of the SNP loci significantly associated with HVP were identified as potential candidate genes based on the Gallus gallus 6.0 reference genome (http://archive.ensembl.org/Gallus_gallus/Info/Index).

The visceral peritoneum tissues were collected from adult chicken with extremely different HVP performance with no peritoneal pigmentation (n = 5) (absent HVP control group) and five HVP-positive (severe HVP) chickens. The total RNA was extracted with RNA TRIzol (Lin et al., 2019) reagents according to standard protocols. An M-MLV RTase cDNA Synthesis Kit was used to prepare cDNA. The visceral peritoneum MFNG, MYH9, RANGAP1, and CYP2D6 gene expression was determined by quantitative real-time polymerase chain reaction (qRT–PCR). qRT–PCR primers were designed using NCBI Primer-Blast. The reactions of qRT–PCR were performed in 20 µL volumes. Parameters of RT-PCR were as follows: 3 min 95°C, followed by 40 cycles of 30 s at 95°C, 30 s at 60°C, and 45 s at 72°C. The relative gene expression levels of the target genes MFNG, MYH9, RANGAP1, and CYP2D6 were calculated using the 2 ^−△△CT method and normalized with housekeeping genes: GAPDH. The differential expression of MFNG, MYH9, RANGAP1, and CYP2D6 in the visceral peritoneum tissue, between the normal and HVP-positive, was determined using a t-test with the software Graphpad Prism 8 (Knapek et al., 2020).

The visualization of population stratification was performed by principal component analysis (PCA) for the first three principal components (PCs) (Figure 3). The first two PCs, gender, and batch were used as covariates in the GWAS to avoid the false positive caused by population stratification. The genomic inflation factor of the logistic mixed model is 1.006, indicating that the rates of false positive errors are well controlled. The decay of LD was the strongest within the first 50 kb (Figure 2B). These SNPs corresponding to the nearest genes within the LD decay range were selected for identification of the candidates.

The HVP heritability was 0.32 in this F₂ population of yellow chicken, which illustrated that the HVP trait could be potentially improved by selective breeding. GWAS performed on HVP showed that 23 consistent SNPs were identified as the most significantly associated loci with HVP performance, which was analyzed by the logistic mixed model (Figure 4). Meanwhile, these SNPs were all from a narrow genomic location on chromosomes 1 (49.2–51.8 Mb). To some extent, the results were partly overlapped with our previous study (Luo et al., 2013).

There were 23 potential candidate genes located in this region, including MFNG, POLDIP3, POLR2F, PICK1, PDXP, SGSM3, RANGAP1, MYH9, RPL3, GALP3, LGALS1, MICALL1, ATF4, and CYP2D6 (Table 1). Of which, rs314284996 and rs317955795 were located in the promoter regions of CYP2D6 (less than 1.3 kb distance from 5′UTR). Interestingly, CYP2D6 was identified as a biomarker associated with increased susceptibility to melanoma, which is coded for another enzyme that may play a role in detoxifying potentially carcinogenic compounds (Bose et al., 2018). A quantile–quantile plot of the logistic mixed model analyzed for HVP traits showed that the observed and expected p-values of the GWAS fitted very well (Figure 5).

There were four haplotype blocks, which spanned 0.80, 0.09, 0.06, and 0.03 kb genomic regions (Table 2), and they were identified to have significant association with HVP traits (Figure 6), with a conservative LD threshold (r ² = 0.8). As shown in Table 3, the haplotypic Block 1 included two significant SNPs, with only two SNP haplotypes constructed. The AG haplotype has a significant minus effect on HVP (p < 0.01), and the GT haplotype has a great plus effect on HVP (p < 0.01). However, a total of two haplotypes were generated in Block 2, and the two haplotypes were significantly associated with HVP. The AT haplotype has a significant minus effect on HVP (p < 0.01), and the GC haplotype has a great plus effect on HVP (p < 0.01). The haplotypic Block 3 included four significant SNPs, with only two SNP haplotypes constructed. The TG haplotype has a significant minus effect on HVP (p < 0.01), and the CA haplotype displayed a plus effect on HVP (p < 0.01). However, a total of four haplotypes were generated in Block 4 with two haplotypes significantly associated with HVP traits. The CT haplotype has a significant minus effect on HVP (p < 0.01), and the TC haplotype has a great plus effect on HVP (p < 0.01).

Exploring expression level changes in candidate genes was proved to be helpful to identify major genes following GWAS analysis. Here, five birds were randomly chosen from both the severe HVP and absent HVP groups and were utilized for detecting transcriptional expression profiles of CYP2D6, MYH9, RANGAP1, and MFNG by qRT-PCR. The primers are listed in Table 4. The results showed that both these functional genes were expressed in the visceral peritoneum tissues (Figure 7). Only the CYP2D6 gene was greatly differentially expressed between the birds from the absent and severe HVP groups (p < 0.05), and CYP2D6 had a much high level in the server individuals, which was consistent with the studies in humans that CYP2D6 genotypes were associated with the hair color, Breslow thickness, and malignant melanoma (Strange et al., 1999). For the MFNG gene, coding an oncogene through Notch-mediated signaling (Kakuda et al., 2020), there was no statistically significant difference in birds with various HVP performances (p > 0.05). For the MYH9 gene, which encodes the non-muscle myosin heavy-chain IIA protein (Pecci et al., 2018), there was no statistically significant difference in birds with various HVP performances (p > 0.05). For the RANGAP1 gene, a crucial component of the RanGTPase system (Lin et al., 2016), there was no statistically significant difference in birds with various HVP performances (p > 0.05).