A Missense Mutation in the MYBPH Gene Is Associated With Abdominal Fat Traits in Meat-Type Chickens

Chicken is an important source of protein for human nutrition and a model system for growth and developmental biology. Although the genetic architecture of quantitative traits in meat-type chickens has been the subject of ongoing investigation, the identification of mutations associated with carcass traits of economic interest remains challenging. Therefore, our aim was to identify predicted deleterious mutation, which potentially affects protein function, and test if they were associated with carcass traits in chickens. For that, we performed a genome-wide association analysis (GWAS) for breast, thigh and drumstick traits in meat-type chickens and detected 19 unique quantitative trait loci (QTL). We then used: (1) the identified windows; (2) QTL for abdominal fat detected in a previous study with the same population and (3) previously obtained whole genome sequence data, to identify 18 predicted deleterious single nucleotide polymorphisms (SNPs) in those QTL for further association with breast, thigh, drumstick and abdominal fat traits. Using the additive model, a predicted deleterious SNP c.482C > T (SIFT score of 0.4) was associated (p-value < 0.05) with abdominal fat weight and percentage. This SNP is in the second exon of the MYBPH gene, and its allele frequency deviates from Hardy–Weinberg equilibrium. In conclusion, our study provides evidence that the c.482C > T SNP in the MYBPH gene is a putative causal mutation for fat deposition in meat-type chickens.

Introduction

Chicken is an important source of protein for human nutrition and a model system for growth and developmental biology (Ellegren, 2005). The genome sequence of the Red Jungle Fowl (Gallus gallus gallus), considered the ancestor of the domestic chicken (G. g. domesticus) (Abplanalp, 1992; Dodgson et al., 2011), and completed in 2004 (Hillier et al., 2004) has allowed the development of new tools for genetc studies. High throughput sequencing of several breeding lines has identified millions of single nucleotide polymorphisms (SNPs) across the chicken genome (Rubin et al., 2010; Boschiero et al., 2018) and these led to the development of high-density SNP panels (Kranis et al., 2013). The most common and frequent DNA variants are SNPs, with a density of 13 SNPs/kb in a Brazilian meat-type chickenline (developed by Embrapa Swine and Poultry) (Boschiero et al., 2018).

High-density SNP panels were previously used in genome wide association studies (GWAS) to identify QTL for body weight (Gu et al., 2011), several measures of fatness (Sun et al., 2013; Moreira et al., 2018b), breast and leg muscle weights, wing weight (Xie et al., 2012), carcass and eviscerated weights (Liu et al., 2013). However, discovering the causative mutations underlying quantitative trait loci (QTL) remains challenging (Ahsan et al., 2013; Moreira et al., 2018a). Linkage disequilibrium between the genetic variant present on a SNP panel and the casual mutation allows QTL detection by genome-wide association study (GWAS), but fine-mapping studies are necessary to identify which sequence mutation within the QTL is the causative mutation responsible for the phenotype of interest. Combining statistical evidence from association studies with functional annotations of the genes or genetic variants, is a helpful approach to identify potential causal mutations (Spain and Barrett, 2015).

Single nucleotide polymorphisms can have a direct impact on coding or indirect impact on regulation of gene expression and thus affect traits of economic interest in animal models and livestock species (Roux et al., 2014). SNPs that occur in coding regions can be classified as missense when triplets code for different amino acids, synonymous when triplets code the same amino acid or nonsense when the mutation results in a premature stop codon. Missense SNPs can be predicted as deleterious or tolerated using the SIFT tool (Sorting Intolerant From Tolerant) (Ng and Henikoff, 2003). Changes at well-conserved positions tend to be deleterious because important amino acids are conserved across a protein family (Ng and Henikoff, 2002, 2003). When a SNP is predicted as a deleterious mutation, it means that a change in the amino acid sequence likely affects the protein structure and function (Ng and Henikoff, 2003; Kumar et al., 2009), and consequently, may alter one or more phenotypes. Previous studies have identified missense SNPs associated with body weight at hatch, semi-eviscerated carcass weight, eviscerated carcass weight, leg muscle weight and carcass weight (Wang et al., 2015), abdominal fat weight, body weight at different ages and body size traits (Han et al., 2012). Phenotype may also be affected by SNPs located on potentially neutral regions (Berulava and Horsthemke, 2010; Jo and Choi, 2015), as these may be regulatory.

Previous whole-genome resequencing studies performed in parental individuals from the meat-type TT chicken reference population studied herein identified several predicted deleterious SNPs (Moreira et al., 2018b). However, the potential role of predicted deleterious SNPs in the regulation of traits of economic interest is still unknown. In this study, we used a GWAS to identify regions associated with carcass traits in a meat-type population, then integrated whole genome sequence data to refine the list of candidate mutations, followed by targetted re-sequencing to identify predicted deleterious SNPs and association study to identify putative causal mutations.

Materials and Methods

Experimental Population

The TT Reference population used in this study was generated from the Brazilian TT broiler line, developed by the Embrapa Swine and Poultry National Research Center. This line has been under multiple-trait index selection since 1992, representing many generations, with the goals of increasing body weight and carcass yield, improving viability, fertility, hatchability, feed conversion, and reducing abdominal fat (Do Rosário et al., 2009). The TT Reference Population was developed from expansion of the TT selection line comprising crossing of 20 males with 92 females in five hatches, yielding 1,430 chickens (Da Cruz et al., 2015; Marchesi et al., 2017).

Phenotype Measurement

