Amino Acid Variation in HLA Class II Proteins Is a Major Determinant of Humoral Response to Common Viruses.

The magnitude of the human antibody response to viral antigens is highly variable. To explore the human genetic contribution to this variability, we performed genome-wide association studies of the immunoglobulin G response to 14 pathogenic viruses in 2,363 immunocompetent adults. Significant associations were observed in the major histocompatibility complex region on chromosome 6 for influenza A virus, Epstein-Barr virus, JC polyomavirus, and Merkel cell polyomavirus. Using local imputation and fine mapping, we identified specific amino acid residues in human leucocyte antigen (HLA) class II proteins as the most probable causal variants underlying these association signals. Common HLA-DRβ1 haplotypes showed virus-specific patterns of humoral-response regulation. We observed an overlap between variants affecting the humoral response to influenza A and EBV and variants previously associated with autoimmune diseases related to these viruses. The results of this study emphasize the central and pathogen-specific role of HLA class II variation in the modulation of humoral immune response to viral antigens in humans.

Main Text

The humoral immune response plays an essential role in the control and prevention of viral infections in humans. It has long been known that concentrations of serum immunoglobulins vary from person to person, and antibody titers against prevalent viruses have been shown to be highly variable in the population.2, 3 A significant fraction of that variation is heritable,1, 3 yet little is known about the human genetic control and regulation of the immunoglobulin response to specific pathogens.

To investigate the impact of common human genetic variation on humoral immunity and to identify pathogen-specific variants associated with antibody response, we measured serum immunoglobulin G (IgG) levels against 14 common viruses (Table 1) in 2,363 immunocompetent adults of European ancestry (Figure S1) with available genome-wide genotype data, comprising 1,147 anonymized blood donors (62.0% male, mean age ± SD = 37.5 ± 13.2) and 1,216 individuals with psychiatric diagnoses (64.9% male, mean age ± SD = 40.6 ± 13.5) who were recruited for the Göttingen Research Association for Schizophrenia (GRAS).5, 6 All study participants provided informed consent, including consent for genetic testing, and the GRAS data collection has been approved by the ethical committee of the Georg-August-Universität Göttingen as well as by the respective local regulatories and ethical committees of all collaborating centers. All subject data were collected in accordance with ethical guidelines and the Helsinki Declaration.

Table 1

List of Analyzed Viruses and Seroprevalence

Family	Species	Name of Test or Epitope	N Seropositive (%)	IgG Levelsa
Herpesviridae	cytomegalovirus	P28	1,065 (45.1)	8.0 ± 1.1
P150	1,049 (44.4)	7.9 ± 1.0
Epstein-Barr virus	EBNA	2,162 (91.5)	8.9 ± 0.5
herpes simplex virus-1	Novagnost HSV1 IgG	1,473 (62.3)	10.4 ± 0.4
varizella zoster virus	Enzygnost Anti-VZV/IgG	2,304 (97.5)	13.7 ± 0.6
Orthomyxoviridae	influenza A virus	Novagnost Influenza A IgG	1,594 (67.5)	10.0 ± 0.4
influenza B virus	Novagnost Influenza B IgG	497 (21.0)	9.5 ± 0.2
Paramyxoviridae	measles virus	Enzygnost Anti-Measles Virus/IgG	2,177 (92.1)	15.1 ± 0.9
mumps virus	Enzygnost Anti-Parotitis Virus/IgG	1,862 (78.8)	14.6 ± 0.7
Parvoviridae	parvovirus B19	Novagnost Parvovirus B 19 IgG	1,664 (70.4)	10.9 ± 0.4
Polyomaviridae	BK polyomavirus	VP1	2,226 (94.2)	8.4 ± 1.0
JC polyomavirus	VP1	1,268 (53.7)	7.1 ± 0.8
Merkel cell polyomavirus	VP1	1,871 (79.2)	8.0 ± 1.0
trichodysplasia spinulosa-associated polyomavirus	VP1	2,020 (85.5)	8.9 ± 1.1
Togaviridae	rubella virus	Enzygnost Anti-Rubella Virus/IgG	2,216 (93.8)	11.4 ± 0.9

Log normalized, mean ± SD.

A list and description of all assays used for determination of IgG levels is provided in Table S1. We used multiplex serology on the Luminex platform, based on glutathione S-transferase (GST) fusion capture immunosorbent assays combined with fluorescent-bead technology, or commercially available ELISA-based Enzygnost or Novagnost assays (Siemens Healthcare Diagnostics). The latter were automatically processed on the BEP III System (Siemens Healthcare Diagnostics) and interpreted (according to the manufacturer’s instructions) as positive, negative, or borderline (the latter of was defined as negative for the purpose of statistics). For the Luminex-based assays, seronegativity was defined as the absence of detectable IgG.

