Skip to content

Advertisement

  • Research article
  • Open Access

A genome-wide association study for natural antibodies measured in blood of Canadian Holstein cows

  • 1,
  • 3, 5,
  • 2,
  • 3, 4,
  • 1 and
  • 3, 5Email author
Contributed equally
BMC Genomics201819:694

https://doi.org/10.1186/s12864-018-5062-6

  • Received: 1 June 2018
  • Accepted: 10 September 2018
  • Published:

Abstract

Background

Natural antibodies (NAb) are an important component of the innate immune system, and fight infections as a part of the first line defence. NAb are poly-reactive and can respond non-specifically to antigens. Therefore, NAb may be a key trait when evaluating an animal’s potential natural disease resistance. Variation in NAb is caused by both genetic and environmental factors. In this study genetic parameters of NAb were estimated and a genome-wide association study (GWAS) was performed to gain further understanding on the genes that are responsible for the observed genetic variation of NAb in Canadian Holsteins.

Results

In total, blood samples of 1327 cows from 64 farms were studied. NAb binding to keyhole limpet hemocyanin (KLH) were determined via indirect ELISA. Immunoglobulin (Ig) isotypes, IgG and IgM, were evaluated. From the sample population, 925 cows were genotyped for 45,187 markers and each individual marker was tested to detect genetic variation in NAb levels. The relationships among animals was accounted for with genomic relationship. Results show heritabilities of 0.27 ± 0.064 (IgG) and 0.31 ± 0.065 (IgM). In total, 23 SNPs were found to be associated with IgG, but no SNPs were associated with IgM (FDR p-value < 0.05). The significant SNPs were located on autosomal chromosomes 1, 20 and 21 of the cow genome. Functional annotation analysis of the positional candidate genes revealed two sets of genes with biologically relevant functions related to NAb. In one set, seven genes with crucial roles in the production of antibody in B cells were associated with the trafficking of vesicles inside the cells between organelles. In the second set, two genes among positional candidate genes were associated with isotype class-switching and somatic hypermutation of B cells.

Conclusions

This study demonstrated the possibility of increasing NAb through selective breeding. In addition, the effects of two candidate pathways are proposed for further investigation of NAb production in Holsteins.

Keywords

  • Genome-wide-association study
  • Natural antibody
  • Dairy cattle
  • Immune system
  • Vesicle trafficking
  • Immunoglobulin class switching

Background

Over the last several decades, breeding in dairy cattle mainly focused on production and fertility traits, with less emphasis on health traits. Health problems, however, can cause substantial economic losses to the dairy industry. The economic losses, together with the rising awareness of animal welfare, increased herd size, and less attention for individual animals, have led to a growing focus on health traits. Health parameters currently used in dairy cattle breeding programs consider milk quality related traits, including somatic cell count or score (SCC or SCS) and clinical mastitis resistance, as well as claw related health parameters, and metabolic disorders, such as ketosis. Somatic cell score is often used as a health indicator in breeding indices and has a heritability between 0.10 and 0.17, when based on lactation means [14]. Somatic cell score was associated with mastitis [5], and is therefore, related to udder-quality rather than the total health status of a cow.

Currently, national animal breeding indices do not select for overall health status based on immune capacity. Therefore, an opportunity exists to find immune response parameters with potential association with overall disease resistance. Natural antibodies (NAb) are a candidate for improving innate host defence in dairy cows. NAb is produced by B1-cells and this type of antibody (Ab) with broad specificity is a component of the innate immune system that can be found in animals without prior exposure to antigens [6]. In this way, NAb can play an important role in the first line of defence against different kinds of pathogens since they are polyreactive Ab that bind different conserved structures including carbohydrate, nucleic acid and phospholipids [7, 8]. Different NAb isotypes exist, with isotype IgM being most common, and isotypes IgG and IgA found to a lesser extent [9]. Additionally, NAb provide a link between the innate and adaptive immune system [10]. In dairy cows, NAb are measurable in both blood- and milk samples [11, 12].

Many studies have shown that it is possible to select for healthier cows based on adaptive immune response traits [1317]. In these studies, the relationships between economically important diseases, including mastitis, metritis, pneumonia, displaced abomasum and retained placenta of dairy cows and their immune responses were examined and it was demonstrated that specific antibodies (SpAb), which are part of the adaptive immune system, are useful biomarkers for disease resistance. Cows were tested for two different adaptive immune responses, antibody-mediated (AMIR) and cell-mediated immune response (CMIR) [17]. A decreased incidence of mastitis was found for cows with a high AMIR, CMIR and overall immune response. Moreover, a higher incidence of metritis was found for cows with a low CMIR and cows with low overall immune response had an increased incidence of retained foetal membrane and displaced abomasum [17]. However, SpAb are only measured after immunization of the cows, whereas NAb can be measured prior to immunization, and therefore may have practical advantages in a breeding program.

Recently, NAb have been studied as a predictor for different diseases in dairy cows. NAb binding KLH and Lipopolysaccharides (LPS) in plasma were found to be negatively associated with metabolic health, but not statistically significant. In contrast, the level of these NAb in milk were positively associated with metabolic health, significantly [18, 19]. Associations between immune response and metabolism are not uncommon [2023]. Van Knegsel et al. found a positive trend between clinical mastitis incidence and NAb binding to the self-antigen, myosin, and a negative trend between mastitis incidence and NAb binding to the self-antigen, transferrin [19]. Banos et al. hypothesized that higher NAb levels (binding to KLH) can be associated with improved capacity of the innate immune system to respond to pathogen challenges. However, they also found that a poorer nutritional status was related with higher NAb levels [24]. Together these studies suggested a potential association between NAb and incidence of infectious and metabolic diseases, but caution must be exercised since some NAb can be cross reactive binding to self antigens.

Variation in NAb stems from both genetic and environmental factors. Heritability estimates of NAb measured in blood or milk of dairy cows, ranges from 0.08–0.45 [11, 2527], where IgM generally has the highest heritability estimates (0.18–0.45) and IgG has lower heritability estimates (0.08–0.31). Additionally, NAb measured in blood serum had higher heritability (0.15–0.25) compared to NAb measured in milk (0.08–0.23) [11]. The relatively high heritability (especially for the IgM isotype) shows potential for genetic selection.

