INTRODUCTION
The increasing proportion of energy contributed by animal foods in daily diet benefits from the improvement in standards of living as well as an increase in awareness of meat nutritive values [1]. Global demand for meat has increased by 70 million tons in 1961 and 356 million tons in 2022. Meat contains about 20% protein and is abundant in essential micronutrients (zinc, iron, selenium, B vitamins, etc.), which make it as an important food for human nutrition [2]. A moderate increase in meat consumption will certainly improve the nutritional adequacy of diets and health outcomes. While several large-scale epidemiological studies have showed that a high consumption of red and especially processed meat with a higher risk of developing various chronic diseases [3]. Poultry meat, high in protein and low in fat, has not been shown to be associated with the above diseases [4]. Nowadays poultry meat consumption is steadily increasing worldwide.
Among poultry meat products, chicken carcasses, cuts, and processed products are the most consumed (~75% of total poultry meat) followed by turkey (~25%) and, to a lesser extent, duck [5]. The heritability of chicken carcass traits was moderate to high and differed between populations and traits, for example, the estimations of heritability for body weight, breast muscle weight, abdominal fat weight ranged from 0.17 to 0.34 [6]. Using genetic markers associated with the genes involved in carcass traits, genetic improvement could be achieved efficiently on these traits. Genome-wide association studies (GWAS) allow the identification of single-nucleotide polymorphisms (SNPs) associated with effective genes that influence these traits, providing a better biological understanding of the trait and a list of candidate genes for fine mapping [7], especially for chickens, the causal genes related to breast muscle, meat flavor, morphological traits, etc. have been revealed [8–10]. However, these causal genes may not be effective in all populations due to differences in selection backgrounds and demographic histories among different chicken breeds. Therefore, under the premise of existing reference genes, it is necessary to perform GWAS for carcass traits in specific chicken populations.
Currently, partial genes have been proved to be effective for carcass traits. For example, the IGF2BP1 gene plays an important role in regulating body size and carcass traits in both duck and chicken [11–13]. After long-term domestication and selection, the genetic variants regulating IGF2BP1 transcription were usually fixed and conserved in commercial chickens, while it’s dispersed in various native chickens. In this study, we aimed to detect genomic peaks and candidate genes associated with carcass traits based on whole-genome sequencing in Sanhuang (SH) chickens. Simultaneously, we compared the expression level of candidate genes in multiple tissues and chicken breeds to refine the functional genes profile, which provided novel insights into genetic mechanism of production traits and facilitated the genomic breeding in the view of biomarker selection associated with carcass traits.
MATERIALS AND METHODS
All procedures followed the guidelines of the Institutional Animal Care and Use Committee at Zhejiang A&F University, China (No. ZAFUAC202434).
Experimental chickens and phenotypic measurement
A total of 565 SH chickens were used in this study and all birds were raised in three-stair step cages (one bird per cage) under the same recommended environmental and nutritional conditions. Basal diets (containing 15.2% crude protein and 18.3 MJ/kg gross energy) were formulated based on the National Resource Council requirements, feed and water were provided ad libitum during this study. Blood samples were obtained using standard venipuncture techniques at 91 days of age. Anticoagulated blood samples were stored at −20°C for genomic DNA extraction. After blood collection, all chickens were euthanized under carbon dioxide anesthesia by severing the carotid artery. After slaughter, slaughter weight (SW), eviscerated weight (EW), and thigh muscle weight (ThW) of each chicken were weighted and recorded. The percentage of these three traits (percentage of slaughter weight [SP], percentage of eviscerated weight; [EP], percentage of thigh muscle weight [ThP]) were calculated according to the standard procedures.
DNA extraction and whole-genome sequencing
DNA was extracted from blood samples using the phenol-chloroform method. A total of 565 DNA samples were used for whole-genome resequencing. The quality and quantity of DNA were examined using a NanoDrop device and by agarose gel electrophoresis. Then, paired-end libraries were generated for each eligible sample using standard procedures. All libraries were sequenced on an Illumina HiSeq X Ten platform and each individual at less with over 10 G raw read sequence.
Genotyping, quality control, and imputation
Filtered FASTQ data (parameters: −q 30; −u 30; −l 150) (https://github.com/OpenGene/fastp) were aligned to the reference genome [Gallus_gallus-6.0 (GCA_000002315.5)] using BWA MEM (v 0.7.10) with default parameters [14]. Mapped reads were converted into BAM files using SAMtools (v 0.1.18), and duplicate reads were removed by SAMtools and Picard MarkDuplicates (http://broadinstitute.github.io/picard). After mapping, SNP calling was performed using exclusively GATK (v 3.5), and the output was further filtered using VCFtools (v 0.1.15) [15,16]. Genotype quality control was conducted using PLINK (v 1.90) [17] to detect and exclude unreliable genotypes. SNPs were selected based on SNP maximum missing rate <5%, minor allele frequency >0.05. Moreover, individuals were excluded due to sample call rate <95%. Totally, 9402908 SNPs on autosomes were retained and missing alleles were imputed using Beagle 5.2 [18].
Genome-wide association study for carcass traits
To minimize false positives, it is essential to account for population structure and family effects. Mixed linear model (MLM) is an effective approach to address population stratification in GWAS. An efficient MLM was used for the association analysis using the GEMMA software [19]. SNP effect was set as fixed effect and kinship matrix as random effect. The GWAS was performed for carcass traits as follows:
where y indicates the phenotypic record of carcass traits; W indicates the covariates matrix, including the top three principal component analysis value; α is vector of the corresponding coefficient (including the intercept); βis the SNP effect; μ~MVNn (0, λτ−1K), μ indicates random polygenic effect, ε is the residual error.
The threshold of genome-wide significance was assessed by simpleM and the Bonferroni correction methods [20]. The significance threshold was 0.05/339733 (−log10 p = 6.83) and the suggestive significance threshold was 1/339733 (−log10 p = 5.53). In this study, SNPs that reached the suggestive significance threshold were identified as significant SNPs associated with traits. Within the peak, the SNP information and candidate genes were presented.
Estimation of heritability
Genomic SNPs information were used to generate genomic relationship matrix (GRM) by the approach from VanRaden [22], and variance component estimates for traits were calculated based on univariate animal model constructed by ASReml v4.1. The univariate model was defined as follows:
where y is the vector of phenotype, b is the vector of fixed effects, including batch effect, a is the vector of random additive genetic effects, e is the vector of random residual effect. X and Z are design matrices relating observation to the corresponding fixed and random effects, respectively. The variance-covariance structure assumed for random effect was following:
where G is the GRM based on SNP data, I is the identity matrix, σa2 is the additive genetic variance, σe2 is the residual variance. The heritability of phenotypes was calculated by: h2 = σa2/(σa2 + σe2).
Estimation of genetic and phenotypic correlation
A bivariate animal model was used to evaluate the genetic correlation and phenotypic correlation between any two traits:
where y, a, b, e, X, and Z were same as above equation.
RNA extraction and mRNA sequencing
Total RNA was extracted from breast muscle, thigh muscle, lung, liver, heart, and fat (n = 5) tissues in SH chickens of 91 d using TRIzol reagent (15596026CN; Thermo Fisher Scientific, Waltham, MA, USA). And breast muscle from commercial broilers and SH chickens (n = 6) of 40 d was also conducted RNA isolation. A total of 42 RNA samples were used to mRNA library construction and sequencing based on Illumina platform (PE150), and >6 G raw reads per sample were produced.
Transcriptomic analysis
The bioinformatic pipeline of mRNA sequencing followed the methodology outlined in previous report [10]. Briefly, the raw reads were firstly trimmed by Fastp (v0.21) software (−q 20 −u 30 −n 5 −l 150) [23], and qualified reads were used to genome alignment (Gallus_gallus-6.0 [GCA_000002315.5]) using HISAT2 (v2.2.0) software [24] with default parameters. After data conversion, quality control, and index, transcript assembling and quantification were implemented by StringTie (v2.1.6) software (−l 150) [25]. Differentially expressed genes were defined when fold change > 1.5 or < 0.67 and p<0.05 using DESeq2 software [26].
Real-time polymerase chain reaction for candidate genes
Gene expression of partial candidate genes (e.g., IGF2BP1, TBXAS1, SKAP1) for carcass traits were examined with real-time polymerase chain reaction (RT-PCR). Total RNA from breast muscle of commercial broilers and SH chicken (n = 6) was isolated using TRIzol reagent (15596026CN, Thermo Fisher Scientific). The measurement of RNA concentration and integrity was conducted with Agilent 2100 Nano and gel electrophoresis. Next, the total RNA was reverse transcribed to cDNA using FastKing RT kit (KR116; TIANGEN, Beijing, China). And the reaction system of RT-PCR was created using SYBR Green Pro Taq HS kit (AG11740; Accurate Biology, Changsha, China). The specific primers were designed based on Oligo 7.0 software (Supplement 1). The condition detail was as follows: 94°C for 10 min, followed by 35 cycles of 94°C for 3 s and an annealing temperature of 32 s. β-actin was regarded as the internal controls. The expression level was compared using the 2−ΔΔCT method.
Statistical analysis
Significance of expression differences between commercial broilers and SH chickens was tested by the Student t-test, and the phenotypic alteration (e.g., EW, SW, ThW) among different chickens with various haplotypes was conducted by Kruskal-Wallis test using the SPSS Version 22.0 (IBM, Armonk, NY, USA). Confidence limits were set at 95% and p<0.05 (*) or p<0.01 (**) was considered significant. Data are presented as the mean±standard deviation.
RESULTS
Descriptive statistics and heritability estimation
Descriptive statistics of the carcass traits (SW, SP, EW, EP, ThW and ThP) were presented in Table 1. The average SW, EW, and ThW in SH population was 1,107.97 g, 923.12 g, and 181.50 g, respectively, the average SP and EP was more than 88% and 73%, respectively, and the ThP ranged from 11.94% to 24.64%. The coefficients of variation (CV) of relative weight (SP, EP, ThP) were maintained below 10%, while the CVs of SW, EW, and ThW were relatively high (>10%). All these traits were followed a normal distribution (Figure 1).
Genetic and phenotypic correlations among carcass traits
Estimates of heritability, genetic and phenotypic correlation based on genomic information were shown in Table 2. The estimates of heritability of six traits were calculated and ranged from 0.08 to 0.50, five of them (besides SP) were classified as traits with moderate to high heritability. High estimates of genetic correlation among SW, EW, and ThW were obtained (ranged from 0.83 to 0.95), and higher results were also found between SW and SP, ThW and ThP, EW and EP. As expected, similar results of phenotypic correlation among these traits were detected. High phenotypic correlation among SW, EW, and ThW was presented (0.72 to 0.93), a moderate to high phenotypic correlation between SW and SP, EW and EP, ThW and ThP was obtained. Additionally, there was a relative low correlations among SP, EP, and ThP.
Genome-wide association studies for slaughter weight and percentage of slaughter weigh
For SW trait, only 3 associated SNPs were detected in GGA1, GGA3, and GGA13 (Figure 2A). The SNP on GGA1 was calculated the high SNP effect (−65.54) and low allele frequency (0.06), while the SNP on GGA13 was identified with a high phenotypic variance explained (PVE) (4.94%) (Table 3). These 3 SNPs were located within RNASEH2B, MMS22L, SPDL1, and DOCK2. For SP trait, a total of 8 significant peaks were found on GGA1, GGA2, GGA3, GGA7, and GGA11 (Figures 2B, 2C). Four of which were located on GGA1, and only one intergenic SNP was identified on GGA2 and GGA7, respectively. Among these SNPs, phenotypes can be explained up to 5.59% by the lead SNP (Table 3), and one major haplotype covering 18 SNPs was formed within the lead peak (Figure 2D), the major haplotype was found a significant enhancing effect on SP (Figure 2E). In different peaks, the genes TEAD4, FOXJ2, CCND2, ZNF384, and others were covered by significant SNPs.
GWASGenome-wide association studies for eviscerated weight and percentage of eviscerated weight
GWAS peaks for EW and EP were presented in Figures 3A, 3B. For EW trait, we found 5 peaks (GGA3, GGA8, GGA9, GGA15, and GGA27) were correlated to EW (Figure 3A). On GGA3, six significant SNPs spanned 5 Mb contributed more than 4.7% of phenotypic variation, and this peak included candidate genes USP45, PNISR, MMS22L, etc. (Table 4). On GGA27, a ~80 kb genomic interval covered 14 SNPs was found significant associations (Figure 3C), and an obvious haplotype including 8 SNPs was detected additive effect on EW (Figures 3D, 3E). IGF2BP1 gene, related to growth traits in farm animals, was annotated by the candidate SNPs. And GIP, SNF8, and others were also included in this peak. Additionally, only 1 to 2 associated SNPs located on GGA8, GGA9, and GGA15. For the EP trait, three prominent peaks were identified on GGA2, GGA13, and GGA27. Additionally, IGF2BP1 showed an association with EP (Supplement 2), and a haplotype located in the upstream of IGF2BP1 exhibited a significant additive correlation with EP (Supplements 2, 3). And the genes BMP6, PTPRM, SGCD, etc. on GGA2 and GGA13 were also correlated to EP. Additionally, single SNP on GGA4 and GGA9 was annotated intergenic variant and contributed more than 4.8% of PVE to EP (Table 4).
Genome-wide association studies for thigh muscle weight and percentage of thigh muscle weight
GWAS peaks for ThW and ThP were presented in Figures 4A, 4B. For ThW trait, five peaks and two isolated SNPs on different chromosomes were identified (Figure 4A). Similar to EW and EP, a peak spanned ~90 kb on GGA27 (including 33 SNPs) was highly correlated to ThW, and IGF2BP1, GIP, SNF8, etc. were covered by these SNPs (Figure 4C). Additionally, three haplotypes were also found important association with ThW (Figures 4D, 4E, Supplements 5, 6). On GGA1, the region from 55,911,047 to 56,322,618, including 71 SNPs, contributed more than 4.5% of PVE to ThW (Supplement 7; Table 5), three haplotypes were revealed in this peak, but a rare haplotype frequency related to enhancing ThW was detected (Supplements 8–11). By annotation, TBXAS1, HIPK2, and ENSGALG00000012796 were located within this region (Supplement 7). On GGA4, peak SNP was detected and contributed 6.85% of PVE to ThW, a rare haplotype was identified within this region and BMP3 gene was found (Supplements 12–14). Additionally, a short peak covering less than 10 SNPs were detected on GGA6 (1), GGA7 (1), GGA13 (7), and GGA15 (6), the genes ARHGAP26, CIT, NR3C1, etc. were overlapped by these candidate regions (Table 5). For ThP, two associations on GGA3 and GGA4 were obtained based on the threshold (Figure 4B), only ENSGALG00000015017 and ATRN genes were annotated. On GGA6 and GGA8, single SNPs were found a significant correlation to ThP, and BEND5 gene was annotated and putatively related to thigh muscle trait (Table 5).
Expression of candidate genes for carcass traits
Based on GWAS summary, a total of 73 protein-coding genes were annotated by candidate SNPs (Supplement 15). Firstly, an expression profile of multi-tissues was founded in SH chicken. And ACSL5, HSD3B1, BMP3, etc. were transcribed specifically in liver, lung, and other tissues (Figure 5A). Next, we compared expression of these genes in muscle of commercial broilers and SH chickens. A total of 17 of 73 genes were defined as differentially expressed gene (DEGs), 15 of which were up-regulated in commercial broilers, while ABI3 and MRPL22 were down-regulated (Figure 5B). Additionally, cell cycle related gene CCND2 was found a trend of increased expression (Fold change = 1.41, p = 0.051). To ascertain the accuracy of sequencing results, we performed quantitative polymerase chain reaction assay for partial genes. Similarly, IGF2BP1, BMP3, ACSL5, and RASGEF1B were highly transcribed in commercial broilers, while MRPL22 and ABI3 were restricted, no difference was detected between two chicken breeds for muscular GIP (Figures 5C–5I). We also confirmed transcriptional difference of CCND2 between two chicken breeds, and significant enhancement was detected in commercial broilers (Figure 5J). All these results proved the accuracy of sequencing results and conclusion.
DISCUSSION
GWAS is a powerful tool for the genetic analysis of important production traits in farm animals [9,27]. Carcass traits are important economic traits, which directly affect the meat production and play an important role in chicken industry. Previously study showed that the heritability of these traits in chicken was moderate to high [6]. Similarly, we also detected a high heritability of EW, SW, ThW, etc., except for SP (h2 = 0.08). Low selective pressure on SP trait may be the reason for low heritability, and limited sample size and absence of other genomic variants (e.g., structural variants) could also be important reasons. High genetic correlations (>0.7) between SW and SP, ThW and ThP were observed, as well as among the SW, EW, and ThW. Which was highly consistent with previous reports. However, it’s moderate between EW and EP (0.24). Totally, mutations associated with these traits can be obtained from GWAS and applied to genomic breeding to improve the performance of meat production.
Usually, GWAS peaks include the causal variants and adjacent SNPs due to the hitchhiking effect, herein, the LD blocks reduce the possibility of causal SNPs and resulted in a high false positive rate in GWAS. In current study, we obtained the effective calculations using simpleM method to avoid the bias from LD, and controlled the false positive rate by Bonferroni correction. Based on that, we found that the significant peaks related SW, EW, and ThW were not completely identical to the SP, EP, and ThP. In an unselected population, the different peak and low significance were often detected for related carcass traits [28], consistent associations with different carcass traits could be revealed in an F2, for example, the effect of terminal region of GGA4 and GGA27 on carcass traits (e.g., body weight, shank length, EW) [29].
For farm animals, carcass traits are typical complex phenotypes, which highly correlated to the muscle development. Genomic peaks on GGA27 have been found significant associations with carcass and growth traits [12,13,29], candidate genes IGF2BP1, GIP, PHOSPHO1, etc. were reported to be correlated to skeletal muscle (breast and thigh muscle) development. Especially for IGF2BP1, the causal role in body size and carcass traits has been proved in chickens and ducks, 15% of body size can be enhanced by regulating this gene. GIP gene is necessary for bone remodeling and bone quality, and activates osteoblasts and promotes bone development [30]. Similarly, BMP3 gene is highly correlated with bone and muscle growth [31], both of them contribute substantially to carcass traits. Cell cycle is basic physiological process affecting cell proliferation, differentiation, and repair after muscle injuries. CCND2 is responsible for the transition of cell cycle from G1 to S phase, and could contribute significantly to the skeletal muscle proliferation and differentiation by subjecting to targeting binding by miR-17 and miR-206 [32,33]. Here, the result that an increasing expression of CCND2 in breast muscle from commercial broilers was consistent with the production performance.
After long-term selection, carcass traits of commercial broilers significantly outperform that of native chickens, including SH chickens, however, this advantage comes at the price of decreased immunity and a set of undesirable performance [34]; additionally, native chickens provide superior meat flavor compared to commercial broilers. Therefore, growth and carcass traits were closely related to immune and metabolism functions [35], which were characterized by a complicated multi-level molecular regulatory network. Usually, carcass performance was negatively correlated with immune performance and lipid metabolism. RASGEF1B gene is involved in immune response during macrophage activation and protection against microbial infections, and could be induced in macrophages on stimulation with toll-like receptor agonists [36,37]. We found RASGEF1B was most affected between commercial broilers and SH chickens, supporting the consistency between carcass traits and immune performance. Interestingly, knockdown of circRASGEF1B could alter the cell cycle and mitotic division progression, which may influence the muscle development [36]. Fatty acid component is tightly correlated with the meat flavor and quality. ACSL5 is a vital factor in fatty acid elongation, degradation, and biosynthesis. Yuan et al. have proved that ACSL5 was associated with meat aldehydes and hexanal components and contributed to flavor formation in native chickens [38]. We have identified the relationship between ACSL5 and ThP, and ACSL5 was more enriched in commercial broilers compared to SH chickens. Demonstrating that lipid metabolism in leg had a significant impact on muscle development. On 42 d, commercial broilers were in a period of rapid muscle development with a high energy demand, while SH chickens were in the early developmental stage rather than the fat deposition stage [39]. Therefore, ACSL5 was more activated in commercial broilers.
In addition to the relationship between genotype to phenotype, we provided a transcriptional profile of candidate genes in part tissues. These evidences refined the range of candidate genes for carcass traits, including IGF2BP1, ACSL5, BMP3, etc., which mainly enriched in muscle growth, metabolism, and immune function. Totally, the findings of this study provide a list of crucial genes (e.g., IGF2BP1, ACSL5, BMP3) related to carcass traits (e.g., SW, EW, ThW) based genomic and transcriptomic evidence. And all these genes could serve as the molecular phenotypes for different carcass traits, and contribute to the development of genomic selection algorithms based on multi-omics technologies [40], which could provide valuable insights into improving carcass traits in chickens.
In SH chickens, the heritability of most carcass traits (except for SP) ranged from moderate to high, and a high genetic correlation was found among carcass traits (except for EW and EP). Based on GWAS summary, a total of 311 SNPs and 73 candidate genes were identified significant associations with six carcass traits (SW, EW, ThW, and relative weight of these traits). And 17 genes (e.g., IGF2BP1, RASGEF1B, ACSL5) were defined as DEGs between commercial broilers and SH chickens. Besides muscle growth and bone formation, the genes involved in lipid metabolism (e.g., ACSL5) and immune function (e.g., RASGEF1B) were also altered with carcass traits.
PDF Links
PubReader
ePub Link
Full text via DOI
Full text via PMC
Download Citation
Supplement1
Print