Genome-wide SNP genotyping was performed with an Axiom myDesign genotyping array (Affymetrix) and was subject to stringent quality control steps as described previously. Imputation of unobserved genotypes was performed with the 1000 Genomes Project phase 1 v.3 haplotypes as a reference panel. Genotypes were pre-phased with MaCH v.1 and subsequently imputed by Minimac. SNPs with a reported r2 quality metric score < 0.8 or a minor allele frequency (MAF) < 5%, as well as reported markers on sex chromosomes, were excluded from downstream analyses. SNPs were also filtered on the basis of “missingness” (excluded if they had < 95% genotyping rate) and marked deviation from the Hardy-Weinberg equilibrium (excluded if p < 5 × 10−7). We then used logistic or linear regression models in PLINK v.1.9 to test for association between ∼six million SNPs and IgG response to 14 common human viruses (Table 1), using both a case-control study design (serostatus: antibody positive versus negative) and a continuous, quantitative approach (log-normalized IgG levels in seropositive samples). The first three principal components, calculated with GCTA (v.1.24), as well as sex and age, which affect humoral response phenotypes (Table S2), were included as covariates in all analyses. In case of significant differences in serostatus or IgG levels between healthy control individuals and individuals affected by neuropsychiatric disease (Table S3), psychiatric diagnosis was included as an additional binary covariate. We observed no evidence of residual inflation in any test statistic (λ = 0.99–1.04, Figure S2).

Correcting for the number of SNPs and viruses tested, we observed genome-wide significant signals (p < 3.57 × 10−9) in the human leucocyte antigen (HLA) class II region of the major histocompatibility complex (MHC) on chromosome 6 for influenza A virus, Epstein-Barr virus (EBV), JC polyomavirus (JCPyV), and Merkel cell polyomavirus (MCPyV) (Table 2, Figure S3). Full summary association results are available for download from Zenodo.

Table 2

Summary of Genome-Wide Significant Association Results

		Study Design: Serostatus			Study Design: IgG Levels
Influenza A	JCPyV	MCPyV	EBV
SNP	ID (coded allele)	rs140012631 (C)	rs9269910 (A)	rs9269268 (C)	rs6927022 (A)
p value	1.06 × 10−14	8.88 × 10−12	2.67 × 10−10	7.35 × 10−26
OR or betaa (95% CI)	2.02 (1.84,2.19)	1.74 (1.58,1.90)	1.53 (1.40,1.66)	0.16 (0.13,0.19)
Classical HLA allele	Allele	HLA-DQB1∗05:01	HLA-DQA1∗01:02	HLA-DQB1∗06:02	HLA-DRB1∗07:01
p value	4.91 × 10−12	4.11 × 10−9	1.35 × 10−9	1.01 × 10−14
OR or betaa (95% CI)	0.53 (0.44,0.63)	0.65 (0.56,0.75)	0.56 (0.46,0.67)	−0.17 (−0.21,−0.12)
Amino acid	Protein (position)	HLA-DRβ1 (96)	HLA-DRβ1 (133)	HLA-DRβ1 (13)	HLA-DRβ1 (11)
Omnibus p value	2.64 × 10−16	1.27 × 10−11	1.87 × 10−10	5.85 × 10−23

The strongest associated GWAS SNP, classical HLA allele, and amino acid position for each genome-wide significant virus is shown, along with corresponding p values and effect estimates. OR, odds ratio; CI, confidence interval.

OR is shown for the serostatus study design, and beta is shown for the IgG levels study design.

To fine map the associated region and pinpoint potentially functional variants, we imputed four-digit classical HLA alleles and variable amino acid positions in the HLA class I and II proteins by using SNP2HLA and the T1DGC Immunochip/HLA reference panel and tested these for association with IgG response (Tables S4 and S5). 101 HLA alleles and 200 amino acids had a MAF > 1% and were included in the analysis (r2 quality metric score: median = 0.99, interquartile range = 0.98–1), and we used a multi-degree-of-freedom omnibus test to test for association at multi-allelic amino acid positions with three or more possible states. For all viruses with genome-wide significant SNPs, we found corresponding associations with HLA class II alleles and amino acids (Table 2), and they explained most of the SNP association signals when included as covariates in conditional analyses (Table S6). The residual SNP significance observed in some of the conditional analyses suggests the existence of additional associated alleles or amino acids that did not reach genome-wide significance, although non-classical MHC associations or HLA imputation errors could also play a role. Among imputed amino acids, the strongest associations were observed for positions in HLA-DRβ1 (Table 2, Table S5). We therefore analyzed the potential impact of common HLA-DRβ1 haplotypes on IgG response to influenza A, EBV, JCPyV, and MCPyV (Figure 1). This analysis showed no consistency in the associations: the size and directionality of the effects were virus specific (Figure 1A). Most prominently, the amino acid haplotypes present in the classical alleles HLA-DRB1 or HLA-DRB1 were associated with influenza A seropositivity and higher anti-EBV IgG levels, but with JCPyV and MCPyV seronegativity. In contrast, the haplotype present in HLA-DRB1 were associated with MCPyV seropositivity and higher anti-EBV IgG levels, but with influenza A seronegativity.

