Glatiramer Acetate modulates ion channels expression and calcium homeostasis in B cell of patients with relapsing-remitting multiple sclerosis

To investigate the effects of Glatiramer Acetate (GA) on B cells by an integrated computational and experimental approach. GA is an immunomodulatory drug approved for the treatment of multiple sclerosis (MS). GA effect on B cells is yet to be fully elucidated. We compared transcriptional profiles of B cells from treatment-naïve relapsing remitting MS patients, treated or not with GA for 6 hours in vitro, and of B cells before and after six months of GA administration in vivo. Microarrays were analyzed with two different computational approaches, one for functional analysis of pathways (Gene Set Enrichment Analysis) and one for the identification of new drug targets (Mode-of-action by Network Analysis). GA modulates the expression of genes involved in immune response and apoptosis. A differential expression of genes encoding ion channels, mostly regulating Ca2+ homeostasis in endoplasmic reticulum (ER) was also observed. Microfluorimetric analysis confirmed this finding, showing a specific GA effect on ER Ca2+ concentration. Our findings unveils a GA regulatory effect on the immune response by influencing B cell phenotype and function. In particular, our results highlight a new functional role for GA in modulating Ca2+ homeostasis in these cells.

Therefore, we decided to investigate GA effects, at the transcriptional level, on treatment-naïve RR-MS patients' B cells by using an integrated computational and experimental approach. Indeed, it has been shown that transcriptional responses of cells treated with small molecules can be used to elucidate their mechanism of action, in the lead optimization phase of drug discovery project and to reveal similarities among drugs, and quickly transfer indications for drug repositioning 7,8 . The Connectivity Map (CMAP), the largest peer-reviewed public database of gene expression profiles following treatment of five human cancer cell lines with 1309 different bioactive small molecules has been extensively used by both the academic and industrial communities 9 .
We hypothesized that measuring the transcriptional response on patients' B cells following drug treatment and comparing this response to those of drugs collected in the CMAP database would help revealing the mechanism of action of GA in these cells. To this end, we applied a bioinformatics tool named Mode of Action by NetwoRk Analysis (MANTRA) (http://mantra.tigem.it) that was recently described in the literature 7,8 . MANTRA is an automatic and robust approach that assesses similarity in gene expression profiles following drug treatment across thousands of small molecules from the CMAP database to predict similarities in drug effect and Mode of Action (MoA). Briefly, each transcriptional response to a drug treatment is represented as a list of genes ranked according to how much they change with respect to untreated cells. Transcriptional similarity among two drugs (i.e. two ranked lists of genes) is quantified by assessing how much two drugs tend to upregulate and downregulate the same genes. MANTRA represents transcriptional similarities among drugs graphically in the form of a network, where two drugs are connected if they are transcriptionally similar according to a metric based on the Gene Set Enrichment Analysis 7 . The position of a drug within the network provides insights about its MoA by exploiting previous knowledge on neighbouring drugs 7,8 . In addition, the network can be partitioned into communities consisting of groups of densely interconnected drugs sharing similar mode of action. The CMAP database is thus represented in MANTRA as a network of 1,309 drugs grouped into 106 communities 7 . MANTRA has been shown to be able to predict drug mechanism of action and functional targets from the analysis of gene expression profiles and for drug repositioning 10 .
Here, by performing a bioinformatics analysis of the CMAP dataset together with gene expression profiles of B cells from patients treated with GA, we revealed a possible role of this drug in modulating immune response, apoptosis, and ion channels expression, in particular voltage gated-and store-operated Ca 2+ channels. To support these findings, we investigated endoplasmic reticulum (ER) Ca 2+ content in B cells from MS patients treated with GA vs healthy subjects, and vs drugs-naïve RR-MS patients or treated with other first line DMTs.

