Immunophenotype and Transcriptome Profile of Patients With Multiple Sclerosis Treated With Fingolimod: Setting Up a Model for Prediction of Response in a 2-Year Translational Study

Background Fingolimod is a functional sphingosine-1-phosphate antagonist approved for the treatment of multiple sclerosis (MS). Fingolimod affects lymphocyte subpopulations and regulates gene expression in the lymphocyte transcriptome. Translational studies are necessary to identify cellular and molecular biomarkers that might be used to predict the clinical response to the drug. In MS patients, we aimed to clarify the differential effects of fingolimod on T, B, and natural killer (NK) cell subsets and to identify differentially expressed genes in responders and non-responders (NRs) to treatment. Materials and methods Samples were obtained from relapsing–remitting multiple sclerosis patients before and 6 months after starting fingolimod. Forty-eight lymphocyte subpopulations were measured by flow cytometry based on surface and intracellular marker analysis. Transcriptome sequencing by next-generation technologies was used to define the gene expression profiling in lymphocytes at the same time points. NEDA-3 (no evidence of disease activity) and NEDA-4 scores were measured for all patients at 1 and 2 years after beginning fingolimod treatment to investigate an association with cellular and molecular characteristics. Results Fingolimod affects practically all lymphocyte subpopulations and exerts a strong effect on genetic transcription switching toward an anti-inflammatory and antioxidant response. Fingolimod induces a differential effect in lymphocyte subpopulations after 6 months of treatment in responder and NR patients. Patients who achieved a good response to the drug compared to NR patients exhibited higher percentages of NK bright cells and plasmablasts, higher levels of FOXP3, glucose phosphate isomerase, lower levels of FCRL1, and lower Expanded Disability Status Scale at baseline. The combination of these possible markers enabled us to build a probabilistic linear model to predict the clinical response to fingolimod. Conclusion MS patients responsive to fingolimod exhibit a recognizable distribution of lymphocyte subpopulations and a different pretreatment gene expression signature that might be useful as a biomarker.

Background: Fingolimod is a functional sphingosine-1-phosphate antagonist approved for the treatment of multiple sclerosis (MS). Fingolimod affects lymphocyte subpopulations and regulates gene expression in the lymphocyte transcriptome. Translational studies are necessary to identify cellular and molecular biomarkers that might be used to predict the clinical response to the drug. In MS patients, we aimed to clarify the differential effects of fingolimod on T, B, and natural killer (NK) cell subsets and to identify differentially expressed genes in responders and non-responders (NRs) to treatment.
Materials and methods: Samples were obtained from relapsing-remitting multiple sclerosis patients before and 6 months after starting fingolimod. Forty-eight lymphocyte subpopulations were measured by flow cytometry based on surface and intracellular marker analysis. Transcriptome sequencing by next-generation technologies was used to define the gene expression profiling in lymphocytes at the same time points. NEDA-3 (no evidence of disease activity) and NEDA-4 scores were measured for all patients at 1 and 2 years after beginning fingolimod treatment to investigate an association with cellular and molecular characteristics.
results: Fingolimod affects practically all lymphocyte subpopulations and exerts a strong effect on genetic transcription switching toward an anti-inflammatory and antioxidant response. Fingolimod induces a differential effect in lymphocyte subpopulations after 6 months of treatment in responder and NR patients. Patients who achieved a promoting Th1 and Th17 differentiation in vitro and in vivo (6). These phenomena could explain the relative increase in Tregs and the reduction of Th1 and Th17 cells in response to fingolimod (7). In addition, an increase in the percentage of B cells with regulatory functions has also been observed in treated patients (8).
The assessment of the effectiveness of treatment for MS is complex and requires the use of clinical and imaging criteria such as the annualized relapse rate (ARR), the cumulative progression of disability, and the number and volume of new T2 or gadoliniumenhanced lesions. NEDA-3 (no evidence of disease activity) is a composite measure of disease activity in MS [including relapses, disability progression, and magnetic resonance imaging (MRI) activity] that is used as a secondary outcome measure in clinical trials for new disease-modifying therapies (DMTs) (9). These measures are focused on the inflammatory process and do not take into account other important aspects of MS pathophysiology such as neurodegeneration. Recently, a fourth measure, brain volume loss, has been suggested to be added to NEDA because of its relation to cognitive decline in MS patients. This new composite measure is called NEDA-4 and has advantages over the conventional NEDA-3 for predicting subsequent disability and structural damage (10). Despite efforts, the use of response biomarkers at the clinical level is limited, and studies designed to identify patients who require an exhaustive follow-up are necessary.
To search for biomarkers of response and to improve our understanding of mechanisms of action of fingolimod, we carried out a translational study of the cellular and transcriptional characteristics of MS patients treated with fingolimod. Patients were classified according to the clinical response which was measured by the NEDA-3 and NEDA-4 status at 1 and 2 years, and samples were analyzed by flow cytometry immunophenotyping and RNA next-generation sequencing before and after 6 months of treatment. The possible markers obtained by a first statistical approach were used to build a model of prediction that might be useful to anticipate the clinical response to fingolimod.