The objective of this study was to identify regions of the bovine genome that are associated with the level of NAb of the IgM and IgG isotypes. Results are expected to provide a further step in unravelling the genetic control of NAb, which is, given the resemblances in NAb between species, not only relevant for dairy cattle but also for many other species.

Results

In total, 925 genotyped animals were included in this analysis. The dataset contained 45,187 SNPs and after removing SNPs on X chromosome and with MAF less than 1%, the number of SNPs used for the final analyses was 43,471. Animals with deviating days in milk (> 500 days in milking) were removed from the total dataset (n = 6). Mean and the corresponding standard deviations for both IgG and IgM isotypes are shown in Table 1. Heritability based on genomic information for isotype IgG was estimated at 0.27 ± 0.064 and for isotype IgM heritability was estimated at 0.31 ± 0.065.
Table 1

Means and standard deviations for the corrected optical destiny of natural antibodies with immunoglobulin isotypes IgG and IgM with and without log-transformation; n =1327

Variable

Mean

Std dev.

Min

Max

IgG

0.07

7.58

0.00

18.45

IgM

1.15

2.01

0.11

9.79

IgG_log10

−1.17

0.88

−4.19

1.27

IgM_log10

0.06

0.30

−0.96

0.99

At a genome-wide FDR of 5%, NAb IgM did not show any significant associations (Fig. 1 and Additional file 1: Figure S1). Significant associations between 23 SNPs and NAb IgG were found (Fig. 2, Additional file 2: Figure S2 and Table 2). The genomic inflation factor for IgG and IgM was 0.99 and .96, respectively. Figures 3 and 4 show Q-Q plots of test statistics. The significant SNPs were located on Bos taurus autosomal chromosome 1, (1 SNP), 20 (1 SNP) and 21 (21 SNPs). The region between the significant SNPs and the adjacent SNP on each side were screened to identify the positional candidate gene. The selected regions around each significant SNPs ranged from 51,008 bp to 276,066 bp with the average of 53,583 bp on each side of the significant SNPs. Thirty-seven genes were found on the bovine genome (UMD3.1, R.91) near the significant SNPs. Two of these genes were located on BTA-1 and the others were located on BTA-21 (Table 3). No gene was found around the significant SNP on BTA-20, despite the large region around this SNP (276,066 bp).
Fig. 1
Fig. 1

Distribution of -Log10 P-values from single SNP analyses for natural antibody isotype IgM binding KLH, for every chromosome. The red line indicates FDR rate of 1% and green line indicates FDR rate of 5%

Fig. 2
Fig. 2

Distribution of -Log10 P-values from single SNP analyses for natural antibody isotype IgG binding KLH, for every chromosome. The red line indicates FDR rate of 1% and green line indicates FDR rate of 5%

Table 2

The list of SNPs that are significantly associated with IgG NAb

No.

SNP ID

BTA:Position

MAF

ASE

FDR Corrected p-value

1

BTA-95285-NO-RS

1:45452593

0.4121

−0.3635

3.93E-02

2

ARS-BFGL-NGS-80731

20:52529599

0.4225

0.358795

1.15E-02

3

ARS-BFGL-NGS-36353

21:56272444

0.3745

0.373862

1.42E-02

4

HAPMAP36811-SCAFFOLD240339_791

21:57920570

0.0753

0.668197

1.16E-02

5

ARS-BFGL-NGS-18270

21:62445122

0.0966

0.669079

1.03E-03

6

ARS-BFGL-NGS-114799

21:66492068

0.4214

−0.44719

6.56E-04

7

ARS-BFGL-NGS-20339

21:67088847

0.3417

−0.37211

1.81E-02

8

ARS-BFGL-NGS-29558

21:67345260

0.2074

−0.52465

6.56E-04

9

ARS-BFGL-NGS-119656

21:67963614

0.1654

0.469889

3.17E-02

10

ARS-BFGL-NGS-21626

21:68105539

0.2091

−0.51097

9.35E-04

11

BTA-03263-RS29011028

21:68356479

0.2298

−0.39036

4.61E-02

12

ARS-BFGL-NGS-86477

21:68399787

0.3504

−0.54776

9.82E-07

13

HAPMAP41289-BTA-53093

21:68751694

0.0464

−0.82127

1.34E-02

14

ARS-BFGL-NGS-39941

21:69279283

0.2102

0.441591

1.25E-02

15

ARS-BFGL-NGS-115062

21:69395154

0.4973

−0.46877

5.27E-05

16

ARS-BFGL-NGS-109976

21:69891992

0.2959

0.422929

3.81E-03

17

ARS-BFGL-NGS-1345

21:69939350

0.4885

−0.36527

1.37E-02

18

ARS-BFGL-NGS-33371

21:69972343

0.4088

−0.56953

8.87E-08

19

ARS-BFGL-NGS-63968

21:70012413

0.0393

1.02312

9.35E-04

20

ARS-BFGL-NGS-10484

21:70050042

0.0486

0.926758

1.03E-03

21

ARS-BFGL-NGS-40895

21:70116946

0.4509

−0.41571

9.13E-04

22

ARS-BFGL-NGS-106472

21:70732448

0.0922

0.67538

1.74E-03

23

BTA-100472-NO-RS

21:70759096

0.185

−0.50167

1.49E-03

BTA Bos taurus autosome, MAF Minor Allele Frequency, ASE Allele Substitution Effect

Fig. 3
Fig. 3

Q-Q plot of –Log10 p-values for natural antibody isotype IgM

Fig. 4
Fig. 4

Q-Q plot of –Log10 p-values for natural antibody isotype IgG

Table 3

The list of positional candidate genes nearby SNPs that are significantly associated with IgG NAb

Gene stable ID

Gene name

Gene start (bp)

Chromosome Number

ENSBTAG00000012140

ADGRG7

45,415,066

1

ENSBTAG00000016363

TFG

45,489,052

1

ENSBTAG00000010416

RIN3

57,859,148

21

ENSBTAG00000018310

SETD3

66,064,391

21

ENSBTAG00000018636

CCNK

