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

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:

$\text{logit OD}=\text{ln}\left(\frac{\text{OD}}{\left({\text{OD}}_{\text{max}}-\text{OD}\right)}\right),$

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:

$\text{titer}=\frac{{\text{logit OD}}_{\text{lpw}}-\left({\text{logit OD}}_{\text{sample}}-\mathrm{\beta }\times {\text{log}}_{\text{2}}\left({\text{dilution}}_{\text{sample}}\right)\right)}{\mathrm{\beta }},$

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).

TABLE 1

Table 1. Descriptive statistics of KLH-binding natural antibody (NAb) titers, and total antibody concentrations (μg/ml) in an adolescent White Leghorn chicken population for KLH-binding IgM NAb (IgM) titer, KLH-binding IgA NAb (IgA) titer, KLH-binding IgG NAb (IgG) titer, total KLH-binding NAb (IgT) titer, tIgM concentration, total IgA (tIgA) concentration, total IgG (tIgG) concentration, and total antibody (tIgT) concentration.

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 ${\text{(χ}}_1^2\text{-value}\ge \text{600)}$ were removed. (4) SNP on the allosome Z with a strong deviation from Hardy Weinberg Equilibrium ${\text{(χ}}_1^2\text{-value}\ge \text{600)}$ 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 (${\mathrm{\sigma }}_{\text{a}}^2$), and the maternal environmental variance (${\mathrm{\sigma }}_{\text{m}}^2$) 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 (${\mathrm{\sigma }}_{\text{a}}^2$), 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

${\text{y}}_{\mathit{ij}}=\mathrm{\mu }+{\text{P}}_i+{\mathrm{\beta }}_{\text{1}}\times {\text{Age}}_{\mathit{ij}}+{\text{id}}_j+{\text{e}}_{\mathit{ij}},$ (1)

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 $~\text{N}\left(\mathrm{0,}\mathbf{\text{A}}{\mathrm{\sigma }}_{\text{a}}^2\right)$, and e_ij is the residual term assumed to be distributed as $~\text{N}\left(\mathrm{0,}\mathbf{\text{I}}{\mathrm{\sigma }}_{\text{e}}^2\right)$. Assumed (co)variance structures of the random model terms are ${\mathbf{\text{A}}\mathrm{\sigma }}_{\text{a}}^2$ and ${\mathbf{\text{I}}\mathrm{\sigma }}_{\text{e}}^2$, in which A is the additive genetic relationship matrix, ${\mathrm{\sigma }}_{\text{a}}^2$ is the additive genetic variance, I is an identity matrix, and ${\mathrm{\sigma }}_{\text{e}}^2$ 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:

${\text{y}}_{\mathit{ijk}}=\mathrm{\mu }+{\text{P}}_i+{\mathrm{\beta }}_{\text{1}}\times {\text{Age}}_{\mathit{ijk}}+{\text{SNP}}_j+{\text{id}}_k+{\text{e}}_{\mathit{ijk}},$ (2a)

and IgM NAb was analyzed using model 2b:

${\text{y}}_{\mathit{ijkl}}=\mathrm{\mu }+{\text{P}}_i+{\mathrm{\beta }}_{\text{1}}\times {\text{Age}}_{\mathit{ijkl}}+{\text{SNP}}_j+{\text{id}}_k+{\text{d}}_l+{\text{e}}_{\mathit{ijkl}},$ (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 $~\text{N}\left(\mathrm{0,}\mathbf{\text{A}}{\mathrm{\sigma }}_{\text{a}}^2\right)$, d_l is the random effect of the lth dam (IgM NAb only) assumed to be distributed as $~\text{N}\left(\mathrm{0,}\mathbf{\text{I}}{\mathrm{\sigma }}_{\text{m}}^2\right)$, and e_ijk and e_ijkl are the residual terms assumed to be distributed as $~\text{N}\left(\mathrm{0,}\mathbf{\text{I}}{\mathrm{\sigma }}_{\text{e}}^2\right)$. Assumed (co)variance structures of the random model terms are ${\mathbf{\text{A}}\mathrm{\sigma }}_{\text{a}}^2$, and ${\mathbf{\text{I}}\mathrm{\sigma }}_{\text{e}}^2$ (for all), and ${\mathbf{\text{I}}\mathrm{\sigma }}_{\text{m}}^2$ (for IgM NAb only), in which A is the additive genetic relationship matrix, ${\mathrm{\sigma }}_{\text{a}}^2$ is the additive genetic variance, I are identity matrices, ${\mathrm{\sigma }}_{\text{e}}^2$ is the residual variance, and ${\mathrm{\sigma }}_{\text{m}}^2$ 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 (${\mathrm{\sigma }}_{\text{a}}^2$) and the maternal environmental variance (${\mathrm{\sigma }}_{\text{m}}^2$) 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 (${\mathrm{\sigma }}_{\text{a}}^2$) 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 ${\mathrm{\sigma }}_{\text{SNP}}^2∕{\mathrm{\sigma }}_{\text{a}}^2\times \mathrm{100\%}$, where ${\mathrm{\sigma }}_{\text{SNP}}^2$ 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 (${\mathrm{\sigma }}_{\text{m}}^2$) 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).