MaTerials anD MeThODs
Patients, control subjects, and samples Enrollment was limited to patients 18-55 years of age with a diagnosis of relapsing-remitting multiple sclerosis (RRMS) according Predicting Response to Fingolimod in Patients With MS Frontiers in Immunology | www.frontiersin.org July 2018 | Volume 9 | Article 1693 to McDonald 2010 criteria (11) with an indication for treatment with fingolimod and an entry score of 0-5.0 in the Expanded Disability Status Scale (EDSS). MS patients categorized as progressive or previously treated with an S1P receptor agonist were excluded. A total of 40 patients were recruited, and all of them signed an informed consent that was previously approved by the Ethics Committee of Puerta de Hierro University Hospital. The patients were treated with fingolimod 0.5 mg once daily, according to the current recommendations of use for fingolimod. The treated patients were subdivided into two subgroups according to the previous use of natalizumab (NTZ): NTZ group (n = 18) and no-NTZ group (n = 22) because switching from NTZ has been associated with a risk of MS reactivation (12). The patients of NTZ-group started fingolimod treatment after a washout period of 1-2 months and in patients of No-NTZ group (4 patients from glatiramer acetate and 18 patients from IFN-β), no washout period was needed according to the current guidelines. Venous blood samples were obtained from patients immediately before starting treatment with fingolimod and 6 months after starting therapy. Samples were also obtained from 10 HC and 10 naïve MS patients (who had never received any disease-modifying treatment) matched for age and sex with treated individuals. The demographic characteristics of the patients and controls and baseline differences between subgroups are shown in Table 1.

Mri Measures and clinical response
Expanded Disability Status Scale, clinical relapses, and ARR were evaluated in all patients at baseline, 12 and 24 months after starting therapy. At the same times, a 1.5 T brain MRI was performed. The measured variables were as follows: number and volume of gadolinium-enhanced T1 lesions (GdE), number of new or enlarged T2-weighted lesions (T2w), total T2 and T1-weighted lesion volume, and annual brain volume loss (BVL), which was obtained using the SIENA method. The percentage of patients with MRI activity measured by the appearance of new GdE and/ or T2w lesions was calculated. Confirmed disability progression (CDP) was defined by an increase of at least 1 EDSS point sustained during 1 year. NEDA-3 was calculated according to published parameters (no MRI activity, no relapses, and no CDP). NEDA-4 was calculated as NEDA-3 plus BVL <0.4%.
The patients who achieved NEDA-4 status were defined as responders at 1 or 2 years and the patient who did not achieve it were defined as non-responders (NRs) at 1 or 2 years.

antibodies and reagents for Flow cytometry
Peripheral blood mononuclear cells (PBMCs) were obtained from all blood samples by centrifugation on a Ficoll-Hypaque gradient and frozen in liquid nitrogen until use. For surface and intracellular staining, PBMCs were incubated using seven combinations of the fluorochrome-conjugated antibodies (FCAs) to define 48 lymphocyte subpopulations of B, T, and natural killer (NK) cells, as shown in Table S1 in Supplementary Material.
For intracellular staining, lymphocytes were stimulated for 4 h with phorbol 12-myristate 13-acetate (100 ng/ml) and ionomycin Io (1 µg/ml) in the presence of brefeldin A (10 µg/ml) from Sigma (St. Louis, MI, USA). The cells were stained with FCAs for surface markers before fixation and permeabilization with the specific INTRACELL-Kit (Immunostep). Fluorochromeconjugated isotypes matched Abs were used as controls. All FCAs were obtained from Miltenyi Biotec (Auburn, CA, USA). The lymphocyte subpopulations and cytokine profiles were analyzed using a MACSQUANT flow cytometer (Miltenyi Biotec) and MACSQuantify 2.5 and FlowJo (tree Star) software. The gating strategy is shown in Figures S1.1-S1.7 in Supplementary Material.