Body weight at 42 days of age (BW42) was measured six hours after fasting, then chickens were euthanized by cervical dislocation followed by bleeding. Blood samples were collected for DNA extraction during bleeding. Feathers were mechanically removed following a hot water bath (60°C for 45 s) and carcass cuts representing breast weight (BTW), thigh weight (THW), drumstick weight (DRW) and abdominal fat weight (ABFW) were individually measured in grams. Drumstick yield (DR%), abdominal fat yield (ABF%), thigh yield (TH%) and breast weight yield (BT%) were estimated as a percentage of live body weight at 42 days of age. More details about the slaughter and phenotype measurements have been previously described (Venturini et al., 2014; Da Cruz et al., 2015).

DNA Extraction and Genotyping

Genomic DNA was extracted from 1,430 blood samples with the PureLink^® Genomic DNA kit (Invitrogen, Carlsbad, CA, United States) and quantified using Qubit^® 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States). For genotyping analyses, we used the 600 K Affymetrix Axiom^TM Genotyping Array (Affymetrix, Inc., Santa Clara, CA, United States). Sample filtering parameters were DishQC ≥ 0.82 performed with Axiom^TM Analysis Suite (Affymetrix^®) software and sample call rate ≥ 90% via PLINK v.1.9 software (Purcell et al., 2007). For loci, the filter parameters of call rate ≥ 98% and minor allele frequency (MAF) ≥ 2% were used. We also excluded SNPs with significant deviations from Hardy–Weinberg equilibrium (HWE) (p-value < 0.000001), those located in the sex chromosomes, and those not annotated in the chicken assembly (Gallus_gallus-5.0, NCBI). More details about the genotyping and filtering steps are in Moreira et al. (2018b). Genomic analysis workflow is presented in Figure 1.

FIGURE 1

Figure 1. Workflow of genomic analysis performed.

Genome Wide Association Analysis

We used GenSel software for GWAS based on a Bayesian methodology approach for genomic prediction. Following Cesar et al. (2014) and Moreira et al. (2018b), a Bayes C model was used to estimate the genetic and residual variances for each trait that were used as priors in fitting a Bayes B model. The mathematical model used was:

$$y=Xb+\sum _{\mathrm{j}=1}^k{a}_{\text{j}}{\mathrm{\beta }}_{\text{j}}{\mathrm{\delta }}_{\text{j}}+e.$$

As described by Moreira et al. (2018b), y represents the vector of phenotypic values (BTW, BT%, THW, TH%, DRW and DR%); X is the incidence matrix for fixed effects; b is the vector of fixed effects; k is the number of SNP loci; a_j is the column vector representing the SNP at locus j as a covariate, coded with the number of B alleles; β_j is the random substitution effect for locus j assumed to be normally distributed N (0, σ²_βj) when δj = 1 but 0 when δj = 0, with δj being a random variable 0/1 indicating the absence (with probability π) or presence (with probability 1-π) of locus j in the model, and e is the residual associated with the analysis. Sex and hatch were included as fixed effects in the model and BW42 (slaughter age) as a fixed covariate for THW, BTW, ABFW and DRW.

We assumed π = 0.9970 in the BayesB models and obtained 41,000 Markov chain Monte Carlo (MCMC) samples with the first 1,000 samples being discarded. A map file was used to position the consecutive markers into 947 non-overlapping 1 Mb windows. So, for each window, the proportion of the genetic variance explained by the QTL was computed among individuals for every 100^th iteration of the MCMC chain based on the marker effects sampled in that iteration as performed by Schurink et al. (2012). The windows that had the marker with higher model frequency in the MCMC iterations had their effect predicted as mentioned by Van Goor et al. (2015). Each window is expected to explain 0.1054% of the genetic variance (100%/947) based on an infinitesimal model (Onteru et al., 2013; Van Goor et al., 2016). Were further considered as significant windows those that explained five times more variation than expected (0.53%) in an infinitesimal model, as used by Moreira et al. (2018b).

SNP Selection and Custom Amplicon Design

Among the 20 founding males of the TT Reference population, 14 were resequenced by our group, resulting in approximately 13X sequencing coverage using a HiScanSQ (Illumina) sequencer and that produced a dataset of good quality SNPs we deposited in the EVA-EMBL database¹. Further details about library preparation, sequencing and filtering are already published (Boschiero et al., 2018; Moreira et al., 2018b). Functional annotation of the SNPs was performed using VEP [Variant Effect Predictor, v.86, (McLaren et al., 2016)] and deleterious predictions were based on SIFT scores (Ng and Henikoff, 2003). To investigate the segregation of the predicted deleterious SNPs detected in the 14 founder males, we selected 237 descendents for targeted resequencing. These offspring represent 14 half-sib and 37 full-sib families. In each of the full-sib families, five to seven animals were re-sequenced. For the genotyping by sequencing methodology, a region of 150 bp centered around each predicted deleterious SNP was used to define a target region for amplicon design (DesignStudio online platform from Illumina).

This study investigated only predicted deleterious SNPs located within significant QTL regions associated with the carcass traits and abdominal fat traits (Moreira et al., 2018b), to narrow the causative mutation search. To avoid overlaps between the amplicons and maximize the number of regions sequenced, we used the Tagger tool in Haploview software (Barrett et al., 2005) to select one SNP per haplotype block of interest. Thus, if adjacent predicted deleterious SNPs exhibited r² > 0.7, just one potential causative SNP was chosen to be genotyped by sequencing in the offspring generation.

Target Sequencing

