Skip to main content

Treatment- and population-specific genetic risk factors for anti-drug antibodies against interferon-beta: a GWAS

Abstract

Background

Upon treatment with biopharmaceuticals, the immune system may produce anti-drug antibodies (ADA) that inhibit the therapy. Up to 40% of multiple sclerosis patients treated with interferon β (IFNβ) develop ADA, for which a genetic predisposition exists. Here, we present a genome-wide association study on ADA and predict the occurrence of antibodies in multiple sclerosis patients treated with different interferon β preparations.

Methods

We analyzed a large sample of 2757 genotyped and imputed patients from two cohorts (Sweden and Germany), split between a discovery and a replication dataset. Binding ADA (bADA) levels were measured by capture-ELISA, neutralizing ADA (nADA) titers using a bioassay. Genome-wide association analyses were conducted stratified by cohort and treatment preparation, followed by fixed-effects meta-analysis.

Results

Binding ADA levels and nADA titers were correlated and showed a significant heritability (47% and 50%, respectively). The risk factors differed strongly by treatment preparation: The top-associated and replicated variants for nADA presence were the HLA-associated variants rs77278603 in IFNβ-1a s.c.- (odds ratio (OR) = 3.55 (95% confidence interval = 2.81–4.48), p = 2.1 × 10−26) and rs28366299 in IFNβ-1b s.c.-treated patients (OR = 3.56 (2.69–4.72), p = 6.6 × 10−19). The rs77278603-correlated HLA haplotype DR15-DQ6 conferred risk specifically for IFNβ-1a s.c. (OR = 2.88 (2.29–3.61), p = 7.4 × 10−20) while DR3-DQ2 was protective (OR = 0.37 (0.27–0.52), p = 3.7 × 10−09). The haplotype DR4-DQ3 was the major risk haplotype for IFNβ-1b s.c. (OR = 7.35 (4.33–12.47), p = 1.5 × 10−13). These haplotypes exhibit large population-specific frequency differences. The best prediction models were achieved for ADA in IFNβ-1a s.c.-treated patients. Here, the prediction in the Swedish cohort showed AUC = 0.91 (0.85–0.95), sensitivity = 0.78, and specificity = 0.90; patients with the top 30% of genetic risk had, compared to patients in the bottom 30%, an OR = 73.9 (11.8–463.6, p = 4.4 × 10−6) of developing nADA. In the German cohort, the AUC of the same model was 0.83 (0.71–0.92), sensitivity = 0.80, specificity = 0.76, with an OR = 13.8 (3.0–63.3, p = 7.5 × 10−4).

Conclusions

We identified several HLA-associated genetic risk factors for ADA against interferon β, which were specific for treatment preparations and population backgrounds. Genetic prediction models could robustly identify patients at risk for developing ADA and might be used for personalized therapy recommendations and stratified ADA screening in clinical practice. These analyses serve as a roadmap for genetic characterizations of ADA against other biopharmaceutical compounds.

Peer Review reports

Background

Interferon β (IFNβ) preparations are a treatment option for multiple sclerosis (MS). IFNβ-1a is produced in Chinese hamster ovary cells and administered either via intramuscular (i.m.) or subcutaneous (s.c.) injection. IFNβ-1b is raised using Escherichia coli and injected subcutaneously. The amino acid sequence of IFNβ-1b differs at two positions from the mammalian protein [1]. Moreover, IFNβ-1b is not glycosylated, which may affect its immunogenicity, e.g., by promoting the formation of protein aggregates [1, 2]. Posttranslational modifications like deamidation, oxidation, and glycation can also occur spontaneously, depending on the manufacturing and processing of biopharmaceuticals [3]. Therefore, also sequence-identical compounds like IFNβ-1a s.c. and i.m. can differ in their immunogenicity [1, 4].

Up to 40% of patients treated with IFNβ develop anti-drug antibodies (ADA) that bind IFNβ (binding ADA, bADA) [1, 5,6,7]. A subset of bADA inhibits the interaction of IFNβ with its receptor and thus neutralizes the drug’s biological activity (neutralizing ADA, nADA) [8, 9]. Previous studies have already identified genetic factors influencing the development of ADA but could not establish a consensus on the human leukocyte antigen (HLA) alleles [10,11,12,13,14,15,16,17] and single nucleotide polymorphisms (SNPs) [14, 15] contributing to ADA development.

The primary aim of the present, retrospective study was to characterize the contribution of genetic risk to ADA development by analyzing a large, cross-sectional sample from two different sites: the Karolinska Institutet Stockholm, Sweden (KI), and the Technical University of Munich, Germany (TUM). In these analyses, it was an objective to establish a consensus on the heterogeneous findings from previous studies, especially regarding the associations of HLA alleles. Both bADA levels and nADA titers were determined in the same patients, allowing for systematic comparisons between the two antibody types. Genome-wide association studies (GWAS) on bADA levels, nADA titers, and nADA presence, as well as analyses of the association of imputed HLA alleles with ADA, were conducted. As primary analyses, results were pooled across treatments; as secondary analyses, treatment-specific results were evaluated. The secondary aim of the study was to use these genetic factors for the prediction of ADA development.

Methods

Sample inclusion criteria

Patient inclusion criteria of this retrospective study were as follows: diagnosis of either clinically isolated syndrome (CIS) or multiple sclerosis (MS), age at first treatment with IFNβ ≥ 18 years, availability of genotype data, and a serum sample fulfilling the sample inclusion criteria. Patients were diagnosed using the current McDonald criteria at the time of diagnosis. The sample inclusion criteria for bADA-/nADA-negative samples were as follows: ≥ 12 months of treatment with IFNβ; if more than one sample was eligible, the first sample available at least 12 months after initiation of treatment with IFNβ was selected; and no previous positive screening for bADA or nADA. The sample inclusion criteria for previously bADA-/nADA-positive samples were as follows: ≥ 6 months of treatment with IFNβ; if previously treated with an IFNβ product, not having been ADA-positive during a previous IFNβ treatment period; and if more than one sample was eligible, the first sample available at least 6 months after initiation of treatment with IFNβ was selected. Based on these criteria, 1810 patients were eligible at KI and 1488 at TUM. The respective local ethics committees approved the study, and all participants provided written informed consent.

Power calculation

In a previous ADA GWAS, Weber et al. identified a genome-wide significant SNP explaining 2.5% of the variance of bADA levels [14]. To have sufficient power for identification of additional associated variants, 2000 patients were assigned to the discovery-stage analyses. In this dataset, 80% of statistical power can be reached for a variant explaining 1.96% of the variance at a p value of 5 × 10−8 (calculated using the R package pwr). Effect sizes in the replication stage are expected to be smaller than in the discovery stage [18]. We thus estimated that at least 682 patients are necessary for replicating up to ten linkage disequilibrium (LD)-independent signals with 80% power, explaining 1.7% of the variance using a one-sided hypothesis. Because of an expected reduction in power due to heterogeneity and an expected decrease in the number of available samples after titration and quality control (QC), we initially selected 800 patients for the replication stage.

Selection of patients

To select approximately 2800 patients for ADA screening and titration, all available previously bADA-/nADA-positive samples (n = 984) were combined with previously ADA-negative samples (n = 2314) best-matching ADA-positive ones (Additional file 1). Propensity score matching was conducted using the R package optmatch [19], based on recruitment site, gender, the age at the blood draw, the IFNβ treatment preparation, the total duration of IFNβ treatment, and eight multi-dimensional scaling (MDS) ancestry components of the genetic identity-by-state (IBS) matrix, calculated from the genotype data to account for population stratification (Additional file 2). From the selected patients, new bADA levels and nADA titers (see below) could be determined for 938 previously bADA-/nADA-positive and 1819 previously ADA-negative samples (Table 1 and Additional file 3). These patients were randomized into a discovery (n = 2000), and a replication (n = 757) set, using adaptive randomization to minimize differences regarding recruitment site, nADA measurement site (Innsbruck or Copenhagen, see below), gender, the age at the blood draw, the IFNβ treatment preparation, and the total duration of IFNβ treatment.

Table 1 Sample characteristics

ADA screening and titration

Binding ADA levels were measured by capture ELISA [20] at a single site (Munich) and were calculated from optical densities using a standard curve (Additional file 2). For the assessment of nADA titers, measured as the inverse of serum dilutions using a luciferase-based bioassay [21], samples were first screened, and titration was only conducted for samples positive during screening [22]. Assessment of nADA titers was conducted at two separate sites (Innsbruck and Copenhagen), to which samples were assigned using adaptive randomization to minimize differences regarding the recruitment site, gender, the age at the blood draw, the IFNβ treatment preparation, and the total duration of IFNβ treatment. We obtained 2748 valid measurements for nADA screening and titers as well as 2752 bADA levels; for 2743 patients; both nADA titers and bADA levels were available (1990 in the discovery and 753 in the replication set). The presence of nADA was defined as samples positive in the screening for nADA and showing a nADA titer ≥ 40 tenfold reduction units per milliliter. Correlations of bADA and nADA were calculated in a combined dataset of all samples. For the estimation of the nADA status from bADA levels, the cutoff was established using nested cross-validation in the discovery dataset (Additional file 2). Sensitivity and specificity were calculated by the application of this cutoff to the replication data.

Genotyping and imputation

SNPs were genotyped on Illumina microarrays, and QC was conducted separately for KI and TUM data in PLINK v1.90b3.44 or higher [23], as described before [24]. Genotype data were imputed to the 1000 Genomes Phase 3 reference panel using SHAPEIT2 and IMPUTE2 [25,26,27]. The resulting datasets contained 9,096,778 and 8,550,834 high-quality variants with a MAF ≥ 1% for KI and TUM, respectively. HLA allele imputation was performed using SNP2HLA v1.0.3/Beagle v3.04 and the Type 1 Diabetes Genetics Consortium imputation panel, as previously described [28,29,30]. The extended haplotypes were determined based on the haplotype phasing estimated in Beagle. An additional file provides further details on QC and imputation (Additional file 2).

Estimation of heritability and GWAS

ADA titers/levels were transformed by rank-based inverse normal transformation before analyses. Sex, age, treatment preparation, treatment duration, titration site, and eight ancestry components were used as covariates in all analyses. The covariate treatment preparation was also used in preparation-specific analyses and controlled, beyond the three preparation types, for (a) whether treatment with IFNβ-1a s.c. had begun before 2008 (change of the formulation [31]); (b) in the TUM cohort, the dose of IFNβ-1a s.c. (22 vs. 44 μg); and (c) in the KI cohort, the IFNβ-1b s.c. brand used.

The SNP heritability and genetic correlations were estimated with GCTA GREML on a combined dataset of KI and TUM genotypes [32,33,34,35], using the covariates mentioned above plus treatment preparation and the recruitment site.

GWAS were conducted separately for the presence of nADA, nADA titers, and bADA levels. ADA titers/levels were analyzed by linear regression models, the presence of nADA by logistic regression. GWAS were run stratified by cohort (KI Sweden and TUM Germany) and by treatment preparation (IFNβ-1a i.m., IFNβ-1a s.c., IFNβ-1b s.c.). For each treatment preparation, samples from Sweden and Germany were analyzed separately in PLINK; GWAS results were pooled per cohort using fixed-effects meta-analysis in METAL [36]. For plots of the ancestry components in both cohorts, see Additional file 4. In the primary analysis (GWAS across treatment preparations), the three treatment groups were subsequently pooled by fixed-effects meta-analysis. The threshold for genome-wide significance was α = 5 × 10−8. For replication, the significance threshold α was corrected for the total number of variants analyzed across all SNP-based analyses in the replication phase (n = 16) using Bonferroni’s method, i.e., α = 0.05/16 = 3 × 10−3. SNPs prioritized for replication had to fulfill the following criteria: (I) genome-wide significance (p < 5 × 10−8) in the discovery-stage GWAS; (II) within each window of 100,000 bp, only the SNP with the lowest p value was selected; and (III) LD with SNPs showing lower p values had to be r2 < 0.2 in each cohort. Although we used a study design involving discovery and replication, for completeness, association results in the pooled complete dataset are reported as secondary results. More details are provided in Additional file 2.