Total rna isolation and cDna library Preparation for Transcriptome sequencing (rna-seq)
Once the NEDA-4 status was determined in all patients, 10 of them were selected according to clinical response, five responders and five NRs at 1 year. In regard to response at 2 years, four of them were responders and six NRs. PBMCs obtained from these patients before and at 6 months of treatment were used to RNA sequencing as well as the samples from five naïve MS patients and five HC. The demographic characteristics of these patients and controls in this subset were similar to the whole sample and are shown in Table 1, as well as the baseline differences between responder and NR patients.
Total RNA was extracted using the Maxwell ® 16 LEV simply RNA Cells kit performed on Promega's robotic platform according to the manufacturer's instructions. The polyA+ fraction was purified and randomly fragmented, converted to doublestranded cDNA and processed through subsequent enzymatic treatments of end-repair, dA-tailing, and ligation to adapters as in Illumina's "TruSeq Stranded mRNA Sample Preparation Part # 15031047 Rev. D" kit. The adapter-ligated library was completed by PCR with Illumina PE primers. The resulting purified cDNA library was applied to an Illumina flow cell for cluster generation and sequenced on an Illumina HiSeq2500 by following the manufacturer's protocols. Image analysis, per-cycle base calling, and quality score assignment were performed with Illumina Real Time Analysis software. Conversion of Illumina BCL files to bam format was performed with the Illumina2bam tool (Wellcome Trust Sanger Institute-NPG). Unaligned single read sequences are provided in BAM files. Read qualities are provided in the standard Qscore Sanger encoding (Illumina 1.9, ASCII offset of 33).

statistic analysis
A logistic regression model was fitted to the lymphocyte subpopulations data to determine the relative importance of each population as a predictor of response for the variables of interest (NEDA3, NEDA-4, MRI activity, clinical activity, and CDP). The model was fitted with the rstanarm package (13). The mean of each posterior distribution was used to create the heatmap with the pheatmap package (14). Four variables from the first approach were selected to first construct a two-dimensional data representation by nonparametric multidimensional scaling (Seuclidean method) and then to visualize differences between

Bioinformatics analysis
To assess the effect of fingolimod on PBMCs from MS patients, the differentially expressed genes (DEGs) of samples before and after 6 months of treatment were analyzed. In samples collected at 6 months, DEGs between responder and NR patients were determined to analyze the differential effect of fingolimod according to the clinical response. To search for possible predictor genes of a clinical response, we analyzed DEGs between responder and NR patients using samples that were collected before treatment. The sequencing quality was analyzed with FastQC (15). Reads were aligned to the human genome (GRch38/hg38) (16) using STAR (17) and Samtools 0.1.19.0 (18), and transcript assembly, abundance estimation and differential expression were calculated with DESeq2 (19). The estimated p-value was corrected to account for multiple hypothesis testing using Benjamini and Hochberg false discovery rate (FDR) adjustment. Genes with an FDR less than or equal to 0.05 were considered differentially expressed.
In addition, we conducted a gene set enrichment analysis (GSEA), a publicly available computational method that employs pre-defined sets of genes to identify statistically significant differences between two biological conditions. This method can detect the upregulation or downregulation of recognized cellular pathways and biological processes under experimental conditions. GSEA was performed using GSEA software v2.2.1 obtained from the Broad Institute (20) and following the developer's protocol (21). The previously obtained DEGs were ranked according to their t-statistics. This ranked file was used as input for the enrichment analysis. All basic and advanced fields were set to default and only those gene sets significantly enriched at an FDR q-values <0.25 were considered. To investigate how selected pathways were modulated by fingolimod therapy, we used the models from the wikipaths platform (22) and drew the differential regulation between responder and NR patients using the PathVisio program (23).

Quantitative Pcr (qPcr) Validation
Based on the RNA-seq analysis, we selected 10 significant DEGs for validation using qPCR analysis. Each RNA sample (500 ng) was reverse-transcribed using the NZY ® First-Strand cDNA Synthesis Kit. qPCR was performed in an LC480-Roche. The reaction conditions consisted of 6 µl cDNA and 1 µl of a 5 µM forward and reverse primer mix (Table S2 in Supplementary Material) and 5 µl of Light Cycler ® SYBR green I Master Mix.
Each sample was analyzed in duplicate for each individual reaction. The PCR program consisted of denaturation at 95°C for 5 min, 45 cycles of denaturing for 10 s at 95°C, annealing for 10 s at 60°C and extension for 10 s. A melting curve was generated by denaturation for 5 s at 95°C, annealing for 1 min at 65°C, and a ramp rate of 0.11°C/s up to 97°C. The gene expression of three endogenous genes (Table S2 in Supplementary Material) was used to normalize the target gene expression, and relative quantification was established using the 2 −ΔΔCt method (24).

