Ethics statement
The work was approved by the Animal Management Committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS, Beijing, China). Ethical approval regarding animal survival was granted by the animal ethics committee of IAS-CAAS (approval number: IASCAAS-AE20140615).
Phenotypes and breeding value estimation
In chickens at 42 days of age, 10 µL of blood was freshly obtained from the wing. Slides were air-dried and dyed with May–Grunwald–Giemsa stain, and the H/L ratio were calculated. The basic data have been used in recent GWAS analysis, and details regarding the genotyping and phenotyping have been reported[14]. In this study, a conventional pedigree-based best linear unbiased prediction model (BLUP) was used to predict breeding values for each trait.
The pedigree-based BLUP model[15] is
y = Xb + Za + e
where y is the vector of the phenotypic records of the trait (H/L), b is the vector of the fixed effect (Batch and Sex), X is the incidence matrix linking b to y, a is a vector of additive breeding values to be estimated, Z is the incidence matrices linking a to y and e is the vector of residuals. We assumed that var(a)=Aσ2a, where A is the pedigree-based genetic relationship matrix. The prediction of breeding values(ebv) was performed the Asreml package[16]. The ebv calculated in this part is used for the next GWAS analysis
Genotyping and quality control
Briefly, 1,317 birds were genotyped with Axiom® SNP arrays[17] with the Affymetrix GeneTitan® system according to the procedure described in the Axiom 2.0 Target Prep 384 Samples Protocol (Affymetrix; Santa Clara, CA, USA) (https://www.thermofisher.com/). Stringent quality control (QC) procedures were performed for the genotype data by using PLINK[18] version 1.09(http://zzz.bwh.harvard.edu/plink/), After QC[call rate >95%, minor allele frequency >0.05, and extreme deviation from Hardy-Weinberg proportions(P > 1E-6)], 1,212 animals and 36,800 SNPs, located on 28 autosomes and in the Z-chromosome, were retained. Slight differences in the number of individuals and SNPs across the recent GWAS analyses are attributed to phenotypic editing.
Single SNP-based association analysis
The basic MLM model was used for association analysis in GAPIT[19]. The P-value was corrected with a strict Bonferroni adjustment based on LD pruning[20]. The sums of the independent LD blocks plus singleton markers were used to define the number of independent statistical comparisons [21]. Finally, 8,068 independent SNPs were used to determine the P-value thresholds, including 5% genome-wide significance (6.80E-6,0.05/8,068), and suggestive association (1.24E-4, 1/8,068). The Quantile-quantile (Q-Q) plot and Manhattan plots of GWAS for H/L were produced with the “qqman” packages[22] in R. The most significant SNP genotypes were added to univariate models to analyze the phenotypic differences among individuals with different SNP genotypes.
Pathway-based association analysis
The pathway-based analysis workflow is represented in Figure 1. In brief, for each trait, nominal P-values < 0.05 from the GWAS analyses were used to identify significant SNPs. With VEP software (http://asia.ensembl.org/Tools/VEP), SNPs were assigned to genes if they were within the genomic sequence of a gene or within a flanking region of 5 kb up- or downstream of a gene, to include SNPs located in regulatory regions. The Ensembl GRCg6a database was used as a reference (http://asia.ensembl.org/Gallus_gallus/Info/Index). The background SNPs represent all SNPs tested in the GWAS analyses, whereas the background genes were the genes associated with those SNPs. For assignment of the genes to functional categories, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway[23] databases were used. Pathway enrichment analysis was performed with the OmicShare tools, a free online platform for data analysis (www.omicshare.com/tools). Significantly enriched pathways in significant genes to compared with the background genes were defined with a hypergeometric test.
Proportion of phenotypic variance explained by SNPs (PVE)
PVE estimation is equivalent to calculating the heritability of SNPs,
The GBLUP model was used to estimate the phenotypic variance explained by SNPs of all genotyped individuals[24]:
y = Xb + Zg + e
where the definitions of y, X, Z and e are the same as those in the BLUP model, g is the vector of genomic breeding values to be estimated, following a normal distribution of N(0, Gσg2), in which σg2 is the variance of additive genetic effect, and G is the marker-based genomic relationship matrix. Z is the incidence matrix linking g to y, and e is the vector of residuals. e is the vector of random errors, following a normal distribution of N (0, Iσe2), in which σe2 is the residual variance. The estimation of variance components was performed with the Asreml package[16].