Permutation analyses

All replicated associations from hypothesis-free linear regression analyses were validated using permutation analyses. In these analyses, the null distribution of test statistics was empirically determined by repeating regression analyses either 200 million or 1 million times with random sampling of phenotype data. To calculate a p value, the number of tests was counted where a model with a random genotype-phenotype association showed the same or a more extreme p value than the correct, non-randomized model; this number was divided by the total number of tests (200 or 1 million). Permutation-based p values were pooled per cohort and treatment using Stouffer’s Z-score method [37]. For GWAS variants, 200 million permutations per dataset (discovery/replication), cohort, and treatment preparation were carried out (allowing for p values down to 1 × 10−8); for stepwise conditional models and HLA alleles, the default was 1 million permutations per group (for p values down to 1 × 10−6). If these permutation p values were < 1× 10−6, 200 million permutations were conducted. If the permutation p values were < 1 × 10−8, they were set to 1 × 10−8.

EQTL analyses

The significant cis-expression quantitative trait loci (eQTLs) in whole blood were looked up in the GTEx v8 database (https://gtexportal.org/) downloaded on April 1, 2020 (dbGaP accession number phs000424.v8.p2) [38].

Gene-set analyses

Gene-set analyses were conducted with MAGMA v1.07b [39]. First, SNPs within gene boundaries were annotated to RefSeq genes (0 bp window). Second, gene analysis was performed on the pooled GWAS summary statistics, based on LD information from the 1000 Genomes EUR reference panel, using both mean- and top-SNP gene models. Third, gene-level analyses used a combination of the curated 186 KEGG and 1499 Reactome pathways from the MSigDB 7.0 database gene sets [40].

HLA and stepwise conditional analyses

The association of HLA alleles was analyzed in R v3.3 or higher. As in the GWAS, sex, age, treatment preparation and duration, titration site, and eight ancestry components were used as covariates. Separate regression models were run per cohort and treatment preparation, followed by a two-level meta-analysis: results were combined using fixed-effects meta-analysis first by cohort and then by preparation. For assessment of significance, we applied Bonferroni correction for testing 131 alleles and extended haplotypes [41] (rounded down to α = 3 × 10−4). In the replication phase, we corrected for multiple testing of 41 HLA alleles and haplotypes prioritized across all analyses (rounded down to α = 1 × 10−3). Note that the associations of all HLA alleles and haplotypes presented in this study also reached genome-wide significance (p < 5 × 108) in the pooled analyses of discovery and replication samples, except for the super-extended haplotypes C7-DQ6 and A3-DQ6, which reached a p < 10 × 10−8.

Stepwise conditional regression was conducted, as previously described [42, 43], first only for HLA alleles and then for a joint dataset of HLA alleles and SNPs mapping to the extended MHC region. In brief, the association of all alleles/SNPs was first tested in separate regression models. The top-associated allele/SNP was then added as a covariate to the regression model, and the analysis was repeated for all remaining alleles/SNPs. This addition of top-associated alleles/SNPs as covariates was repeated until no allele/SNP was significant anymore after correction for multiple testing.

Polygenic risk scores (PRS) and prediction of nADA

PRS were calculated in R v3.33 using imputed genetic data, as described previously [44, 45]. For each PRS, the effect sizes of variants from the discovery-stage analyses (training data), below a selected discovery-stage p value threshold, were multiplied by the imputed SNP dosage in the replication-stage test data and then summed to produce a single PRS per threshold. For each analysis group, eight PRS based on different GWAS p value thresholds were calculated on the discovery data. More details are provided in Additional file 2.

For the prediction of the presence of nADA in the replication dataset, logistic regression of the eight PRS, the top single GWAS variant, and the top HLA allele from the discovery stage was conducted using the GWAS models. The area under the receiver operating characteristic curve (AUC) was calculated using the R package pROC, its 95% confidence interval (CI) with the function ci.auc (2000 stratified bootstrap replicates). At this stage, we adapted the significance threshold for ten tests using the Bonferroni correction (Fig. 3a, b). For each treatment preparation, the model with the highest AUC was selected. The performance of all top models was subsequently compared, with a significance threshold adapted for 160 comparisons (α = 3.13 × 10−4, Fig. 3c).

The sensitivity and specificity of the predictions were calculated using the package OptimalCutpoints, maximizing both measures (MaxSpSe). To avoid overfitting, the cutpoint was selected using nested cross-validation with three outer and four inner folds. In each outer cross-validation instance, the cutoff producing the maximum sensitivity across three inner cross-validation folds was tested on the remaining fold. Nested cross-validation was repeated 100 times, and the mean cutoff of the 100 repetitions was used as the final cutoff. Nagelkerke’s pseudo-R2 was calculated using the package fmsb. For a comparison of patients either at low or high genetic risk, patients within the lower 30% of genetic risk were compared to the patients in the upper 30%. We initially selected a 10% cutoff for this contrast and increased it in 10% steps until the sample size in the replication dataset sufficed for the stable convergence of regression models.

Results

From 2757 MS patients recruited in Sweden and Germany and treated with three different IFNβ preparations (Table 1), bADA levels were measured by capture ELISA [20] and nADA titers using a luciferase-based bioassay [21, 22] (Additional files 1-3). The bADA levels were correlated with nADA presence (Spearman ρ = 0.66) and nADA titers (ρ = 0.71). Compared to the presence of nADA determined via screening and titration, estimation of the nADA status from bADA levels had a sensitivity = 0.85 and a specificity = 0.84 (Fig. 1).

Fig. 1
figure1

Comparison of bADA levels and nADA titers. The orange dashed line indicates the cutoff at a log10 bADA level of 1.442256, which optimized the maximum sensitivity and specificity in the discovery data. This cutoff had a sensitivity = 0.83 and a specificity = 0.82 in the discovery dataset; a sensitivity = 0.85 and a specificity = 0.84 in the replication dataset; and a sensitivity = 0.84 and a specificity = 0.82 in the combined dataset. a Density plot showing log10 bADA levels in the combined discovery and replication dataset stratified by nADA presence. b Density plot showing log10 bADA levels in the replication data stratified by nADA presence. c Density plot showing rank-based inverse-normal transformed bADA levels in the combined dataset stratified by nADA presence. d Comparison of log10 bADA levels to log10 nADA titers in the combined dataset, colored by nADA presence

SNP heritability and genetic correlations

The SNP-based heritability estimated from the genotype data was h2g = 0.47 (standard error (SE) = 0.15, p = 1.4 × 10−4) for the inverse-normal transformed bADA levels and h2g = 0.50 (SE = 0.15, p = 2.9 × 10−4) for the transformed nADA titers. The SNP heritability of the presence of nADA was h2go = 0.48 on the observed scale (SE = 0.15, p = 4.9 × 10−4) and, assuming an incidence of 0.35 for ADA, h2gl = 0.79 (SE = 0.25) on a liability scale. Genetic correlations of bADA levels with nADA presence (rg = 0.89, SE = 0.14, p = 1.2 × 10−3) and titers (rg = 0.95, SE = 0.11, p = 7.0 × 10−4) were very high.

Outline of the genetic association analyses

To improve control for type I errors, patients were randomized a priori into a discovery (n = 2000) and a replication (n = 757) set for all genetic association analyses (Additional file 3). The sizes of both datasets were guided using a power calculation (see the “Methods” section). We conducted three separate analyses: first, a pooled analysis of all three treatment preparation groups; second, an analysis of patients treated with IFNβ-1a s.c.; and third, an analysis of patients treated with IFNβ-1b s.c. For each of these three analysis levels, we conducted separate GWAS of nADA presence, nADA titers, and bADA levels. Because of the small number of IFNβ-1a i.m.-treated patients with ADA in the present study, we did not analyze this treatment preparation on its own. In addition to the GWAS, we analyzed imputed HLA alleles in the same manner. To estimate the number of independent association signals in the major histocompatibility complex (MHC) region, we carried out conditional analyses in a combined dataset of GWAS variants and HLA alleles.

In all GWAS, only variants within the MHC region were significant on a genome-wide scale (p < 5 × 10−8) in the discovery-stage analyses and replicated (Additional files 5-10). There was no indication for systematic inflation of test statistics; all genomic inflation factors λ were in the expected range (Additional file 11).

GWAS across IFNβ preparations

In the discovery GWAS of nADA presence across all three treatment preparation groups, the strongest association was observed for the insertion TTTTTTT of the variant rs9281971, which was associated with decreased risk for nADA (Table 2 and Additional files 12 and 13). This insertion had a frequency of 36.0% in Swedish and 38.6% in German patients. No other genome-wide significant variant with linkage disequilibrium (LD) r2 < 0.2 with the top signal was identified. The insertion replicated at genome-wide significance and was also the top association signal when pooling discovery and replication data (discovery: odds ratio (OR) = 0.59 (95% confidence interval (CI) = 0.50–0.69), p = 1.9 × 10−11; replication: OR = 0.44 (0.32–0.59), pone-sided = 2.4 × 10−08; discovery + replication: OR = 0.55 (0.48–0.63), p = 2.3 × 10−17). Inversely, for each copy lacking the insertion TTTTTTT at this site, the OR associated with risk for nADA was thus 1.82 (1.58–2.08). This association was supported in all three treatment groups but was strongest in IFNβ-1a s.c.-treated patients (Additional file 14). Because the MHC region shows long-range LD patterns, we conducted stepwise conditional regression analyses in the pooled dataset (discovery + replication) to estimate the number of independently associated signals. No variant except rs9281971-TTTTTTT was significantly associated with nADA presence on a genome-wide scale (Additional file 15).

Table 2 Genome-wide significant variants from GWAS across IFNβ preparations

In addition to the analysis of the dichotomous nADA presence, we also conducted GWAS for the quantitative measures nADA titers and bADA levels. In both cases, the top-associated and replicated SNP was rs9271377 (Table 2). As a secondary analysis, we pooled the discovery and replication GWAS. Here, variant rs9281971-TTTTTTT was the strongest association for nADA titers (9558 base pairs (bp) downstream from rs9271377; LD r2 = 0.39 in Swedish and r2 = 0.51 in German patients). For bADA levels, SNP rs9272071 (2759 bp downstream from rs9271377; r2 = 0.77 in Swedish and r2 = 0.81 in German patients) was the top pooled association. These associations were supported in all three treatment groups, but to a lesser degree in IFNβ-1a i.m.-treated patients (Additional file 14). Because of possible deviations from normality, the associations of all replicated ADA variants were confirmed using nonparametric permutation analyses (Additional file 12). In stepwise conditional regression analyses in the pooled dataset, rs9281971-TTTTTTT was the only significant variant for nADA titers, while four variants reached significance for bADA levels (rs9272071, rs28746882, rs1265086, and HLA-DRB1*04:04 (Additional file 13)).

The three top-associated variants map directly upstream of the gene HLA-DQA1 (rs9271377 18.0 kbp, rs9281971 8.5 kbp, and rs9272071 5.7 kbp upstream (Additional file 13)). The variants were all in weak to moderate LD with the HLA allele DQA1*05:01 (LD range: 0.31 ≥ r2 ≤ 0.50) and part of cis-expression quantitative trait loci (eQTLs) with HLA-DRB5 (GTEx v8, see Additional file 16) [38]. In gene-set analyses using KEGG and Reactome gene sets [39, 40], several immune-related pathways were significant after correction for multiple testing, e.g., “antigen processing and presentation,” “Translocation of ZAP-70 to Immunological synapse,” and “PD-1 signaling” (Additional file 17).

Treatment-specific GWAS: IFNβ-1a s.c.

Forest plots of effect sizes in treatment preparation subgroups suggested that preparation-specific genetic risk factors may exist (Additional file 10). Therefore, we also conducted analyses separately for the two main treatment preparations used in our cohorts (IFNβ-1a s.c. and IFNβ-1b s.c.). We did not conduct hypothesis-free preparation-specific analyses for IFNβ-1a i.m., due to its much lower number of ADA-positive patients (Table 1).

In GWAS of nADA presence and nADA titers in IFNβ-1a s.c.-treated patients, variant rs77278603 was genome-wide significant in the discovery stage and confirmed in the replication dataset (OR = 3.55 (2.81–4.48), p = 2.1 × 10−26; Table 3 and Additional file 12). The variant maps downstream of HLA-DRB5 (Additional file 18). In the secondary meta-analysis of discovery and replication GWAS, two different but correlated SNPs upstream of HLA-DQA1 were the most strongly associated signals, rs9271700 for nADA presence and rs9271673 for nADA titers (Table 3 and Additional file 12). Both variants were significant eQTLs with an HLA-DRB5 transcript (Additional file 16). All three SNPs associated with nADA in IFNβ-1a s.c.-treated patients were in LD with the HLA allele HLA-DRB1*15:01 (r2 ≥ 0.71).

Table 3 Genome-wide significant variants from treatment-specific GWAS for IFNβ-1a s.c

For bADA levels in IFNβ-1a s.c.-treated patients, rs9281962 was the top-associated variant in both the discovery and pooled analysis (Table 3 and Additional file 12), which was in very high LD with the nADA-associated SNPs rs9271700 and rs9271673 (r2 ≥ 0.93). In stepwise conditional regression analyses in the pooled dataset, none but the respective top-associated variants were significantly associated with IFNβ-1a s.c.-induced ADA (Additional file 15).

Notably, none of the variants associated at genome-wide significance with ADA measurements in IFNβ-1a s.c.-treated patients showed statistical support for an association in IFNβ-1b s.c.-treated patients with p < 0.001 in the discovery stage (Additional files 12 and 19). When analyzing both IFNβ-1a preparations together, results were highly similar to when analyzing IFNβ-1a s.c.-treated patients alone (Additional file 12).

Treatment-specific GWAS: IFNβ-1b s.c.

In IFNβ-1b s.c.-treated patients, SNP rs28366299 was significantly associated with nADA presence in both the discovery and pooled analysis (OR = 3.56 (2.69–4.72), p = 6.6 × 10−19; Table 4 and Additional file 12). It maps upstream of HLA-DRB1, which is correlated with HLA-DQA1*03:01 (r2 ≥ 0.46) and an eQTL with an HLA-DQA2 transcript (Additional files 16 and 20). We confirmed the association of this SNP with nADA presence in a published, independent study on 941 IFNβ-1b s.c.-treated patients [15], where it replicated robustly (OR 2.37 (1.81–3.08), one-sided p = 9.88 × 10−11; meta-analysis with the present study: OR 2.87 (2.37–3.48), p = 7.74 × 10−27). The same variant was also associated with nADA titers (Table 4 and Additional file 12). SNP rs9272775, intronic in HLA-DQA1 and correlated with rs28366299 (r2 ≥ 0.79), was the top variant in the pooled analysis of nADA titers. In the independent study [15], variant rs9272775 replicated with one-sided p = 6.05 × 10−17 (meta-analysis p = 7.62 × 10−40).

Table 4 Genome-wide significant variants from treatment-specific GWAS for IFNβ-1b s.c

The analysis of bADA levels in IFNβ-1b s.c.-treated patients produced similar results (Table 4 and Additional file 12). In the discovery stage, rs78279385 (LD with rs9272775 r2 ≥ 0.88) was the variant showing the most robust support for an association. In the meta-analysis of discovery and replication GWAS, rs9272775 was the top-associated SNP. None of the variants associated at genome-wide significance with any ADA measurement in IFNβ-1b s.c.-treated patients showed statistical support for an association in IFNβ-1a s.c.-treated patients with p < 0.009 in the discovery stage (Additional files 12 and 21).

In stepwise conditional regression analyses, rs559242105, in LD with HLA-DPB1*03:01 (r2 ≥ 0.78), was identified as a secondary signal for nADA presence and titers when conditioning for the respective top SNP (Additional file 15). In stepwise conditional analyses of bADA levels, the allele HLA-DRB1*04:01 reached a lower p value (β = 0.61, p = 8.3 × 10−20) than the top GWAS variant rs9272775 (β = 0.41, p = 1.6 × 10−19). Note that this constitutes the only analysis where an HLA allele reached a lower p value than the best available SNP. In this analysis, variant rs17205731 was identified as a significant secondary signal. The two secondary variants rs559242105 (nADA presence/titers) and rs17205731 (bADA levels) are not in LD with each other; both map to HLA-associated loci further downstream than the primary association signal (Additional file 20). The variant rs559242105 was independent of HLA-DRB1*04:01 (nADA titers without conditioning for HLA-DRB1*04:01: p = 5.2 × 10−10, conditioned for HLA-DRB1*04:01: 1.1 × 10−08). Possibly, the secondary signal of rs559242105 corresponds to the protective association observed for HLA-DPB1*03:01 (conditional model of HLA-DPB1*03:01 OR = 0.49 (0.35–0.69), p = 4.7 × 10−05) (Additional file 15).

Analysis of HLA variants across IFNβ preparations

Most previous studies have not conducted GWAS but instead analyzed the association of HLA alleles with ADA. Therefore, we also conducted a dedicated association analysis of imputed HLA alleles with ADA (Table 5). In this secondary analysis, the discovery-stage significance threshold was set to α = 3 × 10−4, corresponding to the Bonferroni correction for 131 analyzed alleles and extended haplotypes. The full list of significantly associated HLA alleles and haplotypes is shown in (Additional file 22). In analyses of nADA presence and nADA titers across all treatment preparations, no HLA allele was significant after correction for multiple testing and replicated.

Table 5 Selected significant HLA alleles and haplotypes

For bADA levels, the top-associated and replicated HLA risk allele across preparations was HLA-DQA1*01:02 (p = 4.79 × 10−11). This allele is part of the extended ancestral haplotype B7-DQ6 (the combined presence of HLA-B*07:02, HLA-DRB1*15:01, HLA-DQA1*01:02, and HLA-DQB1*06:02 on the same chromosome), which was also associated (p = 2.70 × 10−8). However, conditional analyses indicated an independent effect of HLA-DQA1*01:02 from the extended haplotype (Additional file 22). The extended ancestral haplotype DR3-DQ2 (the combined presence of HLA-DRB1*03:01, HLA-DQA1*05:01, and HLA-DQB1*02:01 on the same chromosome) was the top protective association for bADA levels across preparations (p = 4.97 × 10−10; Table 5 and Additional file 22). In conditional analyses, none of these HLA associations was independent of the top bADA SNP rs9271377 (Additional file 22).

Treatment-specific HLA analyses

For many identified HLA alleles, support for an association was predominantly observed in patients treated with either IFNβ-1a s.c. or -1b s.c. but not in both groups simultaneously (Additional file 14).

When analyzing IFNβ-1a s.c.-treated patients, allele HLA-DQB1*06:02 and the ancestral haplotype DR15-DQ6, both smaller subsets of B7-DQ6, were the HLA risk variants showing the most robust support for an association in each of the three ADA measurements (DR15-DQ6 nADA presence: OR = 2.88 (2.29–3.61), p = 7.4 × 10−20; Fig. 2a, Table 5 and Additional file 22). All other risk alleles associated with IFNβ-1a s.c.-induced ADA were part of this extended haplotype and did not remain significant when conditioning for DR15-DQ6. None of these variants passed the discovery-stage significance threshold in patients receiving IFNβ-1b s.c. (Additional file 19). When conditioning the IFNβ-1a s.c.-associated risk HLA haplotype and alleles for the top GWAS SNP rs77278603, none of the HLA signals remained significant (Additional file 22). The DR15-DQ6 risk association from the HLA analysis thus likely corresponds to the GWAS association of rs77278603 and correlated SNPs.

Fig. 2
figure2

Treatment-specific HLA haplotypes. a, b The association of nADA titers for selected extended haplotypes showing strong treatment-specific effects. For association statistics, see Table 5 and Additional file 22. Disc. = discovery, Repl. = replication. a The association of the DR15-DQ6 haplotype with nADA titers is specific for IFNβ-1a s.c. b The association of the DR4-DQ3 haplotype with nADA titers is specific for IFNβ-1b s.c

The ancestral haplotype DR3-DQ2 and its allele HLA-DQB1*02:01 were the protective alleles showing the lowest p values in IFNβ-1a s.c.-treated patients (DR3-DQ2 nADA presence: OR = 0.37 (0.27–0.52), p = 3.7 × 10−9; Table 5 and Additional file 22), with all other protective alleles being part of this extended haplotype. No other allele remained significant when conditioning for DR3-DQ2. None of these variants were significantly associated in patients treated with IFNβ-1b s.c. (e.g., DR3-DQ2, nADA presence p = 0.27) or IFNβ-1a i.m. after correction for multiple testing (Additional files 19 and 22).

In IFNβ-1b s.c.-treated patients, HLA-DRB1*04:01 was the risk allele that showed the most robust support for an association for all three ADA measurements, and all associated alleles were part of the haplotype DR4-DQ3 (the combined presence of HLA-DRB1*04:01, HLA-DQA1*03:01, and HLA-DQB1*03:02 on the same chromosome). The pooled association strength of DR4-DQ3 for nADA presence was OR = 7.35 (4.32–12.47), p = 1.5 × 10−13 in patients treated with IFNβ-1b s.c. (Fig. 2b, Table 5, and Additional file 22). Of note, when conditioning for DR4-DQ3, the association of HLA-DRB1*04:01 remained significant, suggesting it to constitute the primary signal (Additional file 22). However, statistical power for fine-mapping this signal was limited because of the low allele frequencies of the allele and the haplotype (Table 3). These alleles and the haplotype were not significantly associated in IFNβ-1a s.c.-treated patients (e.g., DR4-DQ3, nADA presence for IFNβ-1a s.c., p = 0.22 (Additional files 21 and 22)). In conditional analyses of the IFNβ-1b s.c.-associated risk HLA haplotype and alleles for the top GWAS SNP rs28366299, HLA-DRB1*04:01 and DR4-DQ3 remained significant (Additional files 15 and 22) and vice versa (Additional file 8), indicating that rs28366299 represents an independent signal from the HLA association. There were no significant protective alleles for patients receiving IFNβ-1b s.c.

Analysis of candidate variants

As a secondary analysis, we analyzed the association of SNPs and HLA alleles previously published to be associated with ADA (Table 6 and Additional file 23). Here, we used the Bonferroni correction for 20 SNPs and alleles, corresponding to a significance threshold of α = 2.5 × 10−3. We found support for an association of 15 of the 20 candidate variants (Table 6), but not for SNP rs4961252 (which does not map to the MHC region), the MHC class I alleles, or the HLA-DRB1*11 alleles. Notably, HLA-DRB1*04:08 was only associated with bADA levels and nADA titers but not with nADA presence, HLA-DRB1*16:01 only with bADA levels. In follow-up analyses, we observed that these two HLA alleles were also associated with nADA titers in nADA-positive patients (Table 6 and Additional file 23).

Table 6 Association of candidate variants

Prediction of ADA

We predicted the occurrence of nADA in the replication dataset using the genetic models derived in the discovery dataset. For each treatment preparation, we analyzed eight polygenic risk scores (PRS), the top single GWAS variant, and the top HLA allele from the discovery stage (Fig. 3a, b and Additional files 24 and 25). Based on the AUC, the best predictions were achieved in models either featuring only the top variant or by PRS consisting of variants showing strong support for an association (Fig. 3). We thus did not observe evidence for a highly polygenic inheritance of nADA development. Interestingly, only treatment-specific predictions were significant; IFNβ-1a s.c.-specific models could not predict nADA in IFNβ-1b s.c.-treated patients and vice versa (Fig. 3c and Additional file 26). Prediction models containing either the top PRS or SNP showed distinctly increased AUCs, Nagelkerke’s pseudo-R2, sensitivities, and specificities over models containing only the covariates (Table 7).

Fig. 3
figure3

Prediction of nADA. Treatment-specific prediction of the presence of nADA in the replication data. Eight PRS, the top single GWAS variant, and the top HLA allele from the discovery stage were analyzed, with sex, age, treatment preparation and duration, titration site, and ancestry components as covariates. The plots show the area under the receiver operating characteristic curve (AUC) and its 95% confidence interval (CI) calculated using bootstrapping. Bonferroni = significant after Bonferroni correction for multiple testing; nominal = nominally significant (p < 0.05); n.s. = not significant. a, b AUC of all ten prediction models for a IFNβ-1a s.c. and b IFNβ-1b s.c.; Bonferroni correction for ten tests. c The model with the highest AUC for each treatment preparation, Bonferroni correction for 160 tests (α = 3.13 × 10−4). Boxes show the prediction groups (replication data) and columns within each box the training data groups (discovery data)

Table 7 Treatment-specific prediction of the presence of nADA in the replication data

Finally, patients with a high and low genetic risk burden were contrasted [46]. To this end, patients within the upper 30% of the top-associated PRS were compared to the patients within the lower 30%. This specific threshold was set to allow for a large enough sample size in the replication dataset for the stable convergence of regression models when conducting cross-validation. In addition, the respective top SNP was analyzed using dominant coding, thereby comparing no copy of the risk allele to any copy. Here, the best prediction was achieved for patients treated with IFNβ-1a s.c. (Table 7). In the Swedish cohort, it had an AUC = 0.91 (CI = 0.85–0.95), pseudo-R2 = 0.59, sensitivity = 0.78, and specificity = 0.90. Patients with the top 30% of genetic risk had, compared to patients in the bottom 30%, an OR = 73.9 (CI = 11.8–463.6, p = 4.4 × 10−06) of developing nADA. In the German cohort, the same model had an AUC = 0.83 (0.71–0.92), pseudo-R2 = 0.38, sensitivity = 0.80, and specificity = 0.76; the OR of patients with the top 30% of genetic risk was 13.8 (3.0–63.3, p = 7.5 × 10−4).

Discussion

Several studies have previously assessed genetic risk factors for ADA. They were limited by much smaller sample sizes, only analyzed a single population, focused mostly on HLA alleles, or did not consistently assess ADA with sensitive and validated methods. The latter is also reflected in the increased number of nADA-positive samples in the measurements conducted for the present study, compared to previous results (Additional files 1 and 3). Most importantly, the existing studies neither reached a consensus on genetic risk factors nor could they delineate a robust prediction model for ADA. To our knowledge, the present study constitutes the most extensive genetic characterization of ADA risk to date, is the first systematic comparison of the genetics of different ADA types, and includes the first genetic prediction model for ADA against IFNβ.

Previously reported genetic risk factors

All genetic variants robustly associated with ADA in the present study map to the MHC region and are linked to the expression of HLA genes or amino acid changes in the peptide-binding groove of HLA molecules. SNP rs9272105, mapping to the MHC region and previously identified in a study conducted on a subset of the patients analyzed here [14], was significantly associated across treatment preparations in the present study. However, we found no support for an association of variant rs4961252 on chromosome 8, identified in the same study [14], which confirms a previously failed replication attempt [15]. Both variants already identified in an independent study of IFNβ-1b s.c.-treated patients [15] replicated only in individuals treated with IFNβ-1b s.c.

Previous studies prioritized sixteen different HLA alleles as potentially associated with nADA presence, nADA titers, or bADA levels [10,11,12,13, 15, 16]. Of these, eleven were significantly associated with an ADA measurement for any treatment preparation in our study (Table 6), and five were not (Additional file 23). Importantly, the present study does not constitute a formal replication for many of the candidate HLA alleles because of the extensive sample overlap with previous Swedish and German studies. The HLA alleles that did not replicate had low frequencies, with a maximum AF = 0.06, and showed only weak support in the original studies. Notably, in some previous studies, highly correlated alleles were analyzed as if they were independent variants, and some studies failed to correct them appropriately for multiple testing. Both factors may have led to an overestimation of the number of associated HLA alleles in previous studies.

One study analyzed Spanish patients [16], a population whose allele frequencies and linkage patterns differ from the individuals studied here, and whose results may thus not be fully comparable to the present study. Three MHC class I HLA alleles reported to be associated with ADA by Núñez et al. [16] did not replicate in the present study (Additional file 23). All three alleles are more frequent in Spain than in Germany and Sweden, with, e.g., HLA-B*14:02 showing a frequency of 1% in Sweden [47], 2% in Germany, and 4% in Spain [48]. To reliably assess whether the associations of these alleles are specific to Spanish populations or whether the lack of correction for multiple testing led to type I errors in the original study, independent replication studies on Spanish patients are required. Next to population-specific effects, joint analyses of patients receiving different proportions of IFNβ treatment preparations constituted a source of heterogeneity in previous studies. In our comprehensive analyses, we could now consolidate several HLA alleles published in previous studies into few treatment-specific haplotypes.

Treatment preparation-specific risk

Before our study, it was unclear whether treatment preparation-independent or preparation-specific genetic risk factors dominate ADA risk. The associations of the top GWAS SNPs identified in the analyses across all treatment preparations were mostly supported in all treatment preparations (Additional file 14). Nevertheless, we observed lower p values and larger effect sizes in the preparation-specific analyses than in the models combining patients across treatments. The combination of patients receiving different treatment preparations thus created heterogeneity that decreased statistical power. This hypothesis was further supported in stepwise conditional analyses. Here, we observed more evidence for the existence of independent risk loci in analyses across preparations than was the case in treatment-specific analyses. Likely, such presumably independent loci in the combined analysis reflect treatment preparation-specific effects. These findings thus argue in favor of conducting treatment-specific rather than cross-treatment analyses. In future studies of ADA against biopharmaceuticals, analyses of preparation-specific risk factors should, therefore, be prioritized.

Although differences in the antigenic potential of the various IFNβ preparations are known [1], the extent of preparation-specific genetic risk observed in the present study is striking (Fig. 2). There are several plausible explanations for why the preparations might be processed differently by the immune system. While the amino acid sequence of IFNβ-1a is identical to natural human IFNβ, IFNβ-1b diverges at two positions: IFNβ-1b lacks the N-terminal methionine, and a cysteine at position 17 is substituted by serine. Furthermore, the products are raised in different cell types, prokaryotic E. coli and eukaryotic Chinese hamster ovary cells, leading to different post-translational modifications, especially glycosylation [1].

Lack of glycosylation facilitates the formation of protein aggregates, increasing the immunogenicity of IFNβ-1b [1, 2]. Previous research demonstrated that, among the three preparations analyzed in the present study, IFNβ-1b shows the highest tendency to aggregate [49]. IFNβ-1a i.m., which does not contain human serum albumin, forms the fewest aggregates and shows the lowest rate of ADA. Furthermore, aggregates observed with IFNβ-1a s.c. preparations are mainly formed by human serum albumin [49]. Differences in IFNβ protein aggregation might, in addition to increased presentation of peptides by dendritic cells and, thus, increased T cell activation [50, 51], also contribute to the diversification of genetic risk factors. When taken up by antigen-presenting cells, e.g., dendritic cells, IFNβ oligomers are likely degraded differently from monomers. Such differences in processing could produce diverse peptides, which may be presented by different MHC class II molecules [50].

Post-translational modifications not only affect aggregate formation but, together with differences in the amino acid sequences, also alter the biochemical properties of IFNβ-1a and IFNβ-1b. Thereby, both post-translational modifications and differences in the amino acid sequence may contribute to the preparation-specific associations with HLA alleles [4]. For example, altering epitopes by glycosylation strongly affects antigen recognition [52]. Possibly, glycosylated IFNβ-1a peptides are thus preferentially recognized by different peptide-binding grooves of MHC molecules than IFNβ-1b-derived epitopes are. Similarly, also the amino acid changes may alter the binding of IFNβ peptides to MHC molecules and T cell recognition [53].

Additional factors in the processing of treatment preparations can influence how the immune system recognizes them. Spontaneously occurring modifications like deamidation, oxidation, and glycation alter the surface and chemical properties of proteins. These modifications even diverge between preparations sharing the identical amino acid sequence, e.g., IFNβ-1a s.c. and i.m., by differential production, processing, or storing of the biopharmaceuticals [3]. Other chemical alterations of amino acids like phosphorylation, PEGylation, methylation, or acetylation can be applied during the manufacturing of drugs, e.g., to alter their stability, and also change epitopes, leading to differential binding to allelic variants of HLA heterodimers [54, 55]. Importantly, these modifications also happen after administration of the product in vivo, and glycosylation may well affect the likelihood of them taking place.

In summary, diverging post-translational modifications may contribute to the observed differences in preparation-specific genetic risk factors. Notably, the MHC class II peptide-binding groove is formed by heterodimers of two HLA proteins, likely contributing to the association of haplotypes spanning HLA α and β chain genes, like HLA-DQA1 and HLA-DQB1, with IFNβ ADA. However, it is unlikely that preparation-specific risk can entirely be attributed to genetic factors. For example, the dosage and injection frequency of preparations may affect the likelihood of developing ADA, independently of genetic risk [8, 56, 57]. Nevertheless, most patients develop ADA within the first months of IFNβ treatment [58], arguing against pronounced long-term dosage-specific effects and underlining the importance of genetic risk.

Next to having to rely on imputed HLA alleles, the low number of available patients that developed ADA under treatment with IFNβ-1a i.m., rendering IFNβ-1a i.m.-specific analyses unfeasible, constitutes a limitation of the present study. We expect genetic risk factors for IFNβ-1a i.m.-induced ADA to exist, but whether these are independent of IFNβ-1a s.c.-associated risk remains to be shown.

The complexity of the genetic risk landscape

Using conditional analyses, we did not find evidence for more than one genetic risk locus for IFNβ-1a s.c.-induced ADA. Results from previous studies can thus, at least for Swedish and German patients, be consolidated to the extended haplotype DR15-DQ6. In the present dataset, it is impossible to assess whether the combined DR15-DQ6 haplotype constitutes the real risk factor for IFNβ-1a-s.c. or whether any of the single alleles HLA-DQB1*06:02 or HLA-DQA1*01:02 convey this risk, with the haplotype showing an association merely because of LD. DR15-DQ6 (MAFKI = 0.34, MAFTUM = 0.29) is less common than the two single alleles, especially compared to HLA-DQA1*01:02 (MAFKI = 0.42, MAFTUM = 0.36). Because statistical power is dependent on the AF, the slightly lower statistical support for the association of DR15-DQ6, compared to the single alleles, likely reflects these differences in AF and power. We thus hypothesize that the combined haplotype DR15-DQ6 constitutes the primary signal. Nevertheless, such fine-mapping and the differentiation between the correlated alleles is irrelevant for risk predictions. Because of the strong correlation of alleles observed within the extended haplotype, any of these alleles can reliably be used as a proxy for the others in prediction models.

Similarly, conditional analyses support the association of the extended haplotype DR3-DQ2 as the primary protective genetic signal for IFNβ-1a s.c., without evidence for secondary signals. However, in the present sample, the association of this haplotype cannot be separated from HLA-DQB1*02:01. By contrast, genetic risk for IFNβ-1b s.c.-induced ADA appears to be more complicated. The association of the haplotype DR4-DQ3 could not fully explain the signal of its allele HLA-DRB1*04:01. Moreover, we found evidence for a secondary signal in stepwise conditional regression analyses. Notably, the prediction models for IFNβ-1b s.c. did not perform as well as the prediction for IFNβ-1a s.c.-induced ADA did, which possibly reflects this more complex risk landscape. To truly unravel an additional potential polygenic contribution to ADA risk, the current study still lacked the sample size necessary for reliably detecting polygenic variants with their expected small effect sizes [59, 60].

Population-specific risk differences

Alleles from all three associated haplotypes, HLA-DRB1*15:01, HLA-DRB1*04:01, and the protective HLA-DRB1*03:01, are concurrent risk factors for MS [42]. The unfortunate coincidence that HLA-DRB1*15:01 and HLA-DRB1*04:01 are enriched among MS patients and also constitute ADA risk factors likely contributes to the high incidence of IFNβ ADA among MS patients.

Interestingly, these alleles also show substantial population-specific differences [47, 48, 61]: The IFNβ-1a s.c. risk allele HLA-DRB1*15:01 is less frequent in Southern Europe and for Ashkenazi, Southern Hispanic, and African ancestries (e.g., Italy 5.6–6.4%, Southwestern Spain 5.2–8.6%). At the same time, it is especially frequent in individuals with ancestry from other parts of Europe (e.g., Northern Spain 16.7–32.1%, Germany 12.9–17.2%, Sweden 16.1%). Note that being the most important MS risk variant, allele HLA-DRB1*15:01 is more frequent among MS patients than in the respective general population. However, population-specific frequency differences exist on top: In the present study, the frequency of HLA-DRB1*15:01 was markedly higher in Swedish (36.1%) than in German MS patients (30.9%). On average, Swedish patients may thus, in comparison to German patients, be at higher risk of developing nADA against IFNβ-1a s.c. The IFNβ-1b s.c. risk allele HLA-DRB1*04:01 is more frequent in parts of Northwestern, Northern, and Central Europe (e.g., England 12.4–13.5%, Denmark 17.6%, Sweden 13.7%) than in Southern Europe and most other ancestries (e.g., Italy 1.7–4.1%, Spain 2.0–3.8%). Germany lies in between with frequencies of 6.8–9.4%.

While the frequencies of both risk alleles for IFNβ-1a s.c.- and IFNβ-1b s.c.-induced ADA thus roughly decrease along a North-South gradient within Europe, their relative frequencies differ sharply in some ancestries (Additional file 27). For example, in Northern Spain, the major genetic risk factor for IFNβ-1a s.c.-specific ADA occurs > 8 times more often than the one for IFNβ-1b s.c.-induced ADA. Such substantial, population-specific differences in risk allele frequencies likely exist for ADA against any biopharmaceutical. If genetic risk factors for a biopharmaceutical are known, and personalized genotyping data for patients are not available, recommendations for the choice of a specific treatment preparation could thus be made on a population level. Where the availability of genetic testing is limited, patients from populations with higher frequencies of risk alleles could be prioritized for genetic testing, as already practiced for other treatments [62].

Genetic factors contributing to nADA titers

The same or highly correlated risk factors contributed to the presence of nADA and the magnitude of nADA titers and bADA levels. The heritability of nADA, explained by common variants, after correction for confounders like treatment preparation and duration, sex, and age, was very high—h2gl = 0.79 on a liability scale. This result underlines the importance of genetic factors in the occurrence of IFNβ ADA. Although nADA titers need to cross a threshold to become functionally relevant, the associations of genetic risk factors with both nADA presence and titers may indicate that most genetic risk factors mainly influence the likelihood of developing ADA and less the absolute titers. Interestingly, the association of the candidate variants HLA-DRB1*04:08 and HLA-DRB1*16:01 was only significant in analyses of quantitative nADA titers or bADA levels but not for nADA presence. This finding indicates that genetic factors influencing the amount of ADA likely exist. In fact, follow-up analyses found that these two alleles were also associated with nADA titers in nADA-positive patients treated with IFNβ-1b s.c. In the present study, we did not conduct hypothesis-free GWAS of nADA titers in the smaller subsample restricted to nADA-positive patients. To reliably distinguish between influences of genetic variants on either the likelihood of nADA development or the amount of nADA, larger patient samples than analyzed in the present study should be collected for future studies.

Comparison to other MS- and ADA-related analyses

The allele HLA-DRB1*15:01 is associated with the risk of MS [42], earlier age at MS disease onset [24], and developing nADA against IFNβ-1a s.c. Moreover, the same allele and associated haplotypes are also associated with intrathecal immunoglobulin G levels [63] and Epstein Barr viral loads and titers in MS patients [64, 65]. By contrast, a strong negative association between HLA-DRB1*15:01 and JC polyomavirus antibody status was reported [66]. The HLA allele HLA-DRB1*15:01, therefore, constitutes the key genetic risk factor for MS, which also differentially influences gene-by-environment interactions, disease severity, and treatment complications.

We identified HLA-DQA1*05:01 to protect from nADA against IFNβ-1a s.c. Interestingly, the same allele is strongly associated with the risk of ADA against the widely used anti-tumor necrosis factor (TNF) treatments for Crohn’s disease [67]. This association was consistent across the two anti-TNF biopharmaceutical drugs adalimumab and infliximab. The HLA allele HLA-DRB1*03, also protective against IFNβ-1a s.c.-induced nADA, is, together with HLA-DQA1*05:01, part of the haplotype DR3-DQ2. HLA-DRB1*03 was published as a risk factor against adalimumab and infliximab in patients suffering from either inflammatory bowel disease or rheumatoid arthritis [68, 69]. Whether treatment-specific genetic risk factors also exist for anti-TNF biopharmaceuticals and whether the haplotype DR3-DQ2 or one of the single HLA alleles confers the risk for anti-TNF ADA could be an interesting topic of future studies.

Prediction of ADA

The prediction models performed better for IFNβ-1a s.c.-induced than for IFNβ-1b s.c.-induced nADA, and they could predict nADA better in the Swedish KI than in the German TUM cohort. Our results indicate that, compared to IFNβ-1b s.c., genetic risk for IFNβ-1a s.c.-induced ADA is more dominated by a single locus. Overall, more patients receiving IFNβ-1a s.c. than IFNβ-1b s.c. were analyzed (1145 vs. 1010). Both factors likely contributed to better prediction models and performance in IFNβ-1a s.c.-treated patients. The top-associated ADA risk SNP was more frequent in Swedish patients than in German ones (43% vs. 40%), and the Swedish sample contained more patients treated with IFNβ-1a s.c. (Sweden 590, Germany 558), of whom more were nADA-positive (34.6% vs. 33.7%). Although these individual differences were small, they may have contributed to prediction models performing better in the Swedish dataset.

Contrasting the samples in the top and bottom percentiles of polygenic risk score distributions is a common practice to compare individuals carrying a high genetic risk burden to the ones not at risk [46]. The sensitivity and specificity reached in the comparison of individuals in the top 30% nADA risk group compared to the bottom 30% (0.78 and 0.90, respectively, in Swedish IFNβ-1a s.c. patients and 0.80 and 0.76, respectively, in German patients) may still not be sufficient for a routine clinical test. However, these prediction models could certainly be optimized by the inclusion of additional predictive factors, e.g., body mass index [70], not available in the present retrospective setting. The significant predictive improvement of the genetic risk model compared to a model containing only demographic and clinical variables (Table 7) underlines the importance of incorporating genetics in prediction models for ADA. The high odds of patients at genetic risk for nADA (Sweden: OR = 73.9, Germany: OR = 13.8) support the use of genetic stratification as a personalized medicine tool—patients at high genetic risk should either switch to a different drug or be monitored more closely, as suggested for other conditions [71].

Conclusions

We have conducted a comprehensive characterization of genetic risk for IFNβ-induced ADA, consolidating previous research. Next to treatment-specific risk factors, we described ancestry-specific effects relevant for treatment choice in specific populations. Our robust prediction models could be employed for personalized medicine, guiding treatment recommendations, and efficient nADA testing regimes. Importantly, our study can serve as a blueprint for the analysis of genetic factors influencing ADA against other biopharmaceuticals and in the context of further diseases.

Availability of data and materials

The datasets generated and analyzed during the current study are available from the corresponding authors on reasonable request.

Abbreviations

ADA:

Anti-drug antibodies

AF:

Allele frequency

AUC:

Area under the receiver operating characteristic curve

bp:

Base pairs

bADA:

Binding ADA

CI:

95% confidence interval

eQTL:

Expression quantitative trait locus

GWAS:

Genome-wide association study

IFNβ:

Interferon β

i.m. :

Intramuscular

HLA:

Human leukocyte antigen

kbp:

Kilobase pairs

KI:

Karolinska Institutet Stockholm, Sweden

LD:

Linkage disequilibrium

MAD:

Median absolute deviation

MHC:

Major histocompatibility complex

MAF:

Minor allele frequency

MDS:

Multidimensional scaling

MS:

Multiple sclerosis

nADA:

Neutralizing ADA

OR:

Odds ratio

s.c. :

Subcutaneous

PRS:

Polygenic risk scores

SD:

Standard deviation

SE:

Standard error

SNP:

Single nucleotide polymorphism

TUM:

Technical University of Munich, Germany

References

  1. 1.

    Bertolotto A, Deisenhammer F, Gallo P, Sørensen P. Immunogenicity of interferon beta: differences among products. J Neurol. 2004;251(Suppl 2):ii15–24.

    PubMed  Google Scholar 

  2. 2.

    Sethu S, Govindappa K, Alhaidari M, Pirmohamed M, Park K, Sathish J. Immunogenicity to biologics: mechanisms, prediction and reduction. Arch Immunol Ther Exp. 2012;60(5):331–44.

    CAS  Article  Google Scholar 

  3. 3.

    Jawa V, Cousens LP, Awwad M, Wakshull E, Kropshofer H, De Groot AS. T-cell dependent immunogenicity of protein therapeutics: preclinical assessment and mitigation. Clin Immunol. 2013;149(3):534–55.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  4. 4.

    Engelhard VH, Altrich-Vanlith M, Ostankovitch M, Zarling AL. Post-translational modifications of naturally processed MHC-binding epitopes. Curr Opin Immunol. 2006;18(1):92–7.

    CAS  PubMed  Article  Google Scholar 

  5. 5.

    Deisenhammer F. Interferon-beta: neutralizing antibodies, binding antibodies, pharmacokinetics and pharmacodynamics, and clinical outcomes. J Interf Cytokine Res. 2014;34(12):938–45.

    CAS  Article  Google Scholar 

  6. 6.

    Bertolotto A, Gilli F, Sala A, Capobianco M, Malucchi S, Milano E, Melis F, Marnetto F, Lindberg RL, Bottero R, et al. Persistent neutralizing antibodies abolish the interferon beta bioavailability in MS patients. Neurology. 2003;60(4):634–9.

    CAS  PubMed  Article  Google Scholar 

  7. 7.

    Sorensen PS, Ross C, Clemmesen KM, Bendtzen K, Frederiksen JL, Jensen K, Kristensen O, Petersen T, Rasmussen S, Ravnborg M, et al. Clinical importance of neutralising antibodies against interferon beta in patients with relapsing-remitting multiple sclerosis. Lancet. 2003;362(9391):1184–91.

    CAS  PubMed  Article  Google Scholar 

  8. 8.

    Kappos L, Clanet M, Sandberg-Wollheim M, Radue EW, Hartung HP, Hohlfeld R, Xu J, Bennett D, Sandrock A, Goelz S, et al. Neutralizing antibodies and efficacy of interferon beta-1a: a 4-year controlled study. Neurology. 2005;65(1):40–7.

    CAS  PubMed  Article  Google Scholar 

  9. 9.

    Pachner AR, Cadavid D, Wolansky L, Skurnick J. Effect of anti-IFNbeta; antibodies on MRI lesions of MS patients in the BECOME study. Neurology. 2009;73(18):1485–92.

    CAS  PubMed  Article  Google Scholar 

  10. 10.

    Barbosa M, Vielmetter J, Chu S, Smith DD, Jacinto J. Clinical link between MHC class II haplotype and interferon-beta (IFN-β) immunogenicity. Clin Immunol. 2006;118(1):42–50.

    CAS  PubMed  Article  Google Scholar 

  11. 11.

    Buck D, Cepok S, Hoffmann S, Grummel V, Jochim A, Berthele A, Hartung HP, Wassmuth R, Hemmer B. Influence of the HLA-DRB1 genotype on antibody development to interferon beta in multiple sclerosis. Arch Neurol. 2011;68(4):480–7.

    PubMed  Article  Google Scholar 

  12. 12.

    Hoffmann S, Cepok S, Grummel V, Lehmann-Horn K, Hackermuller J, Stadler PF, Hartung HP, Berthele A, Deisenhammer F, Wassmuth R, et al. HLA-DRB1*0401 and HLA-DRB1*0408 are strongly associated with the development of antibodies against interferon-beta therapy in multiple sclerosis. Am J Hum Genet. 2008;83(2):219–27.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  13. 13.

    Link J, Lundkvist Ryner M, Fink K, Hermanrud C, Lima I, Brynedal B, Kockum I, Hillert J, Fogdell-Hahn A. Human leukocyte antigen genes and interferon beta preparations influence risk of developing neutralizing anti-drug antibodies in multiple sclerosis. PLoS One. 2014;9(3):e90479.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  14. 14.

    Weber F, Cepok S, Wolf C, Berthele A, Uhr M, Bettecken T, Buck D, Hartung HP, Holsboer F, Muller-Myhsok B, et al. Single-nucleotide polymorphisms in HLA- and non-HLA genes associated with the development of antibodies to interferon-beta therapy in multiple sclerosis patients. Pharmacogenomics J. 2012;12(3):238–45.

    CAS  PubMed  Article  Google Scholar 

  15. 15.

    Buck D, Andlauer TF, Igl W, Wicklein EM, Muhlau M, Weber F, Kochert K, Pohl C, Arnason B, Comi G, et al. Effect of HLA-DRB1 alleles and genetic variants on the development of neutralizing antibodies to interferon beta in the BEYOND and BENEFIT trials. Mult Scler. 2019;25(4):565–73.

    CAS  PubMed  Article  Google Scholar 

  16. 16.

    Núñez C, Cénit MC, Alvarez-Lafuente R, Río J, Fernández-Arquero M, Arroyo R, Montalbán X, Fernández O, Oliver-Martos B, Leyva L, et al. HLA alleles as biomarkers of high-titre neutralising antibodies to interferon-β therapy in multiple sclerosis. J Med Genet. 2014;51(6):395–400.

    PubMed  Article  CAS  Google Scholar 

  17. 17.

    Stickler M, Valdes AM, Gebel W, Razo OJ, Faravashi N, Chin R, Rochanayon N, Harding FA. The HLA-DR2 haplotype is associated with an increased proliferative response to the immunodominant CD4+ T-cell epitope in human interferon-β. Genes Immun. 2004;5(1):1–7.

    CAS  PubMed  Article  Google Scholar 

  18. 18.

    Kraft P. Curses—winner’s and otherwise—in genetic epidemiology. Epidemiology. 2008;19(5):649–51.

    PubMed  Article  Google Scholar 

  19. 19.

    Hansen BB, Klopfer SO. Optimal full matching and related designs via network flows. J Comput Graph Stat. 2006;15(3):609–27.

    Article  Google Scholar 

  20. 20.

    Pachner AR. An improved ELISA for screening for neutralizing anti-IFN-beta antibodies in MS patients. Neurology. 2003;61(10):1444–6.

    CAS  PubMed  Article  Google Scholar 

  21. 21.

    Hermanrud C, Ryner M, Luft T, Jensen PE, Ingenhoven K, Rat D, Deisenhammer F, Sorensen PS, Pallardy M, Sikkema D, et al. Development and validation of cell-based luciferase reporter gene assays for measuring neutralizing anti-drug antibodies against interferon beta. J Immunol Methods. 2016;430:1–9.

    CAS  PubMed  Article  Google Scholar 

  22. 22.

    Jensen PEH, Warnke C, Ingenhoven K, Piccoli L, Gasis M, Hermanrud C, Fernandez-Rodriguez BM, Ryner M, Kramer D, Link J, et al. Detection and kinetics of persistent neutralizing anti-interferon-beta antibodies in patients with multiple sclerosis. Results from the ABIRISK prospective cohort study. J Neuroimmunol. 2019;326:19–27.

    CAS  PubMed  Article  Google Scholar 

  23. 23.

    Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015;4(1):7.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  24. 24.

    Andlauer TF, Buck D, Antony G, Bayas A, Bechmann L, Berthele A, Chan A, Gasperi C, Gold R, Graetz C, et al. Novel multiple sclerosis susceptibility loci implicated in epigenetic regulation. Sci Adv. 2016;2(6):e1501678.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  25. 25.

    Delaneau O, Zagury JF, Marchini J. Improved whole-chromosome phasing for disease and population genetic studies. Nat Methods. 2013;10(1):5–6.

    CAS  PubMed  Article  Google Scholar 

  26. 26.

    Howie B, Fuchsberger C, Stephens M, Marchini J, Abecasis GR. Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. Nat Genet. 2012;44(8):955.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  27. 27.

    Howie BN, Donnelly P, Marchini J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet. 2009;5(6):e1000529.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  28. 28.

    International HIV Controllers Study, Pereyra F, Jia X, McLaren PJ, Telenti A, Bakker PIWd, Walker BD, Ripke S, Brumme CJ, Pulit SL , et al. The major genetic determinants of HIV-1 control affect HLA class I peptide presentation. Science 2010, 330(6010):1551–1557.

  29. 29.

    Jia X, Han B, Onengut-Gumuscu S, Chen WM, Concannon PJ, Rich SS, Raychaudhuri S, de Bakker PI. Imputing amino acid polymorphisms in human leukocyte antigens. PLoS One. 2013;8(6):e64683.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Browning BL, Browning SR. A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. Am J Hum Genet. 2009;84(2):210–23.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  31. 31.

    Giovannoni G, Barbarash O, Casset-Semanaz F, Jaber A, King J, Metz L, Pardo G, Simsarian J, Sorensen PS, Stubinski B, et al. Immunogenicity and tolerability of an investigational formulation of interferon-beta1a: 24- and 48-week interim analyses of a 2-year, single-arm, historically controlled, phase IIIb study in adults with multiple sclerosis. Clin Ther. 2007;29(6):1128–45.

    CAS  PubMed  Article  Google Scholar 

  32. 32.

    Lee SH, Yang J, Goddard ME, Visscher PM, Wray NR. Estimation of pleiotropy between complex diseases using single-nucleotide polymorphism-derived genomic relationships and restricted maximum likelihood. Bioinformatics. 2012;28(19):2540–2.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  33. 33.

    Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88(1):76–82.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. 34.

    Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, Nyholt DR, Madden PA, Heath AC, Martin NG, Montgomery GW, et al. Common SNPs explain a large proportion of the heritability for human height. Nat Genet. 2010;42(7):565–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  35. 35.

    Lee SH, Wray NR, Goddard ME, Visscher PM. Estimating missing heritability for disease from genome-wide association studies. Am J Hum Genet. 2011;88(3):294–305.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  36. 36.

    Willer CJ, Li Y, Abecasis GR. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics. 2010;26(17):2190–1.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  37. 37.

    Stouffer SA, Suchman EA, Devinney LC, Star SA, Jr RMW: The American soldier: adjustment during army life. Princeton: Princeton University Press; 1949.

  38. 38.

    Mele M, Ferreira PG, Reverter F, DeLuca DS, Monlong J, Sammeth M, Young TR, Goldmann JM, Pervouchine DD, Sullivan TJ, et al. Human genomics. The human transcriptome across tissues and individuals. Science. 2015;348(6235):660–5.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. 39.

    de Leeuw CA, Mooij JM, Heskes T, Posthuma D. MAGMA: generalized gene-set analysis of GWAS data. PLoS Comput Biol. 2015;11(4):e1004219.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  40. 40.

    Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  41. 41.

    Horton R, Gibson R, Coggill P, Miretti M, Allcock RJ, Almeida J, Forbes S, Gilbert JG, Halls K, Harrow JL, et al. Variation analysis and gene annotation of eight MHC haplotypes: the MHC Haplotype Project. Immunogenetics. 2008;60(1):1–18.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  42. 42.

    Patsopoulos NA, Barcellos LF, Hintzen RQ, Schaefer C, van Duijn CM, Noble JA, Raj T, IMSGC, ANZgene, Gourraud PA, et al. Fine-mapping the genetic association of the major histocompatibility complex in multiple sclerosis: HLA and non-HLA effects. PLoS Genet. 2013;9(11):e1003926.

  43. 43.

    Moutsianas L, Jostins L, Beecham AH, Dilthey AT, Xifara DK, Ban M, Shah TS, Patsopoulos NA, Alfredsson L, Anderson CA, et al. Class II HLA interactions modulate genetic risk for multiple sclerosis. Nat Genet. 2015;47(10):1107–13.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. 44.

    Andlauer TFM, Guzman-Parra J, Streit F, Strohmaier J, Gonzalez MJ, Gil Flores S, Cabaleiro Fabeiro FJ, Del Rio NF, Perez FP, Haro Gonzalez J, et al. Bipolar multiplex families have an increased burden of common risk variants for psychiatric disorders. Mol Psychiatry. 2019. https://doi.org/10.1038/s41380-019-0558-2.

  45. 45.

    Andlauer TFM, Mühleisen TW, Hoffstaedter F, Teumer A, Wittfeld K, Teuber A, Reinbold CS, Bülow R, Caspers S, Herms S, et al. Genetic factors influencing a neurobiological substrate for psychiatric disorders. bioRxiv. 2019:774463. https://doi.org/10.1101/774463.

  46. 46.

    Khera AV, Chaffin M, Aragam KG, Haas ME, Roselli C, Choi SH, Natarajan P, Lander ES, Lubitz SA, Ellinor PT, et al. Genome-wide polygenic scores for common diseases identify individuals with risk equivalent to monogenic mutations. Nat Genet. 2018;50(9):1219–24.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  47. 47.

    Nordin J, Ameur A, Lindblad-Toh K, Gyllensten U, Meadows JRS. SweHLA: the high confidence HLA typing bio-resource drawn from 1000 Swedish genomes. Eur J Hum Genet. 2020;28(5):627–35.

    CAS  PubMed  Article  Google Scholar 

  48. 48.

    Faviel, McCabe A, Eduardo, Jones J, Takeshita L, Nestor, Glenda, Ramsbottom K, Ghattaoraya G, Alfirevic A, et al. Allele frequency net database (AFND) 2020 update: gold-standard data classification, open access genotype data and new query tools. Nucleic Acids Res. 2020;48(D1):D783–8.

  49. 49.

    Barnard JG, Babcock K, Carpenter JF. Characterization and quantitation of aggregates and particles in interferon-beta products: potential links between product quality attributes and immunogenicity. J Pharm Sci. 2013;102(3):915–28.

    CAS  PubMed  Article  Google Scholar 

  50. 50.

    Rombach-Riegraf V, Karle AC, Wolf B, Sorde L, Koepke S, Gottlieb S, Krieg J, Djidja MC, Baban A, Spindeldreher S, et al. Aggregation of human recombinant monoclonal antibodies influences the capacity of dendritic cells to stimulate adaptive T-cell responses in vitro. PLoS One. 2014;9(1):e86322.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  51. 51.

    Gallais Y, Szely N, Legrand F-X, Leroy A, Pallardy M, Turbica I. Effect of growth hormone and IgG aggregates on dendritic cells activation and T-cells polarization. Immunol Cell Biol. 2017;95(3):306–15.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Wei X, Decker JM, Wang S, Hui H, Kappes JC, Wu X, Salazar-Gonzalez JF, Salazar MG, Kilby JM, Saag MS, et al. Antibody neutralization and escape by HIV-1. Nature. 2003;422(6929):307–12.

    CAS  PubMed  Article  Google Scholar 

  53. 53.

    Kalluri SR, Grummel V, Hracsko Z, Pongratz V, Pernpeintner V, Gasperi C, Buck D, Hemmer B, Consortium A. Interferon-beta specific T cells are associated with the development of neutralizing antibodies in interferon-beta treated multiple sclerosis patients. J Autoimmun. 2018;88:83–90.

    CAS  PubMed  Article  Google Scholar 

  54. 54.

    Schultz HS, Ostergaard S, Sidney J, Lamberth K, Sette A. The effect of acylation with fatty acids and other modifications on HLA class II:peptide binding and T cell stimulation for three model peptides. PLoS One. 2018;13(5):e0197407.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  55. 55.

    De Groot AS, Scott DW. Immunogenicity of protein therapeutics. Trends Immunol. 2007;28(11):482–90.

    PubMed  Article  CAS  Google Scholar 

  56. 56.

    Ross C, Clemmesen KM, Svenson M, Sorensen PS, Koch-Henriksen N, Skovgaard GL, Bendtzen K. Immunogenicity of interferon-beta in multiple sclerosis patients: influence of preparation, dosage, dose frequency, and route of administration. Danish Multiple Sclerosis Study Group. Ann Neurol. 2000;48(5):706–12.

    CAS  PubMed  Article  Google Scholar 

  57. 57.

    Hemmer B, Stuve O, Kieseier B, Schellekens H, Hartung HP. Immune response to immunotherapy: the role of neutralising antibodies to interferon beta in the treatment of multiple sclerosis. Lancet Neurol. 2005;4(7):403–12.

    CAS  PubMed  Article  Google Scholar 

  58. 58.

    Sorensen PS, Koch-Henriksen N, Ross C, Clemmesen KM, Bendtzen K, Danish Multiple Sclerosis Study G: Appearance and disappearance of neutralizing antibodies during interferon-beta therapy. Neurology 2005, 65(1):33–39.

  59. 59.

    Wray NR, Ripke S, Mattheisen M, Trzaskowski M, Byrne EM, Abdellaoui A, Adams MJ, Agerbo E, Air TM, Andlauer TMF, et al. Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression. Nat Genet. 2018;50(5):668–81.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  60. 60.

    International Multiple Sclerosis Genetics Consortium: Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility. Science 2019, 365(6460):eaav7188.

  61. 61.

    Ameur A, Dahlberg J, Olason P, Vezzi F, Karlsson R, Martin M, Viklund J, Kahari AK, Lundin P, Che H, et al. SweGen: a whole-genome data resource of genetic variability in a cross-section of the Swedish population. Eur J Hum Genet. 2017;25(11):1253–60.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  62. 62.

    Ferrell PB Jr, McLeod HL. Carbamazepine, HLA-B*1502 and risk of Stevens-Johnson syndrome and toxic epidermal necrolysis: US FDA recommendations. Pharmacogenomics. 2008;9(10):1543–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  63. 63.

    Gasperi C, Andlauer TFM, Keating A, Knier B, Klein A, Pernpeintner V, Lichtner P, Gold R, Zipp F, Then Bergh F, et al. Genetic determinants of the humoral immune response in MS. Neurol Neuroimmunol Neuroinflamm. 2020;7(5):e827.

  64. 64.

    Agostini S, Mancuso R, Guerini FR, D'Alfonso S, Agliardi C, Hernis A, Zanzottera M, Barizzone N, Leone MA, Caputo D, et al. HLA alleles modulate EBV viral load in multiple sclerosis. J Transl Med. 2018;16(1):80.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  65. 65.

    Jacobs BM, Giovannoni G, Cuzick J, Dobson R. Systematic review and meta-analysis of the association between Epstein-Barr virus, multiple sclerosis and other risk factors. Mult Scler. 2020. https://doi.org/10.1177/1352458520907901.

  66. 66.

    Sundqvist E, Buck D, Warnke C, Albrecht E, Gieger C, Khademi M, Lima Bomfim I, Fogdell-Hahn A, Link J, Alfredsson L, et al. JC polyomavirus infection is strongly controlled by human leucocyte antigen class II variants. PLoS Pathog. 2014;10(4):e1004084.

    PubMed  PubMed Central  Article  Google Scholar 

  67. 67.

    Sazonovs A, Kennedy NA, Moutsianas L, Heap GA, Rice DL, Reppell M, Bewshea CM, Chanchlani N, Walker GJ, Perry MH, et al. HLA-DQA1*05 carriage associated with development of anti-drug antibodies to infliximab and adalimumab in patients with Crohn’s disease. Gastroenterology. 2020;158(1):189–99.

    CAS  PubMed  Article  Google Scholar 

  68. 68.

    Billiet T, Vande Casteele N, Van Stappen T, Princen F, Singh S, Gils A, Ferrante M, Van Assche G, Cleynen I, Vermeire S. Immunogenicity to infliximab is associated with HLA-DRB1. Gut. 2015;64(8):1344–5.

    CAS  PubMed  Article  Google Scholar 

  69. 69.

    Liu M, Degner J, Davis JW, Idler KB, Nader A, Mostafa NM, Waring JF. Identification of HLA-DRB1 association to adalimumab immunogenicity. PLoS One. 2018;13(4):e0195325.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  70. 70.

    Callahan ST, Wolff M, Hill HR, Edwards KM, Keitel W, Atmar R, Patel S, Sahly HE, Munoz F, Paul Glezen W, et al. Impact of body mass index on immunogenicity of pandemic H1N1 vaccine in children and adults. J Infect Dis. 2014;210(8):1270–4.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  71. 71.

    Mavaddat N, Pharoah PD, Michailidou K, Tyrer J, Brook MN, Bolla MK, Wang Q, Dennis J, Dunning AM, Shah M, et al. Prediction of breast cancer risk based on profiling with common genetic variants. J Natl Cancer Inst. 2015;107(5):djv036.

  72. 72.

    Pruim RJ, Welch RP, Sanna S, Teslovich TM, Chines PS, Gliedt TP, Boehnke M, Abecasis GR, Willer CJ. LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics. 2010;26(18):2336–7.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We are thankful for the excellent technical work preparing samples and measuring nADA by the technical staff, Birgit Kassow and Ulla Abildtrup, of DMSC, Department of Neurology, Rigshospitalet, University of Copenhagen, Denmark; Marlies Jank of the Department of Neurology, Medical University of Innsbruck, Austria; and Anna Mattsson, of Karolinska Institutet, Stockholm, Sweden. We thank Niek de Vries, Faculty of Medicine, University of Amsterdam, Netherlands, for helpful discussions on HLA haplotypes.

Funding

This study was supported by the Innovative Medicines Initiative Joint Undertaking under grant agreement no. 115303 (ABIRISK), resources of which are composed of financial contribution from the European Union’s Seventh Framework Program (FP7/2007–2013) and EFPIA companies’ in kind contribution.

TFMA, CG, and BH were supported by the German Federal Ministry of Education and Research (BMBF) through the DIFUTURE consortium of the Medical Informatics Initiative Germany (grant 01ZZ1804A) and by the European Union’s Horizon 2020 Research and Innovation Programme (grant MultipleMS, EU RIA 733161). TFMA and BMM were supported by the BMBF through the Integrated Network IntegraMent, under the auspices of the e:Med Programme (grant 01ZX1614J). LA received funding from a clinician scientist program by the Munich Cluster for Systems Neurology (SyNergy). CG received a research stipend from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation). TO has research funding from the Swedish Research Council, the Swedish Brain Foundation, and Margaretha af Ugglas Foundation.

