Genome-wide association scan for QTL and their positional candidate genes associated with internal organ traits in chickens

Background Poultry breeding programs have been focused on improvement of growth and carcass traits, however, this has resulted in correlated changes in internal organ weights and increased incidence of metabolic disorders. These disorders can affect feed efficiency or even cause death. We used a high density SNP array (600 K, Affymetrix) to estimate genomic heritability, perform genome-wide association analysis, and identify genomic regions and positional candidate genes (PCGs) associated with internal organ traits in an F2 chicken population. We integrated knowledge of haplotype blocks, selection signature regions and sequencing data to refine the list of PCGs. Results Estimated genomic heritability for internal organ traits in chickens ranged from low (LUNGWT, 0.06) to high (GIZZWT, 0.45). A total of 20 unique 1 Mb windows identified on GGA1, 2, 4, 7, 12, 15, 18, 19, 21, 27 and 28 were significantly associated with intestine length, and weights or percentages of liver, gizzard or lungs. Within these windows, 14 PCGs were identified based on their biological functions: TNFSF11, GTF2F2, SPERT, KCTD4, HTR2A, RB1, PCDH7, LCORL, LDB2, NR4A2, GPD2, PTPN11, ITGB4 and SLC6A4. From those genes, two were located within haplotype blocks and three overlapped with selection signature regions. A total of 13,748 annotated sequence SNPs were in the 14 PCGs, including 156 SNPs in coding regions (124 synonymous, 26 non-synonymous, and 6 splice variants). Seven deleterious SNPs were identified in TNFSF11, NR4A2 or ITGB4 genes. Conclusions The results from this study provide novel insights to understand the genetic architecture of internal organ traits in chickens. The QTL detection performed using a high density SNP array covered the whole genome allowing the discovery of novel QTL associated with organ traits. We identified PCGs within the QTL involved in biological processes that may regulate internal organ growth and development. Potential functional genetic variations were identified generating crucial information that, after validation, might be used in poultry breeding programs to reduce the occurrence of metabolic disorders. Electronic supplementary material The online version of this article (10.1186/s12864-019-6040-3) contains supplementary material, which is available to authorized users.


Background
Poultry breeding programs have for many generations been selecting animals mainly based on feed efficiency, growth, performance and carcass traits, achieving significant improvements in broiler chicken production [1][2][3][4]. However, muscle development has been disproportionate to development of internal organs, increasing the susceptibility of improved broiler chickens to metabolic disorders, and death [5][6][7].
Although the weights of internal organs have not been direct targets for selection in poultry breeding programs, feed intake (commonly targetted for selection) has exhibited positive genetic correlations with heart weight (0.77), liver weight (0.73) and intestine weight (0.92) in a meat-type chicken population [8]. Thus, the detection of genomic regions and/or potentially causative mutations associated with internal organ weight traits may provide important information to poultry breeding programs, facilitating selection of chickens with proportional development of internal organs, improved feed efficiency and reduced susceptibility to metabolic disorders.
The pulmonary and cardiac capacities of commercial broilers are similar to those observed in non-selected chickens [9], despite their larger body size, such that the heart and lungs can be overloaded [10,11]. This increases the occurrence of ascites, a metabolic disorder characterized by accumulation of fluid in the abdominal cavity of birds [10,11], hypertension syndrome (PHS) [12,13] or even heart failure and sudden death syndrome (SDS) [14][15][16][17]. Moreover, changes in the weight of the gizzard can influence nutrient availability in the small intestine, and intestine length is associated with feed efficiency and growth in chickens [18,19]. Long-term intense selection for carcass traits can also increase the occurrence of fatty liver syndrome [20]. Thus, the discovery of genomic regions and positional candidate genes (PCGs) for intestine length or weights and percentages of heart, liver, gizzard or lungs, is a first step for understanding their genetic architecture.
The use of a high density SNP panel covering the whole genome combined with genomic prediction methodology can provide more accurate QTL detection and facilitate the identification of positional candidate genes (PCGs). The main goal of this study was to estimate genomic heritability, and detect QTL and PCGs for internal organs weights and percentages in chickens. Knowledge of haplotype blocks, selection signature regions and the use of sequencing data were integrated to refine the list of PCGs.

