Jingxing-yellow chickens, a domestic broiler breed, were obtained from Changping Base of Chinese Academy of Agricultural Sciences, from which a directional selection line was established based on the 56-day-old peripheral blood H/L ratio of Jingxing-yellow chickens. Thirty roosters and 90 female chickens with low H/L phenotype were selected as parents in each generation, based on the principles that one rooster was mated with three hens, roosters and hens in the same family were not full siblings and hens within a family were not half siblings.
The H/L ratio was evaluated at 56 days old in each generation. Fresh blood was extracted from the wing vein and stored in EDTA tubes. One drop of blood was smeared on to a glass slide for each bird. Once air-dried, smears were stained using a Wright’s Giemsa solution (G1020, Solarbio, Beijing) according to manufacturer's instructions for microscopic observation. One hundred leucocytes including lymphocyte, heterophil and monocyte cells were counted in different fields on each slide, using a Leica DM500 microscope with a magnification of 1000× and the H/L ratio was calculated by dividing the number of heterophils by the number of lymphocytes .
Detection of phagocytosis and oxidative burst function of heterophils
Heterophils were isolated from the peripheral blood of chicks at seven days post-hatch using the commercially-available chicken heterophils isolation kit (Solar Technology Co., China). Phagocytic ability was determined by flow cytometry by the detection of fluorescein isothiocyanate dextran (FITC-dextran), 40 kDa (Chondrex, USA) using a previously published method with modification . Approximately 1×106 isolated heterophils were suspended in one mL RPMI 1640 including 10% fetal bovine serum (FBS) and 1% chicken serum and incubated together with one mg/ml FITC-dextran and 100nM lipopolysaccharide (LPS) at 37 °C with 5% CO2 for two hours. Incubation was stopped by mixing with ice-cold phosphate buffered solution (PBS), after which the cells were washed two times with cold PBS and analyzed using a FACS Calibur flow cytometer (FACS Calibur, BD Bioscience, USA). The endocytic activity was expressed as a positive rate of cells with FITC+.
Production of oxide during the oxidative burst was determined by oxidation of intracellular dihydrorhodamine 123 (DHR123) to green, fluorescent rhodamine 123 . Approximately 1×106 isolated heterophils were suspended in one mL RPMI 1640 including 10% FBS and 1% chicken serum. The heterophils were stimulated by one μg/ml phorbol myristate acetate (PMA), after which one μg/ml DHR123 was added to the medium and co-incubated at 37°C with 5% CO2 for two hours. Heterophils were washed twice with cold PBS, resuspended in PBS, and analyzed with a BD FACS flow cytometer (BD Bioscience, USA). Oxidative burst activity was reported as a positive rate of cells with DHR123+. Flow cytometry data were analyzed with FlowJo software version 10 (TreeStar, Ashland, OR) and results per group were compared using Student’s t-test. All data were shown as mean ± standard error mean (SEM).
Whole genome sequencing and indel detection
Whole-genome resequencing of 368 individuals comprising G5 of selection line (G5), G9 of selection line (G9), non-selection line (NS) and F2 of G9 constructed by G9 of the selection line and non-selection line, was performed. Forty-three chickens obtained from G5, 92 chickens obtained from G9 and 92 chickens obtained from NS, were used in this study. A total of 141 F2 chickens were also used in the current study, obtained from the cross of G9 of selection line and non-selection line. In cross design, the F1 population were established by 15 roosters (G9) mated with 30 hens (NS) and the F2 population was established by using the F1 hybrids based on the principle of random mating. The chickens’ genomic DNAs were extracted from blood samples from wing veins using the standard phenol-chloroform extraction procedure. The quality of the DNA was examined using the NanoPhotometer® (Implen, CA, USA), the concentration of the DNA was determined by Qubit® 3.0 Flurometer (Life Technologies, CA, USA) and the integrity of the DNA was evaluated using agarose gel electrophoresis. For genome sequencing, 0.5 μg of genomic DNA from each sample was used to construct paired-end sequencing libraries using an Annoroad® Universal DNA Fragmentase kit V2.0(AN200101-L) and an Annoroad® Universal DNA Library Prep Kit V2.0 (AN200101-L) following the manufacturer’s recommendations and index codes and sequenced on the Illumina HiSeq X Ten Sequencer (Illumina Inc) at a mean depth of ×10.
Raw data were processed with Perl scripts to ensure the quality of data used in further analysis. The adopted filtering criteria were the removal of the adaptor-polluted reads where reads containing more than five adapter-polluted bases were regarded as adaptor-polluted reads and would be filtered out, removal of the low-quality reads, that had a number of low-quality bases and a phred quality value less than 19 meaning that more than 50% of total bases were regarded as low-quality reads and removal of reads with the number of N bases accounting for more than 5 %. For paired-end sequencing data, both reads were filtered out if any read of the paired-end reads was adaptor-polluted. The obtained clean data after filtering was analyzed statistically on its quantity and quality, including Q30, data quantity and base content statistics.
After sequence quality filtering, the Burrows–Wheeler aligner (BWA 0.7.12)  was used to map the clean reads to the chicken reference genome GRCg6a, GCF_000002315.6. Samtools v1.2  were used to sort reads and duplicate reads resulting from polymerase chain reaction (PCR) were removed using Picard tools v1.13 (http://broadinstitute.github.io/picard/). Indels were filtered before the following analysis with GATK VariantFiltration protocol .The filtering settings were as follows: QD<2.0, ReadPosRankSum<-20.0, FS>200.0, QUAL <30.0, DP<4.0. The annotation was performed using ANNOVAR  for all the qualified variants based on the GFF file.
Population structure and phylogenetic tree
Principal component analysis (PCA) was conducted using the PLINK software (Version: 1.90). For the first PCA plot, all 227 chickens were used, with the first three principal components cumulatively explaining 15.41% of the total variance. Thirty individuals were selected in each population for subsequent analysis. Population structure was evaluated using ADMIXTURE (Version 1.3.0) (https://dalexander.github.io/admixture/) and assuming two ancestral populations (K). After converting the PLINK indel matrix to distance matrix, the “meg.” format obtained by the “plink.distance.matrix.to.mega.pl” script, then MEGA7.0  and iTOL software  were used to visualize the phylogenetic trees.
Detection of selective sweeps
All indels passing the quality control were used to detect signature of selection in G9 of selection line and non-selection line within the 184 genomes included 92 individuals in each population, respectively. The Fst fixation and nucleotide diversity ratio (θπ) were calculated within 40 kb windows sliding in 20 kb steps. The θπ ratios of the two population were calculated and then log2-transformed the data. Those windows comprised within the top 5% quantile of all two statistics were considered as candidate outliers and annotated using the ANNOVAR based on the GFF file. Functional enrichment for Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses was assessed with Metascape .
Genome-wide association study
The univariate linear mixed model (LMM) implemented in GEMMA version 0.98.1 software (https://github.com/genetics-statistics/GEMMA/releases) was used for the association analysis .The first three principal component values PCA eigenvectors derived from the whole-genome indels were fitted as fixed effects in the mixed model to correct for population stratification . The parameter ‘--indep-pairwise 25 5 0.2’ was used to infer effective independent tests. Considering the over-conservation of 5% Bonferroni correction method, the suggestive and whole-genome significance cutoff was defined as the Bonferroni threshold, which is 1.00/ independent indels (-log10 [P] =5.01) and 0.05/ independent indels (-log10 [P] =6.31). Manhattan and quantile–quantile (Q-Q) plots were constructed by the “qqman” package (http://cran.r-project.org/package=qqman) in RStudio-1.3.1093.