This research was performed using a reference panel assembled by the Type 1 Diabetes Genetics Consortium (T1DGC), a collaborative clinical study sponsored by the National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK), National Institute of Allergy and Infectious Diseases (NIAID), National Human Genome Research Institute (NHGRI), National Institute of Child Health and Human Development (NICHD), and Juvenile Diabetes Research Foundation International (JDRF). The data were supplied by the NIDDK Central Repositories. The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS.

Author information

Affiliations

Authors

Consortia

Contributions

TFMA drafted the conceptual study design, administered the project, curated data, devised the methodology for data analysis, conducted and interpreted the formal statistical data analysis, generated the data visualizations, and has prepared and revised the original manuscript draft. JL drafted the conceptual study design, curated and provided data on and investigated the Swedish patients, and revised the manuscript. DM drafted the conceptual study design, acquired funding, recruited and investigated the German patients, curated and provided their data, and revised the manuscript. MR curated data on and investigated the Swedish patients and revised the manuscript. CH curated data on and investigated the Swedish patients and reviewed the manuscript. VG performed the ELISA measurements of bADA levels and acquired data. MA contributed to nADA screening and titration, acquired data, and revised the manuscript. HH contributed to nADA screening and titration, acquired data, and revised the manuscript. LA recruited and acquired data on German patients and revised the manuscript. CG recruited and acquired data on German patients and revised the manuscript. BK recruited and acquired data on German patients and reviewed the manuscript. BMM devised the methodology for data analysis, supervised statistical data analyses, and reviewed the manuscript. PEHJ conducted nADA screening and titration, acquired data, and reviewed the manuscript. FS contributed to nADA screening and titration and reviewed the manuscript. IK provided data on the Swedish patients. TO provided data on the Swedish patients and reviewed the manuscript. MP acquired funding, administered the project, and reviewed the manuscript. SeS drafted the conceptual study design, acquired funding, administered the project, and revised the manuscript. FD drafted the conceptual study design, acquired funding, administered the project, devised the methodology for ADA measurement, conducted nADA screening and titration, acquired data, and revised the manuscript. AFH drafted the conceptual study design, acquired funding, administered and supervised the project, devised the methodology for ADA measurement, provided data on Swedish patients, and revised the manuscript. BH drafted the conceptual study design, acquired funding, supervised the project, provided data on German patients, and revised the manuscript. The author(s) read and approved the final manuscript.