Descriptive statistics and genomic heritability
The number of chickens analyzed for each internal organ trait, along with the averages, standard deviation, variance components (genetic, residual and total variances) and genomic heritability are in Table 1. The estimated genomic heritability ranged from low (0.06) to high (0.45).

Genome-wide association study (GWAS)
We estimated the proportion of the genetic variance explained by each one of the 943 non-overlapping 1 Mb windows for all the traits analyzed (Additional file 1). The Manhattan plots of the posterior means of the proportion of genetic variance explained by each SNP window across the 28 autosomal chromosomes for internal organ traits are in Fig. 1.
The characterization of the significant genomic windows is in Table 2. A total of 28 significant 1 Mb windows (with different positions) were identified on GGA1, 2,4,7,12,15,18,19,21,27, and 28 for LIVP, GIZZWT, GIZZP, LUNGP and INTES. The posterior probability of association (PPA) [30] ranged from 0.38 to 0.8 for each significant region. The percentage of the variance explained by the window ranged from 0.53 to 3.26.
Adjacent genomic windows associated with the same trait were considered as the same QTL and their respective proportions of the variance explained were summed. For GIZZWT, two QTL were identified: one on GGA1 (rs15488636-rs13972577) explaining 5.23% of the genetic variance and one on GGA4 (rs317228804 -rs316997412) explaining 2.17% of the genetic variance. For GIZZP two QTL were identified: one on GGA18 (rs315188853 -rs312861933) explaining 2.98% of the genetic variance and one on GGA28 (rs313774457-rs317469562) explaining 2.71% of the genetic variance. For INTES, one QTL was identified on GGA4 (rs317228804 -rs15617120) explaining 4.98% of the genetic variance. Therefore, considering the adjacent windows associated with the same trait, as the same QTL, a total of 20 QTL were detected herein.