66,159,547

21

ENSBTAG00000025181

CCDC85C

66,175,707

21

ENSBTAG00000013491

EML1

66,358,950

21

ENSBTAG00000007233

WDR25

66,917,780

21

ENSBTAG00000017716

BEGAIN

67,072,575

21

ENSBTAG00000038073

 

68,149,406

21

ENSBTAG00000043578

 

68,150,604

21

ENSBTAG00000037977

 

68,151,114

21

ENSBTAG00000020192

PPP2R5C

68,377,171

21

ENSBTAG00000007007

WDR20

68,639,658

21

ENSBTAG00000007008

MOK

68,706,703

21

ENSBTAG00000007009

 

68,751,572

21

ENSBTAG00000044797

U6

68,761,089

21

ENSBTAG00000010993

 

68,778,065

21

ENSBTAG00000010995

CINP

68,789,514

21

ENSBTAG00000043100

U5

68,797,764

21

ENSBTAG00000012143

TECPR2

68,837,087

21

ENSBTAG00000013475

TRAF3

69,228,047

21

ENSBTAG00000015750

AMN

69,266,381

21

ENSBTAG00000016456

CDC42BPB

69,279,222

21

ENSBTAG00000026893

EXOC3L4

69,418,979

21

ENSBTAG00000035995

TNFAIP2

69,455,152

21

ENSBTAG00000021199

APOPT1

69,847,359

21

ENSBTAG00000017299

KLC1

69,896,140

21

ENSBTAG00000009966

XRCC3

69,960,277

21

ENSBTAG00000003556

ZFYVE21

69,978,023

21

ENSBTAG00000003559

PPP1R13B

69,982,291

21

ENSBTAG00000026886

MP68

70,113,045

21

ENSBTAG00000046186

 

70,114,200

21

ENSBTAG00000020402

TDRD9

70,123,917

21

ENSBTAG00000046401

RD3L

70,133,009

21

ENSBTAG00000022775

C21H14orf180

70,707,408

21

ENSBTAG00000006673

TMEM179

70,717,761

21

Candidate genes based on functional annotation analysis are shown in bold

The Over Representation Analysis (ORA) on the Gene Ontology (GO) terms associated with all the positional candidate genes (37 genes) revealed 3 significant over-represented terms (FDR p-value < 0.05): “thyroxine 5-deiodinase activity” (GO:0033798, for molecular function), “hormone biosynthetic process” (GO:0042446, for biological process) and “exocyst” (GO:0000145, for cellular compartment). It also should be noted that “exocytosis” (GO:0006887, for the biological process) was over-represented, but the FDR p-value was just above significant level (p-value = 0.076).

The functions of all 37 genes were individually checked to investigate any possible biological similarity between the candidate genes [2830]. Seven genes (TFG on BTA-1, and EML1, AMN, CDC428BP, EXOC3L4, TNFAIP2, and KLC1 on BTA-21) were found among the candidate genes that are associated with intracellular transportation of vesicles. In addition, two genes (TRAF3 and XRCC3 on BTA-21) were found among the positional candidates which are associated with immunoglobulin class switching and somatic hypermutation of B-cells [3134].

Discussion

In this study no significant associations were found for NAb IgM even though the heritability of this trait was relatively high. This may indicate that the underlying genetic architecture of IgM is more polygenic with no outstanding major genes. For Isotype IgG significant associations were found on chromosome 1, 20 and 21 (BTA-1, 20 and 21). For both traits, several peaks across genomes were observed that were not significant. These suggestive peaks should be further investigated in a bigger sample size.

In order to check for existence of population stratification or family structure, the GWAS study was repeated without including the animal random effect (the results were not shown). Genomic inflation factors (λ) were 1.25 and 1.44 for IgG and IgM, respectively, showing severe inflation in p-values. After properly correcting for population stratification and family structure, λ should be close to 1 [35]. Incorporating the full genomic covariance structure between individuals by fitting the random additive effect in the model resulted in λ close to 1 for both traits (0.99 for IgG and 0.96 for IgM; Figs. 3 and 4).

The present study included NAb of isotypes IgG and IgM, binding KLH. Natural antibodies are considered as the humoral part of the innate immune system [6]. The model antigen KLH is assumed not to be present in the common daily environment of dairy cows, and therefore can be assumed as an antigen that cows are not exposed. Consequently, any antibody detected to KLH is likely NAb or possibly specific antibody cross reacting to antigens to which the individual has been exposed. These antigens may include self antigens and therefore caution must be exercised if selecting for increased amounts of NAb until any cross reacting self antigens can be determined. Generally, anti-KLH antibody is considered Nab [27].

Natural antibodies are thought to be an indicator of innate host defence. Moreover, NAb are variable between individual cows and are under substantial genetic control [11, 26, 27]. These findings, together with the ability of early protection and poly-reactivity of NAb, make them potential traits to aid in genetic selection for disease resistance.

In this GWAS, significant and suggestive associations between markers and NAb IgG were found on chromosome 1, 20 and 21 (BTA-1, 20 and 21). In this study to identify candidate regions, the immediate adjacent SNPs on each side of the significant SNPs were selected, instead of selecting a fixed size. In this approach, fewer genes can be identified compared to the fixed size approach, but it helps to remove noise by excluding genes when their adjacent SNPs are not significant.

