- Research Article
- Open Access
Detection of QTL for traits related to adaptation to sub-optimal climatic conditions in chickens
Genetics Selection Evolution volume 49, Article number: 39 (2017)
Growth traits can be used as indicators of adaptation to sub-optimal conditions. The current study aimed at identifying quantitative trait loci (QTL) that control performance under variable temperature conditions in chickens.
An F2 population was produced by crossing the Taiwan Country chicken L2 line (selected for body weight, comb area, and egg production) with an experimental line of Rhode Island Red layer R- (selected for low residual feed consumption). A total of 844 animals were genotyped with the 60 K Illumina single nucleotide polymorphism (SNP) chip. Whole-genome interval linkage mapping and a genome-wide association study (GWAS) were performed for body weight at 0, 4, 8, 12, and 16 weeks of age, shank length at 8 weeks of age, size of comb area at 16 weeks of age, and antibody response to sheep red blood cells at 11 weeks of age (7 and 14 days after primary immunization). Relevant genes were identified based on functional annotation of candidate genes and potentially relevant SNPs were detected by comparing whole-genome sequences of several birds between the parental lines.
Whole-genome QTL analysis revealed 47 QTL and 714 effects associated with 178 SNPs were identified by GWAS with 5% Bonferroni genome-wide significance. Little overlap was observed between the QTL and GWAS results, with only two chromosomal regions detected by both approaches, i.e. one on GGA24 (GGA for Gallus gallus chromosome) for BW04 and one on GGAZ for six growth-related traits. Based on whole-genome sequence, differences between the parental lines based on several birds were screened in the genome-wide QTL regions and in a region detected by both methods, resulting in the identification of 106 putative candidate genes with a total of 15,443 SNPs, of which 41 were missense and 1698 were not described in the dbSNP archive.
The QTL detected in this study for growth and morphological traits likely influence adaptation of chickens to sub-tropical climate. Using whole-genome sequence data, we identified candidate SNPs for further confirmation of QTL in the F2 design. A strong QTL effect found on GGAZ underlines the importance of sex-linked inheritance for growth traits in chickens.
Detection of quantitative trait loci (QTL) is a classical approach to better understand the genetic architecture of complex traits and unravel genomic regions that control quantitative variation of traits. Many QTL detection studies have been conducted in chickens, as reviewed by Abasht et al. , most of these using commercial or experimental lines, and only a few using local breeds. Since the publication of the Abasht et al.  review, additional QTL studies have been reported on egg laying [2, 3] and comb traits . In general, these QTL studies are conducted under standard management practice. In the current study, we carried out a QTL detection project under the real conditions of a humid sub-tropical climate in Taiwan, using an F2 cross between a locally adapted line, named L2, and an imported selected line, named R-. The L2 line is a closed line that has been selected for egg production as a dam line in the mating design used for the production of Taiwan Country chicken since 1982 . The R- line is a closed line that has been selected in France for adult residual food consumption since 1975 and was imported in Taiwan in 2003 . This work focused on traits that are measured during the growing period and that likely play a role in the adaptation of chickens to sub-tropical climate. Body weight was shown to play an important role in adaptation to high temperatures [7, 8], shank length and comb area are non-feathered areas that contribute largely to heat dissipation , and antibody response to sheep red blood cells is an indicator of immune response . The analysis used two different methods, genome-wide association analysis (GWAS) and QTL mapping, which were applied to genotypes obtained with the high-density 60 K Illumina iSelect chicken array .
A three-generation design was produced by crossing two pedigreed parental lines, L2 and R-, which are maintained at the experimental farm of the National Chung Hsing University (NCHU). Line L2 was established in 1983 from a sample of Taiwan Country chickens and has been maintained as a closed population since then . L2 males were selected on the basis of the size of the comb area (as an indicator of sexual maturity), body weight at 12 or 14 weeks of age, and average egg production to 40 weeks of age of full- and half-sisters. L2 females were selected for egg production to 40 weeks of age based on their own records and those of their full- and half-sisters. Since 1975, a divergent selection was undertaken at the French National Institute for Agricultural Research (INRA) to study the genetic basis of residual feed intake (RFI) in Rhode Island Red (RIR) layers . The R- line has been selected for low RFI (residual food consumption adjusted for body weight (BW), change in body weight (ΔBW) and total egg mass (EM) by multiple linear regression). Individual feed intake, BW and ΔBW, were measured between 33 and 37 weeks of age in birds of both sexes, and EM was recorded over the same period (i.e. 28 consecutive days) in females. R- males were selected on the basis of their own RFI value, whereas R- females were selected in two steps, combining their RFI value and their own egg production at 40 weeks of age. In 2003, a subset of the 28th generation of selection of R- line was shipped to NCHU to study the adaptation of R- chickens to sub-tropical climate . Since then, this line has been maintained as a closed population with a one-year generation interval.
In 2009, 46 F0 parents, born in the same batch, were used to generate the F1 generation by mating six L2 males to 15 R- females (LR) and seven R- males to 18 L2 females (RL). Then, six F1 males (four LR and two RL) were mated to 51 unrelated F1 dams; 32 RL females were mated to four LR males and 19 LR females were mated to two RL males. The family size of dams ranged from 1 to 35 (see Additional file 1: Figure S1). A total of 743 F2 individuals were produced in four successive batches with the following birth dates: 17 Dec. 2010, 18 Jan. 2011, 31 Jan. 2011, and 18 Feb. 2011 (see Additional file 2: Table S1). The experimental period lasted until 16 weeks of age for each batch.
All F2 chickens were reared in floor pens in an open-sided building, with a temporary fence to close the rooms, and additional heating (24 h/day) for the first two weeks (up to 29 °C). The temporary fences were removed at three weeks of age. The mean ambient temperature in the surrounding area ranged from 18 to 26 °C from batch 1 to batch 4, with minimal values ranging from 6.6 to 20 °C and maximal values ranging from 24 to 35.6 °C (see Additional file 3: Figure S2). Chicks were fed according to recommended nutrition standards, with a starter diet (metabolizable energy: 2830 kcal ME/kg and 19.14% crude protein) from hatch to 4 weeks of age and a grower diet (2818 kcal ME/kg and 16.11% crude protein) from 5 to 16 weeks of age. The same vaccination plan was applied to all F2 birds (see Additional file 4: Figure S3). Natural light was supplied during the rearing period, with an average day length close to 12 h.
Body weight was measured with an electronic scale (±1 g) at hatch (BW00) and every four weeks thereafter, i.e. BW04, BW08, BW12 and BW16. Growth rate was estimated by the difference between successive measurements of body weight every 4 weeks, i.e. between 4 and 8 weeks of age (BW0804), between 8 and 12 weeks of age (BW1208), and between 12 and 16 weeks of age (BW612). Length and height of the comb at 16 weeks of age were measured with a ruler (±1 mm), and comb area was computed as the product of comb length by comb height (CA16). Length of the right shank was measured with a Vernier caliper (±1 mm) at 8 weeks of age (SL08).
All F2 chickens were injected with 0.1 mL of a 0.25% suspension of sheep red blood cell (SRBC) at 11 weeks of age. Blood was sampled from the wing vein of each bird just before immunization to determine the baseline antibody level (SRBC00) and antibody levels at 7 and 14 days after primary immunization, i.e. SRBC07 and SRBC14, respectively. Sera were collected after centrifuging the blood samples, and stored at −20 °C until all assays could be run simultaneously. SRBC antibody levels were assayed by the hemagglutination test  and expressed in log2 units. All experiments involving animals in this study were done according to the approved protocol of the Institutional Animal Care and Use Committees of the NCHU (Taichung, Taiwan; IACUC No. 97-99).
Genotyping and quality control
At each of the three generations of the experimental design, blood was collected from the wing vein of each bird and genomic DNA was extracted using a commercial DNA extraction kit (DNeasy® Blood kit) and diluted to 50 ng/μl. After DNA quality control, each chicken was genotyped using the Illumina 60 K Chicken SNP iSelect chip. Of the 846 chickens (F0: 46, F1: 57, and F2: 743), two F2 individuals were excluded because SNPs for these individuals had a low call rate (<95%). The SNP set used in the current study included 57,636 SNPs. Approximately 30% (17,018) of the SNPs did not meet the following criteria and were removed: low call rate (<95%), low minor allele frequency (<0.05), or unknown chromosomal position on the chicken reference genome GalGal4. Pedigree checking was considered necessary before performing linkage analysis. An in-house software developed by D. Boichard (personal communication) was used to detect incompatible genotypes for each sire/dam/progeny triplet. Each incompatibility with either the sire, the dam, or both was recorded. When the total number of incompatible SNPs exceeded 2% of all SNPs, another in-house software, also developed by D. Boichard (personal communication), was used to assign the most probable parent(s) to each progeny. About 10% of pedigree errors were identified and corrected. When the pedigree was confirmed, only a few incompatible SNP genotypes remained for some F2 individuals compared to their F1 parents, and these were considered as genotyping errors rather than as pedigree errors and replaced by missing values for the corresponding F2 individuals. The final data included genotypes for 844 individuals (F0, F1, and F2) and 40,618 SNPs, distributed on 28 autosomes, two linkage groups and the Z sex chromosome. The F st between parental lines was calculated with PLINK (Version 1.9)  using the SNP genotypes from the 46 F0 parents.
The normality of trait distribution in the F2 population was checked by the SAS® UNIVARIATE procedure (Statistical Analysis System, Version 9.3, SAS, Institute Inc., Cary, NC, USA). Trait differences between parental lines were studied by analysis of variance using the GLM procedure in SAS® and including the effects of sex and line. For the F2 population, the GLM procedure was applied to estimate the fixed effects of sex, dam (full-sib family effect), and batch, taking into account BW00 as a covariate for all other body weight measurements, as well as for SL08 and CA16. For the SRBC traits, the model included BW08 instead of BW00 as a covariate, in addition to the fixed effects.
Genome-wide association analysis
Population stratification of the F2 individuals was assessed with multidimensional scaling (MDS) analysis available from PLINK (Version 1.0.7) . The indep-pairwise option with a window size of 25 SNPs, a step of 5 SNPs, and an r2 threshold of 0.2 was used to obtain independent SNPs. Pairwise identity-by-state (IBS) distances were calculated between all individuals using 3266 independent SNPs, and MDS components were estimated by the mds-plot option based on the IBS matrix. Linkage disequilibrium (LD) blocks were defined as a set of contiguous SNPs with pairwise r2 values exceeding 0.4, which resulted in 7063 LD blocks for growth-related traits and 7048 LD blocks for SRBC traits.
GWAS was carried out by using a univariate linear mixed model in GEMMA . The linear mixed model used in this study was applied for each chromosome, with sex, batch, and dam as fixed effects and by incorporating random genetic effects with genomic relationships to correct for genetic structure in the experimental population. The analysis was implemented on all data for the autosomes, while data from males and females were separated for the analyses on the Z chromosome. P value thresholds were computed based on the number of independent SNPs and LD blocks [17, 18], resulting in a P-value threshold of 5% Bonferroni genome-wide significance threshold of 4.84 × 10−6 (0.05/10329) and a suggestive linkage threshold of 9.68 × 10−5 (1/10329). Manhattan plots of the GWAS results for each trait were produced with the qqman package available from R (Version 3.1.2) .
Single QTL detection analyses were carried out with the QTLMap software , which was developed to handle data on F2 individuals for designs from outbred populations. QTLMap is an interval mapping method, in which QTL are detected using a likelihood ratio test (LRT) calculated under the hypothesis of one versus no QTL linked to the given SNPs. The linkage disequilibrium linkage analysis (LDLA) option of QTLMap was applied to single trait . No assumptions were made on the fixation of alleles in the grand-parental lines (L2 and R-) and on the number of segregating alleles at the QTL. The family structure of the F2 population consisted of six half-sib families, and the minimum size for considering dam families was set at 20 progeny, which was the case for 15 of the 51 dam families. Fixed effects of dam, sex, batch, as well as covariates (BW00 or BW08 depending on the trait studied) were taken into account in the model. For the Z chromosome, an interaction between the QTL and the fixed effect of sex was added to the model. The LRT value was calculated by scanning chromosomes for QTL with a step of 1 cM. The method proposed by Churchill and Doerge  was used to estimate empirical thresholds taking into account the family structure (pedigree) and genotype data (SNPs). Chromosome-wide significance of LRT values was obtained by 10,000 simulations under the null hypothesis (no QTL) to calculate 1 and 5% significance thresholds . Genome-wide thresholds were defined by applying the Bonferroni correction across chromosomes:
where n is the number of chromosomes used in the study. Confidence intervals for QTL were estimated by the one LOD drop-off method . The QTL substitution effect was assessed in each sire family at the position of the LRT maximum, and significance was assessed by Student’s t test. The additive value of a QTL effect was calculated as the mean of the absolute value of the sire substitution effects (P value < 0.05). SNP positions and information were annotated using the GalGal4 genome assembly . Information on potential candidate genes in each QTL region was searched for in the NCBI and Ensembl databases [24, 25]. Furthermore, for genes that were associated with significant SNPs, we performed a search of the gene ontology database  to draw hypotheses about the biological processes and molecular functions that likely influence the trait of interest.
In previous research, whole-genome sequence data was produced for one L2 chicken  and for two R- chickens . The sequenced animals were chosen at random in a recent generation for each line and constitute a resource for identifying SNPs that differ between the lines but were not used as F0 animals in the current F2 design.
Sequence alignment between the reference genome assembly GalGal4 and the whole-genome sequences of the L2 and R- individuals was performed to identify SNPs that differ between the two lines. Sequencing fragments with a Phred quality score lower than 15 and the adapter sequences were removed by Trimmomatic  before alignment to the chicken reference genome assembly Galgal4.84 . Trimmed reads that were shorter than 36 bp were then eliminated. Only read pairs for which both forward and reverse reads remained after trimming were mapped to the chicken reference by Burrows–Wheeler Alignment tool (BWA) .
Genome Analysis Toolkit (GATK) was applied to identify variants across the genome [31–33]. Mutations that were located within repeated sequences were marked by Picard  and discarded. In order to increase the accuracy of insertion and deletion (InDel) calling, reads around InDels were locally realigned by the IndelRealigner tool provided by GATK. The HaplotypeCaller tool of GATK was applied to detect variants. A read base with a Phred score less than 10 was not considered in the variant discovery step. Only mutations that were called with a Phred score equal to or higher than 30 were considered as confident mutations and reported. Our search for SNPs that differed between the L2 and R- birds was focused on the QTL regions that were detected by both GWAS and QTL mapping, with the aim to identify potential candidate mutations for further study. Functional effects of the identified mutations were predicted by SnpEff  and with Variant Effect Predictor (VEP) .
Phenotypic means of the F0 and F2 populations
None of the phenotypes recorded in this study significantly deviated from the normal distribution (P > 0.05). Means and standard deviations for the F0 lines and their F2 crosses are in Tables 1 and 2 for males and females, respectively. Highly significant differences were found between the L2 and R- lines for all growth traits (BW00 to BW16) and for CA16. Lines did not differ for immune response to sheep red blood cells, except for SRBC14 which was lower in R- females (P < 0.05) than in L2 females. Shank length was not available for the F0 individuals. Animals from the R- line were lighter than those from the L2 line, and also had a much smaller comb area. In the F2 generation, the fixed effects of sex, dam, and batch were significant, except for sex for BW00 and SRBC00, batch for BW08, and dam for BW1612, SRBC07, and SRBC14. The R- square values for the model ranged from 0.11 (SRBC00) to 0.77 (CA16). Regression coefficients on BW00 or BW08 were not significant for any trait (P > 0.10). Coefficients of variation for body weight traits and shank length ranged from 10 to 20% but were much higher for comb area (about 40 to 50% depending on sex) and for immune response traits (about 50 to 70%).
The F st between the parental lines was equal to 0.31, which shows significant differentiation between the lines. The structure of the F2 population was assessed by MDS analysis by applying an r2 less than 0.2. Plotting of the first two MDS components showed that F2 individuals clustered according to mating type (LR or RL F1 males) and half-sib family (four families from the LR males and two families from the RL males), which showed that the MDS components represented the family structure of the F2 generation (see Additional file 5: Figure S4). As expected, clustering of individuals within sires reflected the within-sire dam family structure (see Additional file 6: Figure S5).
Genome-wide association analysis
A total of 714 SNP effects associated with 178 SNPs were identified for nine traits with 5% Bonferroni genome-wide significance (P value <4.84 × 10−6) (see Additional file 7: Table S2). No SNP effect reached 5% Bonferroni genome-wide significance for BW00, CA16, SRBC00, and SRBC07. BW0804 had the largest number of genome-wide significant SNP effects (136), followed by SL08 (131) and BW08 (126) (Fig. 1). In contrast, a limited number of significant SNP effects was detected for SRBC14 (1), BW1208 (17), and BW1612 (20). In addition, 364 SNP effects associated with 234 SNPs reached suggestive significance (P value <9.68 × 10−5).
Only three chromosomes (GGA14, 24, and Z) carried genome-wide significant SNPs (Figs. 2, 3), whereas suggestive SNP effects were found on 11 chromosomes (GGA1, 2, 4, 7, 9, 14, 15, 22, 24, 26, and Z) (Table 3). GGAZ carried the most genome-wide significant SNP effects, which were detected only in females and for eight traits (BW04, BW08, BW12, BW16, BW0804, BW1208, BW1612, and SL08). The strongest association was found for a large region on GGAZ that spanned 13.1 Mb (between 8.9 and 22.0 Mb) and contained 146 genome-wide significant SNPs associated with eight growth-related traits. The most significant SNPs (3.57 × 10−42 < P value <1.74 × 10−47) were associated with BW08, BW0804, and SL08. This region contained 30 SNPs, harbored three annotated genes and one uncharacterized gene. Among these, the highest significance (P value ≤1.74 × 10−47) was obtained for SNPs in the prostaglandin E receptor 4 (PTGER4) and complement component 7 (C7) genes and in the uncharacterized gene (LOC100857889). Two overlapping chromosomal regions on GGA24 (between 4.1 and 5.3 Mb and between 4.3 and 6.0 Mb) were associated with BW04 and BW08, respectively. These two regions harbored 13 and seven genome-wide significant SNPs, respectively, and have been reported to carry several growth-related QTL . For SRBC14, we detected one genome-wide significant SNP (rs10724420) that was flanked by suggestive SNPs. This SNP was located in the 3′-UTR of the 3-phosphoinositide dependent protein kinase 1 gene (PDPK1), which is involved in transcriptional regulation of proopiomelanocortin in mice, with consequences for food intake, body weight and size . This region on GGA14 was previously reported in studies on body weight and immune related QTL .
The whole-genome single QTL analysis led to the identification of 47 QTL that corresponded to 34 non-overlapping regions on 20 chromosomes (Table 3), since there were several overlapping QTL regions for different traits. Thirteen of these 47 QTL had high genome-wide significance (P < 0.05) and corresponded to 10 non-overlapping regions on GGA1, 2, 4, 9, 14, 27, and Z. BW00 was associated with the largest number of genome-wide QTL regions, which were distributed on five chromosomes (GGA1, 2, 4, 9, and 27), whereas the QTL regions for other traits were distributed on GGA2 (BW12), GGA14 (BW04), and GGAZ (BW04, BW08, BW12, BW16, BW0804, and SL08). The QTL that were the most highly significant at the chromosome-wide level (P value <0.01) affected BW04 (GGA7), BW08 (GGA2 and GGA14), BW12 (GGA3), CA16 (GGA3 and 5), SRBC00 (GGA17), and SRBC14 (GGA4 and 19). All other identified QTL had chromosome-wide significance P values less than 0.05. Estimates of QTL allele substitution effects ranged from 0.18 (SRBC07 on GGA18) to 2.25 (CA16 on GGA26) residual standard deviations, with a mean of 0.58 standard deviations. Among the autosomes, GGA2 carried the largest number of detected regions, with six identified QTL, followed by GGA3, 4, and 17, which each carried three QTL. The average confidence interval spanned 2.4 cM, ranging from 2 to 5 cM.
Three QTL with chromosome-wide significance that affected SRBC traits on GGA18 and 21 were not reported before. The chromosome-wide QTL found for CA16 on GGA3 was not located at the same position as the one previously described on that chromosome for comb mass .
QTL detected by both GWAS and QTL mapping-LDLA
The position of the SNPs that flanked a QTL region was used to establish the correspondence between the results of QTLMap and those of GEMMA, by considering all SNP effects detected at a P value less than 0.05. Regions detected by both methods were found only on two chromosomes: GGA24 for BW04 and GGAZ for seven traits (BW04, BW08, BW12, BW16, BW0804, BW1612, and SL08). Thus, of the 47 QTL detected with QTLMap, eight were also detected with GEMMA. The 39 QTL that were detected only with QTLMap corresponded to 35 chromosome regions that overlapped with previously published QTL, for either growth- or immune-related traits .
The analysis of whole-genome sequence data focused on the 13 genome-wide significant QTL (10 non-overlapping chromosomal regions) detected by QTLMap, and the region identified with GEMMA on GGA24, which contained genome-wide significant SNPs but was only chromosome-wide significant with QTLMap. The L2 bird and the two R- birds differed for 15,443 SNPs in 106 genes (83 known genes, 16 genes coding for uncharacterized proteins, and seven miRNA); 1698 of these SNPs were novel. Among these 106 genes, 21 carried 41 missense mutations, as identified by VEP. The results are summarized according to QTL region in Table S3 (see Additional file 8: Table S3). An in-depth analysis focused on genes with known functions that were the closest to the position showing the highest probability in QTL mapping or GWAS.
Two genome-wide significant QTL were detected on GGA2 for body weight traits. The QTL that affected BW00 (between 41.1 and 42.4 Mb) contained 18 genes, of which four abhydrolase domain containing 5 (ABHD5), calpain 7 ( CAPN7), maturin, neural progenitor differentiation regulator homolog (MTURN), and testis and ovary specific PAZ domain containing 1 (TOPAZ1), are involved in the regulation of body weight, total amount of body fat, and lean body mass in mice . The TOPAZ1 gene carried 207 SNPs that differed between the L2 and the R- birds, of which 80 were novel SNPs, and one missense SNP that resulted in a valine being replaced by an isoleucine at position 431 of the protein. The QTL that affected BW16 (between 111.2 and 112.0 Mb) harbored two genes, cytochrome P450, family 7, subfamily a, polypeptide 1 (CYP7A1) and inositol monophosphatase domain containing 1 (IMPAD1), which are associated with birth weight  and increased body weight  in mice.
The region on GGA14 (between 2.8 and 3.3 Mb) was associated with a strong effect for BW04 and carried the ACTB gene, which is involved in the regulation of body weight and size in mice . This gene contained 75 SNPs that differed between the L2 and the R- birds, of which 12 were novel SNPs. Three QTL on GGAZ that contain several growth-related genes were detected in this study for BW16 (between 11.7 and 13.0 Mb), BW04 (between 13.7 and 14.7 Mb), and BW08, BW0804, BW12, SL08 (between 12.1 and 13.0 Mb). The PTGER4 gene carried 104 SNPs that differed between the L2 and the R- birds, of which 96 were novel SNPs. A previously described missense mutation (rs315266117, 12708649 bp) was detected in the genome of the R- birds; it is located close to the SNP with the lowest P-value in our GWAS analysis and results in the replacement of a glutamic acid by a lysine at position 375 of the protein, which is deleterious according to the SIFT prediction. The OSMR gene, which is involved in muscle organ and skeletal development in chickens , carried 144 SNPs that differed between the L2 and the R- birds, of which 136 were novel SNPs. One missense mutation (rs312389473) that resulted in the replacement of a threonine by a serine was located at position 501 of the protein. Among the other identified genes, PRKAA1 (protein kinase, AMP-activated, alpha 1 catalytic subunit) and RICTOR (RPTOR independent companion of MTOR, complex 2) are related to viscera weight [42, 43] and prostaglandin E receptor 4 (EP4), fibroblast growth factor 10 (FGF10), and leukemia inhibitory factor receptor (LIFR), which are located in the QTL region for BW16 (between 11.7 and 13.0 Mb on GGAZ), have been reported to be associated with body weight related traits [44, 45]. FGF10 and EP4 contained 189 and 106 SNPs that differed between the L2 and R- lines and each of these genes carried one missense SNP.
Furthermore, three QTL for BW00, on GGA4, 9, and 27, harbored four growth-related genes, among which ubiquitin-like modifier activating enzyme 6 (UBA6) on GGA4 and MDS1 and EVI1 complex locus (MECOM) on GGA9 play a critical role in body size/length and weight loss in mice [46, 47]. The QTL region on GGA27 contained two genes that regulate decrease in body weight in mice, i.e., LIM and SH3 protein 1 (LASP1) and phosphatidylinositol-5-phosphate 4-kinase, type II, beta (PIP4K2B) [48, 49].
In this study, we focused on traits that were measured during the growing period and are likely to play an important role in the adaptation of chickens to sub-optimal conditions, particularly to tropical climate. Since the lines used for the F2 design were original experimental lines, our findings may be specific to these genetic backgrounds and/or to the sub-optimal environmental conditions that differ from other studies. Yet, most of the QTL regions detected for growth traits overlapped with previously published QTL, which suggests that these QTL have an effect across a range of environmental conditions and could be particularly useful for selection. Several genome-wide QTL were detected for BW00, a trait that is strongly influenced by maternal effects. Since data on egg weight of the F1 dams were not available at the time of reproduction, this effect could not be separated from the additive dam effect on chick weight.
Surprisingly, no genome-wide significant QTL effects were detected on comb area, a trait that is likely relevant for climate adaptation and that differed significantly between the parental lines. Because of the large sexual dimorphism of this trait, we analyzed the data from F2 males and females separately to search for sex-specific QTL but none were detected. This may be explained by measurement errors, which can be quite important for both sexes, or by the rather large number of QTL that control variation of this trait, each with a small effect. None of the regions previously described for comb mass  were identified in the current study, which could be due to the use of lines with different genetic backgrounds in the current study.
Very little overlap was observed between the QTL regions detected by QTLMap and GEMMA. The differences between GWAS and QTL mapping-LDLA may be due to GWAS detecting associations based on LD across the F2 population, whereas QTL mapping-LDLA also exploits the information provided by within-family segregation, in order to estimate SNP effects on the trait variation. Legarra et al.  compared LDLA and EMMA methods (such as GEMMA) using real and simulated datasets and showed good agreement for the location of QTL between them. In general, these datasets consisted in a large number of small families, except for the sheep dataset, in which five F1 sires were mated to a large number of F1 dams, so that the number of half-sib families was similar to the design in the current study but the proportion of full-sibs was much lower in the sheep dataset. In our case, differences between the GWAS and LDLA results could be due to the fact that within-family segregation provides the most relevant information for QTL detection.
Although GEMMA and QTLMap identified the same QTL region on GGAZ, with strong effects on body weight and shank length, it was significant only in females with GEMMA, while QTLMap offered the possibility to include an interaction between sex and QTL. Such detection of sex-specific QTL is not new and they could contribute to our results, but in the present case, this result could also be due to the segregation mode of Z-linked alleles according to sex. In females, only two alleles are compared since females are hemizygous, so this is a 1:1 comparison of an R- allele to an L2 allele, whereas in males the within-family comparison will differ according to mating type. When F1 females have inherited their Z chromosome from an L2 F0 female (named ZL2), the F2 male progeny can have two possible genotypes: homozygous ZL2ZL2 or heterozygous ZL2ZR-, with some recombination occurring on the Z chromosome inherited from the F1 male. When F1 females inherit their Z chromosome from an R- F0 female (named ZR-), the F2 male progeny scan have two possible genotypes: homozygous ZR–ZR- or heterozygous ZL2ZR-, with recombination occurring on the Z chromosome inherited from the F1 male. Thus, if the QTL effect is not additive, the different F2 families will not provide the same contribution to estimation of the QTL, which may result in the absence of a significant association. Thus, GWAS involving sex-linked inheritance may require specific analysis.
Relevant candidate genes were identified in this study, according to their position and function. Most of them belong to networks of genes that are involved in embryonic, organism and tissue development, according to Ingenuity Pathway Analysis (IPA) . The sequence data available from previous studies provided a list of SNPs that represents only a subset of all SNPs that can differ between the lines. Some of the alleles at these SNPs are expected to be line-specific as a result of the long selection history of each line and the genetic distance between the lines. Since line-specific alleles may arise because of selection response or because of random drift, follow-up studies are necessary to validate whether or not these SNPs contribute to the QTL effect. Since missense mutations are rare, they could be the first choice for validation analyses, but it is also expected than non-coding variants can be quite relevant for quantitative variation. Thus, the choice of SNPs to be validated will be made according to their position on the gene, heterozygosity in the F1 population, and technical factors depending on the genotyping method.
For instance, PTGER4, a candidate gene for the growth QTL on GGAZ, has no known function in chicken but plays an important role in osteoclast differentiation and physiology in mice [52–54]. Although the role of this gene is not well known in chickens, it is a good candidate for further analysis, starting with the genotyping of F1 and F2 animals for SNPs chosen from the sequence data, followed by expression studies in each parental line if a significant effect on growth is identified in the F2 animals. Furthermore, genotyping candidate SNPs in the F2 population would also help to understand the different results obtained for F2 males and females. An obvious candidate gene to be tested for the QTL on GGAZ is the growth hormone receptor gene, which is known to carry mutations that cause sex-linked dwarfism in chickens. This gene lies close to the QTL region but not in it, so it may be worthwhile to include it in a confirmation study.
Our findings can also be useful for the future management of the L2 and R- lines. The frequency and the phenotypic consequences of the validated SNPs should be determined in the current breeding populations of the L2 and R- lines, to investigate their potential use in the management of these lines. Depending on their frequency and effect, some of these SNPs could be useful in the selection process for sub-tropical climate adaptation. Another option would be to set up a dual-purpose cross, taking advantage of the sex-linked QTL on growth. For instance, R- males could be crossed to L2 females to take advantage of heterosis in hybrids by producing light weight F1 females for egg production and F1 males with a greater body weight for local broiler production.
We have identified QTL for growth and morphological traits that may influence adaptation of chickens to varying environmental conditions. The availability of whole-genome sequence data for each parental line was useful to better document the candidate genes that were identified according to their positions and known functions. Finally, the very strong QTL effects found on the Z chromosome for body weight and shank length underlines the importance of sex-linked inheritance for growth-related traits in chickens, which is particularly relevant for crossbreeding in poultry breeding programs.
Abasht B, Dekkers JC, Lamont SJ. Review of quantitative trait loci identified in the chicken. Poult Sci. 2006;85:2079–96.
Wolc A, Arango J, Settar P, Fulton JE, O’Sullivan NP, Preisinger R, et al. Genome-wide association analysis and genetic architecture of egg weight and egg uniformity in layer chickens. Anim Genet. 2012;43(suppl 1):87–96.
Yuan J, Wang K, Yi G, Ma M, Dou T, Sun C, et al. Genome-wide association studies for feed intake and efficiency in two laying periods of chickens. Genet Sel Evol. 2015;47:82.
Wright D, Kerje S, Brändström H, Schütz K, Kindmark A, Andersson L, et al. The genetic architecture of a female sexual ornament. Evolution. 2008;62:86–98.
Lee YP, Yeh LT, Huang HH. A study on the production system of Taiwan Country chicken: growth traits of three-strain cross country chickens. J Chin Soc Anim Sci. 1997;26:271–84.
Chen CF, Huang NZ, Gourichon D, Lee YP, Tixier-Boichard M, Bordas A. Effect of introducing the naked neck gene in a line selected for low residual feed consumption on performance in temperate or subtropical environments. Poult Sci. 2008;87:1320–7.
Yalçin S, Testik A, Ozkan S, Settar P, Celen F, Cahaner A. Performance of naked neck and normal broilers in hot, warm, and temperate climates. Poult Sci. 1997;76:930–7.
Deeb N, Cahaner A. The effects of naked neck genotypes, ambient temperature, and feeding status and their interactions on body temperature and performance of broilers. Poult Sci. 1999;78:1341–6.
Yunis R, Cahaner A. The effects of the naked neck (Na) and frizzle (F) genes on growth and meat yield of broilers and their interactions with ambient temperatures and potential growth rate. Poult Sci. 1999;78:1347–52.
Nelson NA, Lakshmanan N, Lamont SJ. Sheep red blood cell and Brucella abortus antibody responses in chickens selected for multitrait immunocompetence. Poult Sci. 1995;74:1603–9.
Groenen MA, Megens HJ, Zare Y, Warren WC, Hillier LW, Crooijmans RP, Vereijken A, Okimoto R, Muir WM, Cheng HH. The development and characterization of a 60 K SNP chip for chicken. BMC Genomics. 2011;12:274.
Bordas A, Mérat P. Correlated responses in a selection experiment on residual feed intake of adult Rhode-Island Red cocks and hens. Ann Agri Fenn. 1984;23:233–7.
Carter GR. Studies on Pasteurella multocida. I. A hemagglutination test for the identification of serological types. Am J Vet Res. 1955;16:481–4.
Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM4, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 2015;4:7.
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.
Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 2012;44:821–4.
Nicodemus KK, Liu W, Chase GA, Tsai YY, Fallin MD. Comparison of type I error for multiple test corrections in large single-nucleotide polymorphism studies using principal components versus haplotype blocking algorithms. BMC Genet. 2005;6:S78.
Lander E, Kruglyak L. Genetic dissection of complex traits: guidelines for interpreting and reporting linkage results. Nat Genet. 1995;11:241–7.
R Development Core Team. R: A language and environment for statistical computing. 2008. http://www.R-project.org Accessed 1 May 2015.
Filangi O, Moreno C, Gilbert H, Legarra A, Le Roy P, Elsen JM. QTLMap, a software for QTL detection in outbred population. In Proceedings of the 9th world congress on genetics applied to livestock production: 1–6 August 2010, Leipzig; 2010.
Legarra A, Fernando RL. Linear models for joint association and linkage QTL mapping. Genet Sel Evol. 2009;41:43.
Ott J. Analysis of human genetic linkage. 3rd ed. London: John Hopkins University Press; 1999.
Churchill GA, Doerge RW. Empirical threshold values for quantitative trait mapping. Genetics. 1994;138:963–71.
National Center for Biotechnology Information. Rockville Pike: US National Library of Medicine; 2005. http://www.ncbi.nlm.nih.gov. Accessed 1 Feb 2015.
Yates A, Akanni W, Amode MR, Barrell D, Billis K, Carvalho-Silva D, Cummins C, et al. Ensembl 2016. Nucleic Acids Res. 2016;44:D710–6.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.
Fan WL, Ng CS, Chen CF, Lu MY, Chen YH, Liu CJ, et al. Genome-wide patterns of genetic variation in two domestic chickens. Genome Biol Evol. 2013;5:1376–92.
Frésard L, Leroux S, Servin B, Gourichon D, Dehais P, Cristobal MS, et al. Transcriptome-wide investigation of genomic imprinting in chicken. Nucleic Acids Res. 2014;42:3768–82.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25:1754–60.
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.
DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011;43:491–8.
Van der Auwera GA, Carneiro MO1, Hartl C, Poplin R, Del Angel G, Levy-Moonshine A et al. From FastQ data to high confidence variant calls: the genome analysis toolkit best practices pipeline. Curr Protoc Bioinformatics. 2013;43:11.10.1–33.
Picard. Cambridge: UK Broad Institute. https://broadinstitute.github.io/picard. Accessed 1 April 2016.
Cingolani P, Platts A, le Wang L, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012;6:80–92.
McLaren W, Pritchard B, Rios D, Chen Y, Flicek P, Cunningham F. Deriving the consequences of genomic variants with the Ensembl API and SNP effect predictor. Bioinformatics. 2010;26:2069–70.
Hu ZL, Park CA, Reecy JM. Developmental progress and current status of the Animal QTLdb. Nucleic Acids Res. 2016;44:D827–33.
Iskandar K, Cao Y, Hayashi Y, Nakata M, Takano E, Yada T, et al. PDK-1/FoxO1 pathway in POMC neurons regulates Pomc expression and food intake. Am J Physiol Endocrinol Metab. 2010;298:E787–98.
Eppig JT, Blake JA, Bult CJ, Kadin JA, Richardson JE, The Mouse Genome Database Group. The Mouse Genome Database (MGD): facilitating mouse as a model for human biology and disease. Nucleic Acids Res. 2015;2015(43):D726–36.
Ishibashi S, Schwarz M, Frykman PK, Herz J, Russell DW. Disruption of cholesterol 7alpha-hydroxylase gene in mice. I. Postnatal lethality reversed by bile acid and vitamin supplementation. J Biol Chem. 1996;271:18017–23.
Shmerling D, Danzer CP, Mao X, Boisclair J, Haffner M, Lemaistre M, et al. Strong and ubiquitous expression of transgenes targeted into the beta-actin locus by Cre/lox cassette replacement. Genesis. 2005;42:229–35.
Föller M, Sopjani M, Koka S, Gu S, Mahmud H, Wang K, et al. Regulation of erythrocyte survival by AMP-activated protein kinase. FASEB J. 2009;23:1072–80.
Cybulski N, Polak P, Auwerx J, Rüegg MA, Hall MN. mTOR complex 2 in adipose tissue negatively controls whole-body growth. Proc Natl Acad Sci USA. 2009;106:9902–7.
Nguyen M, Camenisch T, Snouwaert JN, Hicks E, Coffman TM, Anderson PA, et al. The prostaglandin receptor EP4 triggers remodelling of the cardiovascular system at birth. Nature. 1997;390:78–81.
Klar J, Blomstrand P, Brunmark C, Badhai J, Håkansson HF, Brange CS, et al. Fibroblast growth factor 10 haploinsufficiency causes chronic obstructive pulmonary disease. J Med Genet. 2011;48:705–9.
Lee PC, Dodart JC, Aron L, Finley LW, Bronson RT, Haigis MC, et al. Altered social behavior and neuronal development in mice lacking the Uba6-Use1 ubiquitin transfer system. Mol Cell. 2013;50:172–84.
Bard-Chapeau EA, Szumska D, Jacob B, Chua BQ, Chatterjee GC, Zhang Y, et al. Mice carrying a hypomorphic Evi1 allele are embryonic viable but exhibit severe congenital heart defects. PLoS One. 2014;9:e89397.
Lamia KA, Peroni OD, Kim YB, Rameh LE, Kahn BB, Cantley LC. Increased insulin sensitivity and reduced adiposity in phosphatidylinositol 5-phosphate 4-kinase beta-/- mice. Mol Cell Biol. 2004;24:5080–7.
Chew CS, Chen X, Bollag RJ, Isales C, Ding KH, Zhang H. Targeted disruption of the Lasp-1 gene is linked to increases in histamine-stimulated gastric HCl secretion. Am J Physiol Gastrointest Liver Physiol. 2008;295:G37–44.
Legarra A, Croiseau P, Sanchez MP, Teyssèdre S, Sallé G, Allais S, et al. A comparison of methods for whole-genome QTL mapping using dense markers in four livestock species. Genet Sel Evol. 2015;47:6.
IPA®, CA QIAGEN, Redwood City. www.qiagen.com/ingenuity. Accessed 1 Jun 2016.
Li M, Healy DR, Li Y, Simmons HA, Crawford DT, Ke HZ, et al. Osteopenia and impaired fracture healing in aged EP4 receptor knockout mice. Bone. 2005;37:46–54.
Choudhary S, Blackwell K, Voznesensky O, Deb Roy A, Pilbeam C. Prostaglandin E2 acts via bone marrow macrophages to block PTH-stimulated osteoblast differentiation in vitro. Bone. 2013;56:31–41.
Gao Q, Zhan P, Alander CB, Kream BE, Hao C, Breyer MD, et al. Effects of global or targeted deletion of the EP4 receptor on the response of osteoblasts to prostaglandin in vitro and on bone histomorphometry in aged mice. Bone. 2009;45:98–103.
CFC, MTB, and CYL managed the project. CYL carried out the variance analysis, QTL mapping, and GWAS analysis. CYL and MTB drafted the manuscript and the interpretation of the data. WFW, CSN, and CYL performed the sequence analysis. SWW and CFC managed the chicken production and phenotype collection. CYL, SWW, and CFC supervised the genotyping. CFC and MTB contributed to the funding and design of the experiment. All authors read and approved the final manuscript.
CYL was supported by a fellowship from the French Institute of Taipei, Ministry of Science and Technology, and Ministry of Education in Taiwan. This work was part of a co-supervised PhD associating NCHU and AgroParisTech. The staff of the experimental farm of the NCHU is gratefully acknowledged. Authors are also grateful to Tatiana Zerjal for her help in IPA, to Bing-Yen Tsai for his help to combine the information of SNP and gene from different resources, and to the members of the PhD committee, Pascale Le Roy, Catherine Larzul, Marie-Hélène Pinard-van der Laan, and Xavier Rognon, for their advice. This work was dedicated to André Bordas, INRA, who selected the R- line and organized the shipment of a subset of the line to NCHU in 2003 and Yen-Pai Lee, who preserved and maintained the L2 line in NCHU.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
All animals used in this study were processed following the approved protocol of the Institutional Animal Care and Use Committees of the NCHU (Taichung, Taiwan; IACUC No. 97-99).
This study was supported by a grant of the Ministry of Science and Technology (Grant Number NSC 99-2321-B-005-009-MY3), Ministry of Education, and French Institute of Taipei, Taiwan.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
- Quantitative Trait Locus
- Quantitative Trait Locus Mapping
- Quantitative Trait Locus Region
- Quantitative Trait Locus Effect
- Quantitative Trait Locus Detection