Genome-Wide Association Study Revealed the Effect of rs312715211 in ZNF652 Gene on Abdominal Fat Percentage of Chickens

Simple Summary With intensive selection in broilers, excessive abdominal fat accumulation is also present and causes economic concerns. Abdominal fat percentage (AFP) is one of the main indices of abdominal fat traits. We identified key SNP and candidate gene affecting AFP by a genome-wide association study (GWAS). Additionally, the main findings show that rs312715211 on the ZNF652 gene can increase body weight (BW), reduce eviscerated carcass weight (ECW), and increase abdominal fat percentage (AFP). Abstract Abdominal fat percentage (AFP) is an important economic trait in chickens. Intensive growth selection has led to the over-deposition of abdominal fat in chickens, but the genetic basis of AFP is not yet clear. Using 520 female individuals from selection and control lines of Jingxing yellow chicken, we investigated the genetic basis of AFP using a genome-wide association study (GWAS) and fixation indices (FST). A 0.15 MB region associated with AFP was located on chromosome 27 and included nine significant single nucleotide polymorphisms (SNPs), which could account for 3.34–5.58% of the phenotypic variation. In addition, the π value, genotype frequency, and dual-luciferase results identified SNP rs312715211 in the intron region of ZNF652 as the key variant. The wild genotype was associated with lower AFP and abdominal fat weight (AFW), but higher body weight (BW). Finally, annotated genes based on the top 1% SNPs were used to investigate the physiological function of ZNF652. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis suggested that ZNF652 may reduce AFW and BW in broilers through the TGF-β1/SMad2/3 and MAPK/FoxO pathways via EGFR and TGFB1. Our findings elucidated the genetic basis of chicken AFP, rs312715211 on the ZNF652 gene, which can affect BW and AFW and was the key variant associated with AFP. These data provide new insight into the genetic mechanism underlying AF deposition in chickens and could be beneficial in breeding chickens for AF.


Introduction
Carcass traits of chickens are economically important, and include abdominal fat (abdominal fat weight, AFW; abdominal fat percentage, AFP), body weight (BW), eviscerated carcass weight (ECW), etc. Broilers typically contain 150-200 g of fat per kg of BW. Abdominal fat accounts for 22% of body fat, but is generally considered waste, due to its low economic value [1]. Excessive deposition of abdominal fat in broilers not only results in low feed conversion rate, fertility, and semen quality, but also affects the economy of the