Design of Predictive Model
From the cytometry and RNA-seq experiments, variables that showed statistically significant differences between responder and NR were selected. To provide a tool for prediction of response to fingolimod, we fitted a probabilistic linear model to selected data and evaluated the resulting predictions defined as follows: The model is based on the Bernoulli distribution or dichotomous distribution, in which the dependent variable can take the value Yi = 1 for the probability of success (responder to 2 years) or Yi = 0 for the probability of failure (NR to 2 years). The distribution of the sample is located in a cloud of points in such a way that their location divides the observations into two groups according to response. The calculated value of the regression line measures the response probability. Values close to 0 have a low probability of occurrence, while values close to 1 indicate that there is a high probability of a good response to treatment.

Validation of Model
The expression levels by qPCR of genes used in the model building were measured on the pretreatment samples obtained from the remaining 30 patients who were not used for the analysis of RNA-seq. All of them have available data on the clinical evolution up to 2 years (NEDA-4) and results of cytometry experiments before treatment. To evaluate out-of-sample prediction, we computed the exact leave-one-out cross-validation (25).

resUlTs clinical and Mri response
In this study, fingolimod significantly reduced the mean ARR in 47.2 and 88.6% at 1 and 2 years, respectively, and the percentage of relapse-free patients was ~75% at the end of the study without differences between the subgroups. The percentage of patients with CDP was less than 10% at 1 and 2 years. Fingolimod reduced the percentage of patients with MRI activity at 2 years. A transitory increase in the ARR and the number of GdE lesions in the NTZ group was observed through the first year, but at 2 years these variables were decreased with no differences between the subgroups. There were no differences in the mean EDSS and number of new T2-w lesions at 1 and 2 years. At 2 years, the mean annual BVL was lower than at 1 year with no significant differences between the subgroups. At 1 year, the percentage of responder patients was 57.5% (n = 23) and 35% (n = 14) according to NEDA-3 and NEDA-4, respectively. At 2 years, it was 42.5% (n = 17) and 27.5% (n = 11) without differences between the subgroups. The results obtained at 1 and 2 years as well as all clinical and MRI measures are shown in Figures 1A-F and Table  S3 in Supplementary Material. The treatment was discontinued in six patients, one because of adverse events and five due to a lack of efficacy. All of these patients were classified as NRs at 2 years in an intention-to-treat analysis.

characterization of lymphocyte subpopulations in Ms Patients Before and after Fingolimod Treatment
In lymphocyte populations, no significant differences between HC and patients were observed before treatment, except for the percentage of LB1 cells composed mostly of CD11b+ cells, which was smaller in patients. At baseline, the only difference between subgroups was a higher percentage of transitional B cells in the NTZ-group compared with the no-NTZ-group (Tables S4 and S5 in Supplementary Material).
Before treatment, a positive correlation between IFN-γproducing CD4 cells and TEM CD4 cells was found, which suggested that IFN-γ production in CD4 lymphocytes was mainly performed by TEM cells. Interestingly, we found that plasmablasts were positively correlated to NK bright cells and negatively to NK dim cells in MS patients before treatment. After 6 months (e) The Wilcoxon signed-rank test was used to compare the mean of the annual brain volume loss (BVL) between 1 and 2 years (measured by SIENA method). (F) Percentage of patients with a good response to fingolimod according NEDA-3 (no evidence of disease activity: no relapses, no any MRI activity and no clinical disease progression) and NEDA-4 (like NEDA but including BVL < 0.4). p < 0.05 was considered statistically significant. The p-values were calculated using the Wilcoxon signed-rank test to compare differences before and after 6 months of therapy and using Mann-Whitney test to compare differences between subgroups. p < 0.001 was considered statistically significant after Bonferroni's correction for multiple tests. selected lymphocyte subpopulations could Predict the Treatment response to Fingolimod As shown in Figure 2, patients who achieved NEDA-3 and NEDA-4 status at 1 year had a significantly higher percentage of NK bright and plasmablast cells and a lower proportion of NK dim and IL-2-producing cells at baseline than NR patients (p < 0.001 for all). These differences were maintained after 6 months of  treatment for NK bright cells. These patients were significantly more resistant to decreases in percentages of NK bright and LB1 cells and showed a greater decline of CD8 naïve T and CD8+ CCR4+ CCR6+ cells than NR patients after 6 months of treatment (p < 0.001 for all). These differences were maintained at 2 years, but they did not achieve statistical significance. The chemokine receptors CCR4 and CCR6 have been identified as markers of human Th17 cells that secrete IL-17 (27). We also found that the frequency of Th17 cells measured by surface markers was higher than their frequency estimated by intracellular staining, but this phenomenon could not be correlated to any of the response parameters. Logistic regression analysis showed a greater number of subpopulations with high regression coefficients (positive or negative) in pretreatment samples than samples obtained at 6 months. In pretreatment sample, the subpopulations that can be used as biomarkers of response to fingolimod for most of the target variables (clinical activity, radiological activity, NEDA-3, and NEDA-4) at 1 and 2 years are NK bright, NK dim, plasmablasts, and IL-2 producing cells. On the contrary, in the sample obtained at 6 months, the regression coefficients only highlight the naive CD8 T population as a possible predictor of response as is shown in Figure S3 in Supplementary Material. According to the values of these four baseline candidates for prediction, each patient was placed in a 2D representation by nonparametric multidimensional scaling, and we found that responder patients at 1 year were located in a well-differentiated region of space, as shown in Figure 3. At 2 years, the differentiation was less evident.

