Genomic Region Containing Toll-Like Receptor Genes Has a Major Impact on Total IgM Antibodies Including KLH-Binding IgM Natural Antibodies in Chickens

Abstract

Natural antibodies (NAb) are antigen binding antibodies present in individuals without a previous exposure to this antigen. Keyhole limpet hemocyanin (KLH)-binding NAb levels were previously associated with survival in chickens. This suggests that selective breeding for KLH-binding NAb may increase survival by means of improved general disease resistance. Genome-wide association studies (GWAS) were performed to identify genes underlying genetic variation in NAb levels. The studied population consisted of 1,628 adolescent layer chickens with observations for titers of KLH-binding NAb of the isotypes IgM, IgA, IgG, the total KLH-binding (IgT) NAb titers, total antibody concentrations of the isotypes IgM, IgA, IgG, and the total antibodies concentration in plasma. GWAS were performed using 57,636 single-nucleotide polymorphisms (SNP). One chromosomal region on chromosome 4 was associated with KLH-binding IgT NAb, and total IgM concentration, and especially with KLH-binding IgM NAb. The region of interest was fine mapped by imputing the region of the study population to whole genome sequence, and subsequently performing an association study using the imputed sequence variants. 16 candidate genes were identified, of which FAM114A1, Toll-like receptor 1 family member B (TLR1B), TLR1A, Krüppel-like factor 3 (KLF3) showed the strongest associations. SNP located in coding regions of the candidate genes were checked for predicted changes in protein functioning. One SNP (at 69,965,939 base pairs) received the maximum impact score from two independent prediction tools, which makes this SNP the most likely causal variant. This SNP is located in TLR1A, which suggests a fundamental role of TLR1A on regulation of IgM levels (i.e., KLH-binding IgM NAb, and total IgM concentration), or B cells biology, or both. This study contributes to increased understanding of (genetic) regulation of KLH-binding NAb levels, and total antibody concentrations.

Introduction

Modern poultry production is facing high impact changes in production systems and management (1, 2). These changes may enhance the risks of infectious disease, because of higher pathogenic pressure. According to FAO statistics, economic losses in poultry production due to diseases were as high as 10–20% of the gross production value in developed countries in 2014.¹ Only few years ago, diseases were controlled by (preventive and abundant) use of antibiotics, but this resulted in an increased risk for antibiotic resistance. New sustainable strategies for obtaining robust populations need to be developed, in addition to vaccination strategies, and feed formulations. Selective breeding for an improved immune system might be an important additional strategy to improve general disease resistance (3).

Selective breeding for an improved immune system requires a trait that is heritable, easy to measure, and preferentially related to general disease resistance (4). Natural antibodies (NAb) might be a good candidate trait. NAb are defined as antigen-binding antibodies present in individuals without a previous exposure to this antigen (5). NAb play an essential role in the innate, and adaptive immunity against self-antigens, and various types of pathogens, and are contributing to healthy aging, and disease resistance (6–8). Different NAb isotypes are found: predominantly IgM, but also IgA, and IgG (8, 9). In earlier studies, high NAb titers (especially IgM) binding the overt antigen keyhole limpet hemocyanin (KLH) at adolescence, but not at later age, were related to lower mortality of commercial layer chickens (10–12). KLH-binding NAb levels at adolescence were previously estimated to have heritabilities between 0.07 and 0.14 in white layer chickens (4). In addition, some genomic regions underlying genetic variation of NAb titers were identified by using dedicated single-nucleotide polymorphism (SNP) sets with 1,022 SNP (13, 14). However, to the best of our knowledge, no genome-wide association studies (GWAS) for NAb levels in chicken have been reported.

This study describes the first GWAS on NAb in chicken. The studied population consisted of 1,628 adolescent purebred White Leghorn chickens with observations for KLH-binding NAb titers of the isotypes IgM, IgA, IgG, and total KLH-binding NAb titers (IgT). In addition, phenotypes on total antibody concentrations of the isotypes IgM [total IgM (tIgM)], IgA (tIgA), IgG (tIgG), and total antibody (tIgT) concentration were collected to investigate similarities/differences between KLH-binding NAb, and total antibody concentration. Identification of associated genomic regions was performed by using 57,636 SNP in a single SNP GWAS. Regions were fine mapped by performing an association study using imputed genotypes based on sequence data with the aim to identify candidate gene(s), and possibly the causal variant(s). The role of the found association(s) in antibody-related immunity in chickens will be discussed.

Materials and Methods

Ethics Approval Statement

Collection of samples and data was done according to Hendrix Genetics (HG) protocols, under the supervision of HG employees. Samples and data were collected as part of routine data collection in a commercial breeding program for layer chickens in The Netherlands. Samples and data were collected on a breeding nucleus of HG for breeding purposes only and are a non-experimental, agricultural practice, regulated by the Act Animals, and the Royal Decree on Procedures. The Dutch Experiments on Animals Act does not apply to non-experimental, agricultural practices. An ethical review by the Statement Animal Experiment Committee was therefore not required. No extra discomfort was caused for sample collection for the purpose of this study.

Study Population

Samples and data were obtained from a purebred White Leghorn chicken line (in other work referred to as “WA”), which is a layer chicken line selected mainly for egg production, but also for other production traits, e.g., traits related to egg quality. The studied chicken population comprised 1,628 chickens, of which 696 males and 932 females. The study population originated from 112 sires and 288 dams. 437 chickens in the current study were also included in the study of Berghof et al. (4). Of these 437 chickens, 222 female chickens were also included in the study of van der Klein et al. (15). The chickens were kept according to standard management of breeding nucleus farms of HG. Further details can be found in van der Klein et al. (15).

The chickens received obligatory vaccinations against Marek’s disease [1 day of age intramuscular (i.m.)], infectious bronchitis (1 day of age, 12–14 days of age, 10 weeks of age, 12 weeks of age via spray; 16 weeks of age i.m.), Newcastle disease (13 day of age, 42 days of age, 12 weeks of age via spray; 16 weeks of age i.m.), infectious bursal disease (25 days of age via spray; 16 weeks of age i.m.), chicken anemia virus (16 weeks of age via water), fowl pox (16 weeks of age by wing web injection), and avian encephalomyelitis (16 weeks of age by wing web injection).

Blood of the study population was collected once between 15 and 22 weeks of age (i.e., adolescence), without anesthesia/analgesia. Samples of the study population were collected in four batches with approximately 1.75 years between the first and last batch. No chickens were killed for sample collection. The blood samples were centrifuged, and plasmas and blood cells were collected separately, and stored at −20°C until use.

KLH-Binding NAb Titers