Corresponding authors

Correspondence to Till F. M. Andlauer or Bernhard Hemmer.

Ethics declarations

Ethics approval and consent to participate

The ethics committees of the Regionala etikprövningsnämnden i Stockholm for Karolinska Institutet, Stockholm, Sweden, and of the Klinikum rechts der Isar, Technical University of Munich, Germany, approved the study, and all participants provided written informed consent.

Consent for publication

Not applicable.

Competing interests

TFMA, JL, DM, MR, VG, LA, CG, PEHJ, IK, and MP have no competing interests to declare.

CH is an employee of Sanofi Genzyme.

MA has received speaker honoraria and/or travel grants from Biogen, Novartis, Merck, and Sanofi Genzyme.

HH has participated in meetings sponsored by and received speaker honoraria or travel funding from Bayer, Biogen, Merck, Novartis, Sanofi-Genzyme, Siemens, and Teva, and received honoraria for acting as consultant for Biogen and Teva.

BK received a research grant and travel compensations from Novartis outside the submitted work.

FS has served on scientific advisory boards, been on the steering committees of clinical trials, served as a consultant, received support for congress participation, received speaker honoraria, or received research support for his laboratory from Biogen, Merck, Novartis, Roche, Sanofi Genzyme, and Teva.

TO has received unrestricted MS research grants, and honoraria for advisory boards/lectures from Biogen, Novartis, Sanofi, Merck, and Roche.