TABLE 2

Table 2. Estimated genetic correlations, and phenotypic correlations of KLH-binding natural antibody (NAb) titers and total antibody concentrations (μg/ml) in an adolescent White Leghorn chicken population.

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

Figure 1. Manhattan plots of genome-wide association studies (GWAS) for keyhole limpet hemocyanin (KLH)-binding natural antibody (NAb) titers, and total antibody concentrations (μg/ml) in an adolescent White Leghorn chicken population. The figures show GWAS results for: (A) KLH-binding IgM NAb (IgM) titer; (B) total IgM (tIgM) concentration; (C) KLH-binding IgA NAb (IgA) titer; (D) total IgA (tIgA) concentration; (E) KLH-binding IgG NAb (IgG) titer; (F) total IgG (tIgG) concentration; (G) total KLH-binding NAb (IgT) titer; and (H) total antibody (tIgT) concentration. tIgM, tIgA, tIgG, and tIgT were log₁₀-tranformed for the analyses. Genomic control was applied to test statistics of IgM NAb, tIgM, tIgA, tIgG, and tIgT. On the x-axis is the physical position of the 15,579 used SNP across the chicken genome (chromosomes 1–28, linkage groups LGE64, chromosome Z, and unplaced represented by –). Positions are based on Gallus_gallus-5.0. On the y-axis is the −log₁₀(p-value_SNP) for association of the SNP with the investigated trait. Each dot represents one SNP. Note the different y-axes for IgM and tIgM. The false discovery rate (FDR) threshold was set at 0.05 for significant SNP (solid horizontal line) and at 0.20 for suggestive SNP (dashed horizontal line). If no SNP reached the FDR threshold, the threshold could not be estimated. For IgM, no −log₁₀(p-value_SNP) reached 0.05 ≤ FDR ≤ 0.20, and therefore, FDR ≤ 0.20 is not shown.

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.

TABLE 3

Table 3. Genomic regions significantly, or suggestively associated with KLH-binding natural antibody (NAb) titers, and total antibody concentrations (μg/ml) in an adolescent White Leghorn chicken population.

TABLE 4

Table 4. Lead SNP (SNP with highest p-value) within genomic regions significantly, or suggestively associated with KLH-binding natural antibody (NAb) titers, and total antibody concentrations (μg/ml) in an adolescent White Leghorn chicken population.

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