Figure 1

Effect of Common HLA-DRβ1 Haplotypes on Virus Serostatus or Antibody Levels

(A) Estimated effects for HLA-DRβ1 haplotypes as defined by the strongest associated amino acid positions for influenza A virus, JCPyV, and MCPyV serostatus, as well as EBV IgG levels (Table 2). These four positions were imputed with SNP2HLA and all showed an imputation accuracy of > 99% in the original publication. The Val-His-Tyr-Arg-encoding haplotype with a frequency of 14.9%, present in classical HLA alleles HLA-DRB1 and HLA-DRB1, was chosen as reference (i.e., given an OR of 1, or a beta of 0, respectively) and not included in the figure. Only common haplotypes with a frequency of > 10% were included in the analysis and accounted for 93.1% of haplotype diversity. Diamonds designate estimated effect sizes; error bars define the 95% confidence interval. OR, odds ratio

(B) 3D model of the HLA-DR αβ heterodimer. The protein is shown in front (left) and side (right) views. The DR α chain is displayed in gray and the β chain in green. Associated amino acid positions, as selected for haplotype analysis, are highlighted. This figure was prepared with UCSF chimera, with Protein Data Bank code PDB: 4MCY.

We next explored the virus-specific associations in more detail. Anti-influenza A IgG antibodies were detected in approximately two thirds of study participants (Table 1), and two SNPs mapping to the HLA-DRB1 (MIM: 142857), DQA1 (MIM: 146880), and DQB1 (MIM: 604305) gene regions were identified as independently associated with influenza A seropositivity in the genome-wide SNP screen (Figure 2A). Amino acid position 96 in HLA-DRβ1 was associated more strongly with influenza A serostatus than any SNP or classical HLA allele (omnibus p = 2.6 × 10−16, Table 2) and explained most of the association signals of both independently associated SNPs when included in a conditional analysis. The residual significance of p = 9.5 × 10−4 for rs140012631 and p = 1.0 × 10−4 for rs9269912 might, however, imply that these SNPs tag additional classical or non-classical HLA effects that did not reach our significance threshold. Individuals carrying a glutamine at position 96 were more likely to have detectable levels of anti-influenza A IgG (odds ratio [OR] = 1.68), whereas the presence of a glutamic acid was associated with seronegativity (OR = 0.51, Figure 2B). These two amino acids are carried by the HLA-DRB1 and the HLA-DRB1 alleles, respectively. They are in almost perfect linkage disequilibrium (LD) with HLA-DQB1 and HLA-DQB1, respectively (Table S7), which were the most strongly associated classical HLA alleles (Figure 2C). This result is consistent with a recent study associating HLA-DRB1 with increased responsiveness to trivalent influenza vaccines in Hispanic children. Interestingly, the HLA-DRB1-HLA-DQB1 haplotype is found in nearly 100% of individuals with narcolepsy (MIM: 161400),17, 18 a condition with strong evidence for an autoimmune basis and with a possible, but still controversial, mechanistic link to influenza infection.19, 20 The same HLA class II haplotype has also been shown to be positively associated with multiple sclerosis (MIM: 126200), systemic lupus erythematosus (MIM: 152700),22, 23 and Goodpasture disease (MIM: 233450).

Figure 2

GWAS and HLA association results for influenza A serostatus

(A) Regional association plot of GWAS results. Values of –log10(p) are plotted by their position in the MHC genomic region on chromosome 6. The most significant association was observed for rs140012631 (p = 1.1 × 10−14, OR = 2.01 for the C allele). Accounting for the effect of rs140012631, we observed an independent association at rs9269912 (p = 8.5 × 10−10, OR = 0.57 for the T allele). A further step of forward conditioning left no genome-wide significant signal. The dashed horizontal line indicates the threshold for genome-wide significance (p = 3.57 × 10−9). The annotated dashed vertical lines indicate the positions of the classical HLA genes.

(B) Effect estimates for amino acid residues at position 96 of HLA-DRβ1 (omnibus p = 2.6 × 10−16). Designated classical HLA alleles contain the respective amino acid residue at the given position. OR, odds ratio; CI, confidence interval.

(C) Association results and conditional regression for classical HLA alleles. Conditioning on HLA-DQB1 did not reveal further independent associations (threshold: p = 3.57 × 10−9). PC, principal component; OR, odds ratio; CI, confidence interval.