alteration of the gene expression Profile by Treatment
A principal component analysis of gene expression profiles of all samples is shown in Figure 4A. Naïve patients and HC were transcriptomically different from patients in our study. No transcriptional differences were found between the NTZ and no-NTZ group before treatment, which allowed the analysis of all patients without taking into account the previous treatment. Fingolimod exerts powerful transcriptional effects on PBMCs of MS patients.  A total of 16,818 filtered probes were used for the differential analyses, resulting in the identification of 3,805 upregulated and 3,741 downregulated genes in response to fingolimod. As expected, genes related to mechanisms of action of fingolimod such as S1PR1, SELL (CD62L), and CCR7, as well as genes implicated in sphingolipid metabolism SPHK1 and SHPK2 were downregulated after treatment. As shown in Figure 4B, fingolimod produced an anti-inflammatory effect by downregulating CD40L, CD40, IRF4, CR2, IL-23A, CXCR5, CD24, CD2, CD27, and CD19, as well as an increase in the expression of the antiinflammatory cytokine genes and their receptors (IL-10, IL10RA, IL10RB, and IL-15) and the TNFRSF1A receptor gene, a regulator of inflammation. We also found an upregulation of antioxidant genes implicated in the control of reactive oxygen species (ROS), such as SOD2, CAT, MT1X, MGST1, and MAOA. All data for the expression levels, fold changes, and statistical significance in the DEG analysis are shown in published repository (26) and all data of RNA-seq are included in SRA database SRP132699 (28).

Differential Modulation of the gene expression Profiles in responder and nr Patients
After stratification based on drug response, we observed that the number of DEGs was higher in responder than NR patients at 1 and 2 years (~2,500 genes vs ~1,500 genes), indicating that the regulation of gene transcription is intimately linked to the treatment response. C3 and IFIT1 genes were upregulated, and genes such as FOXP3 and IL17A were downregulated in response to fingolimod only in responder patients. By contrast, EGR4 and MAPK11 genes were downregulated exclusively in NR patients ( Figure 4C).
In the S1P pathway, we found that fingolimod inhibited antiapoptotic processes through AKT and NOS3 downregulation and could play a role in decreasing stress-fiber formation and T cell maturation through NFKB downregulation. We did not observe significant differences between responder and NR patients in the S1P pathway except for an upregulation in responders and a downregulation in NRs of MAPK3 (at 1 and 2 years). Similarly, in the oxidative stress response pathway (Ontology terms: PW0000378), we observed a downregulation of FOS and an upregulation of HMOX1 in patients with evidence of disease activity (EDA) at 1 and 2 years ( Figure S4 in Supplementary Material).