Genomic DNA of 237 offspring was extracted using the PureLink^® Genomic DNA kit (Invitrogen, Carlsbad, CA, United States) and then quantified using a Qubit^® 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States). DNA integrity was evaluated in 1% agarose gel. Library preparation was performed according to Truseq^® Custom Amplicon Low Input Kit Reference Guide (Illumina Technology). Libraries were quantified with quantitative real time PCR, using KAPA^® Library Quantification kit (KAPA Biosystem) and fragment size was estimated using either Bioanalyzer^® (Agilent Technologies) or a Fragment Analyzer (Advanced Analytical Technologies). Paired-end sequencing with a read length of 150 bp was performed on a MiniSeq^TM (Illumina Technology).

Data Analyses, Variant Calling and Functional Annotation of Sequencing Data

We aligned the raw sequence reads against the chicken reference genome Gallus_gallus5.0 (NCBI) with BWA-MEM (v.0.7.15). For SNP calling, we used SAMtools v.1.3.1 (Li et al., 2009) with mpileup option (Li, 2011), and filtered based on mapping and base qualities (Phred) ≥ 20. The variant calling was performed after sequence reads from all 237 animals were pooled together. After the initial variant identification, we applied the following filtering options: INDEL removal, minor allele frequency (MAF) ≥ 0.05, SNP call rate ≥ 0.7, biallelic locus, sequencing depth ≥ 15 and Phred score quality ≥ 40.

After filtering, we annotated those SNPs that remained using the VEP tool version 91 (McLaren et al., 2016) available from the Ensembl v.91 website (Zerbino et al., 2018), where the SIFT score was predicted. We searched for gene ontology terms in Amigo 2 GO online server². Deviations from Hardy–Weinberg equilibrium were calculated in Haploview software (Barrett et al., 2005) and Bonferroni multiple test correction was applied.

Association Analysis

For association analysis, we tested 18 predicted deleterious SNPs located in QTL regions using the package lmerTest 3.1-0 (Kuznetsova et al., 2017) in R software (v. 3.5.2)³. To investigate the SNP effects, we fitted both additive and non-additive models for each trait (BTW, BT%, THW, TH%, DRW, DR%, ABFW, and AB%). As adopted by Edwards et al. (2015), for the additive models the SNPs were considered as continuous variables (0, 1, or 2 copies of a given SNP) and in the non-additive models the variable genotypes were considered as a factor (with three levels i.e. aa, Aa, AA). Sex and hatch were added in the model as fixed effects, and family (determined by the mother) was fitted as a random effect. For the carcass weight traits (BTW, THW, DRW and ABFW), BW42 was used as a fixed covariate. For each trait, all the SNPs located in the QTL associated with the respective trait were simultaneously fitted in the following model:

$$y=X\mathrm{\beta }+Wa+Zu+e$$

Where y is the vector of observations for the measured phenotype; X is the incidence matrix relating the fixed effects for sex and hatch to y; β is the vector of sex and hatch fixed effects; W is the genotype matrix for the predicted deleterious SNPs located in the associated QTL regions; and a is the vector of fixed effects for the SNPs. The matrix Z is an incidence matrix relating u to y; u is the vector of length 41 representing the random effect of dam family; and e is the vector of residual effects. Associations were considered significant at comparison-wise p-value < 0.05.

Results

Phenotype Measures and Descriptive Statistics

Descriptive statistics such as number of animals, averages and standard errors from traits and animals utilized in GWAS are presented in Table 1. Summary statistics for all traits and animals used in genotyping and single SNP association analysis are presented in Table 2.

TABLE 1

Table 1. Number of animals (N), mean, standard deviation (SD), minimum and maximum values for carcass traits from animals utilized for GWAS analysis in the TT Reference population.

TABLE 2

Table 2. Number of animals (N), mean, standard deviation (SD), minimum and maximum values for carcass traits from animals utilized for single SNP association analysis in the TT Reference population.

Genome-Wide Association Analysis

We conducted a GWAS to identify genomic regions associated with traits of interest. Nineteen unique QTLs for THW, TH%, DRW, DR%, BTW, and BT%, were detected on GGA 1-4, 6, 12, 14, 17, 20, 22, 25, 26 and 28 (Table 3). The QTL explained from 0.54 to 3.2% of the genetic variance. The posterior probability of association (PPA) (Onteru et al., 2013) ranged from 0.69 to 0.99. Characterization of all 947 genomic windows analyzed is available in Supplementary Table 1. For ABFW and ABF%, nine unique QTL had been previously reported (see Moreira et al., 2018b).

TABLE 3

Table 3. Significant genomic windows associated with THW, TH%, DRW, DR%, BTW, and BT% in the TT Reference population.

Windows that explained a high proportion of genetic variance were selected for each carcass trait to further search for predicted deleterious SNPs. The genomic windows selected were: GGA 1 (166 Mb) for drumstick traits; GGA 22 (4 Mb) for thigh traits; GGA 26 (1 Mb) for abdominal fat traits [9]; GGA 25 (2 Mb) and GGA 26 (3 Mb) for breast traits.

SNP Selection and Amplicon Design

We searched for predicted deleterious SNPs located within the QTLs to identify potential causative mutations. We scrutinized SNP data from whole-genome resequencing of 14 founding animals (TT line), approximately 11 million SNPs (Boschiero et al., 2018), concordant with the five selected GWAS windows. This recovered 89 predicted deleterious SNPs. After linkage disequilibrium pruning in the founder’s dataset, 18 uncorrelated predicted deleterious SNPs were selected for amplicon design and association analysis. One predicted deleterious SNP was tested for each of the following regions: GGA 1 (166 Mb) window associated with DRW and DR%, GGA 22 (4 Mb) associated with THW and TH% and GGA 25 (2 Mb) associated with BRW and BR%. For GGA 26 (3 Mb) associated with BRW and BR% and GGA 26 (1 Mb) associated with ABFW and ABF%, seven and eight predicted deleterious SNPs that were not in linkage disequilibrium were tested, respectively.