SeS is a former employee and has stocks and/or stock options in Novartis.

FD has participated in meetings sponsored by or received honoraria for acting as an advisor/speaker for Almirall, Alexion, Biogen, Celgene, Genzyme-Sanofi, Merck, Novartis Pharma, Roche, and TEVA ratiopharm. His institution has received research grants from Biogen and Genzyme Sanofi. He is section editor of the MSARD Journal (Multiple Sclerosis and Related Disorders).

AFH has received unrestricted research grants from Merck-Serono and BiogenIdec, served as consultant for Johnson & Johnson, and received honoraria for lectures by BiogenIdec and Sanofi-Aventis.

During the last 2 years, BH has served on scientific advisory boards for Novartis; he has served as DMSC member for AllergyCare, Polpharma, and TG therapeutics; he or his institution have received speaker honoraria from Desitin; his institution received research grants from Regeneron for MS research; he holds part of a patent for the detection of antibodies against KIR4.1 in a subpopulation of patients with MS. None of these conflicts are relevant to the topic of the study.

BMM and BH hold parts of a patent for genetic determinants of neutralizing antibodies to interferon.

Additional information

Publisher’s Note

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

Supplementary information

Additional file 1.

Previous measurements and design of new ADA measurements. Previous ADA measurements in the Swedish KI and German TUM cohorts per treatment preparation and distribution of samples for the new ADA measurements. For part of the TUM patients, only previous bADA measurements were available.

Additional file 2.

Additional details supporting the Methods section

Additional file 3.

New ADA measurements and design of the datasets for analyses. New ADA measurements in the Swedish KI and German TUM cohorts per treatment preparation and assignments of samples into the discovery and replication datasets. In the discovery and replication datasets, the first number indicates nADA and the second number bADA measurements. The distinction into negative and positive patients was made using nADA measurements.

Additional file 4.

Visualization and analysis of population stratification. For detailed figure legends, see the file.

Additional file 5.

Manhattan plots of the GWAS across IFNβ preparations. Manhattan plots of the (A-C) discovery-stage, (D-F) replication-stage, and (G-I) pooled discovery + replication GWAS. The red line between -log10p = 7 and -log10p = 8 indicates genome-wide significance; the top genome-wide significant variant is labeled with a red diamond.

Additional file 6.

Manhattan plots of the GWAS on patients treated with IFNβ-1a s.c. Manhattan plots of the (A-C) discovery-stage, (D-F) replication-stage, and (G-I) pooled discovery + replication GWAS. The red line between -log10p = 7 and -log10p = 8 indicates genome-wide significance; the top genome-wide significant variant is labeled with a red diamond.

Additional file 7.

Manhattan plots of the GWAS on patients treated with IFNβ-1b s.c. Manhattan plots of the (A-C) discovery-stage, (D-F) replication-stage, and (G-I) pooled discovery + replication GWAS. The red line between -log10p = 7 and -log10p = 8 indicates genome-wide significance; the top genome-wide significant variant is labeled with a red diamond.

Additional file 8.

Manhattan plots of the MHC region of the GWAS across IFNβ preparations. Manhattan plots of the (A-C) discovery-stage, (D-F) replication-stage, and (G-I) pooled discovery + replication GWAS, showing only the MHC region. The red line between -log10p = 7 and -log10p = 8 indicates genome-wide significance. For (A-C) discovery-stage plots, the prioritized variants are labeled with red diamonds for (D-F) replication-stage plots, the top genome-wide significant variant is labeled with a red diamond, and for (G-I) pooled discovery + replication plots, the replicated variants are labeled with red diamonds, and the top pooled variant is labeled in magenta.

Additional file 9.

Manhattan plots of the MHC region of the GWAS on patients treated with IFNβ-1a s.c. Manhattan plots of the (A-C) discovery-stage, (D-F) replication-stage, and (G-I) pooled discovery + replication GWAS, showing only the MHC region. The red line between -log10p = 7 and -log10p = 8 indicates genome-wide significance. For (A-C) discovery-stage plots, the prioritized variants are labeled with red diamonds, for (D-F) replication-stage plots, the top genome-wide significant variant is labeled with a red diamond, and for (G-I) pooled discovery + replication plots, the replicated variants are labeled with red diamonds, and the top variant from the pooled analysis is labeled in magenta.

Additional file 10.

Manhattan plots of the MHC region of the GWAS on patients treated with IFNβ-1b s.c. Manhattan plots of the (A-C) discovery-stage, (D-F) replication-stage, and (G-I) pooled discovery + replication GWAS, showing only the MHC region. The red line between -log10p = 7 and -log10p = 8 indicates genome-wide significance. For (A-C) discovery-stage plots, the prioritized variants are labeled with red diamonds, for (D-F) replication-stage plots, the top genome-wide significant variant is labeled with a red diamond, and for (G-I) pooled discovery + replication plots, the replicated variants are labeled with red diamonds, and the top variant from the pooled analysis is labeled in magenta.

Additional file 11.

Genomic inflation factors for all GWAS. Lambda = Median genomic inflation factor.

Additional file 12.

Table of the top GWAS associations. Variants prioritized in the discovery GWAS (bold font if replicated) and top variants from the pooled analysis of discovery + replication data. All effect sizes are relative to the minor allele. Bp = base pairs, MAF = minor allele frequency, beta = regression effect size, SE = standard error, P = p-value, cond. = conditional analysis, R2 = linkage disequilibrium r2.

Additional file 13.

Regional association plots of the top GWAS variants in the analysis across IFNβ preparations. Regional association plots of variants from the GWAS generated using LocusZoom v1.4 and the 1000 Genomes 1000G_Nov2014 EUR reference panel [72]. The color of dots indicates LD with the lead variant (pink). Gray dots represent signals with missing LD r2 values. If no LD information was present in the database on the top variant, LD with the variant showing the second-lowest p-value is indicated. The gray line indicates genome-wide significance. cM: centimorgan, chr: chromosome, Mb: mega base pairs.

Additional file 14.

Forest plots of the top GWAS variants and HLA alleles in the analysis across IFNβ preparations. Green: IFNβ-1a s.c., blue: IFNβ-1a i.m., orange: IFNβ-1b s.c., magenta: pooled discovery−/replication-stage analyses. D. = discovery, R. = replication, P. = pooled discovery + replication.

Additional file 15.

Results from stepwise conditional analyses. Results from stepwise conditional analyses on the pooled discovery and replication data. Genome-wide significant p-values are labeled in bold font. All effect sizes are relative to the minor allele. Bp = base pairs, MAF = minor allele frequency, beta = regression effect size, SE = standard error, P = p-value, R2 = linkage disequilibrium r2.

Additional file 16.

Results from eQTL analyses. Summary statistics as downloaded from GTEx v8 (https://gtexportal.org/). Significance thresholds are shown in the column Gene-level P threshold.

Additional file 17.

Results from MAGMA gene set analyses. FDR = 5% false discovery rate.

Additional file 18.

Regional association plots of the top GWAS variants in the analysis of IFNβ-1a s.c.-treated patients. Regional association plots of variants from the GWAS generated using LocusZoom v1.4 and the 1000 Genomes 1000G_Nov2014 EUR reference panel [72]. The color of dots indicates LD with the lead variant (pink). Gray dots represent signals with missing LD r2 values. If no LD information was present in the database on the top variant, LD with the variant showing the second-lowest p-value is indicated. The gray line indicates genome-wide significance. cM: centimorgan, chr: chromosome, Mb: mega base pairs.

Additional file 19.

Forest plots of the top GWAS variants and HLA alleles in the analysis of IFNβ-1a s.c.-treated patients. Green: IFNβ-1a s.c., blue: IFNβ-1a i.m., orange: IFNβ-1b s.c., magenta: pooled discovery−/replication-stage analyses. D. = discovery, R. = replication, P. = pooled discovery + replication.

Additional file 20.

Regional association plots of the top GWAS variants in the analysis of IFNβ-1b s.c.-treated patients. Regional association plots of variants from the GWAS generated using LocusZoom v1.4 and the 1000 Genomes 1000G_Nov2014 EUR reference panel [72]. The color of dots indicates LD with the lead variant (pink). Gray dots represent signals with missing LD r2 values. If no LD information was present in the database on the top variant, LD with the variant showing the second-lowest p-value is indicated. The gray line indicates genome-wide significance. cM: centimorgan, chr: chromosome, Mb: mega base pairs.

Additional file 21.

Forest plots of the top GWAS variants and HLA alleles in the analysis of IFNβ-1b s.c.-treated patients. Green: IFNβ-1a s.c., blue: IFNβ-1a i.m., orange: IFNβ-1b s.c., magenta: pooled discovery−/replication-stage analyses. D. = discovery, R. = replication, P. = pooled discovery + replication.

Additional file 22.

Association statistics of all replicated HLA alleles. AF = allele frequency, beta = regression effect size, SE = standard error, P = p-value, cond. = conditional analysis.

Additional file 23.

Results from analyses of previously published candidate SNPs and HLA alleles. Variants, alleles, and the respective p-values are labeled in bold font if they showed a one-sided p < 2.5 × 10− 3 (Bonferroni correction for 20 tests) either in the pooled analysis across treatment preparations or in the analysis of IFNβ-1a s.c.-treated or IFNβ-1b s.c.- treated patients. “nADA titer in nADA-positive” refers to an analysis of nADA titers restricted to nADA-positive patients. All effect sizes are relative to the minor allele. Bp = base pairs, MAF = minor allele frequency, AF = allele frequency, beta = regression effect size, SE = standard error, P = p-value, cond. = conditional analysis.

Additional file 24.

Treatment preparation-specific prediction of the presence of nADA in the replication data. Eight PRS, the top single GWAS variant, and the top HLA allele from the discovery stage were analyzed in the replication data using the covariates sex, age, treatment preparation, treatment duration, titration site, and eight ancestry components. Upper table: Prediction of the presence of nADA in IFNβ-1a s.c.-treated patients from the replication data using all ten prediction models based on analyses for IFNβ-1a s.c. in the discovery data. Lower table: Prediction of the presence of nADA in IFNβ-1b s.c.-treated patients from the replication data using all ten prediction models based on analyses for IFNβ-1b s.c. in the discovery data. Beta = regression effect size, SE = standard error, P = p-value.

Additional file 25.

Treatment preparation-specific prediction of the presence of nADA in the replication data: performance of single models. Eight PRS, the top single GWAS variant, and the top HLA allele from the discovery stage. Covariates: sex, age, treatment preparation, treatment duration, titration site, and ancestry components. The plots show the area under the receiver operating characteristic curve (AUC) and its 95% confidence interval (CI). Bonferroni = significant after Bonferroni correction for multiple testing; nominal = nominally significant (p < 0.05); n.s. = not significant

Additional file 26.

Treatment preparation-specific prediction of the presence of nADA in the replication data: comparison of top models. For each top model, the plots show either the AUC and its 95% CI or Nagelkerke’s pseudo-R2 and its 95% CI. Boxes show the prediction groups (replication data) and columns within each box the training data groups (discovery data). Bonferroni = significant after Bonferroni correction for multiple testing; nominal = nominally significant (p < 0.05); n.s. = not significant.

Additional file 27.

Comparison of allele frequencies for HLA-DRB1*15:01 and HLA-DRB1*04:01. The allele frequencies (AF) were queried from allelefrequencies.net on July 27th 2020 [48]. All four-digit European Silver and Gold populations with data on both HLA-DRB1*15:01 and HLA-DRB1*04:01 were used and populations with relative differences for both alleles are shown (i.e., with an AF above or below the average for one allele without the other allele being in the same group). Populations with an AF below the average for HLA-DRB1*15:01 and above the average for HLA-DRB1*04:01 are colored in green. Populations with an AF above the average for HLA-DRB1*15:01 and below the average for HLA-DRB1*04:01 are colored in magenta. In addition, the European populations with the highest or lowest AF for the respective allele (if not already present) as well as the largest German population and the Swedish SweHLA sample [47] are shown in gray.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Andlauer, T.F.M., Link, J., Martin, D. et al. Treatment- and population-specific genetic risk factors for anti-drug antibodies against interferon-beta: a GWAS. BMC Med 18, 298 (2020). https://doi.org/10.1186/s12916-020-01769-6

Download citation

Keywords

  • Multiple sclerosis
  • Interferon beta
  • Anti-drug antibodies
  • Human leukocyte antigen (HLA) system
  • Genetics
  • Genome-wide association study
  • Prediction