Baseline signature in responder Patients
The gene expression profile before treatment was compared among responder and NR patients. Before treatment, 127 genes were differentially expressed when the response was measured at 1 year and 36 genes at 2 years (Figures 5A,B). In responder patients, overexpression of genes related to cytokine secretion (IL-21, IFNG, IL-17A, and CXCL9), activation of the T cell response (TNFSF8 and TNFRS9), activation of regulatory mechanisms (FOXP3 and SPP1), immunoglobulin secretion (GPI), apoptosis (AATF) and processing of class I MHC (PSMD9), as well as a down-expression of genes related to modulation of the innate immune response (LEP and B3GAT1), regulation of the Th1 response (UBASH3A and MAPK8IP1), immunoglobulin receptors (FCRL1, FCRL2, and JAK3), and B lymphocyte response (BCL3) before treatment was observed. In addition, 10 DEGs selected for validation by qPCR showed consistent expression between the qPCR and RNA-seq results.
When the gene expression profile was compared at 6 months between subgroups to determine the differential effect of fingolimod according to the clinical response, we found 6 DEGs in responder vs NR patients at 1 year and 51 DEGs at 2 years. The good responders showed a downregulation of specific markers of activated Tregs such as LRRC32, an upregulation of the cellular growth factor ADA-2, apoptotic genes such as SGK1 and BCL2L13 and genes related to the antigenic presentation by MHC-II such as CD74, HLA-DRB6, HLA-DRA, and HLA-DRB1. The most representative DEGs are shown in Figures 5C,D. A summary of the expression levels of the most representative genes and PCR validation are shown in Table S7 in Supplementary Material.

gsea of Degs
The most important genes related to immunity and MS that achieved significant enrichment (positive or negative) after treatment are shown in Figure 6. Fingolimod was observed to exert an important transcriptional effect on PBMCs of MS patients, approaching a profile more similar to the HC profile. The gene sets containing genes related to S1P signaling such as MAPK/ERK, Rho kinases, RAC1, Ras, and sphingolipid metabolism, as well as the genes participating in IL-10 production and innate immune activation such as NK cells, Fc-mediated phagocytosis, CCR3, CCR5, CXCR4, the complement cascade, the innate immune system, and apoptosis were positively enriched in response to fingolimod. In our study, fingolimod downregulated gene sets related to cytokine production, such as the IL12 pathway, IL-17 pathway, and cytokine pathway, and of gene sets implicated in Th1, Th2, and Th17 response, such as Th1/Th2 pathway and GATA3. After stratification based on response to the drug, we observed that antigen processing cross-presentation and oxidative phosphorylation pathways before treatment were positively enriched exclusively in responder patients.
The normalized enrichment scores, p-values, and FDRs of significantly enriched gene sets implicated in the immune response and MS are summarized in Table S8 in Supplementary Material.

Predictive Model of the Fingolimod response
From the baseline cytometry and RNA-seq variables that showed statistically significant differences between responder and NR patients, we selected two lymphocyte subpopulations (NK bright and plasmablasts) and three genes (FOXP3, FCRL1, and GPI).
With these five variables the model was built (Figure 7A) and the individual contribution of each variable to the fitted model was evaluated (Figure 7B). An in-sample posterior prediction check indicated that the model was correctly fitted to the provided data ( Figure 7C). The exact leave-one-out cross-validation showed that the model could be expected to perform well for out-ofsample predictions (Figure 7D).

