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].
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, treatmentspecific results were evaluated. The secondary aim of the study was to use these genetic factors for the prediction of ADA development.

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 ADApositive 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-/nADApositive 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.

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

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 × 10 8 ) in the pooled analyses of discovery and replication samples, except for the superextended 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 crossvalidation was repeated 100 times, and the mean cutoff of the 100 repetitions was used as the final cutoff.
Nagelkerke's pseudo-R 2 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).

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 2 < 0.2 with the top signal was identified. The insertion replicated at genomewide 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 one-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 + nADA pres. 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). 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 2 = 0.39 in Swedish and r 2 = 0.51 in German patients). For bADA levels, SNP rs9272071 (2759 bp downstream from rs9271377; r 2 = 0.77 in Swedish and r 2 = 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)).
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 2 ≥ 0.71).
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 2 ≥ 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 genomewide 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).
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 2 ≥ 0.88) was the variant showing the most robust support for an association. In the meta-

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 Selected four-digit HLA alleles and extended haplotypes that were significantly associated (p < 3 × 10 −4 ) with an ADA measurement and replicated (p one-sided < 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 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.

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 nADApositive patients (Table 6 and Additional file 23).

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.   IFNb-1a  s.c. replication dataset) a b 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) specificities over models containing only the covariates (Table 7). 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 2 = 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 2 = 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 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 2 Nagelkerke's pseudo-R 2 ; KI Karolinska Institutet, Sweden; TUM Technical University of Munich, Germany 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 populationspecific 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 Nterminal 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 posttranslational 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 posttranslational 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 peptidebinding grooves of MHC molecules than IFNβ-1bderived 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 (MAF KI = 0.34, MAF TUM = 0.29) is less common than the two single alleles, especially compared to HLA-DQA1*01:02 (MAF KI = 0.42, MAF TUM = 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].
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 2 gl = 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. 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.