Results
Generation of gene expression profiles. Ten consecutive treatment-naïve patients with RR-MS were enrolled between July 2013 and March 2014. One patient was successively excluded from the study for pregnancy. All patients were females. Mean age at enrolment ± SD was 31.5 ± 2.07. Mean age at MS onset was 29.16 ± 4. Mean EDSS ± SD was 2.5 ± 0.5. Specifically, as shown in Fig. 1A,B cells of each patient were isolated at baseline (treatment free cells, before starting GA) and after six months of 20 mg/mL GA in vivo administration (in vivo treated cells). At baseline, B cells were also cultured and treated or not in vitro with 100 µg/ml GA for six hours Fig. 1A. Acute (6 hours, in vitro) and chronic (6 months, in vivo) responses to GA treatment in B cells were then analyzed at the transcriptional level by Affymetrix microarrays (GSE110023). RNA quality of three patients was poor therefore microarray data were analyzed only in 6 patients (60%).
Functional pathway and MANTRA analysis of acute (in vitro) transcriptional response to GA treatment. Gene expression profiles of B cells following 6 hour GA treatment in vitro were compared to those of untreated cells in order to detect changes in gene expression caused by the drug treatment across the 6 patients. Genes were ordered according to their fold-change in treated versus untreated samples. To assess the molecular pathways modulated in vitro by GA, we performed Gene Set Enrichment Analysis (GSEA). GSEA uses one or more databases of set of genes (i.e. pathways) to identify those gene-sets which are significantly modulated by the treatment. In this study, we selected as pathway databases those of gene expression signatures of immune system cells (C7) and of Gene Ontology (GO) including biological processes, cellular components and molecular function (C5). Pathway enrichments were evaluated by their normalized enrichment score (NES), nominal P value, and FDR. The most significant GO pathways identified by GSEA included those involved in translation, which were up-regulated, and those involved in ion channels expression, including calcium, which were down-regulated ( Table 1). GSEA of immune system signatures highlighted up-regulation of genes specific to naïve T and B cells (Table 1).
We then analysed the ranked list of genes following acute treatment of B cells with GA to identify drugs inducing a significantly similar transcriptional response, and hence which could share a common mechanism of action with GA by means of the MANTRA online tool 7,8 . MANTRA identified 42 drugs eliciting a significantly similar transcriptional response to the in vitro GA treatment ( Fig. 1B and in Supplementary Table S1). The 42 drugs were part of 17 communities sharing similar mode of action 7 . Interestingly, among the 42 drugs we found agents with antiinfiammatory effects (estradiol, corticosterone etc.), widely used for the treatment of relapses in MS and drugs such as vigabatrin an antiepilectic able to suppress voltage-sensitive sodium channels 11 , and naltrexone, an opioid antagonist proposed against spasticity, pain and fatigue in MS 12 , which has been shown to prevent relapses in MS and to reduce the disease 13 .
Functional pathway and MANTRA analysis of chronic (in vivo) transcriptional response to GA treatment. Gene expression profiles of patients' B cells after six months of GA administration (in vivo) were compared to untreated cells of each patient at baseline (before starting GA treatment). We again performed GSEA and found a significant up-regulation of genes encoding ion channels (Table 1), such as genes belonging to the calcium voltage-gated channel subunit (CACN) and to the transient receptor potential (TRP) channel www.nature.com/scientificreports www.nature.com/scientificreports/ (TRPC1-5, TRM2, TRPV2) families. TRPC1, TRPC3 and TRPC4 are molecular component of store-operated Ca 2+ entry (SOCE), an ubiquitous Ca 2+ entry pathway that is activated in response to depletion of ER-Ca 2+ stores 14,15 . Ca 2+ signal, initiated by Ca 2+ release from ER is evoked by tyrosine kinase receptors phosphatase C activation 16 . Interestingly, GSEA shows also the induction of genes involved in protein tyrosine kinase activity, amine-binding function and in nervous system development and function.
A general down-regulation of pathways involved in negative regulation of cell death and cellular biosynthetic process was also identified by GSEA ( Table 1). As in the case of acute response, GSEA for pathways related to the immune system identified a significant up-regulation of genes linked to a less mature B cell phenotype. On the contrary, genes specific for lymphocytes vs others cells were down-regulated (naïve _bcell_vs_neutrophil_up, bcell_vs_monocyte_up), while in the acute response these were up-regulated.
We then performed MANTRA analysis of the B cells in vivo transcriptional response to GA. MANTRA identified a total of 49 drugs grouped into 18 communities (Fig. 1C and Supplementary Table S1). The communities with the largest number of drugs included community n. 14 (cyclin-dependent kinases 2 (CDK2) and Topoisomerase II inhibitors) with 6 drugs, community n. 89 (hemostatic agents) with 7 drugs, and community 90 (mainly anti-inflammatory, antibacterial and antiretroviral compounds) with 14 drugs.