Amplicon Sequencing, Variant Calling and Functional Annotation

Libraries sequenced using MiniSeq produced on average 298,791 raw reads per amplicon. The average overall mapping rate of the raw reads against the Gallus_gallus5.0 (NCBI) genome assembly was 99.74%. After variant calling, quality control and functional annotation, the 18 predicted deleterious SNPs were genotyped by sequencing with an average coverage of 7,000X.

Detailed information of the 18 predicted deleterious SNPs (genome position, SNP ID, located gene, allele and genotype frequencies, HWE test and SIFT score) is in Table 4. The 18 predicted deleterious SNPs were detected in the offspring, validating our whole genome sequence data and SNP calling. Two SNPs are novel so did not have rs IDs, and four SNPs were annotated in novel genes. Seven SNPs did not have any animal genotyped as homozygous for the alternative allele, and five SNPs exhibited significant deviation from HWE.

TABLE 4

Table 4. Deleterious SNPs selected for the association analyses with carcass traits.

Association Analysis for Additive and Dominance Models

We used both additive and non-additive models to verify if any of the deleterious SNP identified in each of the candidate genes were associated with the phenotype detected in the GWAS analysis. Detailed results for association tests (p-values and effects) between the SNPs and the eight traits are presented in Table 5 for the test for additive effects and Table 6 for the test for dominance. It is important to keep in mind that, for each trait, the association test was performed only with the SNPs located in the QTL region chosen for the respective trait. Based on the additive model, from the 18 predicted deleterious SNPs tested, c.482C > T was the only locus significantly associated with ABFW and ABF% (Table 5). In the analysis testing for dominance using a non-additive model, there were no SNPs with significant effects for any of the eight traits (Table 6).

TABLE 5

Table 5. Association analysis results, p-values (upper value in cell) and effects (lower value in cell) for an additive model applied to carcass traits in the TT Reference population.

TABLE 6

Table 6. Association analysis results (p-values) for dominance effects based on a non-additive model applied to carcass traits.

The SNP c.482C > T is in the Myosin Binding Protein H (MYBPH) gene and within the QTL region for ABFW and ABF% traits in GGA26. This polymorphism is a C > T change, resulting in the amino acid change of threonine to methionine at amino acid position 161 (T161M). In the additive model, the effects of the tested predicted deleterious SNP c.482C > T were −0.17% and −3.43g for ABF% and ABFW traits, respectively. Functional annotation analysis performed in Amigo 2 GO of this gene reported gene ontology (GO) terms related to protein binding (GO:0005515), regulation of striated muscle contraction (GO: 0006942), cell adhesion (GO:0007155), structural constituent of muscle (GO:0008307) and myosin filament (GO:0032982).

Discussion

The economic importance of carcass weight and yield motivated the search for genomic regions and candidate genes associated with these traits (Coutinho et al., 2010; Boschiero et al., 2013; Van Goor et al., 2015; Derks et al., 2018). Combining functional studies and association analyses helped us to focus on positionally supported polymorphisms with the main goal of finding putative causal mutations for carcass traits. Here we present a report with genomic association analysis combining information of genomic regions and whole genome sequence to support the selection of predicted deleterious SNPs for association analysis with carcass traits in chickens.

From GWAS analysis, 19 unique QTL for THW, TH%, DRW, DR%, BTW and BT%, were detected, cumulatively explaining 0.54, 1.19, 8.54, 8.11, 2.64, and 4.97% of the genetic variance, respectively, with PPA > 0.69. These QTL are useful for further applications such as positional candidate genes search, gene network analysis and discovery of putative candidate mutations, as we investigated herein.

From the 19 QTL identified, five genomic windows were used here (four identified in this study and one characterized by Moreira et al., 2018b). Whitin these regions, we selected SNPs that had SIFT scores ranging from 0.0 to 0.04. Lower scores represents higher confidence and better sensitivity for detecting deleterious SNPs (Ng and Henikoff, 2003). All SNPs had an amino acid change prediction considered as moderate impact; no substitution resulted in a premature stop codon.

Deviation from HWE is an indication that allele frequency is not consistent with genotype frequency. Several evolutionary forces can affect HWE and a deleterious SNP could drive a deviation from HWE. Four predicted deleterious polymorphisms exhibited significant departures from HWE (p-value < 0.05 after bonferroni correction for multiple testing). Deviations from HWE may suggest that selection has affected the allele frequency.

It is pertinent to note that even though rs738655377 was not significantly associated with BTW or BT%, none of the animals genotyped were homozygous for the alternative allele, and the test for departure from HWE was significant. This variant is within a novel gene (ENSGALG00000028858) that presents gene ontology terms related to oxidoreductase activity. Out of the 237 animals genotyped, the expectation based on the observed allele frequency was that 4.45 animals would be homozygous for the alternative allele. The lack of homozygous animals for the alternative allele provides suggestive evidence of a lethal polymorphism when in homozygosity. However, further genotyping studies should be performed to confirm these findings based on matings between heterozygotes. Previous studies in cattle and pigs have identified lethal mutation by searching for lack of homozygous alleles (VanRaden et al., 2011; Kuehn et al., 2013; Derks et al., 2017; Jenko et al., 2019).

