Pharmacogenomics of clinical response to Natalizumab in multiple sclerosis: a genome-wide multi-centric association study

Abstract

Background: Inter-individual differences in treatment response are marked in multiple sclerosis (MS). This is true for Natalizumab (NTZ), to which a subset of patients displays sub-optimal treatment response. We conducted a multi-centric genome-wide association study (GWAS), with additional pathway and network analysis to identify genetic predictors of response to NTZ.

Methods: MS patients from three different centers were included. Response to NTZ was dichotomized, nominating responders (R) relapse-free patients and non-responders (NR) all the others, over a follow-up of 4 years. Association analysis on ~ 4.7 M imputed autosomal common single-nucleotide polymorphisms (SNPs) was performed fitting logistic regression models, adjusted for baseline covariates, followed by meta-analysis at SNP and gene level. Finally, these signals were projected onto STRING interactome, to elicit modules and hub genes linked to response.

Results: Overall, 1834 patients were included: 119 from Italy (R = 94, NR = 25), 81 from Germany (R = 61, NR = 20), and 1634 from Sweden (R = 1349, NR = 285). The top-associated variant was rs11132400_T (p = 1.33 × 10^-6, OR = 0.58), affecting expression of several genes in the locus, like KLKB1. The interactome analysis implicated a module of 135 genes, with over-representation of terms like canonical WNT signaling pathway (p_adjust = 7.08 × 10^-6). Response-associated genes like GRB2 and LRP6, already implicated in MS pathogenesis, were topologically prioritized within the module.

Conclusion: This GWAS, the largest pharmacogenomic study of response to NTZ, suggested MS-implicated genes and Wnt/β-catenin signaling pathway, an essential component for blood-brain barrier formation and maintenance, to be related to treatment response.

Keywords: GRB2; LRP6; Multiple sclerosis; Natalizumab; Pharmacogenomics.

Fig. 1

Patients’ workflow. For each participating center, the number of patients for each subsequent application of filtering criteria is reported. The sample size for final study cohort is reported in bold character

Fig. 2

Manhattan plot. The Manhattan plot of − log₁₀(p) of associations from fixed-effect meta-analysis. The genome-wide significance level is set at p = 5 × 10^–8 (blue line), whereas suggestive significance threshold at p = 5 × 10^–5 (red line). rsIDs of the top five associated SNPs are marked

Fig. 3

Regional plot of top-associated locus from fixed-effects meta-analysis. The plot shows the genomic context of associated signal (lead variant: rs11132400_T, p = 1.33 × 10^–6, OR = 0.58) mapping to intronic region of F11-AS1 gene. The plot was generated with LocusZoom tool (http://locuszoom.sph.umich.edu). The − log₁₀(p) of associations from fixed-effect meta-analysis is reported on left y-axis, and the recombination rate on right y-axis, over the genomic position (hg19). Each symbol represents one SNP, with the most associated SNP marked in purple and shading of the other points based on the linkage disequilibrium metrics with the top SNP. Positions of genes are shown below the plot

Fig. 4

Detected network module. The final subnetwork resulting from the merge of modules, residing in top 1% of graph scores assigned by dmGWAS search algorithm, associated with response to NTZ. Color coding of nodes represents meta-analytic gene-level p values from MAGMA analysis, as indicated by the legend. Nodes/genes with association p > 0.05 were left white. Two genes which showed high level of centrality metrics and have already been implicated with MS (GRB2 and LRP6) are highlighted

Fig. 5

Enriched terms from over-representation analysis of genes in the detected network module against Gene Ontology Biological Process database. a Dotplot displaying the first 30 associated GO terms. p values were calculated from hypergeometric test, with adjustment for multiple testing with Benjamini–Hochberg procedure at FDR < 5%. The Count parameter in the legend illustrates the number of genes annotated to GO term and belonging to module. On x-axis, Gene Ratio reports the ratio between the number of genes in the module annotated to the term and the overall number of genes in the module (N = 135). b Enrichment map, reporting a graph-based representation of semantic similarity measures between GO terms enriched at FDR < 5% (N = 75, see Supplementary Table 3). Terms with high similarity tend to cluster together: the stronger the similarity, the shorter and thicker the edges. The color of nodes is coded according to p.adjust from hypergeometric test, as reported in the legend. Similarity between terms was computed with Jaccard correlation coefficient