Effect of GA on [Ca 2+ ] i homeostasis in B cells. In order to better investigate the computational results
hinting at a role of GA in modulating Ca 2+ channels expression, we analyzed the effect of GA treatment on ER [Ca 2+ ] i homeostasis. We recorded PBMCs form three healthy subjects and from 13 RR-MS patients: three disease-treatment free, four treated with GA, six with INFβ-1a and 1b (Table 2). Confocal immunofluorescence analysis showed that the anti-CD19 antibodies clearly recognized B cells both in PBMCs isolated from healthy subjects and GA-treated patients ( Fig. 2A,a,b). In this latter group, CD19 + cells appeared surrounded by several vesicles when observed in brightfield microscopy. (Fig. 2A,c,d).
Microfluorimetric experiments showed that the ER Ca 2+ content was significantly higher in PBMC isolated from GA-treated patients than in those isolated from healthy subjects, disease-treatment free, or INFβ1a and 1b treated patients (p = 0.02; Fig. 2B). To assess whether the observed changes in Ca 2+ content occurred in B cells, we recorded the ER Ca 2+ content in both FITC-CD19 + and FITC-CD19 − cells. Interestingly, microfluorimetric analysis revealed that the ER Ca 2+ release elicited by ATP+ thapsigargin was significantly higher in FITC-CD19 + cells than in FITC-CD19 − cells from patients treated with GA (p = 0.01), whereas the ER Ca 2+ content was similar in FITC-CD19 + and FITC-CD19 − cells obtained from healthy subjects (Fig. 2C,D).

Discussion
Our functional analyses of transcriptional profiles after both acute and chronic GA treatment revealed two main perturbated biological functions: (1) immune response and (2) ion channel activity.
The first GA action on immune response seems to be on B phenotype switch from a more to a less mature one. Indeed the acute response shows a relative decrease in plasma cells or memory B cells with concomitant expansion of B-cell precursors and/or naïve B cells. These data are in line with the literature 5,6 . Given the short time of GA exposure for acute treatment, we may only appreciate the initiation of the cells phenotypic switch phenomenon.
In chronic response, we were able to confirm that GA targets immune response by: (1)   www.nature.com/scientificreports www.nature.com/scientificreports/ The apoptosis induction, revealed by GSEA analysis, may be consequence of GA modulation on immune system rather than a direct effect.
Ion channels genes are down-regulated in acute treatment (in vitro) but up-regulated in chronic GA treatment (in vivo).
A possible explanation of this opposite response to GA may be a compensatory effect where GA treatment acutely down-regulates ion channels genes in B cells, which in turn causes a chronic upregulation of these genes to maintain homeostasis.  www.nature.com/scientificreports www.nature.com/scientificreports/ Ion channels genes, whose expression changes following GA treatment, belong the repertoire of ion-conducting proteins like potassium and chloride, Ca 2+ release-activated Ca 2+ (CRAC) and TRP channels. These last channels participate to ER Ca 2+ refilling 17,18 which represents the major Ca 2+ influx mechanism in B lymphocytes controlling antigen-mediated lymphocyte activation, cytokine/chemokine production, exocytosis, enzyme control, proliferation and apoptosis [19][20][21] . In line with the bioinformatic findings, intracellular Ca 2+ concentration was reduced in vitro by GA but increased ex vivo.
A role of GA on intracellular Ca 2+ has been already reported in platelets where, similar to our findings, intracellular Ca 2+ was reduced after 100 µg of GA in vitro 22 . However, if ion channels expression represents a trigger or a consequence of Ca 2+ homeostasis modifications is still to be elucidated.
After acute treatment, according to MANTRA, Brefeldin-A (BFA) and Ambroxol had the most similar gene expression profile to GA. BFA is an inhibitor of intracellular protein transport that leads to blockade the forward transport between the ER and the Golgi apparatus 23 . When anterograde transport is suppressed, retrograde transport from the Golgi to the ER becomes apparent. This BFA-induced retrograde transport depends on Ca 2+ sequestered into intracellular stores 24 . Therefore, it is possible that, due to GA modulation on Ca 2+ homeostasis, and through second messengers and transcription factors, the mode of action of the two drugs could be similar at transcriptional level.
Ambroxol is a mucolytic agent acting by inhibiting Na + channel and by reducing nitric oxide levels 25 . Ambroxol has also antiinflammatory properties including inhibition of oxidative stress, increased of local defense and reduction of pro-inflammatory cytokines 26 . By blocking voltage gated Na+ channels and Ca 2+ channels in animal models Ambroxol promotes axon regeneration in the injured adult mouse central nervous system 27 .
The compounds most transcriptionally similar to chronic GA treatment belonged to the Topoisomerase II and CDK2 inhibitors drug class, agents widely used for human cancer treatment that prevent unregulated cancer cell proliferation leading to apoptosis. Considering these proapoptotic drugs, MANTRA results are in line with GSEA analysis in vivo, that showed a downregulation of pathways involved in negative regulation of apoptosis (Table 1). Among these compounds, MANTRA identified Mitoxantrone (Fig. 1C), currently used for the treatment of secondary progressive or worsening RR-MS 28 . This is noteworthy, since it proves that our protocol was able to identify a previously proven drug against MS.
In conclusion, our combined computational and experimental study hints at a possible role of GA in B cells as a modulator of Ca 2+ homeostasis. We suggest that GA interferes with B cell activity by tuning ion channels expressions and Ca 2+ influx, crucial for the regulation of B cell differentiation, cytokine production and apoptosis. In a second phase, to ascertain computational results, we enrolled new healthy controls and patients treated with GA, with INFβ-1a and 1b, and treatment naïve.