More than 90% of study participants had detectable IgG against the EBV-encoded nuclear antigen-1 (EBNA-1), an antigen expressed in both latent and lytic modes of infection and that is essential for efficient EBV genome replication, persistence, and transcription in dividing cells. A large number of HLA SNPs were significantly associated with anti-EBNA-1 IgG levels (Figure S4A), consistent with the results of a previous study. Among amino acids, positions 11 and 26 of HLA-DRβ1 showed strong, independent associations (Table 2, Figure S4B), which together explained the top SNP association from the genome-wide association studies (GWASs) (conditional p = 0.04, Table S6). Among classical alleles, HLA-DRB1 showed the strongest association result (Figure S4C), and there was an independent association of HLA-DRB1. However, these two alleles could not completely explain the top GWAS hit (residual p = 2.6 × 10−11, Table S6), pointing to a possible relevance of additional classical HLA alleles that did not reach the significance threshold in our study. EBV is known or suspected to play a role in multiple sclerosis and systemic autoimmune diseases, and we observed an overlap between the genetic determinants of IgG response to EBV and the HLA class II variants that have been reported to associate with autoimmunity. Carriage of a glycine at HLA-DRβ1 position 11 was associated with decreased antibody levels (Figure S4B) and a lower risk of seropositive and seronegative rheumatoid arthritis (MIM: 180300). On the contrary, IgG levels were higher in individuals carrying a proline at position 11 (Figure S4B), which is in strong LD with alanine at position 71 (r2 = 0.83, Table S8) and carried by the HLA-DRB1 classical allele. HLA-DRβ1 position 71 and HLA-DRB1 have been reported as the strongest genetic risk factors for multiple sclerosis at the level of amino acids and classical HLA alleles, respectively. However, the directionality of immune response and autoimmune disease risk associations is not always consistent, and larger study cohorts will be required to determine the role and extent of the contribution of EBV immune response to autoimmune disease susceptibility.

We also observed significant SNP associations in the HLA class II region for JCPyV and MCPyV serostatus and for MCPyV IgG levels (Table 2, Figures S5A, S6A, and S7). HLA-DRβ1 positions 133 and 13 showed the strongest amino acid associations with JCPyV and MCPyV serostatus, respectively (Figures S5B and S6B), whereas HLA-DQA1 and HLA-DQB1 were associated with JCPyV and MCPyV seronegativity, respectively (Figures S5C and S6C). These two alleles are in strong LD with HLA-DRB1 (Table S7), which was found to be associated with JCPyV seronegativity previously. No classical HLA alleles or amino acids reached the significance threshold for MCPyV IgG levels.

Genome-wide significant HLA class II amino acids explained between 2.5% (JCPyV serostatus) and 5.4% (EBV anti-EBNA-1 IgG levels) of the phenotypic variance (Table S9). With our sample size, we had > 80% power to detect variants with a MAF of 10% that (1) explain at least 2.5% of the trait variance (IgG levels) in the case of 80% seropositivity or (2) have an OR of 1.5 for a discrete trait with 67% seropositivity.

Although the continuous phenotype—IgG levels in seropositive individuals—is a direct indicator of the intensity of the humoral immune response, the interpretation of the dichotomous phenotype—serostatus—is not as straightforward. Seronegativity can indeed reflect both non-exposure and IgG levels under the detection limit. A correlation analysis of ORs of the dichotomous trait and the respective betas of the continuous trait showed highly significant results (Table S10), suggesting that seronegativity results from very low IgG levels rather than from lack of exposure in a subset of the study population. Thus, we sought to maximize power by using the largest possible sample (i.e., by including seronegative individuals in the dichotomous analyses), particularly for viruses with a large seronegative fraction.

In summary, we here present a cross-pathogen, genome-wide investigation of the role of host genetics in modulating the individual IgG response to common viral antigens. We identified a strong association between HLA class II variation, in particular HLA-DRβ1 amino acid variants, and IgG response to several prevalent viruses. Our study illustrates the value of GWASs in infectious disease research and emphasizes the importance of fine mapping SNP association signals in the MHC via imputation of potentially functional HLA variation; doing so offers the unique possibility for a straightforward investigation of putative causal variants. The pathogen-specific patterns of association that we observed suggest that small differences in the capacity of HLA-DRβ1 proteins to bind specific viral peptides and present them efficiently to CD4+ T cells can have a measurable impact on downstream antibody production. The concordance in associated genetic loci between IgG response to several common viral antigens and antibody-dependent autoimmune diseases observed in this study is intriguing. Although our results do not allow for the conclusion of a direct causative role of these viruses in autoimmune pathogenesis, they strongly encourage further genetic and functional work. Understanding how genetic differences affect human response to infection might ultimately inform drug and vaccine development and facilitate personalized preventive and therapeutic approaches.