Overlap with previously reported QTL
Sixteen QTL identified in this study have previously been annotated/converted in Gallus_gallus-5.0 assembly, available in the Chicken QTL databaserelease 35 [29]. Three detected QTL overlapped with known QTL mapped for internal organ traits (Table 3) in different populations; one detected in a Brazilian paternal broiler population through SNP association analysis (QTL #137389) [23], and two mapped in a meat-type chicken population derived from a F2 cross through linkage analyses by interval mapping method (QTL #96657, QTL #96660) [21].
Four of the 20 QTL overlapped with QTL previously reported for the same population for fatness traits using the same methods applied in this study [32] (Table 4).

Positional candidate genes
A total of 756 genes have been annotated in the regions defined by the 20 QTL (Additional file 2). Among these, 14 PCG (Table 5) were identified as promising candidates based on their Gene Ontology (GO) terms and information available in the literature. The characterization of selection signature regions overlapping with the associated genomic windows and the characterization of haplotype blocks harboring the SNPs with the highest model frequency in each window are available in Additional files 3 and 4, respectively.

SNPs located in positional candidate genes
A total of 13,748 SNPs were annotated in the 14 PCGs. From this total, approximately 2% exhibited multiple annotations, being classified in different annotation categories as reported by Moreira et al. [34]. The number of SNPs annotated and their respective locations in each PCG is in Table 6.
A total of 156 SNPs were annotated in coding regions (124 synonymous, 26 non-synonymous and six splice site variants) and for each of those, the SIFT score was used to predict deleterious mutations. Seven predicted deleterious SNPs were identified: one in the TNFSF11 gene, three in the NR4A2 gene, and three in the ITGB4 gene (Table 7).

Descriptive statistics and genomic heritability
The means and standard deviations observed herein for HWT and LIVWT are in agreement with the values observed in a previous GWAS study for organ traits performed by Dou et al. [31], using an F2 population derived from reciprocal crosses between a commercial broiler line and an indigenous strain [31], but differ for GIZZWT. The means and standard deviations of the internal organ weights and percentages reported herein differ from those reported by Venturini et al. [36] and Gaya et al. [8], however, they used meat-type chickens from broiler lines selected over many generations for rapid growth and meat production [8,36,37]. Those   Fig. 1 Manhattan plots of the posterior means of the percentage of genetic variance explained by each 1 Mb SNP window across the 28 autosomal chromosomes for heart weight (a), heart weight as a percentage of body weight (b), liver weight (c), liver weight as a percentage (d), gizzard weight (e), gizzard weight as a percentage (f), lung weight (g), lung weight as a percentage (h), and intestine length (i). The X-axis represents the chromosomes, and Y-axis shows the proportion of genetic variance explained by each window from Bayes B analysis. Red lines indicate the threshold to deem significant SNP windows  differences are expected since the F2 population used herein was founded from reciprocal crosses between grandparental broiler and layer lines, that exhibit different phenotypic abilities [6,[38][39][40].
Our genomic heritability estimates ranged from low to high for the traits analyzed in the Embrapa F 2 Chicken Resource population (Table 1). Previous studies reported heritability estimates using pedigree information for internal organ traits in meat-type chickens [8,36,41]. The heritability estimates for HWT, GIZZWT, and LIVWT reported in one study were 0.38, 0.39 and 0.25 [8], and in another 0.27, 0.44 and 0.33 [36]. Thus, the heritability estimates for GIZZWT reported in different populations (0.39 and 0.44) corroborates our finding (0.45). In addition to the expected differences in the heritabilities due to the methods used to estimate, it is important to consider that we used an F2 population obtained from reciprocal crossings between two distinct lines, thus, differences from the heritability estimated reported in meat-type populations would not be surprising.
A previous study using a larger dataset from the same chicken population comprising 2063 F 2 -TCTC population, reported pedigree-based heritability for HWT, LIVWT, GIZZWT, LUNGWT and INTES equal to 0.19 ± 0.04, 0.14 ± 0.05, 0.70 ± 0.05, 0.60 ± 0.05 and 0.04 ± 0.02, respectively [42]. The genomic heritability estimates for HWT and LIVWT were similar to the heritability reported herein and, the differences observed in GIZZWT, LUNGWT and INTES estimates could be due the reduced sample size used herein, which may not represent the genetic variance for these traits in the whole population.
In summary, our genomic heritability estimates ranged from low (0.06) to high (0.45) ( Table 1). Based on the high genomic heritability estimated for GIZZWT, GIZZP and INTES, a response would be expected for selection focused on these traits.

GWAS and QTL discovery
We detected 28 significant 1 Mb QTL associated with LIVP, GIZZWT, GIZZP, LUNGP or INTES explaining from 0.53 to 3.26% of the genetic variance, with the PPA ranging from 0.38 to 0.80 (Table 2).
Considering the adjacent associated QTL and their respective cumulative proportions of the variance explained, we detected QTL explaining 5.23% (on GGA1) and 2.17% (on GGA4) of the genetic variance for GIZZWT, 2.98% (on GGA18) and 2.71% (on GGA28) for GIZZP, and 4.98% (on GGA4) for INTES. A total of HWT: heart weight; GIZZWT: gizzard weight; INTES: intestine length 1 Genomic positions based on Gallus_gallus-5.0 assembly (NCBI) 2 Chicken QTLdb ID numbers databaserelease 35 [29] [43], also explaining these greater number of QTL detected for the traits with high heritability estimates. The low PPA values observed in our study can be due to the reduced sample size used herein, that can reduce the power to detect QTL [43]. However, the overlaps observed with known QTL available in the Chicken QTL databaserelease 35 [29], and also the detection of PCGs within the QTL associated with internal organ weights and percentage corroborate our findings.
Only 16 known QTL mapped for the internal organs traits evaluated in this study are well annotated (have genome position available) in the Chicken QTL database release 35 [29]: four for GIZZWT, one for GIZZP, two for HWT, one for HEARTP, one for LIVWT, one for LUNGWT and six for INTES. Thus, not many overlaps between the 20 QTL detected herein with the 16 previously published QTL should be expected. Nevertheless, one QTL mapped for GIZZWT, and one mapped for INTES overlapped with known QTL mapped for internal organ traits (HWT and GIZZWT) in chickens, and one QTL for GIZZWT detected herein overlapped with one known QTL mapped for the same trait suggesting our findings are unlikely to represent false positive discoveries ( Table 3). The QTL detected on GGA1 and GG4 explaining 5.23 and 2.17% of the genetic variance for GIZZWT, respectively, also overlapped with QTL detected in a GWAS performed using the high density SNP array (600 K, Affymetrix) and layer chickens from a Chinese F2 population [31], validating our findings.
We did not observe overlaps with the QTL detected by Moura et al. [5] in a mapping study performed with the same F2 population, using a different number of animals, low density of markers (approximately 127 microsatellite markers) and an interval mapping method. It is important to highlight that in the study by Moura et al. [5], only 11 significant QTL were detected, and from those, six showed similar effects across sexes [5] for HWT, GIZZWT, GIZZP and LIVP. Additionally, different QTL can be detected  Downstream  57  94  104  127  70  129  54  80  60  16  127  80  83  88   Intron  290  1431 332  436  150  1849  13  -404  1789  811  3138 14  490 Non-synonymous -  Deleterious low confidence: little sequence diversity in this position affecting the substitution model, and consequently the averages of conservation value [35] depending on the method and the density of markers used. In this context, in a study performed to identify QTL that control performance under variable temperature conditions in chickens, Lien et al. [44] compared QTL detected in an F2 chicken population by GWAS using univariate linear mixed model (GEMMA software) and using interval mapping method (QTLMap software). This study showed that eight QTL were detected by QTLMap, 47 were detected by GEMMA and from those, only two were common to both approaches [44]. The detected QTL for GIZZWT on GGA21 (1 Mb) overlapped with the known QTL #96660 ( Table 3) that had been localized to a wider region of 3.4 Mb. The use of a high density of SNP and GWAS approach provided us with better resolution in the detection of this QTL, reducing approximately 70% in the size of the associated region.
Four QTL mapped for GIZZWT, GIZZZP and INTES overlapped with QTL previously mapped for fatness, using the method and the exact same birds from Embrapa F2 Chicken Resource population (Table 4). These QTL detected herein may exhibit pleotropic effects with abdominal fat and carcass fat content traits, thus, selection based on genotype for these internal organ traits might affect fat deposition. Moreover, changes in the gizzard weights can influence nutrient availability in the small intestine; intestine lengths are associated with feed efficiency and growth in chickens [18,19].
We confirmed and refined several known QTL, and we discovered 17 novel QTL associated with LIVP, GIZZWT, GIZZP, LUNG and INTES; from the 20 QTL detected, only three QTL overlapped with known QTL mapped for internal organ traits in chickens (Table 3). These QTL can help to understand the genetic architecture of metabolic important organs that potentially affects feed efficiency and the occurrence of metabolic disorders, such as pulmonary hypertension [5].
The Ligand dependent nuclear receptor corepressor like (LCORL) gene is located within the QTL on GGA4 (76 Mb) associated with GIZZWT, and this region overlapped with a QTL previously associated with the same trait in an F2 Chicken population [31]. This gene has been reported as a candidate gene for body size in horses [45] and for performance and carcass traits in cattle [46], suggesting pleotropic effects. Additionally, LCORL gene is located within a haplotype block that harbors the SNP with the highest model frequency, indicating that genetic variants in this gene are in linkage disequilibrium with the QTL. These evidences point to LCORL as a candidate gene for GIZZWT in chickens, however, further validation is needed to confirm our findings.
The Protein tyrosine phosphatase, non-receptor type 11 gene (PTPN11) is located within the QTL on GGA15 (6 Mb), explaining 0.7% of the genetic variance for LIVP and exhibited GO term annotated for organ growth (GO:0035265). This gene encodes a tyrosine phosphatase (SHP2) that plays a role in biological processes such as cell growth and differentiation [47]. SHP2 affects many biological functions in the gastrointestinal tract [47], and in a study with mice, animals knocked-out for SHP2 exhibited an attenuated hepatocyte proliferation [48].
Within the QTL on GGA1 (rs15488636 -rs13551715) explaining 5.23% of the genetic variance for GIZZWT, two PCGs were identified: RB1 and TNFSF11. The Retinoblastoma 1 (RB1) gene encodes proteins from retinoblastoma (RB) family [49] and exhibited GO term annotated for digestive tract development (GO:0048565). In mammals, this family plays a role in the regulation of cell cycle [49] and in the proliferation of pre-adipocytes in chickens [50]. A previous study performed with mouse reported that RB1 expression can affect the mass of internal organs such as liver, spleen, lungs and heart during postnatal development [51]. The TNF superfamily member 11 (TNFSF11) gene encodes members of the tumor necrosis factor (TNF) family, and this family plays a role in different biological processes such as cell growth and differentiation [52], bone remodeling process and osteoclast differentiation [23]. In a study performed in a broiler chicken population, associations of the TNFSF11 gene with weight gain, feed conversion and heart weight were identified [23]. Additionally, this gene exhibited GO term for animal organ morphogenesis (GO:0009887). Further studies with this gene may help elucidate its role in gizzard weight regulation.
The integrin subunit beta 4 (ITGB4) gene is located within the QTL on GGA18 (rs315188853 -rs312861933) explaining 2.98% of the genetic variance for GIZZP and encodes ß4 integrin polypeptides [53]. In humans, several studies reported both lethal and nonlethal variants in ITGB4 gene responsible for epidermal and gastrointestinal disorders [54][55][56][57][58]. Genetic variants in the ITGB4 gene were associated with the occurrence of pyloric atresia [54,55], which can lead to a distended stomach in humans [56]. Additionally, this gene exhibited GO term annotated for digestive tract development (GO:0048565).
The Solute carrier family 6 member 4 (SLC6A4) gene is located within the QTL on GGA19 (6 Mb) explaining 1.01% of the genetic variance for INTES and exhibited GO term annotated for negative regulation of organ growth (GO:0046621). This gene belongs to the solute carrier (SLC) family of proteins that encodes membranebound transporters [59]. In an expression profile study of the SLC gene family in intestine from the late embryonic to early post-hatch stages, the SLC6A4 gene exhibited a transient peak of expression during embryonic day of hatch [60]. In humans, SLC6A4, also known as SERT or 5-HT gene, is expressed across many regions of the intestine [61]. This gene is expressed in human intestinal epithelial cells and changes in its expression pattern can lead to intestinal disorders [62]. In chickens, intestinal disorders can affect feed efficiency, and consequently, animal performance [18,19], and this is consistent with the observation that feed intake exhibited a positive genetic correlation with intestine weight (0.92) in a meattype chicken population [8]. Thus, the SLC6A4 gene should be thoroughly investigated as a candidate for regulation of intestine length in chickens.
Additionally, three PCGs (PCDH7, NR4A2 and GPD2) overlapped with selection signature regions identified in the founder lines of Embrapa F2 Chicken Resource Population [33], suggesting that they are under positive selection, affecting internal organ weights in at least one of the lines (broiler or layer line).
The Protocadherin 7 (PCDH7) gene is located within the QTL on GGA4 (72 Mb) associated with INTES and GIZZWT. This gene is a member of protocadherin gene family and encodes an integral membrane protein that plays a role in cell-cell recognition and adhesion [63]. In a study with humans, PCDH7 was reported as a hepatic stellate cell surface marker [63] and hepatic stellate cells play important role in the occurrence of liver fibrosis [64]. Thus, the PCDH7 gene can be responsible for liver disorders. Chen et al. [65] found PCDH7 with a high expression in human gastric mucosa tissues. This gene is located within a haplotype block that harbor the SNP with the highest model frequency indicating that genetic variants in this gene are in linkage disequilibrium with the QTL. This is the first report showing some evidence to support PCDH7 as a candidate for regulation of INTES and GIZZWT in chickens. Further functional validation studies are necessary to confirm our findings.
The Nuclear receptor subfamily 4 group A member 2 (NR4A2) and Glycerol-3-phosphate dehydrogenase 2 (GPD2) genes are located within the QTL on GGA7 (36 Mb) associated with INTES. Both genes play a role in glucose homeostasis [66,67], and consequently, they are involved in lipid metabolism. In mice, changes in glucose levels can lead to a high susceptibility to intestinal barrier dysfunction and enteric infection [68]. These disorders can affect intestinal permeability [69], affecting animal health and performance [70]. Thus, both genes are candidates for INTES regulation in chickens.
Based on the GO terms and the literature information, the 14 genes discussed here are potential candidates for regulation of the weights of internal organs or their proportionate growth, which influences their percentages. Potentially functional genetic variants in these genes may affect expression or can also affect protein function, may leading to changes in the phenotype. Thus, the identification of genetic variants can help to elucidate and support the role of the candidate genes in the regulation of internal organs traits in chickens. Consequently, this knowledge can help to mitigate the occurrence of metabolic disorders that leads to poultry production losses.

Discovery of potential causative SNP in PCGs
Approximately 81% of the 13,748 SNPs in PCGs were annotated in potentially neutral regions (introns), while 19% were annotated in potentially functional regions (Table 6). Even considered as a potentially neutral region, SNPs in introns can be associated with the phenotype, playing a role in gene expression and alternative splicing regulation as well, affecting the mRNA transport, providing important information to understand the genetic architecture of the trait [71][72][73]. However, SNPs located in introns were not deeply investigated in this study.
Several studies have reported genetic variants in potentially functional regions (3'and 5′-UTR, down and up-stream from the gene and in exonic regions classified as synonymous and non-synonymous) associated with important traits in chickens [22,40,[74][75][76][77]. However, only a few studies have focused on the investigation of deleterious and high impact variants associated with important traits in chickens [34,[78][79][80] that can cause loss of protein function or produce a truncated protein.
Thus, in the search of potentially functional mutations we focus on predicted deleterious and high impact variants annotated in our PCGs.
Seven predicted deleterious SNPs were identified in three candidate genes (Table 7), all annotated as missense variants; variants that results in amino acid changes but the length is preserved [81]. One SNP (rs312442403) was annotated in the TNFSF11 gene, a positional candidate for GIZZWT and this genetic variant can be responsible for changes in the tumor necrosis factor (TNF) family affecting cell growth and differentiation.
Three SNPs (g.36224286 > C/T, g.36225242 > G/T and g.36225278 > C/T) were annotated in the NR4A2 gene, a positional candidate gene for INTES and these genetic variants might be responsible for changes in glucose levels driving the occurrence of intestinal disorders. Three other SNPs (rs313279811, rs316340790 and rs732847450) were annotated in the ITGB4 gene, positional candidate for GIZZP, and these genetic variants can result in epidermal and gastrointestinal disorders. It is important to highlight that these mutations are potentially functional and further association and functional studies should be performed to validate their role in phenotype regulation. No high impact variants were detected in our PCGs.
From all the PCGs, four (RB1, HTR2A, NR4A2 and GPD2) were already selected as candidates for regulating fat deposition in the same population in a study by Moreira et al. [32], suggesting that they may exhibit pleiotropic effects with fat deposition traits.

Conclusions
High genomic heritability estimates were reported for some traits indicating that they can respond successfully to selection based on genotypes. QTL were detected for LIVP, GIZZWT, GIZZP, LUNGP and INTES, across the whole-genome. Within these QTL, PCGs involved in biological processes that might regulate internal organ growth and morphogenesis, digestive tract development, and metabolic disorders and diseases were identified. These biological processes may influence feed efficiency and also the occurrence of metabolic disorders leading to losses in poultry production. Potentially functional genetic variants annotated in PCGs were identified and these mutations can responsible for phenotypic changes. Further validation studies can elucidate their role in internal organ traits regulation, allowing its possible application in poultry breeding programs aiming to select chickens with higher feed efficiency and lower susceptibility for metabolic disorders, thus improving chicken production.

Methods
All experimental protocols related to animal experimentation in this study were performed in agreement with resolution number 011/2010 approved by the Embrapa Swine and Poultry Ethics Committee on Animal Utilization (CEUA) in Concordia, Santa Catarina State -South of Brazil. These protocols are in agreement with the rules of the National Council of Animal Experimentation Control (CONCEA) to ensure compliance with international guidelines for animal welfare. The population, genotypes and methods described are the same as those adopted by Moreira et al. [32] in a previous GWAS performed by our group for fat deposition traits.

Chicken population
Embrapa F 2 Chicken Resource Population used in this study was developed by the Embrapa Swine and Poultry National Research Center, for QTL mapping and genomic studies [6,38], from reciprocal crosses between one layer line (known as CC) and one broiler line (known as TT). The grandparental lines had previously undergone multi-trait selection over many generations. The CC line had been selected for improved egg production, egg weight, feed conversion, viability, sexual maturity, fertility, hatchability, egg quality and for low body weight, while the TT line had been selected for improved body weight, feed conversion, carcass and breast yield, viability, fertility, hatchability, and for reduced abdominal fat and metabolic syndromes. More details about the Embrapa F 2 Chicken Resource Population are in Nones et al. [6] and Rosário et al. [38]. A total of 28 parental (from the broiler and layer lines), five F 1 and 496 F 2 -TCTC chickens were genotyped with a 600 K Axiom Chicken Genotyping Array, which contains SNPs segregating for different chicken lines [27].
The chickens and genotypes are the same as those adopted by Moreira et al. [32] in a previous GWAS performed in the same population, that was capable to detect reliable QTLs for fat deposition traits. Moreover, this sample size is in agreement with other previous GWAS that utilized 600 K Axiom Chicken Genotyping Array and also detected reliable QTLs for production traits [82,83].

Phenotypic measurement
After 6 h of fasting, chickens at 42 days of age had their live body weight measured (BW42) and were then euthanized by cervical dislocation. A blood sample from each chicken was immediately collected for subsequent DNA extraction, then the carcasses were eviscerated, the internal organs (heart, lungs, liver and gizzard) were weighted, the intestine length was measured [5] and the carcasses were cooled. The percentage traits for each of the organs were calculated by dividing their weight by BW42 and multiplying by 100. More details about measurement of internal organs in the F 2 population are in Moura et al. [5].

DNA extraction, genotyping and quality control
The genomic DNA from 529 chickens was extracted from blood samples with DNAzol® protocol. The DNA integrity was evaluated in agarose gel (1%), quantified in a NanoDrop® spectrophotometer (Thermo Fisher Scientific) then diluted to the final concentration of 20 ng/μL. All subsequent DNA preparation steps for genotyping were performed following recommended Affymetrix protocols. The samples were genotyped with the 600 K Axiom Chicken Genotyping Array [27].
Samples and genotypes quality control were performed using Affymetrix Power Tools v1.17.0 (APT) software. Samples with DishQC ≤0.82 and call rate ≤ 90% as well as SNPs with call rate ≤ 98%, minor allele frequency (MAF) ≤ 2% and those not annotated in autosomal chromosomes, were removed from the dataset. The genome position of SNPs was taken to be that on the Gallus_gallus-5.0 chicken assembly (NCBI). At any particular locus, the small number of missing genotypes were replaced by their average covariate value representing the number of copies of one of the alleles at that locus across all the individuals with an observed genotype at that locus [84].

Descriptive statistics and heritability
Descriptive statistics (average and standard deviation) for each trait were calculated using basic functions in R software (http://www.r-project.org/). The variance components reported (genetic, residual and total variance) were obtained from the means of the posterior distributions of those components, obtained by fitting a Bayes B model in GenSel software [85]. The posterior distribution for heritability (the ratio of genetic variance divided by total variance) was constructed from each sample of the posterior distributions of the variance components, and the posterior mean is reported as the estimate of the genomic heritability.

Genome-wide association analysis
Informative and reliable SNPs (after quality control filters) were kept for GWAS analysis with a Bayesian approach in GenSel software [85]. Initially, Bayes C was used to estimate the variance components and then, genetic and residual variances were used as priors to run a Bayes B model for inference about genomic effects [78,84]. The mathematical model was: where, y represents the vector of phenotypic values for one trait, X is the incidence matrix for fixed effects, b is the vector of fixed effects, K is the number of SNPs, a j is the column vector representing the SNP as a covariate for locus j coded with the number of B alleles, βj is the random substitution effect for locus j, assumed to be normally distributed N (0, σ 2 βi ) when δj = 1 but βj = 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 equation and e is the residual associated with the analysis. In the model, sex and hatch were included as fixed effects and BW42 as a fixed covariate for weight traits. We tested the significance of the effects using a linear regression model in R software (http://www.r-project.org/).
Parameterization of the BayesB model assumed π = 0.9988 and involved obtaining 41,000 Markov Chain Monte Carlo (MCMC) samples from which the first 1000 samples were discarded. The genome was partitioned into 943 non-overlapping 1 Mb windows, based on genomic locations of SNP provided in a map file. In order to select associated windows, we assumed that each window was expected to explain approximately 0.1060% of the genetic variance (100%/943) based on an infinitesimal model [30,83]. Windows explaining five-fold more than expected (0.53%) were considered biologically significant and only those windows were used as the basis for identification of PCGs.

Overlap with known QTL
We checked the overlap of associated genomic windows with known QTL for internal organs available in the Chicken QTL databaserelease 35 [29], accessed in July, 2018. The overlaps were checked using QTL coordinates according to the Gallus_gallus-5.0 (NCBI) chicken genome assembly, available in the BED file. We excluded those QTL without genomic positions and also those not properly annotated (provisionally annotated between 0 and 100 bp), leaving 16 well annotated QTL for overlap analysis. Genomic windows that did not overlap with known QTL were considered novel.
We checked the overlap between our associated genomic windows and the QTL for internal organ weights recently reported from the analysis of a Chinese F2 population [31]. Additionally, we investigated the overlap between the QTL detected herein and QTL detected for fatness traits, using the same method and population [32].

Identification of positional candidate genes
We considered the flanking positions from the genomic windows associated with internal organ traits in order to search for PCGs. An entire gene list was obtained using Ensembl BioMart tool [86,87]. In order to refine that gene list and select PCGs, we adopted two criteria: (1) genes located within haplotype blocks harboring the SNP with the highest model frequency in each associated genomic window (the haplotype blocks detection was performed using PLINK v.1.9 [88] software with default parameters); and (2) genes that overlapped with selection signature regions identified in a previous study [33], which evaluated 28 grandparental chickens from the two lines that generated the F 2 population analyzed in this study.
The previous study [33] used whole genome sequences to identify genetic variants and applied Fst method [33,89] to estimate the divergence between populations to identify regions under selection (TT vs. CC lines). We used the CrossMap tool (http://crossmap.sourceforge.net/) to convert selection signature region coordinates from Gallus_gallus-4.0 to the Gallus_gallus-5.0 chicken genome assembly (NCBI).
Additionally, we investigated genes that exhibit biological processes from GO terms directly related to internal organs growth and development. Further, the genes selected in each step were investigated in different databases (NCBI, OMIM), and in the literature, to support/refute their selection as PCGs.

Detection of potential functional mutations
We searched for deleterious and high impact SNPs annotated in PCGs to identify potentially functional mutations. A dataset of approximately 13 million reliable SNPs were initially identified in the whole genome as a result of sequencing 28 grandparental chickens from the Embrapa F 2 Chicken Resource Population, and after filtration, approximately 14 million SNPs were further utilized for the search of SNPs annotated in PCGs. A detailed description of the sequencing, SNP calling, and filtering criteria are available in Boschiero et al. [43].
The filtered SNPs dataset was annotated using Variant Effect Predictor (VEP) tool [90]. For the variants located in coding regions, the VEP tool was also used to calculate the SIFT (sorting intolerant from tolerant) score. This score is an assessment of the level of conservation in homologous protein sequences, used to predict whether amino acid changes caused by SNPs are tolerant nor not (may affect the function of the gene product). Coding SNPs with SIFT scores < 0.05 are predicted as deleterious, whereas coding SNPs with SIFT scores ≥0.05 are predicted as tolerated [91].
Prediction of the putative variant impact was also performed using the VEP tool [90]. All variants annotated as transcript ablation, splice acceptor, splice donor, stop gained, frameshift, stop loss, start lost and transcript amplification which may cause protein truncation, loss of function or trigger nonsense mediated decay were predicted as high impact variants (https://www.ensembl. org/info/genome/variation/prediction/predicted_data. html) [78].
Availability of data and material All SNPs reported (identified by sequencing) were submitted to dbSNP-NCBI with the handle "LBA_ESALQ" (https://www.ncbi.nlm.nih.gov/SNP/snp_ viewBatch.cgi?sbid=1062434). The datasets used and/or analyzed during the current study (genotypes from SNP array and phenotypes) are available from the corresponding author on reasonable request. Ethics approval and consent to participate All experimental protocols related to animal experimentation in this study were performed in agreement with the resolution number 010/2012 approved by the Embrapa Swine and Poultry Ethics Committee on Animal Utilization to ensure compliance with international guidelines for animal welfare.