B cells isolation and treatment.
Patients' peripheral blood samples were collected at baseline and after six months of GA in vivo administration.
Peripheral blood mononuclear cells (PBMC) were separated using Ficoll-Paque centrifugation (GE Healthcare); B cells were purified by autoMACS Pro Separator (Miltenyi biotec) using B cell Isolation kit II (Miltenyi Biotec, Bergisch Gladbach, Germany). The isolated cells, were analyzed using a FACS CANTO II (BD Biosciences, CA, USA) to confirm the purity of B cells (>98%). Cells at baseline were divided in two group. One underwent RNA extraction like B cells successively isolated from the same MS patient after six months of GA. The second group was seeded (1 × 10 6 cells/ml) with RPMI-1640 medium supplemented with 100 UI/ml penicillin, 100 μg/ml streptomycin (Life Technologies) and 5% autologous serum. Cultured cells were treated in vitro or not with 100 μg/ml GA for six hours in a humidified atmosphere containing 5% CO 2 at 37 C°. All experiments were carried out in duplicate.

Microarray experiments.
Microarray hybridization experiments on total RNA extracted with TRIzol reagent (Ambion), were performed at Coriell Insitute for Medical Research, Camden, NJ. The Affymetrix Gene-Chip Human Genome HG-U133A_2 hybridization platform (Affymetrix, Santa Clara, CA) was performed as described previously 30 . Prior to hybridization, RNA was reverted to cDNA using the NuGen Ovation RNA-Seq System V2 (NuGen Technologies, San Carlos, CA).
Bioinformatics analysis. Microarray analysis: GSEA and MANTRA. Microarray data were pre-processed using the Bioconductor package Affy and normalized with the RMA method 31 . Differentially expressed genes (DEGs) between conditions (GA-treated versus untreated B cells in vitro and in vivo) were identified using a Bayesian t-test 32 . For each P value, the Benjamini-Hochberg procedure was used to calculate the false discovery rate (FDR) to correct for multiple testing.
The Microarray data have been deposited in NCBIs Gene Expression Omnibus (GEO) and are accessible through GEOSeries accession number GSE110023 33 . www.nature.com/scientificreports www.nature.com/scientificreports/ Gene set enrichment analysis (GSEA) was then performed using the freely available software GSEA v2.0 from the Broad Institute 34 . To run the GSEA algorithm, RMA-normalized microarray data were used. The GSEA algorithm collapsed the probe sets into gene symbols (~13,300 genes) and ranked the genes on the basis of the fold change after GA treatment (in vitro or in vivo). An enrichment score was calculated for each gene-set by applying the Kolmogorov-Smirnov statistics as described in the GSEA algorithm. Gene-sets were collected from the Molecular signatures database (MSigDB) available as part of the GSEA software by restricting the output to the two MSigDB collections C5 (e.g. GO gene sets) and C7 (e.g. immunologic signatures) 34,35 . The significant gene sets were obtained by using p value < 0.01 and FDR < 0.1 as threshold. Ranking of the genes set was done using GSEA v2.0.
Microarrays before and after acute (six hours in vitro) and chronic (six months in vivo) treatment with GA were also analyzed with the Mode-of-action by Network Analysis (MANTRA) online tool [http:/mantra.tigem.it] 7 . MANTRA's output is a network where each drug is a node and drugs are connected to each other if they elicit a similar transcriptional response 7 . The 'distance' between the connected drugs is a measure of the similarity of their gene expression profiles 8 . The transcriptional distance threshold of ≤0.8 from GA treatment was established to obtain the significant compounds.
[Ca 2+ ] i Measurement. PBMCs of three healthy subjects and of four patients treated with GA, three with INFβ-1a, three with INFβ-1b, and three drug naïve patients, were isolated by Ficoll-Paque density gradients centrifugation.

Data Availability
All data generated or analysed during this study are included in this published article (and its Supplementary Information files).