Figure 2. Manhattan plots of association studies of the selected region on chromosome 4 for keyhole limpet hemocyanin (KLH)-binding IgM NAb (IgM) titer in red and total IgM (tIgM) concentration (log₁₀-tranformed for the analyses) in blue in an adolescent White Leghorn chicken population. On the y-axes is the −log₁₀(p-value_SNP) for the association of the investigated SNP with IgM, and tIgM after genomic control was applied. Each dot represents one SNP. The figures show: (A) on the x-axis is the physical position of the approximately 40,000 SNP used on chromosome 4 for the selected region; (B) zoom of peak-region in figure A [−log₁₀(p-value_SNP) ≥ 5]; (C) zoom of peak-region in figure B [−log₁₀(p-value_SNP) ≥ 10 for IgM and −log₁₀(p-value_SNP) ≥ 6 for tIgM], including candidate gene overview (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.

TABLE 5

Table 5. Significance, and effects of the predicted causal variant (SNP on chromosome 4 position 69,965,939 bp, based on Gallus_gallus-5.0) of KLH-binding IgM natural antibody (NAb) (IgM) titer, and total IgM (tIgM) concentration (μg/ml) on all measured KLH-binding NAb titers, and total antibody concentrations (μg/ml) in an adolescent White Leghorn chicken population.

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.

Ethics 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.

Author Contributions

TB performed the research, analyzed, and interpreted the data, and wrote the manuscript. MV and JA analyzed and interpreted the data, and revised the manuscript. HP, JP, and HB designed the research, interpreted the data, and critically revised the manuscript. AV performed the research and analyzed and interpreted the data.

Conflict of Interest Statement

This study was in kind financed by Hendrix Genetics. In addition, AV is employed by Hendrix Genetics Research. AV was involved in genotyping and sequencing of the chickens and imputation of genomic data to the described genomic data sets. These genomic data sets are also of interest for Hendrix Genetics’ commercial targets, but this interest did not influence the results in this manuscript in any matter. Except for the delivered data, and the results reported in this project, or for other projects, no other shared interests (e.g., employment, consultancy, patents, products) exist between Hendrix Genetics, and Wageningen University & Research. All other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Ger de Vries Reilingh (Adaptation Physiology) and Tessa Brinker and Alex Hulzebosch (Animal Breeding and Genomics) for technical assistance.

Funding

This work is in kind financed and intellectually supported by Hendrix Genetics. This work is part of the research program “Divergent selection for natural antibodies in poultry” with project number 12208, which is financed by the Netherlands Organisation for Scientific Research (NWO).

Supplementary Material

The Supplementary Material for this article can be found online at http://www.frontiersin.org/articles/10.3389/fimmu.2017.01879/full#supplementary-material.

Footnotes

1. ^http://www.fao.org/ag/againfo/themes/en/poultry/animal_health.html.
2. ^See also http://www.antibodyregistry.org.

References

1. Hodges J. Emerging boundaries for poultry production: challenges, dangers and opportunities. Worlds Poult Sci J (2009) 65:5–22. doi:10.1017/S0043933909000014

CrossRef Full Text | Google Scholar

2. Neeteson-Van Nieuwenhoven A-M, Appleby MC, Hogarth G. Making a resilient poultry industry in Europe. In: Burton E, Gatcliffe J, O’Neill HM, Scholey D, editors. Sustainable Poultry Production in Europe. Oxfordshire, UK: CABI (2016). p. 3–24.

Google Scholar

3. Cheng HH, Kaiser P, Lamont SJ. Integrated genomic approaches to enhance genetic resistance in chickens. Annu Rev Anim Biosci (2013) 1:239–60. doi:10.1146/annurev-animal-031412-103701

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Berghof TVL, Van Der Klein SAS, Arts JAJ, Parmentier HK, Van Der Poel JJ, Bovenhuis H. Genetic and non-genetic inheritance of natural antibodies binding keyhole limpet hemocyanin in a purebred layer chicken line. PLoS One (2015) 10:e0131088. doi:10.1371/journal.pone.0131088

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Baumgarth N. The double life of a B-1 cell: self-reactivity selects for protective effector functions. Nat Rev Immunol (2011) 11:34–46. doi:10.1038/nri2901

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Zhou Z-H, Zhang Y, Hu Y-F, Wahl LM, Cisar JO, Notkins AL. The broad antibacterial activity of the natural antibody repertoire is due to polyreactive antibodies. Cell Host Microbe (2007) 1:51–61. doi:10.1016/j.chom.2007.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Panda S, Ding JL. Natural antibodies bridge innate and adaptive immunity. J Immunol (2015) 194:13–20. doi:10.4049/jimmunol.1400844

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Rothstein TL. Natural antibodies as rheostats for susceptibility to chronic diseases in the aged. Front Immunol (2016) 7:127. doi:10.3389/fimmu.2016.00127

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Baumgarth N, Waffarn EE, Nguyen TTT. Natural and induced B-1 cell immunity to infections raises questions of nature versus nurture. Ann N Y Acad Sci (2015) 1362:188–99. doi:10.1111/nyas.12804

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Star L, Frankena K, Kemp B, Nieuwland MGB, Parmentier HK. Natural humoral immune competence and survival in layers. Poult Sci (2007) 86:1090–9. doi:10.1093/ps/86.6.1090

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Sun Y, Parmentier HK, Frankena K, Van Der Poel JJ. Natural antibody isotypes as predictors of survival in laying hens. Poult Sci (2011) 90:2263–74. doi:10.3382/ps.2011-01613

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Wondmeneh E, Van Arendonk JAM, Van Der Waaij EH, Ducro BJ, Parmentier HK. High natural antibody titers of indigenous chickens are related with increased hazard in confinement. Poult Sci (2015) 94:1493–8. doi:10.3382/ps/pev107

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Biscarini F, Bovenhuis H, Van Arendonk JAM, Parmentier HK, Jungerius AP, Van Der Poel JJ. Across-line SNP association study of innate and adaptive immune response in laying hens. Anim Genet (2010) 41:26–38. doi:10.1111/j.1365-2052.2009.01960.x

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Sun Y, Biscarini F, Bovenhuis H, Parmentier HK, Van Der Poel JJ. Genetic parameters and across-line SNP associations differ for natural antibody isotypes IgM and IgG in laying hens. Anim Genet (2013) 44:413–24. doi:10.1111/age.12014

PubMed Abstract | CrossRef Full Text | Google Scholar

15. van der Klein SAS, Berghof TVL, Arts JAJ, Parmentier HK, Van Der Poel JJ, Bovenhuis H. Genetic relations between natural antibodies binding keyhole limpet hemocyanin and production traits in a purebred layer chicken line. Poult Sci (2015) 94:875–82. doi:10.3382/ps/pev052

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

17. De Koning DB, Damen EPCW, Nieuwland MGB, Van Grevenhof EM, Hazeleger W, Kemp B, et al. Association of natural (auto-) antibodies in young gilts with osteochondrosis at slaughter. Livest Sci (2015) 176:152–60. doi:10.1016/j.livsci.2015.03.017

CrossRef Full Text | Google Scholar

18. Browning SR, Browning BL. 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. doi:10.1086/521987

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hayes BJ, Bowman PJ, Chamberlain AJ, Goddard ME. Invited review: genomic selection in dairy cattle: progress and challenges. J Dairy Sci (2009) 92:433–43. doi:10.3168/jds.2008-1646

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Stear M, Nikbakht G, Matthews L, Jonsson N. Breeding for disease resistance in livestock and fish. CAB Rev Perspect Agric Vet Sci Nutr Nat Resour (2012) 7:1–10. doi:10.1079/PAVSNNR20127007

CrossRef Full Text | Google Scholar

21. Gilmour AR, Gogel BJ, Cullis BR, Thompson R. ASReml User Guide, Release 4.1. Hemel Hempstead, UK: VSN International Ltd (2014).

Google Scholar

22. Devlin B, Roeder K. Genomic control for association studies. Biometrics (1999) 55:997–1004. doi:10.1111/j.0006-341X.1999.00997.x

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Lopes MS, Bastiaansen JWM, Harlizius B, Knol EF, Bovenhuis H. A genome-wide association study reveals dominance effects on number of teats in pigs. PLoS One (2014) 9:e105867. doi:10.1371/journal.pone.0105867

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Aulchenko YS, Ripke S, Isaacs A, Van Duijn CM. GenABEL: an R package for genome-wide association analysis. Bioinformatics (2007) 23:1294–6. doi:10.1093/bioinformatics/btm108

CrossRef Full Text | Google Scholar

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. doi:10.1038/nature05911

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Duijvesteijn N, Knol EF, Bijma P. 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. doi:10.2527/jas.2014-7863

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Storey JD, Bass AJ, Dabney A, Robinson D. qvalue: Q-Value Estimation for False Discovery Rate Control. R package version 2.4.2 (2015).

Google Scholar

28. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The ensembl variant effect predictor. Genome Biol (2016) 17:122. doi:10.1186/s13059-016-0974-4

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Sim N-L, Kumar P, Hu J, Henikoff S, Schneider G, Ng PC. SIFT web server: predicting effects of amino acid substitutions on proteins. Nucleic Acids Res (2012) 40:W452–7. doi:10.1093/nar/gks539

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, et al. A method and server for predicting damaging missense mutations. Nat Methods (2010) 7:248–9. doi:10.1038/nmeth0410-248

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yates A, Akanni W, Amode MR, Barrell D, Billis K, Carvalho-Silva D, et al. Ensembl 2016. Nucleic Acids Res (2015) 44:D710–6. doi:10.1093/nar/gkv1157

CrossRef Full Text | Google Scholar

33. Keestra AM, De Zoete MR, Van Aubel RAMH, Van Putten JPM. 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. doi:10.4049/jimmunol.178.11.7110

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Schmid M, Smith J, Burt DW, Aken BL, Antin PB, Archibald AL, et al. Third report on chicken genes and chromosomes 2015. Cytogenet Genome Res (2015) 145:78–179. doi:10.1159/000430927

CrossRef Full Text | Google Scholar

35. Zhang L, Li P, Liu R, Zheng M, Sun Y, Wu D, et al. The identification of loci for immune traits in chickens using a genome-wide association study. PLoS One (2015) 10:e0117269. doi:10.1371/journal.pone.0117269

CrossRef Full Text | Google Scholar

36. Megens H-J, Crooijmans R, Bastiaansen J, Kerstens H, Coster A, Jalving R, et al. Comparison of linkage disequilibrium and haplotype diversity on macro- and microchromosomes in chicken. BMC Genet (2009) 10:86. doi:10.1186/1471-2156-10-86

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Vu TT, Gatto D, Turner V, Funnell APW, Mak KS, Norton LJ, et al. Impaired B cell development in the absence of Krüppel-like factor 3. J Immunol (2011) 187:5032–42. doi:10.4049/jimmunol.1101450

CrossRef Full Text | Google Scholar

38. Huang Y, Temperley ND, Ren L, Smith J, Li N, Burt DW. 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. doi:10.1186/1471-2148-11-149

CrossRef Full Text | Google Scholar

39. Ray A, Biswas T. 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. doi:10.1111/j.1365-2567.2004.02002.x

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Alugupalli KR, Akira S, Lien E, Leong JM. MyD88- and Bruton’s tyrosine kinase-mediated signals are essential for T cell-independent pathogen-specific IgM responses. J Immunol (2007) 178:3740–9. doi:10.4049/jimmunol.178.6.3740

CrossRef Full Text | Google Scholar

41. Bekeredjian-Ding I, Jego G. Toll-like receptors – sentries in the B-cell response. Immunology (2009) 128:311–23. doi:10.1111/j.1365-2567.2009.03173.x

CrossRef Full Text | Google Scholar

42. Bicocca VT, Chang BH, Masouleh BK, Muschen M, Loriaux MM, Druker BJ, et al. Crosstalk between ROR1 and the pre-B cell receptor promotes survival of t(1;19) acute lymphoblastic leukemia. Cancer Cell (2012) 22:656–67. doi:10.1016/j.ccr.2012.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Temperley ND, Berlin S, Paton IR, Griffin DK, Burt DW. Evolution of the chicken toll-like receptor gene family: a story of gene gain and gene loss. BMC Genomics (2008) 9:62. doi:10.1186/1471-2164-9-62

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Brownlie R, Allan B. Avian toll-like receptors. Cell Tissue Res (2011) 343:121–30. doi:10.1007/s00441-010-1026-0

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Leulier F, Lemaitre B. Toll-like receptors - taking an evolutionary approach. Nat Rev Genet (2008) 9:165–78. doi:10.1038/nrg2303

CrossRef Full Text | Google Scholar

46. Higuchi M, Matsuo A, Shingai M, Shida K, Ishii A, Funami K, et al. Combinational recognition of bacterial lipoproteins and peptidoglycan by chicken toll-like receptor 2 subfamily. Dev Comp Immunol (2008) 32:147–55. doi:10.1016/j.dci.2007.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Gautam JK, Ashish, Comeau LD, Krueger JK, Smith MF. 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. doi:10.1074/jbc.M602057200

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Jin MS, Kim SE, Heo JY, Lee ME, Kim HM, Paik S-G, et al. Crystal structure of the TLR1-TLR2 heterodimer induced by binding of a tri-acylated lipopeptide. Cell (2007) 130:1071–82. doi:10.1016/j.cell.2007.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Van Bergenhenegouwen J, Plantinga TS, Joosten LAB, Netea MG, Folkerts G, Kraneveld AD, et al. TLR2 & Co: a critical analysis of the complex interactions between TLR2 and coreceptors. J Leukoc Biol (2013) 94:885–902. doi:10.1189/jlb.0113003

CrossRef Full Text | Google Scholar

50. Ratcliffe MJH, Härtle S, Kaspers B, Kaiser P. Chapter 4 – B cells, the bursa of Fabricius and the generation of antibody repertoires. 2nd ed. In: Schat KA, Kaspers B, Kaiser P, editors. Avian Immunology. Boston, MA: Academic Press (2014). p. 65–89.

Google Scholar

51. Iqbal M, Philbin VJ, Smith AL. Expression patterns of chicken toll-like receptor mRNA in tissues, immune cell subsets and cell lines. Vet Immunol Immunopathol (2005) 104:117–27. doi:10.1016/j.vetimm.2004.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Yilmaz A, Shen S, Adelson DL, Xavier S, Zhu JJ. Identification and sequence analysis of chicken toll-like receptors. Immunogenetics (2005) 56:743–53. doi:10.1007/s00251-004-0740-8

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Hu Z-L, Park CA, Reecy JM. Developmental progress and current status of the animal QTLdb. Nucleic Acids Res (2016) 44:D827–33. doi:10.1093/nar/gkv1233

CrossRef Full Text | Google Scholar

54. Heifetz EM, Fulton JE, O’Sullivan NP, Arthur JA, Wang J, Dekkers JCM, et al. Mapping quantitative trait loci affecting susceptibility to Marek’s disease virus in a backcross population of layer chickens. Genetics (2007) 177:2417–31. doi:10.1534/genetics.107.080002

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Hu X, Zou H, Qin A, Qian K, Shao H, Ye J. Activation of toll-like receptor 3 inhibits Marek’s disease virus infection in chicken embryo fibroblast cells. Arch Virol (2016) 161:521–8. doi:10.1007/s00705-015-2674-x

CrossRef Full Text | Google Scholar

56. Podisi BK, Knott SA, Dunn IC, Law AS, Burt DW, Hocking PM. 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. doi:10.1530/rep-10-0276

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Sasaki O, Odawara S, Takahashi H, Nirasawa K, Oyamada Y, Yamamoto R, et al. Genetic mapping of quantitative trait loci affecting body weight, egg character and egg production in F2 intercross chickens. Anim Genet (2004) 35:188–94. doi:10.1111/j.1365-2052.2004.01133.x

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Goto T, Ishikawa A, Onitsuka S, Goto N, Fujikawa Y, Umino T, et al. Mapping quantitative trait loci for egg production traits in an F2 intercross of Oh-Shamo and White Leghorn chickens. Anim Genet (2011) 42:634–41. doi:10.1111/j.1365-2052.2011.02190.x

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Wolc A, Arango J, Jankowski T, Dunn I, Settar P, Fulton JE, et al. Genome-wide association study for egg production and quality in layer chickens. J Anim Breed Genet (2014) 131:173–82. doi:10.1111/jbg.12086

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Tuiskula-Haavisto M, Honkatukia M, Vilkki J, De Koning DJ, Schulman NF, Maki-Tanila A. Mapping of quantitative trait loci affecting quality and production traits in egg layers. Poult Sci (2002) 81:919–27. doi:10.1111/j.1365-2052.2009.01914.x

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Schreiweis MA, Hester PY, Settar P, Moody DE. 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. doi:10.1111/j.1365-2052.2005.01394.x

CrossRef Full Text | Google Scholar

62. Carlborg O, Hocking PM, Burt DW, Haley CS. Simultaneous mapping of epistatic QTL in chickens reveals clusters of QTL pairs with similar genetic effects on growth. Genet Res (2004) 83:197–209. doi:10.1017/S0016672304006779

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Zhou H, Deeb N, Evock-Clover CM, Ashwell CM, Lamont SJ. 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. doi:10.1093/ps/85.10.1700

CrossRef Full Text | Google Scholar

64. Nadaf J, Pitel F, Gilbert H, Duclos MJ, Vignoles F, Beaumont C, et al. QTL 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. doi:10.1152/physiolgenomics.90384.2008

CrossRef Full Text | Google Scholar

65. Ankra-Badu GA, Bihan-Duval EL, Mignon-Grasteau S, Pitel F, Beaumont C, Duclos MJ, et al. Mapping QTL for growth and shank traits in chickens divergently selected for high or low body weight. Anim Genet (2010) 41:400–5. doi:10.1111/j.1365-2052.2009.02017.x

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Goraga ZS, Nassar MK, Brockmann GA. Quantitative trait loci segregating in crosses between New Hampshire and White Leghorn chicken lines: I. Egg production traits. Anim Genet (2012) 43:183–9. doi:10.1111/j.1365-2052.2011.02233.x

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Podisi BK, Knott SA, Burt DW, Hocking PM. 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. doi:10.1186/1471-2156-14-22

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Nassar MK, Goraga ZS, Brockmann GA. Quantitative trait loci segregating in crosses between New Hampshire and White Leghorn chicken lines: IV. Growth performance. Anim Genet (2015) 46:441–6. doi:10.1111/age.12298

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Navarro P, Visscher PM, Knott SA, Burt DW, Hocking PM, Haley CS. Mapping of quantitative trait loci affecting organ weights and blood variables in a broiler layer cross. Br Poult Sci (2005) 46:430–42. doi:10.1080/00071660500158055

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Zhou H, Deeb N, Evock-Clover CM, Ashwell CM, Lamont SJ. Genome-wide linkage analysis to identify chromosomal regions affecting phenotypic traits in the chicken. II. Body composition. Poult Sci (2006) 85:1712–21. doi:10.1093/ps/85.10.1712

CrossRef Full Text | Google Scholar

71. Van Egmond M, Damen CA, Van Spriel AB, Vidarsson G, Van Garderen E, Van De Winkel JGJ. IgA and the IgA Fc receptor. Trends Immunol (2001) 22:205–11. doi:10.1016/S1471-4906(01)01873-7

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Cerutti A. The regulation of IgA class switching. Nat Rev Immunol (2008) 8:421–34. doi:10.1038/nri2322

CrossRef Full Text | Google Scholar

73. Lebacq-Verheyden A-M, Vaerman JP, Heremans JF. Quantification and distribution of chicken immunoglobulins IgA, IgM and IgG in serum and secretions. Immunology (1974) 27:683–92.

Google Scholar

74. Rombout JHWM, Sijtsma SR, West CE, Karabinis Y, Sijtsma OKW, Van Der Zijpp AJ, et al. Effect of vitamin A deficiency and Newcastle disease virus infection on IgA and IgM secretion in chickens. Br J Nutr (1992) 68:753–63. doi:10.1079/BJN19920131

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Harris JR, Markl J. Keyhole limpet hemocyanin (KLH): a biomedical review. Micron (1999) 30:597–623. doi:10.1016/S0968-4328(99)00036-0

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Geyer H, Wuhrer M, Resemann A, Geyer R. Identification and characterization of keyhole limpet hemocyanin N-glycans mediating cross-reactivity with Schistosoma mansoni. J Biol Chem (2005) 280:40731–48. doi:10.1074/jbc.M505985200

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Baumgarth N. B-1 cell heterogeneity and the regulation of natural and antigen-induced IgM production. Front Immunol (2016) 7:324. doi:10.3389/fimmu.2016.00324

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Aguirre-Gamboa R, Joosten I, Urbano PCM, Van Der Molen RG, Van Rijssen E, Van Cranenbroek B, et al. Differential effects of environmental and genetic factors on T and B cell immune traits. Cell Rep (2016) 17:2474–87. doi:10.1016/j.celrep.2016.10.053

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Li Y, Oosting M, Smeekens SP, Jaeger M, Aguirre-Gamboa R, Le KTT, et al. A functional genomics approach to understand variation in cytokine production in humans. Cell (2016) 167:1099–110.e1014. doi:10.1016/j.cell.2016.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar