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

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.

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–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-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

Treatment preparation	IFNβ-1a i.m.		IFNβ-1a s.c.		IFNβ-1b s.c.
Cohort	KI Sweden	TUM Germany	KI Sweden	TUM Germany	KI Sweden	TUM Germany
N (%)	345 (24.7)	251 (18.4)	590 (42.3)	558 (40.9)	459 (32.9)	554 (40.6)
Mean age (SD)	46.7 (9.9)	40.1 (9.6)	44.1 (9.9)	38.9 (9.6)	45.4 (10.4)	41.4 (10.4)
Female sex (%)	216 (62.6)	191 (76.1)	440 (74.6)	406 (72.8)	334 (72.8)	390 (70.4)
Median treatment duration in months (MAD)	21.0 (8.1)	40.0 (20.4)	30.0 (15.0)	55.2 (23.1)	24.9 (12.9)	46.9 (23.4)
Progressive MS (%)	60 (17.4)	34 (13.5)	120 (20.3)	90 (16.1)	130 (28.3)	128 (23.1)
nADA positive (%)	45 (13.0)	41 (16.3)	204 (34.6)	188 (33.7)	245 (53.4)	255 (46.0)
Median nADA titer (MAD) nADA-positive samples	320 (280)	320 (280)	640 (600)	1280 (1240)	320 (280)	320 (280)
Median bADA level (MAD) all samples	13.9 (6.5)	9.6 (5.7)	23.2 (13.9)	16.3 (10.5)	35.8 (19.9)	29.4 (20.0)
Median bADA level (MAD) nADA-positive samples	63.8 (42.6)	25.8 (22.3)	109.0 (84.8)	115.0 (104.0)	69.0 (37.0)	73.8 (50.7)

N (%) refers to the entire cohort, the other percentages to the respective column. The nADA and bADA measurements shown here were obtained within the present study. Non-parametric summary statistics are provided for variables that were not normally distributed. Progressive MS = patients with a primary or secondary progressive disease course, as opposed to clinically isolated syndrome and relapsing-remitting MS. The dataset contained 1.6% primary progressive, 0.6% progressive-relapsing, and 18.2% secondary progressive MS patients. The frequency of nADA did not differ between progressive (35.2%) and other (35.5%) MS patients. Patients were diagnosed using the current McDonald criteria at the time of diagnosis. KI Karolinska Institutet, Sweden; TUM Technical University of Munich, Germany; SD standard deviation; MAD median absolute deviation

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

Fig. 3

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)

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-R 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

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 h = 0.47 (standard error (SE) = 0.15, p = 1.4 × 10−4) for the inverse-normal transformed bADA levels and h = 0.50 (SE = 0.15, p = 2.9 × 10−4) for the transformed nADA titers. The SNP heritability of the presence of nADA was h = 0.48 on the observed scale (SE = 0.15, p = 4.9 × 10−4) and, assuming an incidence of 0.35 for ADA, h = 0.79 (SE = 0.25) on a liability scale. Genetic correlations of bADA levels with nADA presence (r = 0.89, SE = 0.14, p = 1.2 × 10−3) and titers (r = 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) r < 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), p = 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

ADA	Variant	Chr.	Pos. (bp)	MA	MAF (KI)	MAF (TUM)	OR/β (Disc.)	p (Disc.)	OR/β (Repl.)	p(one-sided) (Repl.)	OR/β (Pool.)	p (Pool.)
nADA pres.	rs9281971	6	32,596,722	(T)7	0.36	0.39	0.59OR	1.9 × 10− 11	0.44 OR	2.4 × 10−08	0.55 OR	2.3 × 10−17
nADA titer	rs9271377	6	32,587,165	G	0.30	0.33	− 0.15β	4.0 × 10−11	− 0.17β	4.0 × 10−05	− 0.16β	1.5 × 10−14
nADA titer	rs9281971	6	32,596,722	(T)7	0.36	0.39	− 0.14β	1.4 × 10−10	− 0.23β	4.6 × 10−09	−0.16β	2.5 × 10−17
bADA level	rs9271377	6	32,587,165	G	0.30	0.33	− 0.23β	1.3 × 10−13	− 0.20β	8.5 × 10−05	−0.22β	1.2 × 10−16
bADA level	rs9272071	6	32,599,487	C	0.32	0.38	− 0.21β	1.8 × 10−11	− 0.28β	2.0 × 10−08	−0.23β	4.5 × 10−18

The top GWAS association signals that showed genome-wide significance in the discovery-stage analysis (α = 5 × 10−8) and replicated (α = 3 × 10−3) in the analysis across all three treatment preparations. For nADA presence, odds ratios are provided (marked by OR), and for quantitative ADA measures effect sizes (marked by β). For detailed association statistics, including conditional analyses, correlated HLA alleles, nearby genes, eQTL results, permutation p values, and preparation-specific association results, see Additional file 12. For locus-specific Manhattan plots of each locus, see Additional file 13. For forest plots of each association, including treatment preparation-specific effects, see Additional file 14. Abbreviations: Chr. chromosome; Pos. position in base pairs (build hg19); MA minor and effect allele; MAF minor allele frequency; KI Karolinska Institutet, Sweden; TUM Technical University of Munich, Germany; OR odds ratio; β effect size; Disc discovery; Repl. replication; Pool. pooled; pres. presence; (T) TTTTTTT

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 r = 0.39 in Swedish and r = 0.51 in German patients). For bADA levels, SNP rs9272071 (2759 bp downstream from rs9271377; r = 0.77 in Swedish and r = 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 ≥ r ≤ 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 (r ≥ 0.71).

Table 3

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

ADA	Variant	Chr.	Pos. (bp)	MA	MAF (KI)	MAF (TUM)	OR/β (Disc.)	p (Disc.)	OR/β (Repl.)	p(one-sided) (Repl.)	OR/β (Pool.)	p (Pool.)
nADA pres.	rs77278603	6	32,469,421	A	0.43	0.40	3.33OR	5.3 × 10−19	4.43OR	1.9 × 10−09	3.55OR	2.1 × 10−26
nADA pres.	rs9271700	6	32,593,198	G	0.42	0.39	3.16OR	8.6 × 10−19	5.08OR	1.2 × 10−10	3.48OR	5.4 × 10−27
nADA titer	rs77278603	6	32,469,421	A	0.43	0.40	0.38β	2.2 × 10−19	0.36β	8.1 × 10−10	0.37β	2.4 × 10−27
nADA titer	rs9271673	6	32,592,833	C	0.41	0.39	0.37β	3.0 × 10−19	0.38β	5.2 × 10−11	0.37β	2.1 × 10−28
bADA level	rs9281962	6	32,594,597	T	0.44	0.43	0.51β	6.0 × 10−22	0.51β	5.5 × 10− 12	0.51β	4.6 × 10−32

The top GWAS association signals that showed genome-wide significance in the discovery-stage analysis (α = 5 × 10−8) and replicated (α = 3 × 10−3) in the analysis of IFNβ-1a s.c.-treated patients. For nADA presence, odds ratios are provided (marked by OR), and for quantitative ADA measures effect sizes (marked by β). For detailed association statistics, including conditional analyses, correlated HLA alleles, nearby genes, eQTL results, permutation p values, and preparation-specific association results, see Additional file 12. For locus-specific Manhattan plots of each locus, see Additional file 18. For forest plots of each association, including treatment preparation-specific effects, see Additional file 19. Abbreviations: Chr. chromosome; Pos. position in base pairs (build hg19); MA minor and effect allele; MAF minor allele frequency; KI Karolinska Institutet, Sweden; TUM Technical University of Munich, Germany; OR odds ratio; β effect size; Disc discovery; Repl. replication; Pool. pooled; pres. presence

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 (r ≥ 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 (r ≥ 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 (r ≥ 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

ADA	Variant	Chr.	Pos. (bp)	MA	MAF (KI)	MAF (TUM)	OR/β (Disc.)	p (Disc.)	OR/β (Repl.)	p(one-sided) (Repl.)	OR/β (Pool.)	p (Pool.)
nADA pres.	rs28366299	6	32,560,870	A	0.20	0.19	3.11OR	2.1 × 10−12	5.84OR	5.2 × 10−09	3.56OR	6.6 × 10−19
nADA titer	rs28366299	6	32,560,870	A	0.20	0.19	0.40β	1.8 × 10−16	0.52β	1.1 × 10−09	0.43β	5.3 × 10−24
nADA titer	rs9272775	6	32,610,257	C	0.22	0.23	0.38β	2.5 × 10−16	0.51β	3.2 × 10−10	0.41β	2.5 × 10−24
bADA level	rs78279385	6	32,451,758	A	0.23	0.25	0.39β	2.7 × 10−14	0.39β	2.0 × 10−06	0.39β	5.5 × 10−19
bADA level	rs9272775	6	32,610,257	C	0.22	0.23	0.39β	1.0 × 10−13	0.46β	1.1 × 10−07	0.41β	1.6 × 10−19

The top GWAS association signals that showed genome-wide significance in the discovery-stage analysis (α = 5 × 10−8) and replicated (α = 3 × 10−3) in the analysis of IFNβ-1b s.c.-treated patients. For nADA presence, odds ratios are provided (marked by OR), and for quantitative ADA measures effect sizes (marked by β). For detailed association statistics, including conditional analyses, correlated HLA alleles, nearby genes, eQTL results, permutation p values, and preparation-specific association results, see Additional file 12. For locus-specific Manhattan plots of each locus, see Additional file 20. For forest plots of each association, including treatment preparation-specific effects, see Additional file 21. Abbreviations: Chr. chromosome; Pos. position in base pairs (build hg19); MA minor and effect allele; MAF minor allele frequency; KI Karolinska Institutet, Sweden; TUM Technical University of Munich, Germany; OR odds ratio; β effect size; Disc discovery; Repl. replication; Pool. pooled; pres. presence

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 r ≥ 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 (r ≥ 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

ADA	Allele/HT	AF (KI)	AF (TUM)	OR/β (Disc.)	p (Disc.)	OR/β (Repl.)	p(one-sided) (Repl.)	OR/β (Pool.)	p (Pool.)
All preparations: risk alleles
bADA levels	HLA-DQA1*01:02	0.42	0.36	0.18β	1.6 × 10−08	0.17β	3.7 × 10−04	0.17β	4.8 × 10−11
bADA levels	B7-DQ6	0.18	0.17	0.17β	8.6 × 10−06	0.21β	3.8 × 10−04	0.18β	2.7 × 10−08
All preparations: protective alleles
bADA levels	DR3-DQ2	0.13	0.12	− 0.23β	9.6 × 10−08	− 0.21β	6.7 × 10−04	− 0.23β	5.0 × 10−10
IFNβ-1a s.c.: risk alleles
nADA pres.	DR15-DQ6	0.34	0.30	2.73OR	3.1 × 10−14	3.41OR	1.5 × 10−07	2.89OR	7.4 × 10−20
nADA titer	DR15-DQ6	0.34	0.30	0.36β	1.3 × 10−15	0.32β	2.3 × 10−07	0.35β	3.6 × 10−21
bADA levels	DR15-DQ6	0.34	0.30	0.44β	6.8 × 10−15	0.45β	4.4 × 10−08	0.44β	3.5 × 10−21
IFNβ-1a s.c.: protective alleles
nADA pres.	DR3-DQ2	0.13	0.12	0.40OR	9.1 × 10−07	0.29OR	4.0 × 10−04	0.37OR	3.7 × 10−09
nADA titer	DR3-DQ2	0.13	0.12	− 0.31β	8.6 × 10−08	− 0.28β	5.0 × 10−04	− 0.30β	3.4 × 10−10
bADA levels	DR3-DQ2	0.13	0.12	− 0.41β	2.1 × 10−08	− 0.40β	1.5 × 10−04	− 0.41β	2.5 × 10−11
IFNβ-1b s.c.: risk alleles
nADA pres.	HLA-DRB1*04:01	0.09	0.07	6.82OR	3.4 × 10−13	14.7OR	1.7 × 10−07	7.95OR	1.4 × 10−18
nADA pres.	DR4-DQ3	0.07	0.05	6.23OR	1.2 × 10−09	14.7OR	6.1 × 10−06	7.35OR	1.5 × 10−13
nADA titer	HLA-DRB1*04:01	0.09	0.07	0.56β	4.1 × 10−15	0.56β	9.2 × 10−06	0.56β	3.7 × 10−19
nADA titer	DR4-DQ3	0.07	0.05	0.50β	1.6 × 10−09	0.54β	1.3 × 10−04	0.51β	1.9 × 10−12
bADA levels	HLA-DRB1*04:01	0.09	0.07	0.62β	4.6 × 10−15	0.61β	1.8 × 10−06	0.62β	8.7 × 10−20
bADA levels	DR4-DQ3	0.07	0.05	0.56β	4.9 × 10−10	0.53β	2.2 × 10−04	0.56β	9.0 × 10−13

Selected four-digit HLA alleles and extended haplotypes that were significantly associated (p < 3 × 10−4) with an ADA measurement and replicated (p < 1 × 10−3). Alleles that are part of one of the listed extended haplotypes and showed a similar or weaker association than the haplotypes and which did not remain significant when conditioning for the haplotypes are not displayed separately. For nADA presence, odds ratios are provided (marked by OR), and for quantitative ADA measures effect sizes (marked by β). For a detailed table of all results, see Additional file 22. For forest plots of each association, including treatment preparation-specific effects, see Additional files 14, 19, and 21. HT haplotype; AF allele frequency; KI Karolinska Institutet, Sweden; TUM Technical University of Munich, Germany; β effect size; OR odds ratio; Disc discovery; Repl. replication; Pool. pooled; pres.presence. Abbreviations of the haplotypes: B7-DQ6, HLA-B*07:02 + HLA-DRB1*15:01 + HLA-DQA1*01:02 + HLA-DQB1*06:02; DR15-DQ6, HLA-DRB1*15:01 + HLA-DQA1*01:02 + HLA-DQB1*06:02; DR3-DQ2, HLA-DRB1*03:01 + HLA-DQA1*05:01 + HLA-DQB1*02:01; DR4-DQ3, HLA-DRB1*04:01 + HLA-DQA1*03:01 + HLA-DQB1*03:02

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

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

Measurement	SNP/HLA allele	First publication	Published effect direction	OR/β (pooled)	p(one-sided) (all prep.)	p(one-sided) (IFNβ-1a s.c.)	p(one-sided) (IFNβ-1b s.c.)
nADA presence	rs2454138-A	[15]	risk	1.40OR	1.59 × 10−05	9.65 × 10−01	2.41 × 10−16
nADA presence	rs522308-T	[15]	risk	1.39OR	2.06 × 10−05	9.80 × 10−01	8.76 × 10−17
nADA presence	rs9272105-A	[14]	protective	0.62OR	1.21 × 10−13	2.76 × 10−12	1.14 × 10−05
nADA presence	HLA-DQA1*02:01	[16]	risk	1.27OR	2.12 × 10−02	5.74 × 10−01	2.99 × 10−04
nADA presence	HLA-DRB1*04:01	[12]	risk	1.76OR	9.46 × 10−06	8.71 × 10−01	1.62 × 10−19
nADA titer	HLA-DRB1*04:08	[12]	risk	1.14β	7.50 × 10−09	1.56 × 10−05	6.21 × 10−05
nADA titer in nADA-positive	HLA-DRB1*04:08	[12]	risk	0.48β	2.32 × 10−03	2.45 × 10−02	3.59 × 10−04
nADA presence	HLA-DRB1*07:01	[10]	risk	1.28OR	1.90 × 10−02	5.47 × 10−01	3.03 × 10−04
nADA presence	HLA-DRB1*08:01	[13]	risk	1.38OR	2.52 × 10−02	1.66 × 10−05	9.14 × 10−01
nADA presence	HLA-DRB1*15:01	[13]	risk	1.32OR	6.27 × 10−05	4.24 × 10−19	1.00 × 100
bADA levels	HLA-DRB1*16:01	[11]	risk	0.54β	7.10 × 10−06	1.12 × 10−02	3.34 × 10−05
nADA titer in nADA-positive	HLA-DRB1*16:01	[11]	risk	0.62β	2.11 × 10−04	2.93 × 10−01	5.07 × 10−05
nADA presence	HLA-DQA1*05:01	[13] / [16]	protective	0.68OR	3.09 × 10−07	3.15 × 10−09	4.59 × 10−03
nADA presence	HLA-DQB1*02:01	[16]	protective	0.65OR	6.24 × 10−06	1.46 × 10−09	1.94 × 10−01
nADA presence	HLA-DRB1*03:01	[11] / [16]	protective	0.65OR	5.43 × 10−06	1.55 × 10−09	1.83 × 10−01
nADA presence	HLA-DRB1*04:04	[11]	protective	0.56OR	2.29 × 10−03	6.91 × 10−03	1.03 × 10−01

The table shows previously published SNPs and HLA alleles that 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. P values below this threshold are labeled in bold font. For nADA presence, odds ratios are provided (marked by OR), and for quantitative ADA measures effect sizes (marked by β), both refer to the pooled analysis across treatment preparations. For a detailed table of all results, see Additional file 23. β effect size; OR odds ratio; prep treatment preparations. “nADA titer in nADA-positive” refers to an analysis of nADA titers restricted to nADA-positive patients

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-R, sensitivities, and specificities over models containing only the covariates (Table 7).

Table 7

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

Preparation	Model	Cohort	OR	95% CI	p	AUC	R2	Sensitivity	Specificity
IFNβ-1a s.c.	Without genetics	KI				0.65	0.08	0.66	0.52
IFNβ-1a s.c.	Without genetics	TUM				0.60	0.06	0.71	0.42
IFNβ-1a s.c.	PRS nADA presence 5 × 10−08	KI	3.89	2.35–6.45	1.44 × 10−07	0.85	0.42	0.78	0.78
IFNβ-1a s.c.	PRS nADA presence 5 × 10−08	TUM	2.56	1.56–4.21	2.11 × 10−04	0.76	0.24	0.68	0.65
IFNβ-1a s.c.	PRS top 30% vs. bottom 30%	KI	73.86	11.77–463.61	4.42 × 10−06	0.91	0.59	0.78	0.90
IFNβ-1a s.c.	PRS top 30% vs. bottom 30%	TUM	13.78	3.00–63.28	7.45 × 10−04	0.83	0.38	0.80	0.76
IFNβ-1a s.c.	SNP rs77278603 additive coding	KI	4.49	2.41–8.36	2.14 × 10−06	0.82	0.36	0.74	0.74
IFNβ-1a s.c.	SNP rs77278603 additive coding	TUM	3.88	1.78–8.47	6.67 × 10−04	0.73	0.21	0.63	0.74
IFNβ-1a s.c.	SNP rs77278603-A dominant coding	KI	9.16	2.48–33.79	8.79 × 10−04	0.78	0.31	0.57	0.76
IFNβ-1a s.c.	SNP rs77278603-A dominant coding	TUM	3.85	1.10–13.49	3.51 × 10−02	0.72	0.20	0.56	0.68
IFNβ-1b s.c.	Without genetics	KI				0.70	0.15	0.48	0.82
IFNβ-1b s.c.	Without genetics	TUM				0.58	0.02	0.82	0.20
IFNβ-1b s.c.	PRS nADA presence 1 × 10−06	KI	2.40	1.45–3.97	6.46 × 10−04	0.78	0.33	0.57	0.85
IFNβ-1b s.c.	PRS nADA presence 1 × 10−06	TUM	2.15	1.43–3.23	2.28 × 10−04	0.73	0.22	0.73	0.58
IFNβ-1b s.c.	PRS top 30% vs. bottom 30%	KI	10.16	2.30–44.95	2.25 × 10−03	0.83	0.46	0.58	0.87
IFNβ-1b s.c.	PRS top 30% vs. bottom 30%	TUM	5.97	2.03–17.52	1.14 × 10−03	0.78	0.33	0.69	0.75
IFNβ-1b s.c.	SNP rs28366299 additive coding	KI	4.51	1.72–11.80	2.14 × 10−03	0.77	0.31	0.57	0.82
IFNβ-1b s.c.	SNP rs28366299 additive coding	TUM	6.91	3.18–15.03	1.07 × 10−06	0.79	0.32	0.74	0.61
IFNβ-1b s.c.	SNP rs28366299-A dominant coding	KI	9.78	2.68–35.74	5.62 × 10−04	0.83	0.40	0.62	0.83
IFNβ-1b s.c.	SNP rs28366299-A dominant coding	TUM	7.56	3.01–19.02	1.71 × 10−05	0.80	0.33	0.77	0.57

Predictors in the model without genetics: sex, age, treatment duration, and titration site. The genetic models contained the same base model plus the indicated genetic factors and ancestry components. The top models are indicated in bold font. OR odds ratio; CI 95% confidence interval; p p value of the genetic component; AUC area under the receiver operating characteristic curve; R Nagelkerke's pseudo-R; KI Karolinska Institutet, Sweden; TUM Technical University of Munich, Germany

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-R = 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-R = 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–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—h = 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.

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.

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

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

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 r 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 r 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-R 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.