The significantly associated polymorphism, c.482C > T, located in the MYBPH gene is a novel SNP that caused a threonine-to-methionine substitution at amino acid position 161 (T161M). Threonine amino acid is polar (uncharged), while methionine is a non-polar amino acid. In Myosin-binding protein H, the substitution is located in the Fibronectin type-II domain (UniProt accession number Q05623 and PROSITE accession number PS50853) and this domain is the most common of the fibronectin subdomains (Hynes, 1990). Fibronectin is a dimeric glycoprotein involved in cell adhesion, cell morphology, cell migration, and embryonic differentiation.

The MYBPH presents gene ontology terms related to muscle growth and development such as: protein binding, regulation of striated muscle contraction, cell adhesion, structural constituent of muscle and myosin filament. In chickens, an increase in the amount of MYBPH protein in breast muscle is associated with muscle development at 56 days of age (Liu et al., 2016). In cattle, the MYBPH gene was previously shown to be associated with intramuscular fat deposition (Poleti et al., 2018; Bazile et al., 2019). A proteomic study performed in Longissimus dorsi samples in cattle found that MYBPH protein was significantly more abundant (p-value < 0.05) in a group of cattle exhibiting low intramuscular fat content (Poleti et al., 2018). Another bovine proteomic study found a significantly lower abundance of MYBPH protein in the high adiposity group (Bazile et al., 2019). Other studies with Nellore cattle found MYBPH gene differentially expressed when comparing groups with high and low intramuscular fat content (Cesar et al., 2015; dos Santos Silva et al., 2019). In chickens, the MYBPH was expressed in skeletal muscle tissue, heart, spleen, brain, kidney and female gonad (UniProt Databse, Bateman et al., 2021). Based on the association results in our population, the allele substitution effect of the predicted deleterious SNP was a reduction in 3.43 g of ABWF and 0.17% of ABF% for each copy of the alternative allele in the genotype.

It is important to mention that the polymorphism c.482C > T herein might not be the causal mutation because it was in strong linkage disequilibrium (LD block GGA26 1008728 Mb – 1015596 Mb with 0.92 r²) with two other polymorphisms in the MYBPH gene, rs313948592 (A/C, located in the intron) and rs312491203 (G/A change, located upstream of the gene). We are aware that polymorphism in intron or upstream regions can have regulatory impact and we cannot rule predicted deleterious SNP as causative mutations (Jo and Choi, 2015; Li et al., 2020). Nevertheless, this study indicates that MYBPH gene could have an impact on fat deposition in chickens. Further functional studies may help to elucidate the role of MYBPH and its polymorphisms in regulating fatness in meat-type chickens. Also, the other polymorphisms identified in the target regions are available at European Variation Archive (EVA) – EMBL-EBI (accession PRJEB25004) for further studies of causal mutations in these regions.

Our strategy has provided novel evidence in the identification of a putative causative mutation associated with abdominal fatness in meat-type chickens. We report one predicted deleterious SNP associated with abdominal fat traits located in the MYBPH gene that is a candidate gene for fat deposition. The main limitation of our study is not being able to discriminate whether the identified mutation is the causative mutation, or it is simply in strong linkage disequilibrium with the real causal mutation. Nevertheless, our findings can be applied in poultry breeding programs focused on carcass traits improvement.

Data Availability Statement

All data generated or analyzed during this study are public and included in this published article. The SNP report (identified by sequencing) was submitted to European Variation Archive (EVA) – EMBL-EBI, accession number PRJEB25004. The datasets used and/or analyzed during the current study (genotypes and phenotypes) are available from the corresponding author on reasonable request.

Ethics Statement

In this study, all experimental protocols that used animals were performed in agreement with resolution number 010/2010 approved by the Embrapa Swine and Poultry Ethics Committee on Animal Utilization (CEUA) in Concordia, Santa Catarina State – South of Brazil, in agreement with the rules of National Council of Animal Experimentation Control (CONCEA) to ensure compliance with international guidelines for animal welfare.

Author Contributions

PT, GCM, CB, ML, GBM, and LC conceived the idea of this research and participated in the experimental design. ML and LC provided the experimental environment, phenotype, and data analysis support. PT, GCM, CB, JP, AC, GRM, and DG performed data analysis. PT drafted the manuscript. PT, GCM, CB, AC, JP, GRM, ML, GBM, DG, and LC collaborated with interpretation, discussion, and writing of the manuscript. All authors have read and approved the final manuscript.

Funding

The authors declare that this study received funding from São Paulo Research Foundation (FAPESP) thematic grant 14/08704-0 and Brazilian Agricultural Research Corporation – Embrapa (project number 01.11.07.002.04.02), both founders contributed to experimental environment, phenotype measure, sequencing, genotyping. The TT Reference Population was subsidized by the National Council of Scientific and Technological Development (CNPq) grant number 481755/2007-1 from the Brazilian Government. PT received a fellowship from FAPESP, grant 2016/13589-0. GCM received fellowships from FAPESP, grants 14/21380–9 (in cooperation agreement with CAPES), and 16/00569–1. CB received a fellowship from the program Science without Borders from CNPq grant 370620/2013–5. ML is employee of the Brazilian Agricultural Research Corporation – Embrapa. L.L.C. GRM, ML, and GBM are recipient of productivity fellowship from CNPq.

Conflict of Interest

