The study population was obtained from a Chinese indigenous chicken breed that has been maintained for 25 generations by the Wens Nanfang Poultry Breeding Co. Ltd (Yun Fu, Guangdong, China). We obtained 435 male birds from the 25th generation, which were generated by 30 males and 360 females from the 24th generation and were reared with the same recommended nutritional and environmental conditions. The FCR was calculated according to recorded data during the feeding trial. For more details about this population, please refer to Zhang et al [10].

All of 435 male chickens from the 25th generation and 15 sires of the 24th generation were selected for genotyping. The genomic DNA from the 450 birds was extracted from blood samples using the NRBC Blood DNA Kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer’s instructions, and DNA samples were analyzed for genome DNA concentration and integrity. DNA samples of the appropriate quality were genotyped using 600K Affymetrix Axiom HD chicken genotyping array, which contains 580,961 SNPs distributed on 28 autosomes, two linkage groups (LGE64 and LGE22C19W28_E50C23), and two sex chromosomes. Then, twenty-four key individuals were selected for sequencing from the 435 birds and their 15 sires based on a strategy of maximizing the expected genetic relationship, while maximizing the proportion of unique genomes sequenced in the population. In this process, we used G matrix to calculate the genetic relationship between key individuals and the remaining population, details of the selection of key individuals in our study were described by Ye et al [11,12]. Finally, 21 males and 3 sires were selected as key individuals and re-sequenced with 150-bp paired-end reads on the Illumina HiSeq 3000 platform with the average sequence depth of 14.62. Sequencing reads were aligned to the Gallus gallus 4.0 genome using the Burrows-Wheeler Alignment tool [13]. Duplicated reads were removed using Picard release 1.119 (http://sourceforge.net/projects/picard/files/picard20tools/1.119/). GATK software [14] was used for SNP calling with the default parameters. Finally, 11,645,758 SNPs remained for further analysis and the concordance between chip data and sequence data was 98.0%. After that, the SNP chip was imputed to sequence data using Fimpute software [15] with key sequenced individuals as the reference. For more details about the genotype imputation, please refer to Ye et al [12]. Quality control both for the chip and imputed sequences genotypes were conducted use PLINK 1.90 software [16] with the same criteria: SNP call rates >0.97, minor allele frequencies >0.05, obeys Hardy-Weinberg equilibrium (p>0.00001), and individual call rate >0.95. After quality control, 447,766 SNPs from SNP chip and 8,626,020 SNPs from imputed sequences remained for subsequent analysis. These SNPs distributed on autosome with the average distance between adjacent SNPs ranged from 92.63 to 168.27 bp among different autosomes.

Mixed linear model was used for the GWAS, and the analysis was performed with GEMMA software (Zhou and Stephens [17]). The statistical model is described as follows:

where Y is the vector of phenotypic values for all individuals, b is a vector of fixed effects including the batch effect, which has three levels, a is the substitution effect of the SNP under consideration, u is the random additive genetic effect; X and Z are the design matrices for b and u, respectively; S is the design vector for a; e is the vector of random residuals. In this model, u and e were assumed to have the structure u~N(0,Gσa2) and e~N(0,Iσe2), where G is the additive genomic relationship matrix which was derived from the sequence data, G=1p∑i=1p(Xi-InX¯i)(Xi-InX¯i)T, X is the n×p matrix of genotypes of n individuals’ p markers, X_i as its ith column representing genotypes of ith SNP X̄_i, as the sample mean and I_n as a n×1 vector of 1’s, σa2 is the polygenic additive variance, I is an identity matrix and σe2 is the residual variance. To compare the GWAS using the sequence data with chip data, the GWAS was performed using the chip data with the same model. DMU software [18] was used to estimate the heritability explained by the sequence data with the same model as the one described above.

To visualize the result from the GWAS, Manhattan plots and quantile-quantile plots were drawn by the qqman package [19] in R for each trait. In this study, the significant level was adjusted based on the effective number of independent tests, which was calculated by the simpleM software [20].

Candidate genes which were located within or nearby the significant SNPs (within 0.3 Mb) were identified by Ensembl (http://www.ensembl.org/index.html) and NCBI (https://www.ncbi.nlm.nih.gov/) annotation of the Gallus gallus 4.0 genome version.

SNPs with p-value <1.0×10⁻⁵ were compared with all reported QTLs that associated with FCR. The QTLs were obtained from the Chicken QTLdb (https://www.animalgenome.org/cgi-bin/QTLdb/GG/index). Based on the physical location, QTLs close to the significant SNPs were extracted, with a maximum distance of 5 Mb between the QTL and the SNP to be compared.

Descriptive statistics for trait FCR are presented in Table 1. We found that the phenotypes of FCR were normally distributed (Shapiro-Wilk test, p>0.05), and its heritability was 0.33.

The effective number of independent tests was 631,181 as calculated by the simple method. Hence, the threshold p-value was adjusted to 7.92×10⁻⁸ (0.05/631,181) for a genome-wide significance level, and to 1.58×10⁻⁶ (1.00/631,181) for a genome-wide suggestive significance level.

In the case of using imputed whole sequence, one genome- wide significant SNP on chromosome 8 locating at 5th intron of zinc finger and BTB domain containing 41 (ZBTB41) was found to be associated with FCR, and 9 SNPs reached the suggestive significance level (Table 2). Among these 9 SNPs, 3 of them were located between 45.43 Mb and 45.61 Mb of chromosome 1 with the nearest genes being ubiquitin specific peptidase 44 (USP44), leukotriene A4 hydrolase (LTA4H), and ETS transcription factor (ELK3), respectively. One SNP was located on chromosome 2 and 60 Kb distance away from its nearest R-spondin 2 (RSPO2) gene. Two SNPs were located on a narrow region of chromosome 4 (from 15.59 Mb to 16.04 Mb), and the nearest genes were inhibitor of apoptosis protein 3 (IAP3) and sosondowah ankyrin repeat domain family member D (SOWAHD). Another two SNPs were located on the positions of 1.39 Mb and 2.71 Mb of chromosome 8 respectively, with the first SNP locating on the 14th intron of calmodulin regulated spectrin associated protein family member 2 (CAMSAP2) and the other locating on the 4th intron of potassium sodium-activated channel subfamily T member 2 (KCNT2). The remaining 1 SNP was located on the position of 3.20 Mb of chromosome 26 and member of RAS oncogene family (RAP1A) was its nearest gene. In the case of using chip data, only one SNP reach a suggestive significance level, and the SNP was also detected by using imputed whole sequence (Figure 1).

In our study, p-value of 31 SNPs was smaller than 1.0×10⁻⁵. They were compared with the reported QTLs collected from Chicken QTLdb. The results are shown in Supplementary Table S1. A total of 32 QTLs on autosomes were reported to be associated with FCR, and 8 of those QTLs were found near the 31 SNPs (Supplementary Table S1). Besides, seven SNPs located on chromosome 1 and 2 located far from reported QTL. Three out of the 7 SNPs located on a narrow region from 102.70 Mb to 106.46 Mb. The distribution of chicken FCR related QTLs as compared with SNPs with p-value lower than 1×10⁻⁵ are graphically displayed in Figure 2.