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).

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).

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.

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:

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).

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.

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).

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.

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:

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.

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.

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).

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.

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.

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.

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).

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).