Optical density (OD) of KLH-binding immunoglobulins of the isotypes IgM, IgA, and IgG, and the total KLH-binding (IgT) immunoglobulins were determined in individual plasma samples by an indirect two-step ELISA as described by van der Klein et al. (15), and Berghof et al. (4). Briefly, plasma samples were 1:10 prediluted (for IgM, IgG, and IgT analyses) or were 1:5 prediluted (for IgA analysis) with dilution buffer [PBS (10.26 g/L Na₂HPO₄⋅H₂O, 2.36 g/L KH₂PO₄, and 4.50 g/L NaCl; pH 7.2) containing 0.5% normal horse serum, and 0.05% Tween^® 20]. Predilutions were stored at 4°C until use the next day, or were stored at −20°C until use (with a maximum storage time of 3 months). Flat-bottomed, 96-well medium binding plates (Greiner Bio-One, Alphen a/d Rijn, The Netherlands) were coated with 2 μg/mL KLH in 100 μL coating buffer (5.3 g/L Na₂CO₃, and 4.2 g/L NaHCO₃; pH 9.6) per well, and incubated at 4°C overnight. After washing for 6 s with tap water containing Tween^® 20, plates were tapped dry. The 1:10 predilution of the samples were further diluted in the KLH-coated plates with dilution buffer to 1:40, 1:160, 1:640, and 1:2,560 test dilutions for IgM, IgG, and IgT, 1:10, 1:20, 1:40, and 1:80 for IgA. Duplicate standard positive plasma samples (a pool of male plasmas) were stepwise 1:1 diluted with dilution buffer, and pipetted into the KLH-coated plates. A minimum of five samples per plate was maintained to ensure proper adjustment of the titers for plate effects in the statistical analyses. The plates were incubated for 1.5 h at room temperature (20–25°C). After washing, plates were incubated with 1:20,000-diluted goat-antichicken IgM labeled with peroxidase (PO) (Cat# A30-102P, RRID:AB_66857),² 1:7,500-diluted goat-antichicken IgA labeled with PO (Cat# A30-103P, RRID:AB_66833) (see text footnote 2), 1:40,000-diluted goat-antichicken IgG(Fc) labeled with PO (Cat# A30-104P, RRID:AB_66843) (see text footnote 2), or 1:20,000-diluted rabbit-antichicken IgG heavy and light chain (IgT) labeled with horse radish PO (Cat# A30-107P, RRID:AB_67386) (see text footnote 2), (all polyclonal antibodies from Bethyl Laboratories, Montgomery, TX, USA), and incubated for 1.5 h at room temperature. After washing, binding of the antibodies to KLH was visualized by adding 100 μL substrate buffer [containing reverse osmosis purified water, 10% tetramethylbenzidine buffer (15.0 g/L sodium acetate and 1.43 g/L urea hydrogen peroxide; pH 5.5), and 1% tetramethylbenzidine (8 g/L TMB in DMSO)] at room temperature. After approximately 15 min, the reaction was stopped with 50 μL of 1.25 M H₂SO₄. OD was measured with a Multiskan Go (Thermo scientific, Breda, The Netherlands) at 450 nm.

Antibody titers were calculated as described by Frankena (16) [taken from De Koning et al. (17)]. Briefly, OD of the duplicate standard positive plasma samples were averaged for each plate. Logit values of OD per plate were calculated using: where OD is the OD of a well, and OD_max is the maximum averaged OD of the duplicate standard positive plasma samples. The last positive well (lpw) of the averaged duplicate standard positive plasma sample was set to the sixth dilution. A linear regression line of the logit OD against the respective log₂-dilution values of the averaged duplicate standard positive plasma samples was determined, which resulted in regression coefficient β. Titers of plasma samples per plate were calculated using: were logit OD_lpw is the estimated logit OD at the lpw calculated with the estimated linear regression function using the log₂-dilution value of that well, logit OD_sample is the logit OD calculated of the OD closest to 50% of OD_max for a plasma sample of an individual (OD_sample), β is the regression coefficient of the estimated linear regression function of the averaged duplicate standard positive plasma samples, and log₂ (dilution_sample) is the log₂-dilution value at which OD_sample occurred, as described by De Koning et al. (17).

The total number of observations was 1,627 for IgM NAb, 1,608 for IgA NAb, 1,623 for IgG NAb, and 1,625 for IgT NAb (see Table 1).

Total Antibody Concentrations

Total concentrations (μg/ml) of the immunoglobulin isotypes IgM (tIgM), IgA (tIgA), and IgG (tIgG) were determined in individual plasma samples by an indirect sandwich ELISA according to manufacturer’s protocols (IgM: Cat# E30-102; IgA: Cat# E30-103; IgG: Cat# E30-104; all from Bethyl Laboratories) with minor additions: (1) plasma samples were 1:5 prediluted, and step-wise diluted to 1:500 (for tIgA), 1:5,000 (for tIgM), and 1:50,000 (for tIgG) with sample/conjugate diluent. (2) Flat-bottomed, 96-well medium binding plates (Greiner Bio-One) were used for ELISA. (3) Washing procedure consisted of washing each well two times with 200 μl and subsequently three times with 100 μl wash solution. Plates were emptied, and tapped dry in between, and at the end. (4) 100 μl ELISA blocking solution was added to each well. (5) A minimum of five samples per plate was maintained to ensure proper adjustment of antibody concentrations for plate effects in the statistical analyses. (6) A duplicate of a standard positive plasma sample (obtained from Bethyl Laboratories) was used per plate. (7) Enzyme substrate (TMB), and ELISA stop solution were made and applied as specified for “KLH-binding natural antibody titers” (see above). (8) OD were measured with a Multiskan Go. (9) All samples on a plate were corrected for the average background OD of the two blanks. (10) SoftMax^® Pro 7.0 build 226962 (Free Trial) was used for generating a 4-parameter logistic curve fit based on the averages of the two standard ranges, and calculating isotype concentrations per sample.

The total concentration (μg/ml) of antibodies (tIgT) was calculated by summing the concentrations for tIgM, tIgA, and tIgG. tIgM, tIgA, tIgG, and tIgT were log₁₀-tranformed to normalize the data for statistical analyses. The total number of observations was 1,619 for tIgM, 1,586 for tIgA, 1,615 for tIgG, and 1,573 for tIgT (see Table 1).

Genotypes

DNA was extracted from the collected blood cells. The study population was genotyped for a 2,740 SNP set (n = 488) (Illumina, San Diego, CA, USA), or for a 52,232 SNP set from which 11,173 SNP were used (n = 1,140) (Illumina). Both sets were imputed with Beagle 4.0 (18) (accuracy ≥ 97%) to a 57,636 SNP set, based on approximately 120 key ancestors of this chicken line. The SNP distribution over the chicken chromosomes (based on Gallus_gallus-5.0) can be found in Table S1 in Supplementary Material.

Quality control was applied in four steps: (1) All monomorphic SNP were removed. (2) SNP with at least one genotype class with less than nine chickens with phenotypic observations (~0.5% of the population) were removed. (3) SNP on autosomes with a strong deviation from Hardy Weinberg Equilibrium were removed. (4) SNP on the allosome Z with a strong deviation from Hardy Weinberg Equilibrium in males were removed for the whole study population. These criteria were chosen to exclude SNP which are likely to have a high rate of genotyping error and to exclude very low frequency SNP from the data set (19).

Genetic Parameters

The additive genetic variances (), and the maternal environmental variance () for KLH-binding NAb titers were previously estimated in a similar, but larger population of this chicken line (4) (see Table 1). The additive genetic variances (), and heritabilities for tIgM, tIgA, tIgG, and tIgT concentrations were estimated with model 1, in a similar approach as for the NAb titers described by Berghof et al. (4).

Model 1 was where y_ij is the tIgM concentration, tIgA concentration, tIgG concentration, or tIgT concentration, μ is the overall mean, P_i is the fixed effect of plate on which a sample was analyzed (i = 1–26 for tIgM and tIgA, or i = 1–27 for tIgG and tIgT), Age_ij is the covariate describing the effect of age at sampling (15–22 weeks of age) with regression coefficient β₁, id_j is the random additive genetic effect of the jth chicken assumed to be distributed as , and e_ij is the residual term assumed to be distributed as . Assumed (co)variance structures of the random model terms are and , in which A is the additive genetic relationship matrix, is the additive genetic variance, I is an identity matrix, and is the residual variance. The pedigree used to construct A consisted of 2,537 individuals, and was based on a minimum of three generations of ancestors for each individual in the study population. The estimated plate effect also corrects for (confounded) effects on the samples, such as sex and storage. The estimated regression coefficient of age also corrects (partly) for (confounded) effects of vaccinations at 16 weeks of age.

Model 1 extended with a random maternal environmental effect [for IgM NAb only, see Berghof et al. (4)], were also used to estimate genetic, and phenotypic correlations with bivariate analyses, as described by Berghof et al. (4). The genetic correlations give estimates on the extent to which two traits are influenced by the same genetic variations (i.e., genes). The phenotypic correlations give estimates on the extent to which two traits are influenced by the same genetic variation, and the same environmental variation (e.g., position of cage in the stable, temperature, age). See also Stear et al. (20) for a detailed explanation.

The statistical analyses were performed using ASReml^® 4.1 (21).

Genome-Wide Association Studies

A single-SNP linear animal model was used for studying KLH-binding IgM, IgA, IgG, and IgT NAb titers and tIgM, tIgA, tIgG, and tIgT concentrations associations.

All traits, except IgM NAb, were analyzed using model 2a: and IgM NAb was analyzed using model 2b: where y_ijk is the IgA NAb titer, IgG NAb titer, IgT NAb titer, tIgM concentration, tIgA concentration, tIgG concentration, or tIgT concentration, y_ijkl is the IgM NAb titer, μ is the overall mean, P_i is the fixed effect of plate on which a sample was analyzed (i = 1–91 for IgM NAb, IgG NAb, and IgT NAb; i = 1–100 for IgA NAb; i = 1–26 for tIgM; and tIgA, or i = 1–27 for tIgG, and tIgT), Age_ijk and Age_ijkl are the covariates describing the effect of age at sampling (15–22 weeks of age) with regression coefficient β₁, SNP_j is the fixed effect of genotype class j (j = AA, AB, or BB) of the analyzed SNP, id_k is the random additive genetic effect of the kth chicken assumed to be distributed as , d_l is the random effect of the lth dam (IgM NAb only) assumed to be distributed as , and e_ijk and e_ijkl are the residual terms assumed to be distributed as . Assumed (co)variance structures of the random model terms are , and (for all), and (for IgM NAb only), in which A is the additive genetic relationship matrix, is the additive genetic variance, I are identity matrices, is the residual variance, and is the maternal environmental variance (for IgM NAb only).

The variances for random factors were fixed to values obtained from analyses from models without a SNP effect: the additive genetic variances () and the maternal environmental variance () for NAb titers were fixed to values estimated by Berghof et al. (4) in a similar population of this chicken line. These estimates are more accurate, because the studied population was bigger. The additive genetic variances () for tIgM, tIgA, tIgG, and tIgT, concentrations were estimated with model 1.

The pedigree used to construct A consisted of 2,537 individuals, and was based on a minimum of three generations of ancestors for each individual in the study population. The estimated plate effect also corrects for (confounded) effects on the samples, such as sex and storage.

The statistical analyses were performed using ASReml^® 4.1 (21).

p-values of SNP effects were checked for inflation (22), which was expressed as inflation factor λ (23). λ was estimated using the “estlambda()” function of the R package “GenABEL” (24). Genomic control was applied when λ ≥ 1.1 [Wellcome Trust Case Control Consortium (25), taken from Duijvesteijn et al. (26)] by dividing the test statistics (F-values) by λ, and subsequently recalculating the p-values. A genome-wide false discovery rate (FDR) was estimated based on observed p-values using the R package “qvalue” (27). A FDR ≤ 0.05 was used to indicate significant associations, and a FDR ≤ 0.20 was used to indicate suggestive associations. The lead SNP [SNP with highest −log₁₀(p-value)] variance as a percentage of the additive genetic variance within an associated region was estimated as , where was calculated based on estimated SNP effects from model 2a for all, except IgM NAb, or model 2b for IgM NAb and the genotypic frequencies. Significant differences between genotype classes were tested with a T-test using the “!TDIFF” statement in ASReml^® 4.1 (21). A p-value ≤0.05 was considered to be significant, and a p-value ≤0.10 was considered to be suggestive.

Fine Mapping, and Identification of Possible Causal Variant

The chickens were imputated to full genome sequence in two steps based on the full genome sequences of 70 key ancestors of this chicken line. Imputation was done with Beagle 4.0 (18) (accuracy ≥ 88%) by imputing the 57,636 SNP set (for the 2,740 SNP set), or the 52,232 SNP set to full genome sequence with an average sequence depth of 12.4 (and SD of 2.1).

Associated genomic regions were selected for fine mapping, if at least 1 SNP was significant (FDR ≤ 0.05). Fine mapping regions were defined based on the last significant SNP on the border of significantly associated regions. The fine mapping regions consisted of the first SNP outside the significantly associated regions −5 and +5 Mbp.

Quality control of the imputed SNP that were used in the fine mapping studies, was applied in four steps as described above. Statistical analyses for the selected regions were performed as described above. Genomic control was applied using λ values from previous analyses, to be able to compare between analyses.

Single-nucleotide polymorphism of candidate genes in the fine mapped regions were checked for amino acid substitutions and predicted changes in protein functioning as a result of the nucleotide change in coding regions by using the Ensembl Variant Effect Predictor (28), including the integrated SIFT option (29), and by using PolyPhen-2 (30).

Results

In total 1,628 chickens, of which 696 males, and 932 females, were in the study population. Descriptive statistics and genetic parameters are shown in Table 1. Mean titers for KLH-binding NAb were 5.8 for IgM, 5.4 for IgA, 5.5 for IgG, and 5.6 for total KLH-binding NAb (IgT). Mean concentrations for total antibody concentrations were 355 μg/ml for tIgM, 314 μg/ml for tIgA, 7,856 μg/ml for tIgG, and 8,484 μg/ml for total antibody concentration (tIgT). Heritabilities for total antibody concentrations were estimated to be 0.23 for tIgM, 0.22 for tIgA, 0.06 for tIgG, and 0.08 for tIgT. Maternal environmental effects were not significant for any of the total antibody concentrations. Age was not significant for tIgA, and tIgT but was kept in the model to be consistent with the other analyses. Heritabilities and maternal environmental variances () for NAb were previously estimated in a similar, but larger population of this chicken line (see Table 1) (4).

Genetic, and phenotypic correlations between the corresponding types of KLH-binding NAb titers and total antibody concentrations are shown in Table 2. Genetic correlations were: 0.91 for IgM NAb/tIgM, 0.38 for IgA NAb/tIgA, −0.61 for IgG NAb/tIgG, and −0.27 for IgT NAb/tIgT. Phenotypic correlations were: 0.41 for IgM NAb/tIgM, 0.26 for IgA NAb/tIgA, 0.08 for IgG NAb/tIgG, and 0.03 for IgT NAb/tIgT. In addition, genetic and phenotypic correlations between all total antibody concentrations were estimated (see Table S2 in Supplementary Material).

Quality control of the SNP resulted in removal of 37,053 SNP, because they were monomorphic (step 1), and an additional removal of 5,005 SNP, because of a low number of observations per genotype class or a strong deviation from Hardy-Weinberg equilibrium (steps 2–4) (see Table S1 in Supplementary Material). In total, 15,579 SNP were used for GWAS, though the actual number of used SNP per trait varied between 15,431 and 15,578 SNP due to missing phenotypes. Figure 1 shows the “Manhattan plots” of the GWAS. Results shown for KLH-binding IgM NAb, and all total antibody concentrations (tIgM, tIgA, tIgG, and tIgT) were adjusted for inflation factor λ, which was estimated to be 1.67 for IgM NAb, 1.23 for tIgM, 1.27 for tIgA, 1.13 for tIgG, and 1.14 for tIgT. This indicates a deviation of the p-values from their expected distribution under the null hypotheses of no genetic association. After correcting the test statistics, and subsequently recalculating the p-values, all new λ were below 1. Inflation factors for IgA NAb, IgG NAb, and IgT NAb were below 1.1, and required no adjustment. For KLH-binding IgA NAb, and KLH-binding IgT NAb, none of the SNP were significantly associated (FDR ≤ 0.05), but some regions showed suggestive associations (FDR ≤ 0.20) (see Table S3 in Supplementary Material). For KLH-binding IgM NAb, and tIgM several SNP were significantly associated (see Table S3 in Supplementary Material). For KLH-binding IgG NAb, tIgA, tIgG, and tIgT no significant or suggestive associations were detected.

Figure 1

Tables 3 and 4 characterize the significantly, and suggestively associated genomic regions based on the GWAS results. The associated regions, and additional information can be found in Table 3. The lead SNP and additional information can be found in Table 4. Table S3 in Supplementary Material gives an overview of associated genotyped and imputed SNP.

One genomic region showed an association with KLH-binding IgT NAb, and tIgM and showed the strongest association for KLH-binding IgM NAb. The genomic region was located on chromosome 4 around 70 M base pairs (bp), and consisted of 3 suggestive SNP for IgT NAb, 4 significant, and 1 suggestive SNP for tIgM, and 35 significant SNP for IgM NAb. The genomic region on chromosome 4 had its lead SNP for IgT NAb on 69,587,709 bp (rs313004783), for tIgM on 69,814,286 bp (rs313437715), for IgM NAb on 69,702,481 bp (rs15614874). Remarkably, the heterozygous genotype class was not significantly different from one of the homozygous genotype classes, indicating full dominance. The lead SNP variance as a percentage of the additive genetic variance was approximately 16.0% for KLH-binding IgT NAb, 14.1% for tIgM, and 57.6% for KLH-binding IgM NAb, illustrating that the genomic region on chromosome 4 has a major effect on the genetic variation of KLH-binding IgT NAb, tIgM, and especially on KLH-binding IgM NAb.

Two genomic regions were detected for KLH-binding IgA NAb. The first region was located on chromosome 9 around 12 Mbp and consisted of two suggestive SNP. The lead SNP was located at 12,261,658 bp (rs15969591). The lead SNP variance was approximately 13.5% of the additive genetic variance for IgA NAb. The second region was located on chromosome 18 around 10 Mbp. This region consisted of three suggestive neighboring SNP with equal p-values. The region showed complete dominance. The variance of the region was approximately 0.2% of the additive genetic variance for IgA NAb.

The significantly associated region on chromosome 4 was imputed to whole genome sequence. The region selected for further investigation ranged from approximately 64.7to 74.8 Mbp. After quality control 43,675 SNP remained for KLH-binding IgM NAb, and 39,792 SNP for tIgM, deviating because of missing phenotypes for tIgM. Figure 2 shows the results of the association study for IgM NAb, and tIgM of the selected region on chromosome 4. The highest signal for IgM NAb, and tIgM was detected in the region 69.5–70.2 Mbp with 3,153 SNP (see Figure 2C). Candidate genes (5′–3′ direction) in this region are: PDS5 cohesin-associated factor A (PDS5A), ubiquitin conjugating enzyme E2 K (UBE2K), small integral membrane protein 14 (SMIM14), UDP-glucose 6-dehydrogenase (UGDH), lipoic acid synthetase (LIAS), ribosomal protein L9 (RPL9), klotho beta (KLB), WD repeat domain 19 (WDR19), replication factor C subunit 1 (RFC1), kelch like family member 5 (KLHL5), transmembrane protein 156 (TMEM156), family with sequence similarity 114 member A1 (FAM114A1), Toll-like receptor 1 family member B (TLR1B), Toll-like receptor 1 family member A (TLR1A), Kruppel like factor 3 (KLF3), and TBC1 domain family member 1 (TBC1D1) [from NCBI (31), and Ensembl (Release 87) (32)] (see Figure 2C). The strongest associated subregion contained FAM114A1, TLR1B, TLR1A, and KLF3.

Figure 2

32). Positions are based on Gallus_gallus-5.0. Arrows (“ <” and “ >”) indicate direction of the genes. Abbreviations genes: PDS5A, PDS5 cohesin-associated factor A; UBE2K, ubiquitin conjugating enzyme E2 K; SMIM14, small integral membrane protein 14; UGDH, UDP-glucose 6-dehydrogenase; LIAS, lipoic acid synthetase; RPL9, ribosomal protein L9; KLB, klotho beta; WDR19, WD repeat domain 19; RFC1, replication factor C subunit 1; KLHL5, kelch like family member 5; TMEM156, transmembrane protein 156; FAM114A1, family with sequence similarity 114 member A1; TLR1B, Toll-like receptor 1 family member B; TLR1A, Toll-like receptor 1 family member A; KLF3, Kruppel-like factor 3; TBC1D1, TBC1 domain family member 1 [from NCBI (31) and Ensembl (Release 87) (32)].

The 3,153 SNP in this region were checked for amino acid substitutions in coding regions, and consequential changes in protein functioning as a result of the nucleotide change. The Ensembl Variant Effect Predictor predicted 28 SNP to result in different amino acids. All 28 SNP could be predicted with SIFT, PolyPhen-2 or both. SIFT predicted amino acid substitutions to be “tolerated” for 11 SNP, and “deleterious” for 2 SNP. The remaining 15 SNP could not be predicted. PolyPhen-2 predicted amino acid substitution to be “benign” for 22 SNP, “possibly damaging” for 1 SNP, and “probably damaging” for 1 SNP. SIFT and PolyPhen-2 gave both the maximum impact score (“deleterious,” and “probably damaging”) to 1 SNP (69,965,939 bp), which is considered the most likely candidate as causal variant for the major effects on KLH-binding IgT NAb, tIgM, and especially on KLH-binding IgM NAb. This SNP consists of a cytosine (C)/guanine (G) polymorphism, and results in a phenylalanine (F)/leucine (L) amino acid substitution in TLR1A at protein position 126. The C-variant had an allele frequency of 0.45, and the G-variant had an allele frequency of 0.55. Based on Figure 1 from Keestra et al. (33), the SNP is located in leucine-rich repeat (LRR) 4.

Significance and effects of the predicted causal variant on all measured traits can be found in Table 5. KLH-binding IgM NAb, KLH-binding IgG NAb, KLH-binding IgT NAb, and tIgM showed a significant effect of the predicted causal variant. For these traits (except IgG NAb), the GG genotype class had significantly lower phenotypic values than the CG and CC genotype classes. The CG and CC genotype classes were not significantly different from each other, showing full dominance. The estimated genetic variance explained by the predicted causal variant was between 7.3 and 15.8% for KLH-binding IgG NAb, KLH-binding IgT NAb, and tIgM, and was 63.5% for KLH-binding IgM NAb. KLH-binding IgA NAb, tIgA, tIgG, and tIgT were not significantly associated.

Discussion

Relatively few GWAS have focused on detecting genomic regions associated with immunity in chicken (3, 34, 35). Though two KLH-binding NAb association studies with dedicated SNP sets were conducted across chicken lines (13, 14). To the best of our knowledge, this study is the first genome-wide association study on (KLH-binding) NAb levels in chickens. None of previously associated genes (13, 14) could be confirmed in this study, likely because of the high genetic uniformity in our data set as a result of a within line study instead of an across line study. No other study found the associated region on chromosome 4 reported in this study, which might indicate that this association is unique for this chicken line. This is also the first study that investigated heritabilities of total antibody concentrations or performed GWAS on all total antibody isotype concentrations in chickens.

57,636 SNP were used to perform the GWAS. 73.0% of the SNP were non-informative, mainly because the SNP were fixed (62.2%). This is in agreement with previously reported genetic uniformity in White Leghorn layer chicken lines (36). In total 15,580 SNP passed quality control. For five out of eight analyzed traits population stratification was found to be present, which suggests family structures within the studied population that are not captured by the pedigree information. Other factors, like heterogeneous variance, or model misspecification, likely also contributed to this, since inflation factors were not of equal size among the traits. Especially for KLH-binding IgM NAb, the population stratification was high. After correcting the test statistics for the inflation factor λ, the new inflation factor indicated that KLH-binding IgM NAb was overcorrected (λ = 0.88). Thereby being possibly too conservative on the associations with KLH-binding IgM NAb.

The heritability of tIgM was low, but twice the heritability of KLH-binding IgM NAb (4). KLH-binding IgM NAb have been found to be influenced by maternal environmental effects (4). In case maternal effects are not accounted for in genetic analyses, heritabilities will be overestimated as was seen for KLH-binding IgM NAb (4). The size of the present data set did not allow detection of maternal environmental effect for KLH-binding IgM NAb. The estimated heritability (without maternal environmental effect) for IgM NAb in this study population was 0.31 (data not shown), which is similar to the heritability of tIgM. In addition, the genetic correlation between KLH-binding IgM, and tIgM was very high. This suggests that maternal environmental effects for tIgM are present, but could not be estimated, and resulted in overestimation of tIgM heritability.

One significant genomic region was identified for KLH-binding IgM NAb, and tIgM on chromosome 4 around 70 Mbp. When correcting KLH-binding IgM NAb, and tIgM data for the estimated effects, the genomic region disappeared (data not shown). This shows that only one genomic region (i.e., one gene with one causal variant) is responsible for the observed effects.

All of the candidate genes in the IgM-associated region on chromosome 4 have been related to organism development, cell cycle, cell proliferation, metabolism, cancer development, viral infections, or a combination (mostly based on homologues in human and mouse studies), except for SMIM14, TMEM156, FAM114A1 (no literature available) (31). This implies that all candidate genes have a potential to influence IgM antibodies. However, some genes are reported in association with antibodies, or B cells, which makes them more likely candidates: KLF3 knockout mice were reported to have impaired B cell differentiation, including NAb-producing B1 B cells (37). TLR, including TLR1/TLR6 [homolog/paralog to TLR1A/TLR1B (38)], are well-known for their essential role in antibody production after stimulation of B cells, including B1 B cells (39–41). TBC1D1 has been reported to have a potential downstream role of silencing of pre-BCR signaling in acute lymphoblastic leukemia in human (42). However, none of these studies reported an effect on IgM without deliberate stimulation of B cells in healthy individuals, as is the case in this study.

Fine mapping of the associated region allowed to predict changes in protein functioning as a result of nucleotide changes in coding regions. A nucleotide change within TLR1A was predicted to have severe impact on protein folding/functioning, making this SNP the most likely candidate for the major effects on KLH-binding IgM NAb, tIgM, and KLH-binding IgT NAb. TLR1A, also known as TLR1-like a (TLR1La), TLR1.1, TLR1/6/10, and TLR16 (43), is one of the 10 known chicken TLR [see Brownlie and Allan (44) for a review on avian TLR]. TLR are among the most studied immune receptors. They are a family of transmembrane proteins that recognize conserved molecular patterns (pathogen-/microbe-associated molecular patterns; PAMP/MAMP), and are conserved in evolution (45). TLR1A dimerizes with TLR2 (either TLR2A or TLR2B), and recognizes peptidoglycans and related structures, including FSL-1, and PAM₃CSK₄ (33, 46), though some anomalies exist between studies (44). TLR1A contains 19 LRRs (33), of which LRR6-16 are required for ligand specificity (33). The region close to the C-terminal end of TLR1A is required for dimerization with TLR2 (in humans) (47, 48). The candidate SNP is located in LRR4 near the N-terminal end, suggesting not a direct influence on ligand recognition, or dimerization. Instead, it is likely influencing the tertiary protein structure, or mediation of coreceptors with TLR [see Van Bergenhenegouwen et al. (49) for a review], or both. Also, only coding regions were checked for amino acid substitutions, but not, for example, promotor regions influencing expression of genes. Based on the data in this study, no conclusion can be drawn on whether the most likely causal variant results in a complete loss of function, or in reduced functioning. Nevertheless, full dominance was observed for this TLR1A variant, suggesting that one functional copy of the gene results in a sufficient expression of functional TLR1A heterodimer and enhancement by the coreceptors for activation. Future studies should confirm the predicted TLR1A variant and its functionality.

The most likely causal variant has an effect on KLH-binding IgM NAb, and tIgM. The same region was found for IgM binding the self-antigen cardiolipin, and IgM binding the self-antigen ovalbumin (Bao et al., submitted). It is tempting to speculate that this TLR1A variant has a role on maturation/proliferation/survival of naive (NAb-producing) B cells, (KLH-binding) IgM production, or a combination of these. The most likely causal variant can affect the three stages of the (avian) B cell development: prebursal (only before hatch), bursal, and postbursal (50). Considering TLR1A as the most likely candidate, the effect on B cells, and IgM antibody levels is expected after exposure to PAMP, i.e., after hatch. TLR1A is expressed on various types of cells, including macrophages, and B cells, and in various organs, including thymus, spleen, liver, and bursa of Fabricius (44, 51). The effect of the TLR1A variant on IgM (NAb) can therefore be either directly on B cells or indirectly on B cells via for example macrophages, or both. An effect of the TLR1A variant is also expected on other cells of the immune system, both of the innate and adaptive part. This TLR1A variant is of particular interest, because it might explain part of the association of especially IgM NAb at adolescence with survival in chicken (11, 12). Yilmaz et al. (52) suggested that the TLR1 and TLR2 subfamilies are under positive selection (favoring polymorphisms) in chicken, while all other TLR are under negative/purifying selection (favoring conserved structures) (52). Huang et al. (38) suggested that only certain sites of TLR1 might be under positive selection, but that TLR1 is mostly under purifying selection (38). Nevertheless, the results imply that TLR1A is not essential for healthy chickens, but that it does influence IgM levels.

The identified region on chromosome 4 has previously been associated with many other characteristics in chickens (layer, broiler, or both) [based on Animal QTLdb (53)]. Only one of these are directly related to disease susceptibility: for susceptibility to Marek’s disease virus (54). Though the association with Marek’s susceptibility was stronger around TLR3 (also located on chromosome 4, around 62 Mbp), which is involved in the recognition of viruses, including Marek’s disease virus infections (55). Other characteristics associated with these region were for example age at first egg (56), eggshell breaking strength (57), eggshell color (58, 59), egg weight (57, 60, 61), feed intake (60), growth/body weight at several ages (56, 62–68), liver weight (69), and spleen weight (70). To investigate these associations, production traits of 222 females of the study population were previously collected [described in van der Klein et al. (15)], and could therefore be tested for their association with the most likely causal variant. Of the investigated traits, only the number of laid eggs between 17 and 24 weeks of age (i.e., adolescence) tended (p = 0.09) to a significant association with the TLR1A variant, where the genotype associated with lower IgM levels was associated with higher egg production (data not shown). Such a beneficial effect of the TLR1A variant on one of the most important production traits in layer chickens, in combination with a reduced natural selection pressure (due to high sanitary barriers on the breeding nucleus), might explain why the variant is present in the study population, even though it might impair (humoral) immunity.

The heritability of tIgA was twice as large as the heritability of KLH-binding IgA NAb. Interestingly, similar concentrations for tIgA, and tIgM were found, which is the first time to be reported. Genetic, and phenotypic correlations between KLH-binding IgA NAb and tIgA were positive, and low to moderate. This suggests that variation in KLH-binding IgA NAb, and tIgA are mainly regulated by different mechanisms. The role of blood IgA antibodies is not well understood so far. It has been suggested to be a response isotype with a proinflammatory function, similar to IgG antibodies (41, 71, 72). However, in chicken most of the intestinal IgA is secreted via the bile (73, 74), which means that intestinal IgA first has to be transported through the blood. IgA in blood might therefore represent the intestinal health status in chickens.

Two suggestively associated regions were identified for KLH-binding IgA NAb. These regions were not confirmed by tIgA concentration, even though a positive genetic correlation was estimated. This is the first time an association study on IgA (either NAb or total antibody concentration) in healthy individuals has been done. Observed results on KLH-binding IgA NAb were suggestive, and therefore these results first need to be confirmed in independent studies before further investigation of these associations.

Heritabilities of tIgG, and tIgT antibody concentrations were low, but comparable to heritabilities of KLH-binding IgG NAb titers, and KLH-binding IgT NAb titers (4). This indicates a strong environmental influence on IgG antibodies, and consequently on IgT, since IgT mainly consists of IgG antibodies: all genetic, and phenotypic correlations were higher than 0.8 [this study and Berghof et al. (4)]. However, KLH-binding IgT NAb, but not tIgT, was suggestively associated with the same region of KLH-binding IgM NAb, and tIgM. This confirms the high genetic correlation between KLH-binding IgT NAb, and KLH-binding IgM NAb (4), and the moderate genetic correlation between tIgT, and tIgM (this study).

Remarkably the genetic correlations between KLH-binding IgG NAb, and tIgG, and KLH-binding IgT NAb, and tIgT were negative. This suggest that KLH-binding IgG NAb, and tIgG, and KLH-binding IgT NAb, and tIgT are genetically regulated in a similar way, but with opposing effects. However, phenotypically there were very weak correlations between KLH-binding IgG NAb, and tIgG, and KLH-binding IgT NAb, and tIgT. In addition KLH-binding IgG NAb, and KLH-binding IgT NAb titers were significantly associated with the TLR1A variant, but tIgG, and tIgT concentrations were not. Meaning that a genetic predisposition for low (or high) levels of KLH-binding IgG NAb, and KLH-binding IgT NAb, or different genotypes for the TLR1A variant do not negatively (or positively) affect the tIgG antibody concentrations, and tIgT antibody concentrations. This implies that the humoral adaptive immune response, and memory formation are not negatively influenced by low NAb levels, or the TLR1A variant.

Natural antibodies are defined as antigen binding antibodies present in individuals without a previous exposure to this antigen (5). Layer chickens receive an extensive number of obligatory vaccinations, that may affect antibody levels in a (non-specific) fashion. In addition, although thought to be unlikely, KLH-binding antibodies could have been generated due to cross-reactivity with a component(s) of the vaccination(s). Also, due to KLH’s molecular weight (75), cross-reactivity of KLH-binding NAb with other antigens cannot be excluded. For example, cross reactivity of antibodies binding KLH with antigens of the tropical parasitic worm Schistosoma mansoni has been suggested (76). This might suggest that NAb can be made in the presence of exogenous antigens, which is currently under investigation [e.g., Baumgarth (77)]. Nevertheless, the whole study population received the same vaccinations, and a possible effect of moment of blood sampling (between 15 and 22 weeks of age) was corrected by the regression coefficient of “age” in the statistical model. Therefore, the vaccinations likely do not underlie the observed results.

In this study, we show the potential of using a large database with defined immunological traits to discover underlying genomic regions, and new fundamental components in the development, and functioning of the immune system. Recently in humans, similar approaches have been used to detect genetic variation for several immune traits, not only measuring baseline immunity but also functional immunological assays (78, 79). GWAS for defined immunological traits offer exiting, and promising opportunities for unraveling genetic variation, and fundamental mechanisms in the immune system in human, and animal species.

In summary, in this GWAS for KLH-binding NAb, and total concentrations of antibodies, we identified one significantly associated region on chromosome 4 with KLH-binding IgT NAb, total IgM concentration, and especially with KLH-binding IgM NAb. Further investigation predicted the causal variant on chromosome 4 to be located in TLR1A. As far as we know, this is the first time a TLR has been connected to IgM antibody levels in healthy individuals. This suggests a fundamental role of TLR in KLH-binding IgM NAb production, total IgM production, and possibly in (NAb-producing) B cell development, and proliferation. Future research should confirm the role of TLR1A, its role in KLH-binding IgM production, and total IgM production, and its role in general disease resistance in chickens.

References

• 1

  HodgesJ. Emerging boundaries for poultry production: challenges, dangers and opportunities. Worlds Poult Sci J (2009) 65:5–22.10.1017/S0043933909000014

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  Neeteson-Van NieuwenhovenA-MApplebyMCHogarthG. Making a resilient poultry industry in Europe. In: BurtonEGatcliffeJO’NeillHMScholeyD, editors. Sustainable Poultry Production in Europe. Oxfordshire, UK: CABI (2016). p. 3–24.

  • Google Scholar
  • View reference in article

• 3

  ChengHHKaiserPLamontSJ. Integrated genomic approaches to enhance genetic resistance in chickens. Annu Rev Anim Biosci (2013) 1:239–60.10.1146/annurev-animal-031412-103701

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BerghofTVLVan Der KleinSASArtsJAJParmentierHKVan Der PoelJJBovenhuisH. Genetic and non-genetic inheritance of natural antibodies binding keyhole limpet hemocyanin in a purebred layer chicken line. PLoS One (2015) 10:e0131088.10.1371/journal.pone.0131088

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BaumgarthN. The double life of a B-1 cell: self-reactivity selects for protective effector functions. Nat Rev Immunol (2011) 11:34–46.10.1038/nri2901

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  ZhouZ-HZhangYHuY-FWahlLMCisarJONotkinsAL. The broad antibacterial activity of the natural antibody repertoire is due to polyreactive antibodies. Cell Host Microbe (2007) 1:51–61.10.1016/j.chom.2007.01.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  PandaSDingJL. Natural antibodies bridge innate and adaptive immunity. J Immunol (2015) 194:13–20.10.4049/jimmunol.1400844

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  RothsteinTL. Natural antibodies as rheostats for susceptibility to chronic diseases in the aged. Front Immunol (2016) 7:127.10.3389/fimmu.2016.00127

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BaumgarthNWaffarnEENguyenTTT. Natural and induced B-1 cell immunity to infections raises questions of nature versus nurture. Ann N Y Acad Sci (2015) 1362:188–99.10.1111/nyas.12804

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  StarLFrankenaKKempBNieuwlandMGBParmentierHK. Natural humoral immune competence and survival in layers. Poult Sci (2007) 86:1090–9.10.1093/ps/86.6.1090

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  SunYParmentierHKFrankenaKVan Der PoelJJ. Natural antibody isotypes as predictors of survival in laying hens. Poult Sci (2011) 90:2263–74.10.3382/ps.2011-01613

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  WondmenehEVan ArendonkJAMVan Der WaaijEHDucroBJParmentierHK. High natural antibody titers of indigenous chickens are related with increased hazard in confinement. Poult Sci (2015) 94:1493–8.10.3382/ps/pev107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  BiscariniFBovenhuisHVan ArendonkJAMParmentierHKJungeriusAPVan Der PoelJJ. Across-line SNP association study of innate and adaptive immune response in laying hens. Anim Genet (2010) 41:26–38.10.1111/j.1365-2052.2009.01960.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  SunYBiscariniFBovenhuisHParmentierHKVan Der PoelJJ. Genetic parameters and across-line SNP associations differ for natural antibody isotypes IgM and IgG in laying hens. Anim Genet (2013) 44:413–24.10.1111/age.12014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  van der KleinSASBerghofTVLArtsJAJParmentierHKVan Der PoelJJBovenhuisH. Genetic relations between natural antibodies binding keyhole limpet hemocyanin and production traits in a purebred layer chicken line. Poult Sci (2015) 94:875–82.10.3382/ps/pev052

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  FrankenaK. The Interaction between Cooperia spp. and Ostertagia spp. (Nematoda: Trichostrongylidae) in Cattle [PhD Dissertation]. Wageningen (NL): Agricultural University Wageningen (1987).

  • Google Scholar
  • View reference in article

• 17

  De KoningDBDamenEPCWNieuwlandMGBVan GrevenhofEMHazelegerWKempBet alAssociation of natural (auto-) antibodies in young gilts with osteochondrosis at slaughter. Livest Sci (2015) 176:152–60.10.1016/j.livsci.2015.03.017

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  BrowningSRBrowningBL. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am J Hum Genet (2007) 81:1084–97.10.1086/521987

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  HayesBJBowmanPJChamberlainAJGoddardME. Invited review: genomic selection in dairy cattle: progress and challenges. J Dairy Sci (2009) 92:433–43.10.3168/jds.2008-1646

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  StearMNikbakhtGMatthewsLJonssonN. Breeding for disease resistance in livestock and fish. CAB Rev Perspect Agric Vet Sci Nutr Nat Resour (2012) 7:1–10.10.1079/PAVSNNR20127007

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  GilmourARGogelBJCullisBRThompsonR. ASReml User Guide, Release 4.1. Hemel Hempstead, UK: VSN International Ltd (2014).

  • Google Scholar
  • View reference in article

• 22

  DevlinBRoederK. Genomic control for association studies. Biometrics (1999) 55:997–1004.10.1111/j.0006-341X.1999.00997.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  LopesMSBastiaansenJWMHarliziusBKnolEFBovenhuisH. A genome-wide association study reveals dominance effects on number of teats in pigs. PLoS One (2014) 9:e105867.10.1371/journal.pone.0105867

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  AulchenkoYSRipkeSIsaacsAVan DuijnCM. GenABEL: an R package for genome-wide association analysis. Bioinformatics (2007) 23:1294–6.10.1093/bioinformatics/btm108

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  Wellcome Trust Case Control Consortium. Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature (2007) 447:661–78.10.1038/nature05911

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  DuijvesteijnNKnolEFBijmaP. Boar taint in entire male pigs: a genomewide association study for direct and indirect genetic effects on androstenone. J Anim Sci (2014) 92:4319–28.10.2527/jas.2014-7863

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  StoreyJDBassAJDabneyARobinsonD. qvalue: Q-Value Estimation for False Discovery Rate Control. R package version 2.4.2 (2015).

  • Google Scholar
  • View reference in article

• 28

  McLarenWGilLHuntSERiatHSRitchieGRSThormannAet alThe ensembl variant effect predictor. Genome Biol (2016) 17:122.10.1186/s13059-016-0974-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  SimN-LKumarPHuJHenikoffSSchneiderGNgPC. SIFT web server: predicting effects of amino acid substitutions on proteins. Nucleic Acids Res (2012) 40:W452–7.10.1093/nar/gks539

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  AdzhubeiIASchmidtSPeshkinLRamenskyVEGerasimovaABorkPet alA method and server for predicting damaging missense mutations. Nat Methods (2010) 7:248–9.10.1038/nmeth0410-248

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  NCBI Resource Coordinators. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res (2016) 44:D7–19.10.1093/nar/gkv1290

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  YatesAAkanniWAmodeMRBarrellDBillisKCarvalho-SilvaDet alEnsembl 2016. Nucleic Acids Res (2015) 44:D710–6.10.1093/nar/gkv1157

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  KeestraAMDe ZoeteMRVan AubelRAMHVan PuttenJPM. The central leucine-rich repeat region of chicken TLR16 dictates unique ligand specificity and species-specific interaction with TLR2. J Immunol (2007) 178:7110–9.10.4049/jimmunol.178.11.7110

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  SchmidMSmithJBurtDWAkenBLAntinPBArchibaldALet alThird report on chicken genes and chromosomes 2015. Cytogenet Genome Res (2015) 145:78–179.10.1159/000430927

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  ZhangLLiPLiuRZhengMSunYWuDet alThe identification of loci for immune traits in chickens using a genome-wide association study. PLoS One (2015) 10:e0117269.10.1371/journal.pone.0117269

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  MegensH-JCrooijmansRBastiaansenJKerstensHCosterAJalvingRet alComparison of linkage disequilibrium and haplotype diversity on macro- and microchromosomes in chicken. BMC Genet (2009) 10:86.10.1186/1471-2156-10-86

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  VuTTGattoDTurnerVFunnellAPWMakKSNortonLJet alImpaired B cell development in the absence of Krüppel-like factor 3. J Immunol (2011) 187:5032–42.10.4049/jimmunol.1101450

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  HuangYTemperleyNDRenLSmithJLiNBurtDW. Molecular evolution of the vertebrate TLR1 gene family – a complex history of gene duplication, gene conversion, positive selection and co-evolution. BMC Evol Biol (2011) 11:149.10.1186/1471-2148-11-149

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  RayABiswasT. Porin of Shigella dysenteriae enhances toll-like receptors 2 and 6 of mouse peritoneal B-2 cells and induces the expression of immunoglobulin M, immunoglobulin G2a and immunoglobulin A. Immunology (2005) 114:94–100.10.1111/j.1365-2567.2004.02002.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  AlugupalliKRAkiraSLienELeongJM. MyD88- and Bruton’s tyrosine kinase-mediated signals are essential for T cell-independent pathogen-specific IgM responses. J Immunol (2007) 178:3740–9.10.4049/jimmunol.178.6.3740

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  Bekeredjian-DingIJegoG. Toll-like receptors – sentries in the B-cell response. Immunology (2009) 128:311–23.10.1111/j.1365-2567.2009.03173.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  BicoccaVTChangBHMasoulehBKMuschenMLoriauxMMDrukerBJet alCrosstalk between ROR1 and the pre-B cell receptor promotes survival of t(1;19) acute lymphoblastic leukemia. Cancer Cell (2012) 22:656–67.10.1016/j.ccr.2012.08.027

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  TemperleyNDBerlinSPatonIRGriffinDKBurtDW. Evolution of the chicken toll-like receptor gene family: a story of gene gain and gene loss. BMC Genomics (2008) 9:62.10.1186/1471-2164-9-62

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  BrownlieRAllanB. Avian toll-like receptors. Cell Tissue Res (2011) 343:121–30.10.1007/s00441-010-1026-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  LeulierFLemaitreB. Toll-like receptors - taking an evolutionary approach. Nat Rev Genet (2008) 9:165–78.10.1038/nrg2303

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  HiguchiMMatsuoAShingaiMShidaKIshiiAFunamiKet alCombinational recognition of bacterial lipoproteins and peptidoglycan by chicken toll-like receptor 2 subfamily. Dev Comp Immunol (2008) 32:147–55.10.1016/j.dci.2007.05.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  GautamJKAshishComeauLDKruegerJKSmithMF. Structural and functional evidence for the role of the TLR2 DD loop in TLR1/TLR2 heterodimerization and signaling. J Biol Chem (2006) 281:30132–42.10.1074/jbc.M602057200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  JinMSKimSEHeoJYLeeMEKimHMPaikS-Get alCrystal structure of the TLR1-TLR2 heterodimer induced by binding of a tri-acylated lipopeptide. Cell (2007) 130:1071–82.10.1016/j.cell.2007.09.008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  Van BergenhenegouwenJPlantingaTSJoostenLABNeteaMGFolkertsGKraneveldADet alTLR2 & Co: a critical analysis of the complex interactions between TLR2 and coreceptors. J Leukoc Biol (2013) 94:885–902.10.1189/jlb.0113003

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  RatcliffeMJHHärtleSKaspersBKaiserP. Chapter 4 – B cells, the bursa of Fabricius and the generation of antibody repertoires. 2nd ed. In: SchatKAKaspersBKaiserP, editors. Avian Immunology. Boston, MA: Academic Press (2014). p. 65–89.

  • Google Scholar
  • View reference in article

• 51

  IqbalMPhilbinVJSmithAL. Expression patterns of chicken toll-like receptor mRNA in tissues, immune cell subsets and cell lines. Vet Immunol Immunopathol (2005) 104:117–27.10.1016/j.vetimm.2004.11.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  YilmazAShenSAdelsonDLXavierSZhuJJ. Identification and sequence analysis of chicken toll-like receptors. Immunogenetics (2005) 56:743–53.10.1007/s00251-004-0740-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  HuZ-LParkCAReecyJM. Developmental progress and current status of the animal QTLdb. Nucleic Acids Res (2016) 44:D827–33.10.1093/nar/gkv1233

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  HeifetzEMFultonJEO’SullivanNPArthurJAWangJDekkersJCMet alMapping quantitative trait loci affecting susceptibility to Marek’s disease virus in a backcross population of layer chickens. Genetics (2007) 177:2417–31.10.1534/genetics.107.080002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  HuXZouHQinAQianKShaoHYeJ. Activation of toll-like receptor 3 inhibits Marek’s disease virus infection in chicken embryo fibroblast cells. Arch Virol (2016) 161:521–8.10.1007/s00705-015-2674-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  PodisiBKKnottSADunnICLawASBurtDWHockingPM. Overlap of quantitative trait loci for early growth rate, and for body weight and age at onset of sexual maturity in chickens. Reproduction (2011) 141:381–9.10.1530/rep-10-0276

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  SasakiOOdawaraSTakahashiHNirasawaKOyamadaYYamamotoRet alGenetic mapping of quantitative trait loci affecting body weight, egg character and egg production in F2 intercross chickens. Anim Genet (2004) 35:188–94.10.1111/j.1365-2052.2004.01133.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 58

  GotoTIshikawaAOnitsukaSGotoNFujikawaYUminoTet alMapping quantitative trait loci for egg production traits in an F2 intercross of Oh-Shamo and White Leghorn chickens. Anim Genet (2011) 42:634–41.10.1111/j.1365-2052.2011.02190.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 59

  WolcAArangoJJankowskiTDunnISettarPFultonJEet alGenome-wide association study for egg production and quality in layer chickens. J Anim Breed Genet (2014) 131:173–82.10.1111/jbg.12086

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 60

  Tuiskula-HaavistoMHonkatukiaMVilkkiJDe KoningDJSchulmanNFMaki-TanilaA. Mapping of quantitative trait loci affecting quality and production traits in egg layers. Poult Sci (2002) 81:919–27.10.1111/j.1365-2052.2009.01914.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 61

  SchreiweisMAHesterPYSettarPMoodyDE. Identification of quantitative trait loci associated with egg quality, egg production, and body weight in an F2 resource population of chickens1. Anim Genet (2006) 37:106–12.10.1111/j.1365-2052.2005.01394.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 62

  CarlborgOHockingPMBurtDWHaleyCS. Simultaneous mapping of epistatic QTL in chickens reveals clusters of QTL pairs with similar genetic effects on growth. Genet Res (2004) 83:197–209.10.1017/S0016672304006779

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 63

  ZhouHDeebNEvock-CloverCMAshwellCMLamontSJ. Genome-wide linkage analysis to identify chromosomal regions affecting phenotypic traits in the chicken. I. Growth and average daily gain. Poult Sci (2006) 85:1700–11.10.1093/ps/85.10.1700

  • CrossRef
  • Google Scholar
  • View reference in article

• 64

  NadafJPitelFGilbertHDuclosMJVignolesFBeaumontCet alQTL for several metabolic traits map to loci controlling growth and body composition in an F2 intercross between high- and low-growth chicken lines. Physiol Genomics (2009) 38:241–9.10.1152/physiolgenomics.90384.2008

  • CrossRef
  • Google Scholar
  • View reference in article

• 65

  Ankra-BaduGABihan-DuvalELMignon-GrasteauSPitelFBeaumontCDuclosMJet alMapping QTL for growth and shank traits in chickens divergently selected for high or low body weight. Anim Genet (2010) 41:400–5.10.1111/j.1365-2052.2009.02017.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 66

  GoragaZSNassarMKBrockmannGA. Quantitative trait loci segregating in crosses between New Hampshire and White Leghorn chicken lines: I. Egg production traits. Anim Genet (2012) 43:183–9.10.1111/j.1365-2052.2011.02233.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 67

  PodisiBKKnottSABurtDWHockingPM. Comparative analysis of quantitative trait loci for body weight, growth rate and growth curve parameters from 3 to 72 weeks of age in female chickens of a broiler-layer cross. BMC Genet (2013) 14:22.10.1186/1471-2156-14-22

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 68

  NassarMKGoragaZSBrockmannGA. Quantitative trait loci segregating in crosses between New Hampshire and White Leghorn chicken lines: IV. Growth performance. Anim Genet (2015) 46:441–6.10.1111/age.12298

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 69

  NavarroPVisscherPMKnottSABurtDWHockingPMHaleyCS. Mapping of quantitative trait loci affecting organ weights and blood variables in a broiler layer cross. Br Poult Sci (2005) 46:430–42.10.1080/00071660500158055

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 70

  ZhouHDeebNEvock-CloverCMAshwellCMLamontSJ. Genome-wide linkage analysis to identify chromosomal regions affecting phenotypic traits in the chicken. II. Body composition. Poult Sci (2006) 85:1712–21.10.1093/ps/85.10.1712

  • CrossRef
  • Google Scholar

• 71

  Van EgmondMDamenCAVan SprielABVidarssonGVan GarderenEVan De WinkelJGJ. IgA and the IgA Fc receptor. Trends Immunol (2001) 22:205–11.10.1016/S1471-4906(01)01873-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 72

  CeruttiA. The regulation of IgA class switching. Nat Rev Immunol (2008) 8:421–34.10.1038/nri2322

  • CrossRef
  • Google Scholar
  • View reference in article

• 73

  Lebacq-VerheydenA-MVaermanJPHeremansJF. Quantification and distribution of chicken immunoglobulins IgA, IgM and IgG in serum and secretions. Immunology (1974) 27:683–92.

  • Google Scholar
  • View reference in article

• 74

  RomboutJHWMSijtsmaSRWestCEKarabinisYSijtsmaOKWVan Der ZijppAJet alEffect of vitamin A deficiency and Newcastle disease virus infection on IgA and IgM secretion in chickens. Br J Nutr (1992) 68:753–63.10.1079/BJN19920131

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 75

  HarrisJRMarklJ. Keyhole limpet hemocyanin (KLH): a biomedical review. Micron (1999) 30:597–623.10.1016/S0968-4328(99)00036-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 76

  GeyerHWuhrerMResemannAGeyerR. Identification and characterization of keyhole limpet hemocyanin N-glycans mediating cross-reactivity with Schistosoma mansoni. J Biol Chem (2005) 280:40731–48.10.1074/jbc.M505985200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 77

  BaumgarthN. B-1 cell heterogeneity and the regulation of natural and antigen-induced IgM production. Front Immunol (2016) 7:324.10.3389/fimmu.2016.00324

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 78

  Aguirre-GamboaRJoostenIUrbanoPCMVan Der MolenRGVan RijssenEVan CranenbroekBet alDifferential effects of environmental and genetic factors on T and B cell immune traits. Cell Rep (2016) 17:2474–87.10.1016/j.celrep.2016.10.053

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 79

  LiYOostingMSmeekensSPJaegerMAguirre-GamboaRLeKTTet alA functional genomics approach to understand variation in cytokine production in humans. Cell (2016) 167:1099–110.e1014.10.1016/j.cell.2016.10.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article