Genome-Wide Association Study
We applied a linear mixed model (LMM) to conduct a GWAS examining the traits underlying AFP in the selection (n = 252) and control (n = 268) lines, using GEMMA software [19]. The population structure was assessed via principal component analysis (PCA), using PLINK 1.9 to correct for differences between the two lines. The genetic relationship matrix and the first three PCA were added for correction during the LMM calculation based on the model: y = Wα + xβ + u + e, where y represents the phenotypic values of n samples, W represents the fixed effect matrix, α represents the fix effect (population structure, including PCA1, PCA2, and PCA3), x represents the SNPs, β represents the effects of corresponding markers, u represents the random effect vector, and e represents a vector of random residuals. Finally, we analyzed single loci one-by-one and calculated p-values using a derived Wald test. The Bonferroni correction multiple test was used to determine the significance threshold after LD linkage modified SNPs. The sums of the independent LD blocks plus singleton markers were used to define the number of independent statistical comparisons [20]. Finally, 378,446 independent SNPs were used to determine the p-value thresholds, including genome-wide significance (−log10(0.05/378,446)) and suggestive association (−log10 (1/378,446)). The Manhattan plots of GWAS for AFP were produced using the "qqman" package in R (https://www.r-project.org/) (accessed on 10 July 2022).

Fixation Indices (F ST ) and Heterozygosity (π)
Selection signals combined with nucleotide polymorphism analysis can reveal the genetic mechanisms underlying population evolution. Based on the SNPs remaining after quality control, F ST and π analyses were performed on the selected and control lines using Vcftools (V0.1.13) [21]. In F ST analysis, a single SNP was used as the step size for genomewide scanning. The top 5% of F ST values were defined as selected loci between the selection and the control lines. Then, π values were calculated for the significant SNPs screened by F ST . During π analysis, the window 40,000 and step 10,000 were used to calculate the region.

Dual-Luciferase Reporter Assay
To quantify the interaction of rs312715211 and its potential target gene ZNF652, we constructed vectors based on PGL4.18 plasmids. According to the chicken 6.0 reference genome sequence, we constructed a 1500 bp promoter sequence upstream of the transcription start site and designated it PGL4.18-ZNF652-pro. We also constructed wild and mutant luciferase reporter vectors of rs312715211 and designated them PGL4.18-SNP-TT and PGL4.18-SNP-AA. Next, chicken embryo fibroblasts (DF1) cells were cultured on a 24-well plate and were co-transfected with 250 ng of PGL4.18-ZNF652-pro and 250 ng of PGL4.18-SNP-TT or PGL4.18-SNP-AA. After 24 h, 100 µL lysate was added according to the Dual Luciferase Reporter Gene Assay Kit (Promega, Madison, WI, USA) and centrifuged for 15 min. Firefly and Renilla activities were measured using 20 ul supernatant, and each group was replicated three times.

RNA-Seq and Weighted Gene Correlation Network Analysis (WGCNA)
Transcriptomic sequencing was performed on 98-day-old Wenchang chickens' abdominal fat, including 18 samples. After eliminating abnormal individuals, high and low AFP groups were selected. For detailed sequencing steps, please refer to previous studies [22]. Then, R package "WGCNA" was used to analyze the weighted gene correlation network of five transcriptomes from abdominal fat tissue of chickens with high or low AFP. A weighted co-expression network was constructed using β = 7 to calculate the adjacency between genes. In addition, parameters mergeCutHeight = 0.25 and minmodule-size = 30 were selected for calculation, and gene significance, correlation of modules, and gene expression profile were calculated. The RNA-seq data describing the abdominal fat are included in a previous report by our group and are available at https://bigd.big.ac.cn/gsa/ (accessed on 10 July 2022) (accession data code CRA006031).

Kyoto Encyclopedia of Genes and Genomes (KEGG)
The top 1% of SNPs were selected for gene annotation, and these genes were enriched. The genes of the significant modules analyzed from the WGCNA were also annotated and then enriched. The KEGG database (http://www.genome.jp/kegg) (accessed on 10 July 2022) [23] is an important public database used for metabolic analysis and regulatory network research. KEGG enrichment analysis of annotated genes was performed using KOBAS, and KEGG pathways with a Q value (corrected p value) ≤ 0.05 were considered significantly enriched. The results were drawn using the "ggplot2" package in R.

Statistical Analysis
R 4.0.4 and SAS 9.4 (SAS Institute, Cary, NC, USA) were used to generate descriptive statistics and for normal distribution tests of AFP and to test the significance of the differences between the groups using Student's t-test. Confidence limits were set at 95%, and p < 0.05 (*) or < 0.01 (**) was considered significant. Data are presented as mean ± standard error. ASReml 3.0 software was used to estimate genetic variance, genetic correlation, and heritability. Genetic parameter estimation was based on the animal single trait model using restricted maximum likelihood (REML) and the model: Y = Xb + Za + e, where Y represents the observed value of traits, b represents the fixed effect vector, and X represents the incidence matrices of fixed effects. a represents the additive genetic effect vector of the individual, Z represents the incidence matrices of the additive genetic effect of the individual, and e represents the random residual effect vector.
The phenotypic variation explained (PVE) can be estimated using equation [24]: where β represents the effect value for the GWAS result, MAF represents SNP minor allele frequency, and N represents the number of individuals included in the GWAS analysis.

Phenotypic Statistics and Heritability Evaluation
Descriptive statistics for AFP, ECW, and AFW for both the selection and control lines are presented in Table 1 and Figure 1. The results showed that the AFP, AFW, ECW, etc., of the selection line was higher than control line after selection. Genetic parameter analyses revealed that AFP has high heritability-the heritability was 0.64. Moreover, AFP, BW, and IMF showed positive genetic correlation and phenotypic correlation. The AFW and ECW were significantly (p < 0.01) increased in the selection line, compared to the control line, resulting in no significant difference in AFP (AFP correlates with the AFW/ECW ratio) between the selection and control lines. of the selection line was higher than control line after selection. Genetic parameter analyses revealed that AFP has high heritability-the heritability was 0.64. Moreover, AFP, BW, and IMF showed positive genetic correlation and phenotypic correlation. The AFW and ECW were significantly (p < 0.01) increased in the selection line, compared to the control line, resulting in no significant difference in AFP (AFP correlates with the AFW/ECW ratio) between the selection and control lines.

GWAS Identified the Effective Variants and Candidate Genes
To find significant variation at the genomic level, we performed GWAS of AFP. The GWAS results are summarized in Figure 2 and Table 2. The chromosomal significance threshold was −log 10 (0.05/378,446). The potential significance level threshold was −log 10 (1/378,446). These analyses showed that the significant SNPs associated with AFP phenotype were located on chromosomes 1, 14, and 27. An approximately 0.15 MB region on chromosome 27 (Chr27: 5,963,734-6,119,680) was strongly associated with AFP and included nine significant SNPs with AFP, which were annotated on IGF2BP1, ZNF652, GIP, UBE2Z, and ETV4, respectively. The most significant SNP (rs312351828) accounted for 5.35% of the observed phenotypic variance. Taken together, all of the SNPs that reached the GWAS threshold explain 3.34−5.58% of the observed phenotypic variance.

rs312715211 in the Intron Region of ZNF652 Was the Primary Variant Associated with AFP
To further identify candidate SNPs and genes, the F ST and π values of all SNPs located on chromosome 1, 14, and the candidate region (Chr27: 5,963,734-6,119,680) were calculated. F ST is calculated with a single SNP as the step size, and the F ST threshold was 0.1. These data are summarized in Figures 3 and 4. F ST analysis showed that SNPs reaching the GWAS threshold line on chromosomes 1, 4, and 14 were not selected ( Figure 3).
On chromosome 27, F ST analysis showed that 159 SNPs in the target region were selected ( Figure 4 and Table S1), and this region included the ZNF652, IGF2BP1, ETV4, DHX8, GIP, PHOSPHO1, SNF8, ABI3, and UBE2Z genes. However, only rs312715211 reached the threshold for GWAS and F ST ( Figure 4A). We calculated the π values in this region, in order to confirm whether the candidate sites were selected. The π values represent nucleotide polymorphisms, which decreased after selection; therefore, the π value was less in the selection than in the control line ( Figure 4B). The results show that the two populations were strongly selected near the location of chr27:6,000,000, and rs312715211 was in this region. We have calculated the single point π value of rs312715211, as shown in Table 3. Taking all the results together, rs312715211, located on intron2 of ZNF652, reaches the thresholds of GWAS, F ST , and π; therefore, it can be considered a key variant associated with AFP.

rs312715211 in the Intron Region of ZNF652 Was the Primary Variant Associated with AFP
To further identify candidate SNPs and genes, the FST and π values of all SNPs located on chromosome 1, 14, and the candidate region (Chr27: 5,963,734-6,119,680) were calculated. FST is calculated with a single SNP as the step size, and the FST threshold was 0.1. These data are summarized in Figures 3 and 4. FST analysis showed that SNPs reaching the GWAS threshold line on chromosomes 1, 4, and 14 were not selected (Figure 3). On chromosome 27, FST analysis showed that 159 SNPs in the target region were selected ( Figure 4 and Table S1), and this region included the ZNF652, IGF2BP1, ETV4, DHX8, GIP, PHOSPHO1, SNF8, ABI3, and UBE2Z genes. However, only rs312715211 reached the threshold for GWAS and FST ( Figure 4A). We calculated the π values in this region, in order to confirm whether the candidate sites were selected. The π values represent nucleotide polymorphisms, which decreased after selection; therefore, the π value was less in the selection than in the control line ( Figure 4B). The results show that the two populations were strongly selected near the location of chr27:6,000,000, and rs312715211 was in this region. We have calculated the single point π value of rs312715211, as shown in Table 3. Taking all the results together, rs312715211, located on intron2 of ZNF652,

rs312715211 Can Affect the Activity of ZNF652 Promoter
We speculated that rs312715211 affects the expression of ZNF652, so a double-luciferase validation was performed to verify the regulation of rs312715211 on expression of ZNF652. Then, PGL4.18/ZNF652pro/TT/AA dual-luciferase plasmids were transferred into 293T cells. When rs312715211 was wild genotype (TT), ZNF652 promoter activity was Biology 2022, 11, 1849 9 of 16 not changed ( Figure 5A). rs312715211 significantly increased the activity of the ZNF652 promoter, when the rs312715211 wild genotype TT was mutated to AA. The binding of transcription factors, meanwhile, was predicted using the PROMO online website after the rs312715211 mutation and wild-type. We found that most transcription factors were changed after mutation at this site ( Figure 5B), which may affect gene expression. The results show that the rs312715211 region not only acts as an enhancer, but may be the transcription factor's binding region.

rs312715211 Can Affect the Activity of ZNF652 Promoter
We speculated that rs312715211 affects the expression of ZNF652, so a double-luciferase validation was performed to verify the regulation of rs312715211 on expression of ZNF652. Then, PGL4.18/ZNF652pro/TT/AA dual-luciferase plasmids were transferred into 293T cells. When rs312715211 was wild genotype (TT), ZNF652 promoter activity was not changed ( Figure 5A). rs312715211 significantly increased the activity of the ZNF652 promoter, when the rs312715211 wild genotype TT was mutated to AA. The binding of transcription factors, meanwhile, was predicted using the PROMO online website after the rs312715211 mutation and wild-type. We found that most transcription factors were changed after mutation at this site ( Figure 5B), which may affect gene expression. The results show that the rs312715211 region not only acts as an enhancer, but may be the transcription factor's binding region.

rs312715211 Mutation Can Increase AFP and AFW and Decrease ECW
AFP consists of AFW and ECW. We compared the genotype frequency and phenotype (AFP, AFW, and ECW) of individuals carrying wild or mutant genotypes of rs312715211 SNP. According to these results, rs312715211 caused significant differences in the AFP, AFW, and ECW phenotypes between the selected and control lines. We also confirmed that rs312715211 can increase AFW and reduce ECW, thereby affecting AFP ( Figure 6B). Both the genotype frequency and the phenotype associated with this SNP changed significantly after selection. The frequency of rs312715211 in the wild genotype increased after selection, as did AFP. The frequency of rs312715211 in the mutant genotype decreased after selection, while AFP increased ( Figure 6A). rs312715211 SNP. According to these results, rs312715211 caused significant differences in the AFP, AFW, and ECW phenotypes between the selected and control lines. We also confirmed that rs312715211 can increase AFW and reduce ECW, thereby affecting AFP ( Figure 6B). Both the genotype frequency and the phenotype associated with this SNP changed significantly after selection. The frequency of rs312715211 in the wild genotype increased after selection, as did AFP. The frequency of rs312715211 in the mutant genotype decreased after selection, while AFP increased ( Figure 6A).

Identification of Candidate Genes and Pathways Related to AFP
To identify the candidate genes and pathways related to AFP, we widened the field of investigation. The top 1% of SNPs associated with AFP were screened accordingly, based on the GWAS results, and a total of 4736 genes were annotated and subjected to KEGG pathway enrichment analysis. KEGG enrichment analysis indicated that 52 pathways were significantly enriched (p < 0.05) (Figure 7 and Table S2). Among them, the metabolic pathways, MAPK, neuroactive ligand-receptor interaction, and calcium signaling pathways were the most significant. Furthermore, some pathways related to fat or carbohydrate metabolism were significantly enriched, including PPAR, adipocytokine, glycerophospholipid metabolism, glycolysis/gluconeogenesis, and other glycan degradation signaling pathways. Unfortunately, ZNF652 was not enriched in the related pathways. However, studies have shown that TGFB1 and EGFR are the target genes of ZNF652 and can be directly recruited and bound [25]. We also found that EGFR and TGFB1 were the target genes of ZNF652 and significantly enriched in 13 pathways (p < 0.05) ( Table 4). Among these pathways, EGFR and TGFB1 were simultaneously enriched in the MAPK and FoxO signaling pathways.

Identification of Candidate Genes and Pathways Related to AFP
To identify the candidate genes and pathways related to AFP, we widened the field of investigation. The top 1% of SNPs associated with AFP were screened accordingly, based on the GWAS results, and a total of 4736 genes were annotated and subjected to KEGG pathway enrichment analysis. KEGG enrichment analysis indicated that 52 pathways were significantly enriched (p < 0.05) (Figure 7 and Table S2). Among them, the metabolic pathways, MAPK, neuroactive ligand-receptor interaction, and calcium signaling pathways were the most significant. Furthermore, some pathways related to fat or carbohydrate metabolism were significantly enriched, including PPAR, adipocytokine, glycerophospholipid metabolism, glycolysis/gluconeogenesis, and other glycan degradation signaling pathways. Unfortunately, ZNF652 was not enriched in the related pathways. However, studies have shown that TGFB1 and EGFR are the target genes of ZNF652 and can be directly recruited and bound [25]. We also found that EGFR and TGFB1 were the target genes of ZNF652 and significantly enriched in 13 pathways (p < 0.05) ( Table 4). Among these pathways, EGFR and TGFB1 were simultaneously enriched in the MAPK and FoxO signaling pathways.

ZNF652 May Regulate Abdominal Fat and BW through MAPK/FoxO Signaling Pathways
To verify the above results, WGCNA was performed on the transcriptomes of the high and low AFP chickens. We found that ZNF652 was expressed less in high AFP group than in the AFP low group, consistent with the dual-luciferase result ( Figure 8A,B). We also found that AFP, AFW, BW, and ECW were simultaneously mapped in a significant module (darkseagreen2) (p < 0.05) ( Figure 8C). Although ZNF652 was not enriched in this significant module, it was co-expressed with classic lipid metabolism genes, such as LIPN1, LIPN2, GATA3, PPARG, DGKQ, and DGKE. Another KEGG analysis was performed on the enriched genes in the significant module, and a total of 19 significantly pathways were enriched (p < 0.05). Of these, the MAPK and FoxO signaling pathways were consistent with our above results ( Figure 8D). In conclusion, we speculate that ZNF652, similar to the lipid metabolism genes, plays an important role in fatty acid degradation and may intervene in the MAPK and FoxO signaling pathways by binding with target genes, thus affecting body weight by decreasing AFP and AFW.

Discussion
Excessive abdominal fat in broilers not only reduces reproductive performance and causes metabolic diseases, but also reduces meat quality [26,27]. An increasing number of researchers are focusing on the genetic mechanisms underlying abdominal fat deposition. In this study, the heritability of AFP was 0.64, similar to the results reported by Demeure, Duclos et al. [7] and Chabault, Baéza et al. [6].
With the application of GWAS and selection signals, AFP in chickens has been fully utilized and SNPs and QTLs associated with abdominal fat have been identified [28]. According to the database, 292 QTLs and numerous SNPs were located in abdominal fat traits [2]. These QTLs and SNPs were verified in several broiler breeds and were significantly associated with fat deposition [29,30]. We attempted to identify the major genetic markers related to abdominal fat and growth traits in broilers. GWAS, F ST , and π analysis identified one significant SNPs on chromosome 27, and the SNP differed significantly between wild and mutant individuals in phenotype (p < 0.05). This is consistent with Zhang [3] and suggests that the SNPs are important genetic markers for reducing abdominal fat deposition in broilers.
As a result, rs312715211 conformed with our expectation. The SNP was mapped to ZNF652 and may simultaneously regulate abdominal fat deposition and growth development in chickens. rs12715211 can increase AFP, AFW, and BW in broilers after mutant. After selection, AFW significantly increased and BW significantly decreased. Moreover, rs312715211 significantly increased ZNF652 promoter activity after mutant. In order to improve chicken breeding and increase its economic benefits, we should select wild homozygous individuals from rs312715211 for breeding and eliminate mutant homozygous individuals, so as to achieve increased body weight and reduce abdominal fat.
In the target region, ZNF652 and IGF2BP1 may be important for the control of abdominal fat deposition in broilers. IGF2BP1 (insulin-like growth factor 2 mRNA-binding protein 1) can bind IGF2 mRNA and is a member of the single-stranded RNA-binding protein family [31]. IGF2BP1 can function by affecting cell proliferation, migration, and apoptosis [32]. IGF2BP1 knockout mice can reduce Ramp3 mRNA expression; thus, IGF2BP1 indirectly affects glucose and lipid metabolism [33]. In chicken pan-genomic studies, a high expression of IGF2BP1 lead to high weight [34]. IGF2BP1 can increase body weight and affects feed efficiency in ducks [35], promotes adipocyte proliferation and differentiation, and regulates expression of genes related to fatty acid metabolism in broilers [36]. IGF2BP1 is also associated with tumors, cancer, dwarfism, and other diseases, but its mechanism of fat regulation requires further study.
ZNF652 is the classical C2H2 zinc finger DNA binding protein [37] and an inhibitor of gene transcription and plays a primary role in tumor invasion [25]. ZNF652 may be involved in human lipid and carbohydrate metabolism, e.g., sex hormone binding globulin (SHBG) [38] and may be related to body weight and bone growth [39,40]. However, its fat and growth regulatory mechanism in chickens is not currently known. It has been reported that the inhibition of ZNF652 can increase the expression of EGFR and TGFB1, and EGFR can regulate fat deposition by regulating fatty acid synthase [38,41]. The inhibition of EGFR expression can significantly reduce BW and subcutaneous and abdominal fat mass in mice, and more importantly, SREBP-1 and FASN expression decreased after EGFR inhibition [42]. EGFR is an epidermal growth factor receptor that regulates fat metabolism genes, especially PPAR, and the inhibition of EGFR can reduce the expression of adipose synthase [41]. TGFB1 (transforming growth factor β) [43] can bond with the SMAD family. Our studies showed that AFW and BW decreased after ZNF652 expression increased, consistent with the results after EGFR inhibition. Previous reports suggest the inhibition of the TGFB1 gene by PPARG activation, and PPARG can inhibit cell proliferation through the TGF-β1/Smad2/3 signaling pathways [44]. Therefore, we hypothesized that the effect of ZNF652 was associated with EGFR, TGFB, and PPARG by targeting the TGF-β1/SMad2/3 signaling pathways and MAPK/FoxO signaling pathways, resulting in reduced proliferation of adipocytes and accelerated fatty acid oxidation.

Conclusions
Our findings elucidated the genetic basis of chicken AFP. rs312715211 in the ZNF652 gene is the key variant associated with AFP of chickens, and the wild genotype is the favorable genotype to lower AFW and heighten BW. Additionally, ZNF652 is the key gene related to AFP by affecting BW and AFW. These data provide a new insight into our understanding of the genetic mechanisms underlying abdominal fat deposition in chickens and will aid in the breeding of broilers with lower AFP. In the future, the further study of these genetic variations may be applied in marker-assisted selection to reduce abdominal fat deposition in broilers.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/xxx/s1, Table S1: Reach the F ST threshold line of SNPs, Table S2: KEGG of top 1% gene pathway.  Informed Consent Statement: Not applicable.

Data Availability Statement:
The whole-genome resequencing datasets generated in this study were submitted to https://bigd.big.ac.cn/gsa (accessed on 10 July 2022), with IDs CRA002643 and CRA00265. The RNA-resequencing datasets of this study have been deposited at https://bigd.big.ac. cn/gsa/, accessed on 10 July 2022, (accession data code CRA006031).

Conflicts of Interest:
The authors declare that they have no competing interests.