FAPESP, CAPES, CNPq, and the Brazilian Government through Embrapa provided financial support to generate the data; however, they did not participate in the design of the study, sample collection, analysis, data interpretation, and writing of the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as apotential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors are grateful to Ricardo Brassaloti (University of São Paulo) for the help in library preparation and sequencing. The authors are also grateful for the support from the FAPESP, CAPES, and CNPq mentioned in funding statement. The authors would like to acknowledge the collaborative efforts among EMBRAPA, University of São Paulo.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.698163/full#supplementary-material

Supplementary Table 1 | Characterization of all the 947 SNP windows analyzed.

Footnotes

1. ^ https://www.ebi.ac.uk/ena/data/view/PRJEB25004
2. ^ http://amigo.geneontology.org
3. ^ http://www.r-project.org/

References

Abplanalp, H. (1992). Inbred lines as genetic resources of chickens. Poult. Sci. Rev. 4, 29–39.

Google Scholar

Ahsan, M., Li, X., Lundberg, A. E., Kierczak, M., Siegel, P. B., Carlborg, O., et al. (2013). Identification of candidate genes and mutations in QTL regions for chicken growth using bioinformatic analysis of NGS and SNP-chip data. Front. Genet. 4:226. doi: 10.3389/fgene.2013.00226

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrett, J. C., Fry, B., Maller, J., and Daly, M. J. (2005). Haploview: analysis and visualization of LD and haplotype maps. Bioinform. Appl. 21, 263–265. doi: 10.1093/bioinformatics/bth457

PubMed Abstract | CrossRef Full Text | Google Scholar

Bateman, A., Martin, M. J., Orchard, S., Magrane, M., Agivetova, R., Ahmad, S., et al. (2021). UniProt: The universal protein knowledgebase in 2021. Nucleic Acids Res. 49, D480–D489. doi: 10.1093/nar/gkaa1100

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazile, J., Picard, B., Chambon, C., Valais, A., and Bonnet, M. (2019). Pathways and biomarkers of marbling and carcass fat deposition in bovine revealed by a combination of gel-based and gel-free proteomic analyses. Meat Sci. 156, 146–155. doi: 10.1016/J.MEATSCI.2019.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Berulava, T., and Horsthemke, B. (2010). The obesity-associated SNPs in intron 1 of the FTO gene affect primary transcript levels. Eur. J. Hum. Genet. 18, 1054–1056. doi: 10.1038/ejhg.2010.71

PubMed Abstract | CrossRef Full Text | Google Scholar

Boschiero, C., Jorge, E. C., Ninov, K., Nones, K., Do Rosário, M. F., Coutinho, L. L., et al. (2013). Association of IGF1 and KDM5A polymorphisms with performance, fatness and carcass traits in chickens. J. Appl. Genet. 54, 103–112. doi: 10.1007/s13353-012-0129-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Boschiero, C., Moreira, G. C. M., Gheyas, A. A., Godoy, T. F., Gasparin, G., Mariani, P. D. S. C., et al. (2018). Genome-wide characterization of genetic variants and putative regions under selection in meat and egg-type chicken lines. BMC Genomics 19:83. doi: 10.1186/s12864-018-4444-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cesar, A. S. M., Regitano, L. C. A., Koltes, J. E., Fritz-Waters, E. R., Lanna, D. P. D., Gasparin, G., et al. (2015). Putative regulatory factors associated with intramuscular fat content. PLoS One 10:e0128350. doi: 10.1371/journal.pone.0128350

PubMed Abstract | CrossRef Full Text | Google Scholar

Cesar, A. S. M., Regitano, L. C. A., Mourão, G. B., Tullio, R. R., Lanna, D. P. D., Nassu, R. T., et al. (2014). Genome-wide association study for intramuscular fat deposition and composition in Nellore cattle. BMC Genet. 15:39. doi: 10.1186/1471-2156-15-39

PubMed Abstract | CrossRef Full Text | Google Scholar

Coutinho, L. L., Rosário, M. F. D., and Jorge, E. C. (2010). Biotecnologia animal. Butanta: Instituto de Estudos Avançados da Universidade de São Paulo. doi: 10.1590/S0103-40142010000300009

CrossRef Full Text | Google Scholar

Da Cruz, V. A. R., Schenkel, F. S., Savegnago, R. P., Grupioni, N. V., Stafuzza, N. B., Sargolzaei, M., et al. (2015). Association of apolipoprotein b and adiponectin receptor 1 genes with carcass, bone integrity and performance traits in a paternal broiler line. PLoS One 10:e0136824. doi: 10.1371/journal.pone.0136824

PubMed Abstract | CrossRef Full Text | Google Scholar

Derks, M. F. L., Megens, H. J., Bosse, M., Lopes, M. S., Harlizius, B., and Groenen, M. A. M. (2017). A systematic survey to identify lethal recessive variation in highly managed pig populations. BMC Genomics 18:858. doi: 10.1186/s12864-017-4278-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Derks, M. F. L., Megens, H.-J., Bosse, M., Visscher, J., Peeters, K., Bink, M. C. A. M., et al. (2018). A survey of functional genomic variation in domesticated chickens. Genet. Sel. Evol. 50:17. doi: 10.1186/s12711-018-0390-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Do Rosário, M. F., Ledur, M. C., Moura, A. S. A. M. T., Coutinho, L. L., and Garcia, A. A. F. (2009). Genotypic characterization of microsatellite markers in broiler and layer selected chicken lines and their reciprocal F1s. Sci. Agric. 66, 150–158. doi: 10.1590/S0103-90162009000200002