The ORA showed “thyroxine 5-deiodinase activity” and “hormone biosynthetic process” as the most represented terms among the candidate genes. However, closer inspection revealed that both of these terms were overrepresented due to the presence of two genes (DIO3 and IDBG-647912) among the candidate genes. Both genes, along with a microRNA, are located in the vicinity of one SNP (21:68105539). Due to the lack of any supporting data that correlate the production of antibody with the function of these two genes, or the presence of genes with similar activity near other SNPs, this enrichment may be a type II error. However, it is possible that the microRNA (ENSBTAG00000038073) controls the production of these antibodies or there is one or more unknown genes in this region that control the production of antibody. Further investigation on the cellular compartment GO term (“exocyst”) revealed evidence to support the “trafficking of vesicles” as a candidate pathway. TFG, located on BTA-1, is responsible for the trafficking of newly synthesized protein from the endoplasmic reticulum (ER) to golgi apparatus (GA). KLC1, located on BTA-21, is responsible for retro-trafficking of vesicles from GA to ER; and AMN, on BTA-21, is responsible for the trafficking of vesicles from GA to the plasma membrane (PM). In addition, EXOC3L4 and TNFAIP2, are responsible in exocytosis, the process of releasing vesicle content to the extracellular environment. Moreover, EML1 and CDC42BPB play a role in organization of cytoskeleton that provides a pathway for the movement of vesicles inside cells in the process of exocytosis [36]. The molecular pathway of Ab production by B lymphocytes, involves several glycosylation steps that requires trafficking between ER, GA, and PM and the final release of antibody from the cell [37]. This pathway aligns with the function of these seven candidate genes. Moreover, the presence of TRAF3, with a major role in class-switching of immunoglobulin, and XRCC3, with a role in DNA recombination and hypermutation in B cells, align with the isotype of the NAb that was measured in in this study [3134].

These results from the pathway analysis also point out the similarity between the production of SpAb and NAb. It has been previously shown that B1 cells in mice producing IgG and IgA NAb undergo some level of affinity maturation [38, 39]. Given the fact that affinity maturation happens after the presentation of antigens to the B cells, it is possible that B1 cells are also exposed to an antigen to produce IgG. The effect of antigen on B1 lymphocytes and NAb has been shown by the reduced level of IgG NAb, but not IgM, in germ-free antigen-free mice [40], however, the source of the antigenic signal has not yet been found. T cell-independent antigens from microbiota and self-antigens were proposed as the source of antigenic stimulation [38].

Recently, the composition of gut microbiota was found to be correlated with the level of poly-reactive IgA in germ-free mice, but this correlation was caused by the genetic structure of the mice rather than gut microbiota [41]. It is noteworthy to mention that genetic association studies have not found any significant correlation between the SpAb and NAb. The genetic correlation in different studies can be found when the phenotype is based on NAb even against different antigens or isotypes rather than between SpAb and NAb [11, 2527, 42].

The available information from GWAS of NAb in cows is limited, but BTA-21 is in common between this study and the study by Wijga et al. [43]. However, there was no similarity between markers associated with SpAb and NAb. From the immunological perspective, it can be hypothesized that the main signals from antigen presenting cells to B2 lymphocytes probably overcome the effect of vesicle trafficking at the level of Ab production, but in the absence of strong classical signal, vesicle trafficking gains more effect on the production of NAb.

From the breeding perspective, the absence of correlation between NAb and SpAb provides an opportunity to include NAb as a new trait to improve dairy cattle health given the notable heritability of NAb. Despite numerous reports on the beneficial role of IgM NAb and to a lesser extent IgG NAb in mice [38, 39, 44], the beneficial correlations reported between NAb and resistance against infectious diseases are minimal in livestock [27, 45, 46]. Therefore, the role of NAb in immunity of livestock needs further investigation. It should also be noted that a number of studies are proposing self-antigens as the stimulator of B1 cells and reporting both negative and positive roles of self-reactive NAb in autoimmune diseases [39]. In mice, lupus-like autoimmune diseases developed when NAb production was impaired, but mice that could produce only IgM NAb, but not IgG NAb, did not develop the autoimmune disease [47, 48]. Therefore, to develop breeding strategies to increase NAb in animals with relatively long lifespan such as cattle the isotype of NAb and antigen should be cautiously explored and animals monitored for any signs of autoimmunity.

Conclusion

The level of NAb in cattle is moderately heritable. Despite similar heritability of IgG and IgM classes of NAb, the structure of genes that are associated with them seem to be different. The absence of any significant SNPs with IgM, given the sample size of this study, is likely representing the polygenetic control of the production of IgM. In the case of IgG, BTA-1,20 and 21 seem to carry the main quantitative trait loci. Based on the function of positional candidate genes that were found in this study, vesicle trafficking in B1 cells from production to secretion of IgG NAb play an important role on the level of IgG NAb in serum. Self-reactive antigens have been reported as stimulators of B1-NAb-producing cells and therefore caution must be exercised if selection based on NAb is to be developed.

Methods

Animals and phenotypes

Blood serum samples were taken of 1327 cows from 64 privately owned herds in Canada. Cows were on average 161 days in lactation when the blood samples were taken, with a range from 0 to 498, and 13 cows with DIM > 500. Cows had on average a parity of 2.16, ranging from 0 to 12. Number of sampled cows per farm ranged from 2 to 160 cows per farm. The cows were the offspring of 472 unique bulls (on average: 2.77 cows per sire) and 1010 unique dams. The pedigree contained 13,128 animals and was provided by the Canadian Dairy Network (CDN; Guelph, Ontario, Canada). Immune response phenotypes used in this study were based on NAb of the isotypes, IgG and IgM, tested against the model antigen, not previously seen by cattle, keyhole limpet hemocyanin (KLH). To obtain the NAb phenotypes an indirect ELISA procedure was used, as described by Thompson-Crispi (2013a). Briefly, 96 well plates were coated with 5 μg/ml KLH) (MP Biomedicals, Solon, OH) in carbonate-bicarbonate buffer (pH 9.6), and incubated over night at a temperature of 4 °C. Next day, plates were washed 3 times with PBS and 0.05% Tween 20 (Sigma-Aldrich Canada Ltd., Oakville, ON, Canada) (wash buffer pH 7.4) and blocked with PBS, 3% Tween 20, 1,5% BSA and 1.5% FCS for 1 h at room temperature (RT) then washed again 3 times. Four serial dilutions starting with 1/40 of the serum samples in wash buffer were added to the plate and incubated for 2 h at RT. Plates were wash 5 times and the secondary antibodies conjugated to alkaline phosphatase dissolved in Tris-Tween buffer with 0.05% Tween 20 (pH 7.4) were added to the plates: either 1: 10.000 monoclonal anti-bovine IgG from mouse ascites fluid (Sigma-Aldrich, St. Louis. MO, USA) or 1:5000 anti-bovine IgM produced in sheep (Bethyl Laboratories, Montgomery, TX, USA) and incubated for 1 h at RT. All wash steps were performed with ELx405 Auto Plate washer (Biotek Instruments Inc., Winooski, USA). Substrate (p-nitrophenyl phosphate) (Sigma-Aldrich Canada Ltd., Oakville, On, Canada) was added and incubated for about 30–60 min. Optical density (OD) values at 405 nm were obtained using EL808 plate reader (BioTek Instruments Inc., Winooski, VT, USA). Optical density values were corrected to the rolling mean of the positive controls for each plate to account for day and plate variation, as described by Heriazon et al. [49]. The dilutions of the corrected OD values were summed and duplicates averaged for statistical analysis. NAb, both IgG and IgM isotype, were log-transformed to accomplish normality. After removing outliers for both IgG and IgM by Median Absolute Deviation method (considering three equivalent standard deviations). All experimental procedures were approved by the Animal Care Committee of the University of Guelph under guidelines of the Canadian Council of Animal Care.

