INTRODUCTION
The internal organs are involved in many biological occurrences, including oxygen transportation, lipid metabolism, and digestion. Each organ has its own developmental process and gene expression profile [1]. Liver development affects abdominal fat, yolk deposit [2] which are all associated with economic performance, and provides a model of non-alcoholic steatohepatitis [3]. Heart is an organ for carrying oxygen [4], affect well-being of the chicken [5]. Physical and chemical digestion takes place in the proventriculus and gizzard, which regulate feed intake and energy balance [6]. Efficient utilization of feed is essential to reduce feed cost and improve economic benefit. Moreover, internal organ like gizzard and heart provide proximate composition and amino acid same as chicken meat [7]. In order to support an economical viable meat system, it is essential to efficient utilize these byproducts. Hence, improve these internal organs will not only improve well-being of the organism but also take full advantage of byproduct from chicken.
With the development of deep sequencing and statistical methods, it is possible to detect the quantitative trait loci (QTLs) of complex traits. In recent years, great advances have been gained in the selection of economic traits in chickens with molecular biology techniques. Nones et al [8] identified novel QTLs associated with gizzard, liver, heart, lung, and thigh located on Gallus gallus autosome (GGA) 1. Rosário et al [9] found novel QTLs on GGA1, GGA3, and GGA4 related with carcass traits. Boschiero et al [10] identified microsatellite markers associated with performance and carcass traits on GGA1 and GGA13. Most studies on chicken internal organ weight have focused on the early stages of development, which has hindered the progress of similar research in the layer hen, especially for spent layer hens in China.
In this study, we conducted a genome-wide association study (GWAS) on an F2 resource population generated from Dongxiang Blue-Shelled (DX) and White Leghorn (WL) chickens with a 600 K Affymetrix chip. Furthermore, a linkage disequilibrium (LD) analysis, variance estimated percentage, and gene annotation were performed by biotechnology. This research will provide the genetic architecture that underlies the weight of internal organs in layer chickens.
MATERIALS AND METHODS
Ethical statement
Animal procedures and sample collections were carried out following the guidelines established by the Ministry of Agriculture of China. The Institution of Animal Care and Use Committee in the Poultry Institute, Chinese Academy of Agricultural Science, Yangzhou, China approved the procedures.
Phenotype collection and analyses
We constructed an F2 resource population by reciprocal crosses between WL and DX chickens. The WL was from a commercial line and DX from a Chinese indigenous strain, therefore the two lines represent two very different genetic pools, including morphological, physiological, and production traits, such as feed efficiency [11], follicle number [12], and the yolk proportion [13]. The F2 resource population of 1512 birds from 49 half-sib and 590 full-sib families were produced from a WL/DX (25 ♂:407 ♀) and DX/WL (24 ♂:235 ♀) cross in the F1 generation. The F1 population was generated from WL (6 ♂)×DX (133 ♀) and DX (6 ♂)×WL (80 ♀) by initial reciprocal crosses, producing 1,029 and 552 chicks, respectively. The experimental animals were kept indoors under standardized conditions at the research base in the Jiangsu Institute of Poultry Science, Yangzhou, China. Birds had ad libitum access to feed and water that met all NRC requirements and were housed in single-hen cages under a 16 L:8 D lighting regime. The F2 birds were slaughtered at 72 weeks; immediately after the slaughter we measured the weights of heart (HW), liver (LW), proventriculus (PW), and gizzard (GW). DNA was collected by standard venipuncture. Data including a descriptive phenotype and normality test statistics were summarized using R software.
Genotyping and quality control
The standard phenol/chloroform method was used to extract genomic DNA, which was genotyped against a 600 K Affymetrix Axiom Chicken Genotyping Array (Affymetrix, Inc. Santa Clara, CA, USA). Affymetrix Power Tools v1.16.0 (APT) (http://affymetrix.com/) software was used to analyze genotype calling and quality control (QC). A total of 1,512 hens and 532,299 single nucleotide polymorphisms (SNPs) remained valid after the application of APT for QC. An R script was used to calculate the SNP QC metrics and filter out individual SNPs falling below given thresholds. The effectiveness of the detecting quality was tested by PLINK v1.90 [14]. We imputed for some sporadic missing genotypes using the BEAGLE v4.0 package [15]. Briefly, samples with a dish QC of more than 0.82, call rate less than 97%, minor allele frequency (MAF) higher than 5%, Hardy–Weinberg equilibrium test P of more than 1×E10−6 were kept for further analysis. After QC, 1,512 individuals and 435,867 SNPs were used in further GWA analysis.
Genome-wide association analysis
A principal component analysis was used to eliminate spurious associations resulting from the presence of cryptic relatedness or hidden population stratification. Then, the full SNP set to 41,130 independent SNPs were pruned via the –indep-pairwise 25 5 0.2 command (PLINK), and we included five principal components (PCs) as covariates in the mixed model. After simpleM was used to determine the threshold by correcting the number of multiple tests, the genome-wide significant and suggestive p-values were 8.43×10−7 (0.05/59,308) and 1.69×10−5 (1.00/59,308), respectively.
The GEMMA v0.94 package [16], with the exact mixed model approach, was implemented with the valid individuals and SNPs for univariate analysis. The independent SNPs were used to compute the centered relatedness matrix, and the significance p-value level between SNPs and phenotypes was calculated from the derived Wald test. The univariate linear mixed model was denoted as:
Where y represents an n-vector of phenotypic values for n individuals; W is a matrix of covariates (fixed effects with a column of 1s and top five PCs), α is a vector of the corresponding coefficients including the intercept; x is an n-vector of the genotypes of the SNP marker, β is the effect size of the marker; μ is an n-vector of random effects; and ɛ is an n-vector of errors.
The Manhattan plots and quantile–quantile (QQ) plots were drawn by “gap” packages in R software. The genomic inflation factor was calculated by using the GenABEL package in the R project.
Linkage disequilibrium analysis and gene annotation
The LD were implemented by Haploview v4.2 to test the association of significant SNPs with each other. When the distance between two consecutive genome-wide significant SNPs was greater than 10 Mb, they were considered as two separate QTL. The LD analyses would find a haplotype blocks because of these loci located in a strong LD region. But, a genuine causal locus could not be distinguished by present GWAS method. Therefore, a functional gene annotation on these significant SNPs were essential performed to characterize potential candidate genes. In our study, we searched candidate genes based on Galgal5 assembly, using BioMart system supported by Ensemble (http://asia.ensembl.org/Gallus_gallus/Info/Index). The genes nearest or harboring significant SNPs of internal organ weight were chosen as candidate locations.
SNP effects and chromosome heritability
The contributions to phenotypic variance (CPV) explained by the significant association SNPs was calculated with GCTA software [17] using univariate restricted maximum likelihood. The heritability explained by the eligible SNPs (h2snp) for GWAS or the variance contributed by each chromosome were also estimated with GCTA.
RESULTS
Phenotypic statistics and genetic parameters
Descriptive data for each organ are displayed in Table 1. After rank-based inverse normal transformation, all phenotypic values conformed to a normal distribution. The range of all internal organ is large, and with the coefficient of variance more than 19%. Table 2 displayed the genetic analysis, which showed that the heritability of internal organs was moderate to high, and the GW had the highest heritability (0.640). Genetic correlation between each organ weight was higher than 0.5, the highest genetic correlation was PW and GW (0.711).
Genome-wide association study results
The results of the GWAS were as follows. There are 380 significant SNPs on GGA1, 49 SNPs on GGA4 were related to internal organ weights, and six SNPs on GGA2 were related to HW (Table 2, Figure 1, Supplementary Figure S1, S2, S3). The regions ranged from 165.63 to 173.78 Mb on GGA1, and from 71.30 to 77.19 Mb on GGA4. After Venn diagram analysis, 17 and 6 significant SNPs located on GGA1 and GGA4, respectively, were found to affect all internal organ weights (Figure 2).
Estimation of contributions to phenotypic variance and plausible candidate genes
By using gene annotation analysis, it was obtained that the 17 significant SNPs on GGA1 were harbored in nine genes, general transcription factor IIF polypeptide 2 (GTF2F2), potassium channel tetramerization domain containing 4, spermatid associated, succinate-coA ligase ADP-forming β subunit, 5-hydroxytryptamine receptor 2A (HTR2A), WD repeat and FYVE domain containing 2 (WDFY2), asparagine-linked glycosylation 11, NIMA related kinase 3, and cytoskeleton associated protein 2. The SNP on GGA1 possessing a MAF of more than 0.4 and the values of beta on GGA1 were all negative (Table 4). We then analyzed the allelic contribution to phenotypic variation (Table 4), the results show that the smaller the p value, the larger the CPV value. All significant SNPs on GGA1 explained over 5% of CPV for PW and GW, while the CPV for HW and LW was from 2.781% to 4.723%. The CPV of SNPs rs312726815, rs315120631, rs13972990 were highest for HW, LW, and PW/GW. The rs312726815 SNP located in the upstream (16,831 bp) of HTR2A, rs315120631 located on the intron of GTF2F2, rs13972990 was downstream (95,255 bp) of WDFY2.
Table 5 displayed that the SNPs on GGA4 possess a MAF more than 0.05 and less than 0.06, with values of beta being positive. These SNPs located in or near genes ligand dependent nuclear receptor corepressor like, non-SMC condensin I complex subunit G (NCAPG), leucine aminopeptidase 3, LIM domain binding 2, quinoid dihydropteridine reductase, prominin 1, and transmembrane anterior posterior transformation 1. We found one locus rs14491030 that located in the missense of the gene NCAPG and another SNP close to NCAPG.
Table 6 showed that the SNPs on GGA2 association with HW association with HW displayed a MAF of ~0.3, with the values of beta being positive. These SNPs located on the intron of sonic hedgehog (SHH) and close to genes 5-Hydroxytryptamine receptor 5A, insulin induced gene 1, and limb development membrane protein. The SNPs rs316413705 explained the highest CPV (3.109%) for HW.
Linkage disequilibrium analysis
The LD analysis was performed because of the potentially strong LD between neighboring variants. The LD analysis results revealed that significant SNPs on GGA1, GGA4, or GGA2 were all extremely strong in a LD status (Figure 3). Combined the aforementioned CPV calculation analysis, LD analysis and gene annotation, we considered rs312726815, rs315120631, and rs13972990 on GGA1 as primary candidate loci, with correspond to genes HTR2A, GTF2F2, and WDFY2 associated with HW, LW, and PW/GW, respectively. The NCAPG gene harbored the missense locus rs14491030 on GGA4 were considered as a candidate gene for all internal organs, SHH gene with intron harbored the SNP rs14491030 on GGA2 had an association with HW.
Genome partitioning of genetic variation
The genetic architecture of internal organs was further illustrated by partitioning the genetic variation onto chromosome segments with an exploratory analysis. HW and PW in the joint model could not converge due to the relatively small sample size in the F2 population. Only the partitioning spectrum of LW and GW were estimated. The estimates of variance contributed by each chromosome exhibited a medium and strong linear relationship with the length of the chromosome for LW (R2 = 0.493, 0.493 Figure 4A) and GW (R2 = 0.753, 0.493 Figure 4B). For GGA1 this explained 18.29% and 10.03% of the phenotypic variance for GW and LW, respectively. While for GGA4 it explained 11.30% and 4.09% for GW and LW, respectively.
DISCUSSION
To accomplish the goal of “500 eggs in 100 weeks”, it is essential to continuous select egg production trait as well as fitness trait (such as increasing heart capacity [18]) for which populations were not routinely selected. Because of the correlation between egg production and fitness showed positive or negative was not absolute, it is essential to detect molecular markers associated with economic traits as well as fitness traits to improve animals simultaneously. Moreover, some body composition trait like gizzard and proventriculus indirectly reflects feed efficiency [5]. Hence, improving these traits will not only make chicken well-being better but also take advantage of these organs. But, these traits are inconvenient or difficult to measure for animals should be killed. Application of DNA markers could simultaneously select economic trait and these traits.
Our results showed that phenotypic data displayed a large variation for all internal weights, which was probably because these traits were not chosen as selection index. The heritability showed a moderate to high, higher than that the previous reports. For LW, the heritability estimates were 0.31 [19] and 0.33 [20] at 6 weeks, and the heritability of HW and GW was 0.27 and 0.44, respectively, in a report by Venturini et al [20]. The difference of the heritability may be due to different breeds and ages [21]. In addition, we dissected genome partitioning of genetic variation for GW and LW. A positive linear correlation between the variance explained by each chromosome and its length were obtained in our work, which corresponds to previous reports [22]. In particular, GGA1 accounted for the largest genetic variance (18.29%) for GW, which is consistent with the highest heritability (0.640), whereas GGA1 accounted for 10.03% genetic variance for LW, the heritability of which is moderate (0.355). We then considered that with the higher heritability of the trait, the genetic variance of GGA will also be higher. This finding is in line with previous research especially in our resource population [23].
Using GWAS method, we totally found 17 SNPs on GGA1 and 6 SNPs on GGA4 associated with all internal organ weights, and 6 SNPs on GGA2 were associated with HW. Previous studies showed that QTL for internal organ weight, such as HW and LW, at 9 weeks was detected on GGA1 and GGA4, respectively [24]. HW and GW at 6 weeks were mapped to GGA13 [10]. The difference between the present study and previous reports might be due to differences in the age and population of the birds used in different studies. These studies used younger birds, the age of the birds in the present study was 72 weeks. In addition, one trait, such as internal organ weight, was controlled by more than one QTL, which will help us understand the genetic architecture underlying the quantitative traits that are controlled by polygenes.
The most significant loci in the GGA1 spanned from 166.99 to 169.66 Mb, which is also related to comb weight [25], egg weight [23], ovary weight [26], and feed intake [27] in our previous studies and growth traits in other reports [28–30]. The internal organs as fitness traits participate in metabolism that affects growth performance and economic traits, therefore the selection of body weight or egg weight will indirectly result in internal organ weight change. This has important ramifications for understanding the pleiotropic effects of the locus or gene [31].
After rigorous statistic and LD analysis, three SNPs rs3127 26815, rs315120631, and rs13972990, correspond to genes HTR2A, GTF2F2, and WDFY2, respectively, were selected as the candidate loci. HTR2A was considered as a candidate gene for HW. The mRNA of the gene, which encodes a receptor of 5-hydroxytryptamine (serotonin), is expressed within vascular smooth muscle, endothelial cells, and cardiomyocytes, and plays a role in vasoconstriction and cellular proliferation [32]. Also, the function of the HTR2A-encoded protein is associated with blood pressure and heart rate and is a contributing factor to cardiovascular disease in humans [33]. In birds, the size of heart was association with the rate of metabolism and indicates the capacity to move the blood [34]. Therefore, it is reasonable to propose that HTR2A was responsible for heart weight, the detailed information needs to be further investigated. Second gene GTF2F2 nearby SNP rs315120631 was considered as candidate gene for LW. GTF2F2 forms a heteromeric general transcription initiation factor [35]. GTF2F2 widely exists in many tissues and organs, especially in the liver, lung, and kidney. The liver is the main organ regulating body homeostasis [36], where a large amount of active transcription and translation takes place. Reports said that GTF2F2 is association with rat liver regeneration [37]. Liver is responsible for detoxification and energy balance, which involved in various of biological process including transcription and translation. Hence, further studies about the function of GTF2F2 in chicken internal organs are needed. The third SNP, rs13972990, is near gene WDFY2, which is associated with PW and GW. WDFY2 (WD repeat and FYVE domain containing 2) is an endosomal protein, modulating the PI3K/AKT pathway, which is known to be involved in oncogenesis [38]. WDFY2 overexpression can increase adipogenesis, which may play a role in metabolic disorders such as diabetes in humans [39]. Proventriculus and gizzard modulated feed behavior [6] and feed overconsumption [5], indirectly regulating glycogen synthesis. It is tempting to speculate that WDFY2 may participate in proventriculus and gizzard growth though its involvement is not clear.
Moreover, the region on GGA4 was related to egg weight, eggshell weight, oviduct weight in our population [23,40]. NCAPG harbored the missense rs14491030 was a candidate gene for all internal weight. The gene was related to residual feed intake in cattle [41] and withers height in horses [42] in previous reports. Most breeders focus on egg production in layer chicken, resulting indirectly in selection for fitness trait or digestive organs. NCAPG was considered as candidate gene for all internal organ weight though no selection for these traits, pleiotropic effects of the locus on other traits that are under selection, or close linkage and LD with QTL that are under selection [43]. We provide an evidence that quantitative traits that may be controlled by polygenes and a single gene or mutation may lead to a host of alterations in multiple traits.
In addition, GWAS results showed that one QTL on GGA2 associated with HW. The SHH could be considered as a candidate gene. The SHH gene is one of the Hh proteins, which plays a crucial role in the development of all animals and regulates morphogenesis of a variety of tissues and organs in the embryo [44]. Previous studies showed that SHH is necessary for secondary heart field proliferation in humans [45]. Consequently, the function of SHH in the chicken heart may involve heart differentiation.
CONCLUSION
In present study, a GWAS strategy was performed to detect potential QTLs or genes association with internal organ traits. Our study provides evidence that the internal organ weight at 72-weeks old appears to moderate a high heritability and may share the similar genetic mechanisms. Five candidate genes were identified with significant effect on internal weights. These promising loci or genes could be helpful to simultaneous improve layer chicken economic trait and fitness trait to accomplish sustainable use of the chicken.
PDF Links
PubReader
ePub Link
Full text via DOI
Download Citation
Supplement
Print