CrossRef Full Text | Google Scholar

Dodgson, J. B., Delany, M. E., and Cheng, H. H. (2011). Poultry genome sequences: progress and outstanding challenges. Cytogenet. Genome Res. 134, 19–26. doi: 10.1159/000324413

PubMed Abstract | CrossRef Full Text | Google Scholar

dos Santos Silva, D. B., Fonseca, L. F. S., Pinheiro, D. G., Muniz, M. M. M., Magalhães, A. F. B., Baldi, F., et al. (2019). Prediction of hub genes associated with intramuscular fat content in Nelore cattle. BMC Genomics 20:520. doi: 10.1186/s12864-019-5904-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Edwards, H. A., Hajduk, G. K., Durieux, G., Burke, T., and Dugdale, H. L. (2015). No association between personality and candidate gene polymorphisms in a wild bird population. PLoS One 10:e0138439. doi: 10.1371/journal.pone.0138439

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellegren, H. (2005). The avian genome uncovered. Trends Ecol. Evol. 20, 180–186. doi: 10.1016/j.tree.2005.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, X., Feng, C., Ma, L., Song, C., Wang, Y., Da, Y., et al. (2011). Genome-wide association study of body weight in chicken F2 resource population. PLoS One 6:e21872. doi: 10.1371/journal.pone.0021872

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, R., Wei, Y., Kang, X., Chen, H., Sun, G., Li, G., et al. (2012). Novel SNPs in the PRDM16 gene and their associations with performance traits in chickens. Mol. Biol. Rep. 39, 3153–3160. doi: 10.1007/s11033-011-1081-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillier, L. W., Miller, W., Birney, E., Warren, W., Hardison, R. C., Ponting, C. P., et al. (2004). Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution. Nature 432, 695–716. doi: 10.1038/nature03154

PubMed Abstract | CrossRef Full Text | Google Scholar

Hynes, R. O. (1990). Introduction and Historical Overview. Fibronectins 1990, 1–6. doi: 10.1007/978-1-4612-3264-3_1

CrossRef Full Text | Google Scholar

Jenko, J., McClure, M. C., Matthews, D., McClure, J., Johnsson, M., Gorjanc, G., et al. (2019). Analysis of a large dataset reveals haplotypes carrying putatively recessive lethal and semi-lethal alleles with pleiotropic effects on economically important traits in beef cattle. Genet. Sel. Evol. 51:9. doi: 10.1186/s12711-019-0452-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Jo, B.-S., and Choi, S. S. (2015). Introns: the functional benefits of introns in genomes. Genomics Inform. 13, 112–118. doi: 10.5808/GI.2015.13.4.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Kranis, A., Gheyas, A. A., Boschiero, C., Turner, F., Yu, L., Smith, S., et al. (2013). Development of a high density 600K SNP genotyping array for chicken. BMC Genomics 14:59. doi: 10.1186/1471-2164-14-59

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuehn, J. S., Gorden, P. J., Munro, D., Rong, R., Dong, Q., Plummer, P. J., et al. (2013). Bacterial community profiling of milk samples as a means to understand culture-negative bovine clinical mastitis. PLoS One 8:e61959. doi: 10.1371/journal.pone.0061959

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, P., Henikoff, S., and Ng, P. C. (2009). Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat. Protoc. 4, 1073–1081. doi: 10.1038/nprot.2009.86

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuznetsova, A., Brockhoff, P. B., and Christensen, R. H. B. (2017). lmertest package: tests in linear mixed effects models. J. Stat. Softw. 82, 1–26. doi: 10.18637/jss.v082.i13

CrossRef Full Text | Google Scholar

Li, H. (2011). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27, 2987–2993. doi: 10.1093/bioinformatics/btr509

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and samtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Lee, M. O., Davis, B. W., Lamichhaney, S., Dorshorst, B. J., Siegel, P. B., et al. (2020). Mutations upstream of the TBX5 and PITX1 transcription factor genes are associated with feathered legs in the domestic chicken. Mol. Biol. Evol. 37, 2477–2486. doi: 10.1093/molbev/msaa093

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Fu, R., Liu, R., Zhao, G., Zheng, M., Cui, H., et al. (2016). Protein profiles for muscle development and intramuscular fat accumulation at different post-hatching ages in chickens. PLoS One 11:e0159722. doi: 10.1371/journal.pone.0159722

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, R., Sun, Y., Zhao, G., Wang, F., Wu, D., Zheng, M., et al. (2013). Genome-wide association study identifies loci and candidate genes for body composition and meat quality traits in beijing-you chickens. PLoS One 8:e61172. doi: 10.1371/journal.pone.0061172

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchesi, J. A. P., Buzanskas, M. E., Cantão, M. E., Ibelli, A. M. G., Peixoto, J. O., Joaquim, L. B., et al. (2017). Relationship of runs of homozygosity with adaptive and production traits in a paternal broiler line. Animal 12, 1126–1134. doi: 10.1017/S1751731117002671

PubMed Abstract | CrossRef Full Text | Google Scholar

McLaren, W., Gil, L., Hunt, S. E., Riat, H. S., Ritchie, G. R. S., Thormann, A., et al. (2016). The ensembl variant effect predictor. Genome Biol. 17:122. doi: 10.1186/s13059-016-0974-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreira, G. C. M., Boschiero, C., Cesar, A. S. M., Reecy, J. M., Godoy, T. F., Pértille, F., et al. (2018a). Integration of genome wide association studies and whole genome sequencing provides novel insights into fat deposition in chicken. Sci. Rep. 8:16222. doi: 10.1038/s41598-018-34364-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreira, G. C. M., Boschiero, C., Cesar, A. S. M., Reecy, J. M., Godoy, T. F., Trevisoli, P. A., et al. (2018b). A genome-wide association study reveals novel genomic regions and positional candidate genes for fat deposition in broiler chickens. BMC Genomics 19:374. doi: 10.1186/s12864-018-4779-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, P. C., and Henikoff, S. (2002). Accounting for human polymorphisms predicted to affect protein function. Genome Res. 12, 436–446. doi: 10.1101/gr.212802

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, P. C., and Henikoff, S. (2003). SIFT: Predicting amino acid changes that affect protein function. Nucleic Acids Res. 31, 3812–3814. doi: 10.1093/nar/gkg509

PubMed Abstract | CrossRef Full Text | Google Scholar

Onteru, S. K., Gorbach, D. M., Young, J. M., Garrick, D. J., Dekkers, J. C. M., and Rothschild, M. F. (2013). Whole genome association studies of residual feed intake and related traits in the pig. PLoS One 8:e61756. doi: 10.1371/journal.pone.0061756

PubMed Abstract | CrossRef Full Text | Google Scholar

Poleti, M. D., Regitano, L. C. A., Souza, G. H. M. F., Cesar, A. S. M., Simas, R. C., Silva-Vignato, B., et al. (2018). Longissimus dorsi muscle label-free quantitative proteomic reveals biological mechanisms associated with intramuscular fat deposition. J. Proteomics 179, 30–41. doi: 10.1016/J.JPROT.2018.02.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., et al. (2007). PLINK: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/519795

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, P.-F., Boutin, M., Désert, C., Djari, A., Esquerré, D., Klopp, C., et al. (2014). Re-sequencing data for refining candidate genes and polymorphisms in QTL regions affecting adiposity in chicken. PLoS One 9:e111299. doi: 10.1371/journal.pone.0111299

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubin, C.-J., Zody, M. C., Eriksson, J., Meadows, J. R. S., Sherwood, E., Webster, M. T., et al. (2010). Whole-genome resequencing reveals loci under selection during chicken domestication. Nature 464, 587–591. doi: 10.1038/nature08832

PubMed Abstract | CrossRef Full Text | Google Scholar

Schurink, A., Wolc, A., Ducro, B. J., Frankena, K., Garrick, D. J., Dekkers, J. C. M., et al. (2012). Genome-wide association study of insect bite hypersensitivity in two horse populations in the Netherlands. Genet. Sel. Evol. 44, 1–12. doi: 10.1186/1297-9686-44-31

PubMed Abstract | CrossRef Full Text | Google Scholar

Spain, S. L., and Barrett, J. C. (2015). Strategies for fine-mapping complex traits. Hum. Mol. Genet. 24, R111–R119. doi: 10.1093/hmg/ddv260

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Zhao, G., Liu, R., Zheng, M., Hu, Y., Wu, D., et al. (2013). The identification of 14 new genes for meat quality traits in chicken using a genome-wide association study. BMC Genomics 14:458. doi: 10.1186/1471-2164-14-458

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Goor, A., Ashwell, C. M., Persia, M. E., Rothschild, M. F., Schmidt, C. J., and Lamont, S. J. (2016). Quantitative trait loci identified for blood chemistry components of an advanced intercross line of chickens under heat stress. BMC Genomics 17:287. doi: 10.1186/s12864-016-2601-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Goor, A., Bolek, K. J., Ashwell, C. M., Persia, M. E., Rothschild, M. F., Schmidt, C. J., et al. (2015). Identification of quantitative trait loci for body temperature, body weight, breast yield, and digestibility in an advanced intercross line of chickens under heat stress. Genet. Sel. Evol. 47:96. doi: 10.1186/s12711-015-0176-7

PubMed Abstract | CrossRef Full Text | Google Scholar

VanRaden, P. M., Olson, K. M., Null, D. J., and Hutchison, J. L. (2011). Harmful recessive effects on fertility detected by absence of homozygous haplotypes. J. Dairy Sci. 94, 6153–6161. doi: 10.3168/jds.2011-4624

PubMed Abstract | CrossRef Full Text | Google Scholar

Venturini, G. C., Cruz, V. A. R., Rosa, J. O., Baldi, F., Faro, L., Ledur, M. C., et al. (2014). Genetic and phenotypic parameters of carcass and organ traits of broiler chickens. Genet. Mol. Res. 13, 10294–10300. doi: 10.4238/2014.December.4.24

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y. C., Jiang, R. R., Kang, X. T., Li, Z. J., Han, R. L., Geng, J., et al. (2015). Identification of single nucleotide polymorphisms in the ASB15 gene and their associations with chicken growth and carcass traits. Genet. Mol. Res. 14, 11377–11388. doi: 10.4238/2015.September.25.5

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, L., Luo, C., Zhang, C., Zhang, R., Tang, J., Nie, Q., et al. (2012). Genome-wide association study identified a narrow chromosome 1 region associated with chicken growth traits. PLoS One 7:e30910. doi: 10.1371/journal.pone.0030910

PubMed Abstract | CrossRef Full Text | Google Scholar

Zerbino, D. R., Achuthan, P., Akanni, W., Amode, M. R., Barrell, D., Bhai, J., et al. (2018). Ensembl 2018. Nucleic Acids Res. 46, D754–D761. doi: 10.1093/nar/gkx1098

PubMed Abstract | CrossRef Full Text | Google Scholar