Genotyping and quality control

DNA was extracted from hair follicles and genotyping was performed with the Illumina Bovine SNP50 BeadChip by Zoetis Canada (Kirkland, Quebec, Canada). The initial dataset contained 45,187 SNP markers that are used in routine official genomic evaluation in Canada by the Canadian Dairy Network (CDN) (Guelph, ON, Canada). Details of quality control were explained in Wiggans et al. [50]. In the present study SNPs located on the X chromosome were not included and due to relatively small sample size, SNPs with MAF < 1% in the 1716 SNPs were excluded resulting in 43,471 SNPs for GWAS. The sporadic missing genotypes were imputed using 50,000 reference Holsteins from the CDN database by FImpute software [51]. After considering all the quality control measures, 925 genotyped cows were included in the association study.

Statistical analysis

The association of the individual NAb with each individual SNP was estimated following a univariate mixed linear model:
$$ \mathrm{y}=\mathrm{Xb}+\mathrm{c}\beta +\mathrm{Wg}+\mathrm{e} $$

where y is the log-transformed corrected OD for NAb, b is a vector of fixed effects including overall mean, days in milk (classes: 1 = 0-20dim, 2 = 21–105 dim, 3 = 106-235dim, 4 > 235dim), parity (classes: 0 = heifers before calving, 1 = parity 1, 2 = parity 2, 3 = parity 3 and 4 = parity 4 and higher), herd (classes: 1 to 64), β is the gene substitution effect for the SNP, g is the random genetic effects, e is the random residual effects, X is an incidence matrix relating elements of b to y, c is a vector of genotypes for the SNP coded as 0 = BB, 1 = AB, 2 = AA and W is a standardized genotype matrix with element \( {\mathrm{w}}_{ij}=\left({c}_{ij}-2{p}_i\right)/\sqrt{2\sum p\left(1-p\right)} \), where pi is the allele frequency of the ith SNP and cij is genotype of ith SNP of jth individual.

The co-variance matrix for the vector y is:
$$ \mathrm{V}=\mathrm{G}{\sigma}_g^2+\mathrm{I}{\sigma}_e^2 $$
with \( \mathrm{g}\sim N\left(0,\mathrm{G}{\sigma}_g^2\right) \) and\( \mathrm{e}\sim N\left(0,\mathrm{I}{\sigma}_e^2\right) \), where \( {\sigma}_g^2 \) and \( {\sigma}_e^2 \) denote variance of random genetic effects and residual variance, respectively. G is the genomic relationship matrix calculated according to VanRaden (2008) using genome-wide SNP information as G = WW [52].

Additive genetic variance, residual variance and subsequently heritability were estimated with restricted maximum likelihood (REML) method and the average information algorithm [53].

Fitting random animal effect with the use of genomic relationship matrix prevents false-positives association due to population stratification and cryptic relationships between individuals and also increases the power [54]. Therefore, the above model should be proper for the Holstein population that has strong family structure due to the widespread use of few top bulls each year.

Inflation or deflation in p-values due to stratification or family structure was assessed by genomic inflation factor (λ) and also visually inspected by quantile-quantile (Q-Q) plot. λ is calculated as the median of the χ2 test statistics (1 degree of freedom) divided by its theoretical median under the null distribution [35]. In order to adjust for multiple comparisons, false discovery rate was controlled at 1 and 5% genome-wise levels [55]. GWAS and heritability estimation were carried out by snp1101 software [51].

Functional annotation of the positional candidate genes

To identify the positional candidate genes that are associated with NAb, the genomic regions around the significant SNPs (FDR corrected p-value < 0.05) up to the next immediate SNP in both directions (before and after the significant SNP) were selected. These regions were cross-referenced against the cow genome (UMD 3.1, Ensemble Genes Release 91) using the BioMart tool in the Ensemble website to identify genes that are located in the vicinity of the significant SNPs. The extracted genes were then submitted into the innateDB online database to identify the GO terms that are associated with the candidate genes and subsequent ORA [30, 56].

Notes

Abbreviations

Ab: 

Antibody

AMIR: 

Antibody-mediated immune response

BTA: 

Bos taurus autosomal chromosome

CMIR: 

Cell-mediated immune response

FDR: 

False discovery rate

GO: 

Gene ontology

GWAS: 

Genome-wide association study

Ig: 

Immunoglobulin

KLH: 

Keyhole limpet hemocyanin

LPS: 

Lipopolysaccharides

NAb: 

Natural antibodies

ORA: 

Over Representation Analysis

SCC or SCS: 

Somatic cell count or score

SNP: 

Single nucleotide polymorphism

SpAb: 

Specific antibodies

Declarations

Acknowledgements

The authors acknowledge the members of Dr. Mallard’s lab involved in the immune response sampling, Shannon Cartwright and Dr. Julie Schmied for technical laboratory advice, and Dr. Filippo Miglior (CDN) for providing the pedigree. The authors would like to thank Dr. Ducro and Dr. van Arendonk for training Britt de Klerk in her graduate program at Wageningen University.

Funding

This research was funded by grants to B.A. Mallard from the Natural Sciences and Engineering Research Council (NSERC). This paper is also a contribution to the Food from Thought research program supported by the Canada First Research Excellence Fund. The funders did not have any role in the design, collection, analysis, and interpretation of data and in writing the manuscript. ME was supported by Natural Sciences and Engineering Research Council via scholarship.