Validation of Predictive Model
The differential expression of FOXP3, FCRL1, and GPI between responder and NR was validated by qPCR in the remaining 30 patients who were not used for the analysis of RNA-seq and the results are shown in the Table 2.
The model was tested in this cohort for the good response prediction with a sensitivity of 85%, a specificity of 95%, a positive predictive value of 85%, and negative predictive value of 95%. The mean squared error is 0.1111438. The response probability obtained for each sample and the comparison with real clinical response is shown in Figure 8. The mean, lower, and upper posterior predictive for each sample in the validation was summarized in Table S9 in Supplementary Material.
In a previous report, pretreatment differences were found between the responder and NR patients in the transitional B and RTE cells (34,35). Our study differed from those results probably due to differences in patient selection. In our study, naïve patients were excluded and the percentage of patients whose previous treatment was NTZ was 45% while the previous report included 37.5% of naive patients and only 20.2% of patients whose previous treatment was NTZ (n = 7). In our study, the percentage of transitional B cells in the pretreatment samples was higher in NTZ group. That is consistent with the increase in transitional B cells and RTEs induced by NTZ that can remain even months after the withdrawal (36)(37)(38). For the analysis of RTE and transitional B cells as cytometry biomarkers, we think that the stratification of the patients according to the previous treatment is necessary and patients who have received NTZ during the last months should be excluded. In addition, although we did not analyze the RTE cells by cytometry, we did not find differences in the expression levels of the genes coding for their surface markers CCR7, CD45RA, CD31, and PTK7 (CCR7, PTPRC, PECAM1, and PTK7) between responder and NR patients before treatment.
In this study, as in most studies of this type (41)(42)(43)(44), the PBMCs are considered as a valid cellular group since they are able to reflect patterns of expression characteristic of certain diseases and treatment effects. The majority of studies evaluating the transcriptomic effect of a treatment do it on whole PBMCs. In fact, the changes in expression levels on PBMCs in patients with MS Color key annotation indicating the source of each sample (healthy controls, naïve patients, patients at baseline, and patients at 6 months of fingolimod treatment) and NEDA-4 status at 1 and 2 years (responder: NEDA-4, non-responder: EDA-4) for the 10 patients (4, 6, 8, 18, 25, 10, 28, 36, 37, and 39). *Gene sets enriched before treatment only in responder patients at 1 and 2 years. treated with IFN-β are known as IFN signature and it has been proposed as a response marker (44). In addition, a recent study showed that the transcriptional change induced by fingolimod belong almost entirely to CD4 lymphocytes and based on these results, the efficiency of the sorting could be lower than the analysis of the entire population of PBMCs (45). Although in our study a discrimination by lymphocyte subtypes was not conducted, we obtained the same results as the previously described genes when the effect of fingolimod was analyzed.
It is clear that fingolimod exerts a potent anti-inflammatory effect by decreasing the percentage of cytokine-producer T and B cells (7). In this study, we confirmed this effect not only at the cellular level but also at the molecular level, since fingolimod induced a downregulation of genes such as CD40L, CD40, IRF4, CR2, IL-23A, CD2, IL17A, and IL17D. Similarly, the cytokine activity pathways were downregulated after treatment.
Natural killer cells are lymphocytes of the innate immune system that are involved in defense against viruses, bacteria, or  parasites, or malignant transformation. NK bright and invariant NKT cells, a subset of NKT, are often related to regulatory processes (46). In MS patients, a high proportion of NK bright cells has recently been associated with stable MRI (47), and we found that a higher percentage of these cells before treatment was associated with a good response and that they decreased to a lesser degree after treatment in responder compared with NR patients. In our study, fingolimod not only increased the percentage of NK cells but it also upregulated NCAM1 (CD56) and FCGR3A (CD16a) at the transcriptional level; an enrichment of the NK and innate response pathways was also observed as a result of the therapy. The role of innate immunity in MS has recently been appreciated: its cells are capable of producing direct damage to CNS myelin and play an important role as a source of ROS contributing to axonal damage (48). Despite this phenomenon, in our study, fingolimod downregulated genes implicated in oxidative stress and inhibition of stress-fiber formation through the RhoA-ROCK1-NFKB pathway, suggesting that the antioxidant effect of fingolimod was not dependent on innate immune activation. In vitro, fingolimod exerts an antioxidant effect, and in the EAE model, it blocks astrocyte activation and nitric oxide production (49). In MS patients, lower serum levels of ROS were observed in those treated with fingolimod compared with the first line-treated patients, and a correlation between ROS levels and disease duration was also observed (50).
The percentage of CD4 T regulatory cells do not differ between patients and controls, but their function is compromised (maturation and migration) (51) and the transcription factor FoxP3, crucial for the function of regulatory T-cells, has been reported to be diminished in RRMS patients (52). According to previous reports (7), we found that fingolimod induces an increase in the percentage of regulatory T lymphocytes but surprisingly a downregulation of expression levels of FOXP3 was observed. This divergence had already been previously described and commented on in a study by van Pesch et al. (53) and it could be partially explained due to the reduction of the absolute number of Tregs. FoxP3 is not only a transcription factor of Tregs but it is also expressed transiently in effector T cells after its activation (54). The decrease in lymphocyte activation in response to fingolimod could contribute to the decrease in FOXP3 expression levels. On the other hand, the increase in the expression of IL-10, IL10RA, IL10RB, IL-15, and TNFRSF1A reveals the powerful inductive effect of regulatory mechanisms in response to fingolimod observed in this study.
In the present study, we found that higher levels of FOXP3 before treatment were correlated with a good response, raising the possibility that patients with stronger regulatory mechanisms can obtain a greater benefit from immunomodulatory therapies.
FOXP3 and LRRC32, a recently identified specific marker of activated regulatory T-cells (55), were downregulated by fingolimod only in responder patient as well as the ADA2 gene, also known as CECR1, which intervenes in the activation of Tregs through binding to the CD39 receptor (56). That suggested that once the excessive inflammatory response is controlled the compensatory mechanisms are less necessary.
Upregulation of the apoptotic pathway indicated that fingolimodinduced programmed cell death, a mechanism that might contribute to lymphopenia. However, we observed an overexpression of ADA2, SGK1, and BCL2L13 after 6 months of treatment only in responder patients, potentially suggesting a differential effect of fingolimod on the mechanisms of cell proliferation and differentiation.
B cells and B cell activation are recognized to play a key role in MS pathophysiology (57). The tumor necrosis factor receptor superfamily member 13B (TNFSF13) and its ligand B-cell activating factor have also been shown to play an important role in the proliferation and differentiation of B cells in MS. In fact, oligoclonal bands have been associated with increased levels of TNFSF13B in the CSF of MS patients (58). The CXCL13 and CXCR5 molecules have also been related to the recruitment of B cells to the CNS during neuroinflammation (59). We observed an important downregulation of IRF4 (a highly expressed molecule in activated B cells) (60), CD27 (a memory B cell marker) and CXCR5, as well as a decrease in the percentages of memory B cells after treatment. Interestingly, we observed a higher percentage of plasmablasts and an overexpression of CXCL13 before treatment in patients who achieved a NEDA response, suggesting that the modulation of B cell activation plays an important role in the clinical response to fingolimod. The plasmablasts are cells highly differentiated toward the production of antibodies and are related to clinical and radiological activity in EM (61). In this study, there was a higher percentage of plasmablasts at baseline in responder patients to fingolimod, as well as higher transcriptional activity in terms of antibody synthesis. An example of this observation is the basal overexpression of the GPI gene that codes for glucose phosphate isomerase (GPI), a lymphokine that induces the production of immunoglobulins (62). In rheumatoid arthritis, GPI levels correlate with disease activity and it has been shown that fingolimod can prevent GPI-induced outbreaks (63). An opposite behavior was observed in this study in the genes that code for the immunoglobulin receptors FCRL1 and FCRL2 indicating that the responder patients despite exhibiting an increase in the production of immunoglobulins have a lower activity of their receptors. On the other hand, low levels of FCRL1 and 2 have been associated with an increase in radiological activity in MS (64) and this is in agreement with the fact that fingolimod exerts a better effect in patients with a more inflammatory component.
The interleukin-2 receptor α-subunit (IL-2RA) has generated great interest since polymorphisms of this gene have been associated with the risk of developing MS (65) and due to the recent approval of daclizumab, a humanized IgG1 monoclonal antibody directed against the alpha subunit of the high-affinity IL-2 receptor, which has demonstrated high efficacy in the control of MS activity (66). Fingolimod treatment decreases soluble IL-2RA plasma levels (67). In accordance with those results, we detected a decrease in IL-2-producing cells and a downregulation of the IL-2RA gene in PBMCs of treated MS patients. Although we observed a differential decrease in IL-2-producing cells between NEDA and EDA patients, no differences were found at the transcriptional level.
Surprisingly, we found no differences in the transcriptomic profile of PBMCs from patients of NTZ-group and no-NTZ group although previous reports of transcriptional changes induced by NTZ and IFN-β in MS patients have been published (68,69). In the NTZ-group, this could be explained by the absence of long-lasting immunological and transcriptional changes by NTZ after the washout period and we could think that the immune reactivation is related just to the reversibility of the effect of NTZ on the lymphocyte transmigration and not to activation changes. In the no-NTZ group, the absence of an interferon/glatiramer acetate signature could be related to the lack of effectiveness of first-line DMTs, which was the main reason for the switch to NTZ.
In addition to the differential effect of fingolimod between patients, the baseline cellular and genetic characteristics can help us to predict the clinical response to the drug. High percentages of NK bright cells and plasmablasts, high expression levels of FOXP3 and GPI, as well as low expression levels of FCRL1 and low EDSS before treatment could be correlated with a good response, which allowed us to design a prediction model of the fingolimod response that must be validated in an independent cohort. cOnclUsiOn Fingolimod treatment affects practically all lymphocyte subpopulations and exerts a strong effect on genetic transcription. In PBMCs from MS patients, fingolimod induces a decrease in the inflammatory response at cellular and molecular levels, improves regulatory mechanisms, and has important antioxidant effects. The differential fingolimod-induced effects on PBMCs at 6 months as well as the baseline cellular and molecular characteristics can be used to predict the clinical response at 1 and 2 years. We postulate that a combination of cellular (NK bright and plasmablasts), molecular (FOXP3, GPI, and FCRL1), and clinical markers (EDSS and gender) are possible response biomarkers. However, our data are only descriptive, and more studies with a greater number of patients are necessary to validate this model.

eThics sTaTeMenT
This study was carried out in accordance with the recommendations of Spanish Real Decreto 223/2004. The protocol was approved by the Clinical Research Ethics Committee of the Puerta de Hierro Hospital (CEIC). All subjects gave written informed consent in accordance with the Declaration of Helsinki.