Availability of data and materials

The datasets used during the current study are available from the corresponding author on reasonable request. In addition, the list of proposed genes and their position that were found in this study are provided in the Table 3 in this manuscript.

Authors’ contributions

BK measured the NAb in serum samples and drafted the manuscript during her visit of BM’s lab. ME ran the functional annotation analysis and drafted the manuscript. KTC optimized the NAb ELISA, collected serum samples and hair follicles, measured the NAb in serum samples and contributed to the drafting of the manuscript. MS ran the GWAS study and prepared the graphs. JP helped in the experimental design of this project. As the PI, BM designed and spearheaded the project. All authors read and approved the final manuscript.

Ethics approval

All experimental procedures were approved by the Animal Care Committee of the University of Guelph under guidelines of the Canadian Council of Animal Care. Written consents from the owners of participant herds were obtained to use their cows in this study.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Animal Breeding and Genomics Centre, Wageningen University, P.O. Box 338, Wageningen, The Netherlands
(2)
Trouw Nutrition Canada Inc., 150 Research Lane, Suite 200, Guelph, ON, N1G 4T2, Canada
(3)
Centre for Genetic Improvement of Livestock, University of Guelph, Guelph, ON, Canada
(4)
Semex Alliance, Guelph, ON, Canada
(5)
Department of Pathobiology, Ontario Veterinary College, Genetic improvement of livestock, University of Guelph, Guelph, ON, N1G 2W1, Canada

References

  1. Pryce JE, Esslemont RJ, Thompson R, Veerkamp RF, Kossaibati MA, Simm G. Estimation of genetic parameters using health, fertility and production data from a management recording system for dairy cattle. Anim Sci. 1998;66:577–84. https://doi.org/10.1017/S1357729800009152.View ArticleGoogle Scholar
  2. Rupp R, Boichard D. Genetic parameters for clinical mastitis, somatic cell score, production, udder type traits, and milking ease in first lactation Holsteins. J Dairy Sci. 1999;82:2198–204. https://doi.org/10.3168/jds.S0022-0302(99)75465-2.View ArticlePubMedGoogle Scholar
  3. Carlén E, Strandberg E, Roth A. Genetic parameters for clinical mastitis, somatic cell score, and production in the first three lactations of Swedish Holstein cows. J Dairy Sci. 2004;87:3062–70. https://doi.org/10.3168/jds.S0022-0302(04)73439-6.View ArticlePubMedGoogle Scholar
  4. Rupp R, Boichard D. Genetics of resistance to mastitis in dairy cattle. Vet Res. 2003;34:671–88. https://doi.org/10.1051/vetres:2003020.View ArticlePubMedGoogle Scholar
  5. Harmon RJ. Physiology of mastitis and factors affecting somatic cell counts. J Dairy Sci. 1994;77:2103–12. https://doi.org/10.3168/jds.S0022-0302(94)77153-8.View ArticlePubMedGoogle Scholar
  6. Baumgarth N, Tung JW, Herzenberg LA. Inherent specificities in natural antibodies: a key to immune defense against pathogen invasion. Springer Semin Immunopathol. 2005;26:347–62. https://doi.org/10.1007/s00281-004-0182-2.View ArticlePubMedGoogle Scholar
  7. Boes M. Role of natural and immune IgM antibodies in immune responses. Mol Immunol. 2000;37:1141–9. https://doi.org/10.1016/S0161-5890(01)00025-6.View ArticlePubMedGoogle Scholar
  8. Casali P, Schettino EW. Structure and Function of Natural Antibodies. Berlin, Heidelberg: Springer; 1996. p. 167–79. https://doi.org/10.1007/978-3-642-85226-8_17.View ArticleGoogle Scholar
  9. Kohler H, Bayry J, Nicoletti A, Kaveri SV. Natural autoantibodies as tools to predict the outcome of immune response? Scand J Immunol. 2003;58:285–9. https://doi.org/10.1046/j.1365-3083.2003.01314.x.View ArticlePubMedGoogle Scholar
  10. Ochsenbein AF, Zinkernagel RM. Natural antibodies and complement link innate and acquired immunity. Immunol Today. 2000;21:624–30. https://doi.org/10.1016/S0167-5699(00)01754-0.View ArticlePubMedGoogle Scholar
  11. de Klerk B, Ducro BJ, Heuven HCM, den Uyl I, van Arendonk JAM, Parmentier HK, et al. Phenotypic and genetic relationships of bovine natural antibodies binding keyhole limpet hemocyanin in plasma and milk. J Dairy Sci. 2015;98:2746–52. https://doi.org/10.3168/jds.2014-8818.View ArticlePubMedGoogle Scholar
  12. Ploegaert TCW, Tijhaar E, Lam TJGM, Taverne-Thiele A, van der Poel JJ, van Arendonk JAM, et al. Natural antibodies in bovine milk and blood plasma: variability among cows, repeatability within cows, and relation between milk and plasma titers. Vet Immunol Immunopathol. 2011;144:88–94. https://doi.org/10.1016/J.VETIMM.2011.07.008.View ArticlePubMedGoogle Scholar
  13. Wagter LC, Mallard BA, Wilkie BN, Leslie KE, Boettcher PJ, Dekkers JC. A quantitative approach to classifying Holstein cows based on antibody responsiveness and its relationship to peripartum mastitis occurrence. J Dairy Sci. 2000;83:488–98.View ArticleGoogle Scholar
  14. Hernández A, Karrow N, Mallard BA. Evaluation of immune responses of cattle as a means to identify high or low responders and use of a human microarray to differentiate gene expression. Genet Sel Evol. 2003;35(Suppl 1):S67. https://doi.org/10.1186/1297-9686-35-S1-S67.View ArticlePubMedPubMed CentralGoogle Scholar
  15. Mallard BA, Emam M, Paibomesai M, Thompson-Crispi K, Wagter-Lesperance L. Genetic selection of cattle for improved immunity and health. Jpn J Vet Res. 2015;63(Suppl 1):S37–44 http://www.ncbi.nlm.nih.gov/pubmed/25872325.
  16. Thompson-crispi KA, Mallard BA. Type 1 and type 2 immune response profiles of commercial dairy cows in 4 regions across Canada. Can J Vet Res. 2012;76:120-8. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3314434/.
  17. Thompson-Crispi K, Miglior F, Mallard B. Incidence rates of clinical mastitis among Canadian Holsteins classified as high, average, or low immune responders. Clin Vaccine Immunol. 2013;20:106–12. http://www.ncbi.nlm.nih.gov/pubmed/23175290.
  18. van Knegsel ATM, de Vries Reilingh G, Meulenberg S, van den Brand H, Dijkstra J, Kemp B, et al. Natural antibodies related to energy balance in early lactation dairy cows. J Dairy Sci. 2007;90:5490–8. https://doi.org/10.3168/jds.2007-0289.View ArticlePubMedGoogle Scholar
  19. van Knegsel ATM, Hostens M, de Vries Reilingh G, Lammers A, Kemp B, Opsomer G, et al. Natural antibodies related to metabolic and mammary health in dairy cows. Prev Vet Med. 2012;103:287–97. https://doi.org/10.1016/J.PREVETMED.2011.09.006.View ArticlePubMedGoogle Scholar
  20. Ganeshan K, Chawla A. Metabolic regulation of immune responses. Annu Rev Immunol. 2014;32:609–34.View ArticleGoogle Scholar
  21. Loftus RM, Finlay DK. Immunometabolism: cellular metabolism turns immune regulator. J Biol Chem. 2016;291:1–10.View ArticleGoogle Scholar
  22. Martinez J, Verbist K, Wang R, Green DR. The relationship between metabolism and the autophagy machinery during the innate immune response. Cell Metab. 2013;17:895–900.View ArticleGoogle Scholar
  23. Cheng S-C, Joosten LAB, Netea MG. The interplay between central metabolism and innate immune responses. Cytokine Growth Factor Rev. 2014;25:707–13.View ArticleGoogle Scholar
  24. Banos G, Wall E, Coffey MP, Bagnall A, Gillespie S, Russell GC, et al. Identification of immune traits correlated with dairy cow health, reproduction and productivity. PLoS One. 2013;8:e65766. https://doi.org/10.1371/journal.pone.0065766.View ArticlePubMedPubMed CentralGoogle Scholar
  25. Ploegaert TCW, Wijga S, Tijhaar E, van der Poel JJ, Lam TJGM, Savelkoul HFJ, et al. Genetic variation of natural antibodies in milk of Dutch Holstein-Friesian cows. J Dairy Sci. 2010;93:5467–73. https://doi.org/10.3168/jds.2010-3264.View ArticlePubMedGoogle Scholar
  26. Wijga S, Bovenhuis H, Bastiaansen JWM, van Arendonk JAM, Ploegaert TCW, Tijhaar E, et al. Genetic parameters for natural antibody isotype titers in milk of Dutch Holstein-Friesians. Anim Genet. 2013;44:485–92. https://doi.org/10.1111/age.12038.View ArticlePubMedGoogle Scholar
  27. Thompson-Crispi KA, Miglior F, Mallard BA. Genetic parameters for natural antibodies and associations with specific antibody and mastitis in Canadian Holsteins. J Dairy Sci. 2013;96:3965–72. https://doi.org/10.3168/jds.2012-5919.View ArticlePubMedGoogle Scholar
  28. Fabregat A, Jupe S, Matthews L, Sidiropoulos K, Gillespie M, Garapati P, et al. The Reactome pathway knowledgebase. Nucleic Acids Res. 2018;46:D649–55. https://doi.org/10.1093/nar/gkx1132.View ArticlePubMedGoogle Scholar
  29. Milacic M, Haw R, Rothfels K, Wu G, Croft D, Hermjakob H, et al. Annotating cancer variants and anti-cancer therapeutics in reactome. Cancers (Basel). 2012;4:1180–211. https://doi.org/10.3390/cancers4041180.View ArticleGoogle Scholar
  30. Zerbino DR, Achuthan P, Akanni W, Amode MR, Barrell D, Bhai J, et al. Ensembl 2018. Nucleic Acids Res. 2018;46:D754–61. https://doi.org/10.1093/nar/gkx1098.View ArticlePubMedGoogle Scholar
  31. Jabara HH, Weng Y, Sannikova T, Geha RS. TRAF2 and TRAF3 independently mediate Ig class switching driven by CD40. Int Immunol. 2009;21:477–88. https://doi.org/10.1093/intimm/dxp013.View ArticlePubMedPubMed CentralGoogle Scholar
  32. Grammer AC, Swantek JL, McFarland RD, Miura Y, Geppert T, Lipsky PE. TNF receptor-associated factor-3 signaling mediates activation of p38 and Jun N-terminal kinase, cytokine secretion, and Ig production following ligation of CD40 on human B cells. J Immunol. 1998;161:1183–93 http://www.ncbi.nlm.nih.gov/pubmed/9686578.
  33. Sharma SK. New trends in telescopic remote Raman spectroscopic instrumentation. Spectrochim Acta A Mol Biomol Spectrosc. 2007;68:1008–22. https://doi.org/10.1016/j.saa.2007.06.047.View ArticlePubMedGoogle Scholar
  34. Kajita M, Magari M, Todo K, Kanayama N, Ohmori H. Conditional transformation of immunoglobulin mutation pattern from gene conversion into point mutation by controlling XRCC3 expression in the DT40 B cell line. J Biosci Bioeng. 2010;109:407–10. https://doi.org/10.1016/j.jbiosc.2009.09.050.View ArticlePubMedGoogle Scholar
  35. Devlin B, Roeder K. Genomic control for association studies. Biometrics. 1999;55:997–1004. https://doi.org/10.1111/j.0006-341X.1999.00997.x.View ArticlePubMedGoogle Scholar
  36. Porat-Shliom N, Milberg O, Masedunskas A, Weigert R. Multiple roles for the actin cytoskeleton during regulated exocytosis. Cell Mol Life Sci. 2013;70:2099–121. https://doi.org/10.1007/s00018-012-1156-5.View ArticlePubMedGoogle Scholar
  37. Jennewein MF, Alter G. The Immunoregulatory roles of antibody glycosylation. Trends Immunol. 2017;38:358–72. https://doi.org/10.1016/j.it.2017.02.004.View ArticlePubMedGoogle Scholar
  38. Panda S, Ding JL. Natural antibodies bridge innate and adaptive immunity. J Immunol. 2015;194:13–20. https://doi.org/10.4049/jimmunol.1400844.View ArticlePubMedGoogle Scholar
  39. Mannoor K, Xu Y, Chen C. Natural autoantibodies and associated B cells in immunity and autoimmunity. Autoimmunity. 2013;46:138–47. https://doi.org/10.3109/08916934.2012.748753.View ArticlePubMedGoogle Scholar
  40. Hooijkaas H, Benner R, Pleasants JR, Wostmann BS. Isotypes and specificities of immunoglobulins produced by germ-free mice fed chemically defined ultrafiltered &quot;antigen-free&quot; diet. Eur J Immunol. 1984;14:1127–30. https://doi.org/10.1002/eji.1830141212.View ArticlePubMedGoogle Scholar
  41. Fransen F, Zagato E, Mazzini E, Fosso B, Manzari C, El Aidy S, et al. BALB/c and C57BL/6 mice differ in Polyreactive IgA abundance, which impacts the generation of antigen-specific IgA and microbiota diversity. Immunity. 2015;43:527–40. https://doi.org/10.1016/j.immuni.2015.08.011.View ArticlePubMedGoogle Scholar
  42. 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. https://doi.org/10.1371/journal.pone.0131088.View ArticlePubMedPubMed CentralGoogle Scholar
  43. Wijga S, Bastiaansen JWM, Van der Poel JJ, Parmentier HK, van Arendonk JAM, Bovenhuis H. Genome-wide associations for natural antibodies in milk of Holstein-Friesians. (PhD thesis). Wageningen: Wageningen University; 2013. p. 37–64. ISBN: 978-84-6173-727-4Google Scholar
  44. Grönwall C, Silverman GJ. Natural IgM: beneficial autoantibodies for the control of inflammatory and autoimmune disease. J Clin Immunol. 2014;34(Suppl 1):S12–21. https://doi.org/10.1007/s10875-014-0025-4.View ArticlePubMedPubMed CentralGoogle Scholar
  45. van Altena SEC, Peen MA, van der Linden FH, Parmentier HK, Savelkoul HFJ, Tijhaar EJ. Bovine natural antibodies in antibody-dependent bactericidal activity against Escherichia coli and Salmonella typhimurium and risk of mastitis. Vet Immunol Immunopathol. 2016;171:21–7. https://doi.org/10.1016/j.vetimm.2016.01.009.View ArticlePubMedGoogle Scholar
  46. Berghof TVL, Arts JAJ, Bovenhuis H, Lammers A, van der Poel JJ, Parmentier HK. Antigen-dependent effects of divergent selective breeding based on natural antibodies on specific humoral immune responses in chickens. Vaccine. 2018;36:1444–52. https://doi.org/10.1016/j.vaccine.2018.01.063.View ArticlePubMedGoogle Scholar
  47. Jiang C, Foley J, Clayton N, Kissling G, Jokinen M, Herbert R, et al. Abrogation of lupus nephritis in activation-induced deaminase-deficient MRL/lpr mice. J Immunol. 2007;178:7422–31 http://www.ncbi.nlm.nih.gov/pubmed/17513793.
  48. Jiang C, Zhao M-L, Scearce RM, Diaz M. Activation-induced deaminase-deficient MRL/lpr mice secrete high levels of protective antibodies against lupus nephritis. Arthritis Rheum. 2011;63:1086–96. https://doi.org/10.1002/art.30230.View ArticlePubMedPubMed CentralGoogle Scholar
  49. Heriazon A, Thompson KA, Wilkie BN, Mathes-Sears W, Quinton M, Mallard BA. Antibody to ovalbumin and delayed-type hypersensitivity to Candida albicans and mycobacteria in lactating Holstein cows using Quil A or Freund’s complete adjuvant. Vet Immunol Immunopathol. 2009;127:220–7. https://doi.org/10.1016/J.VETIMM.2008.10.002.View ArticlePubMedGoogle Scholar
  50. Wiggans GR, Sonstegard TS, VanRaden PM, Matukumalli LK, Schnabel RD, Taylor JF, et al. Selection of single-nucleotide polymorphisms and quality of genotypes used in genomic evaluation of dairy cattle in the United States and Canada. J Dairy Sci. 2009;92:3431–6. https://doi.org/10.3168/jds.2008-1758.View ArticlePubMedGoogle Scholar
  51. Sargolzaei M, Chesnais JP, Schenkel FS. A new approach for efficient genotype imputation using information from relatives. BMC Genomics. 2014;15:478. https://doi.org/10.1186/1471-2164-15-478.View ArticlePubMedPubMed CentralGoogle Scholar
  52. VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91:4414–23. https://doi.org/10.3168/JDS.2007-0980.View ArticlePubMedGoogle Scholar
  53. Gilmour AR, Thompson R, Cullis BR. Average information REML: an efficient algorithm for variance parameter estimation in linear mixed models. Biometrics. 1995;51:1440. https://doi.org/10.2307/2533274.View ArticleGoogle Scholar
  54. Yang J, Zaitlen NA, Goddard ME, Visscher PM, Price AL. Advantages and pitfalls in the application of mixed-model association methods. Nat Genet. 2014;46:100–6. https://doi.org/10.1038/ng.2876.View ArticlePubMedPubMed CentralGoogle Scholar
  55. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodological). 1995;57:289–300. https://doi.org/10.2307/2346101.View ArticleGoogle Scholar
  56. Breuer K, Foroushani AK, Laird MR, Chen C, Sribnaia A, Lo R, et al. InnateDB: systems biology of innate immunity and beyond--recent updates and continuing curation. Nucleic Acids Res. 2013;41(Database issue):D1228–33. https://doi.org/10.1093/nar/gks1147.View ArticlePubMedGoogle Scholar

Copyright

© The Author(s). 2018

Advertisement