Informed consent was obtained from all participants. The study was approved by the Vall d’Hebron Hospital, Barcelona [EPA(AG)57/2013 (3834)] and the Hospital Ramón y Cajal, Madrid.

The effect of cladribine on gene, miRNA and protein expression was explored in vitro in PBMCs obtained from ten RRMS patients. Seven patients were treatment naïve and three patients received in the past teriflunomide, glatiramer acetate and interferon-beta but were untreated for the last six months before blood collection. Table 1 summarizes the main demographic and baseline clinical characteristics of patients included in the study.

Human PBMCs were isolated using standard Ficoll gradient centrifugation. Cells were cultured in complete medium made with RPMI 1640 growth medium (Life Technologies, CA, USA) supplemented with 1% penicillin, 1% streptomycin, 1% glutamine, and 5% fetal calf serum (FCS; Gibco, Carlsbad, CA. USA). Cladribine (Merck, Darmstadt, Germany) and deoxycytidine (Sigma-Aldrich, St. Louis, MO, USA) were dissolved in DMSO and water respectively and further diluted to the appropriate concentration in complete medium

Isolated PBMCs were cultured in complete medium at 5% CO₂, 37°C, and >95% humidity on 24-well plates, at a concentration of 1x10⁶ cells/ml. Cells were stimulated with 10 ng/ml phorbol 12-myristate 13-acetate (PMA) plus 500 ng/ml ionomycin (Sigma-Aldrich) and simultaneously treated with 0.1 μM cladribine in the presence or absence of 250 μM deoxycytidine (added 30 min prior to the addition of cladribine) for 4h. After in vitro treatment with cladribine, 2.5x10⁶ PBMCs were harvested in separated tubes (for RNA and protein extraction), pelleted (10 min, 400 × g) and frozen as dry pellets on dry ice, and stored at −80°C until all samples were collected. In this context, the following experimental groups were considered in the study: no cladribine, cladribine, cladribine plus deoxycytidine, and only deoxycytidine.

These in vitro culture conditions were based on previous studies performed by our group which demonstrated that untreated PBMCs needed to be stimulated first in order to induce activation and cell proliferation (8). The activation status after PBMC culture was tested by flow cytometry (by means of CD69 expression), and only experiments in which PBMCs showed a percentage of activation higher than 90% compared with the unstimulated condition were used to obtain RNA and protein for further purposes (Figure S1). The concentration of cladribine used in the in vitro experiments (0.1 μM) corresponds to the mean concentration present in plasma of patients with MS under cladribine treatment (15) The use of deoxycytidine, a molecule that competes with cladribine for the phosphorylation via DCK, will allow to investigate whether the inhibition of cladribine activation by DCK produces changes in expression patterns.

Total RNA was extracted from PBMC pellets with the miRNeasy Mini Kit Qiagen (Qiagen, Basel, Switzerland) following the manufacturer’s protocol. RNA quantity and quality were determined with a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Gene expression was assessed using the Clariom™ S array (Affymetrix, Santa Clara, CA, USA), and miRNA expression was determined with the GeneChip^® miRNA 4.0 Array (Affymetrix). Data processing and normalization were carried out using the R statistical language (https://www.r-project.org/) (version 3.6.3) and the libraries created for microarray analysis in the Bioconductor Project (16). In order to get normalized expression values (log2 scale), raw expression values obtained directly from.CEL files were processed using the Robust Multiarray Average (RMA) algorithm in oligo package (17). Linear models for microarray data (LIMMA) R package was used to identify differentially expressed genes and miRNAs between the different groups (18). Genes and miRNAs showing a differential expression by two-sample t-test (p-value <0.05), were considered significant.

PBMC pellets were lysed in RIPA buffer (1% sodium deoxycholate, 1% Triton X-100, 0.15 M NaCl, 0.1% SDS, 0.01 M sodium phosphate, pH 7.2) supplemented with 1 mM Na₃VO₄, 0.5 mM DTT, and protease inhibitor mixture. Protein concentration was determined using Pierce BCA Protein Assay Kit (Thermo-Fisher Scientific, San Jose, CA, USA) as recommended by the manufacturer. Then, samples (10ug) were precipitated with cold acetone and the protein pellets were dissolved in 6M Urea/200mM ammonium bicarbonate. Samples were reduced and alkylated (dithiothreitol and iodoacetamide) prior to tryptic digestion (LysC and Trypsin) and the resulting peptide mixtures were acidified with formic acid and desalted with a MicroSpin C18 column (The Nest Group, Inc; Southborough, MA).

Samples were analyzed by mass spectrometry (Orbitrap Fusion Lumos) coupled to an EASY-nLC 1200 (Thermo-Fisher Scientific). Peptides were separated by reversed-phase chromatography using C18 columns of 50 cm length and 75 μm inner diameter in chromatographic gradients of 90 minutes (water, acetonitrile, 1% formic acid) and a flow rate of 300 nl/min. The mass spectrometer was operated in data-dependent acquisition mode, MS1 scans were acquired in a mass range of 350-1,400 m/z and detected in the Orbitrap at 120.000 resolution. The most intense precursor ions selected by the instrument (“Top Speed” algorithm and 60 seconds dynamic exclusion) were then fragmented by high-energy collision dissociation (HCD) at 28% normalized collision energy and detected in the ion trap. To control the instrument performance (19) and avoid for carryover, digested bovine serum albumin was run between each sample.

Mass spectrometry data was analyzed by Proteome Discoverer v2.0 (Thermo-Fisher Scientific) and Mascot search engine v2.6 (Matrix Science) (20) using the Swiss-Prot human database (February 2020) plus a list of common contaminants and all the corresponding decoy entries (21). Mass tolerances were set to 7 ppm and 0.5 Da for MS1 and MS2 respectively. Tryptic peptides with a maximum of three missed cleavages were permitted. Carbamidomethylation on cysteines was set as a fixed modification and the variable modifications of N-terminal protein acetylation and methionine oxidation were also considered. For peptide identification, we set a maximum of 5% false discovery rate.

Protein Abundances were estimated as the average of the three most abundant peptides and they were normalized by equalizing the median protein abundance in all samples. The raw proteomics data can be found in the repository PRIDE (22) with the identifier PXD040572.

Data obtained from arrays of genes, miRNAs, and mass spectrometry were processed as described above to get the normalized, log-transformed expression values, and then analyzed using two different integrative bioinformatics approaches:

(i) DIABLO (Data Integration Analysis for Biomarker discovery using Latent variable approaches for Omics studies): This algorithm selected the most predictive or discriminative molecules in the datasets that helped to classify the samples (23). Before running this analysis, three types of parameters were tuned. (i) Matrix design: a minimal link of 0.1 between the datasets was used to maximize the discriminative power of the analysis. (ii) Number of components: was set on 2. (iii) Number of features to select on each dataset/component: based on a sPLS-DA (sparse Partial Least Squares-Discriminant Analysis) method, the optimal number of features to select on gene, protein and miRNAs was 40, 20, and 40, respectively. This approach will allow to have an overall multiomics picture of the most relevant molecular patterns associated with the different phenotypes.

(ii) Multi-stage: the individual omics databases are first analyzed separately and the results are then merged and compared (18). For the differential expression analysis, genes and proteins with a p-value <0.5, and miRNAs with a p-value <0.5 and absolute log2Fold-Change > 0.5 in the comparison between cladribine treated and the untreated conditions were selected for further analysis. This approach will allow to focus on individual candidate biomarkers functionally connected in the three datasets that are associated with the effect of cladribine.

The multiMiR Bioconductor’s package was employed to identify validated targets for the miRNAs differentially expressed between the treated and untreated conditions. Such analysis enables retrieval of validated miRNA-gene target interactions from several external databases (miRTarBase, miRecords, and TarBase databases) (24).

Before integration of the three datasets (gene/protein/miRNA), a filter was applied in order to eliminate molecules that either did not differentiate from the rest in the same dataset or were not detected in most of the samples included in the study. In the gene expression dataset, a filtering step was performed based on a variance threshold of 65%, and transcripts with a low variability or without a valid Entrez ID annotation were removed from the analysis. For the miRNA dataset, probes belonging to other species (non-human), control probesets, and duplicates were removed. For the proteomic dataset, proteins with a confidence score < 20, as well as proteins detected in < 25% of the samples in each group were discarded (23).

Gene and miRNA expression levels were determined ex vivo by qPCR in PBMCs from five RRMS patients at baseline and after 3 and 12 months of cladribine treatment onset. Three patients were treatment naïve, and 2 patients were treated before with dimethyl fumarate and one with glatiramer acetate. A summary of demographics and clinical characteristics of the ex vivo cohort is given in Table 1. Total RNA was extracted from PBMCs as previously described for the microarray studies, following a reverse transcription into complementary DNA (cDNA) using the High capacity cDNA Archive kit (Applied BiosystemsTM, Foster City, CA, USA) for gene expression, and the TaqManTM Advanced miRNA cDNA synthesis kit (Applied Biosystems) for miRNA expression experiments. Gene expression levels of NIBAN2 (niban apoptosis regulator 2), NHLRC2 (NHL repeat containing 2), PPIF (Cyclophilin D), and JUN (Transcription factor Jun) transcripts were determined with TaqMan gene expression assays (Hs00382858_m1, Hs00612275_m1, Hs00194847_m1, and Hs01103582_s1, respectively) and TaqMan Fast Advanced Master Mix (Applied Biosystems). The obtained values were normalized according to the level of expression of the housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The miRNAs expression levels of miR-30e-5p, miR-766-3p, miR-21-5p, miR-30b-5p, and miR-484 were determined with TaqMan Advanced miRNA assays (479235_mir, 478342_mir, 477975_mir, 478007_mir, 478308_mir, respectively) and the TaqMan Fast Advanced Master Mix (Applied Biosystems), according the manufacturer. Assay for miR-191-5p (477952_mir) was used as a reference miRNA for normalization purposes. Both assays (gene and miRNA) were run on the ABI PRISM^® 7900HT system (Applied Biosystems) Data was analyzed with the 2-ΔCT method (25). Results were expressed as fold change in gene and miRNA expression at 3 and 12 months of cladribine treatment compared with the baseline (calibrators).

Statistical analysis was conducted using SPSS (Version 20) software (SPSS Inc, Chi-cago, IL, USA) for MS-Windows and GraphPad Prism 6.0 software (GraphPad Prism Inc, La Jolla, CA, USA). Expression values are represented using boxplots. Adjusted p-values for omics data were obtained using the Benjamini & Hochberg method. One-way ANOVA with repeated measures following a Dunnett multiple comparisons test was used, taking the control condition (baseline) as reference.

The data presented in this study are available on request from the corresponding author.

A schematic overview summarizing the main steps of the study design and analysis is represented in Figure 1.

To reveal the molecular expression pattern associated with cladribine treatment, we applied the DIABLO bioinformatics approach which identified a multiomics molecular signature, by finding co-expression between the molecules from the different datasets measured in the same samples that best explained the study phenotypes: no cladribine, cladribine, and cladribine plus deoxycytidine. The only deoxycytidine condition was not included in the analysis, as this group showed a similar expression pattern to the no cladribine condition (data not shown). Figure 2A shows the molecules in the three databases with the highest differences between the groups. Figure 2B shows the unsupervised hierarchical clustering analysis performed with the features presenting the highest differences at the gene, protein, and miRNA levels. Overall, cladribine induced an expression signature characterized by general downregulation of the three different types of molecules, most noticeable at the gene expression level compared with the untreated condition. Interestingly, the cladribine plus deoxycytidine condition was characterized by an intermediate expression signature between the cladribine treated and untreated conditions at the three molecular levels, genes, proteins and miRNA. The lists of features included in the clustered heatmap and loading plots are shown in Supplementary Table S1 (for genes), Supplementary Table S2 (for proteins), and Supplementary Table S3 (for miRNAs).

In order to identify individual biomarkers candidates associated with cladribine treatment, we applied a multi-stage approach to analyze the common genes and proteins differentially expressed in the comparison between the cladribine treated and untreated conditions. The gene expression microarray and mass spectrometry analyses in PBMCs allowed to identified, after applying the corresponding filtering steps, 6,484 genes and 4,056 proteins respectively, of which 1,767 molecules were found in both studies (Figure 3A). Of these, 265 genes and 208 proteins were differentially expressed (p-values < 0.05) between the treated and untreated conditions. In line with the previous integrative analysis (DIABLO), cladribine induced a major downregulation of gene expression with 209 molecules downregulated vs. 56 upregulated (Figure 3A). Cladribine also reduced protein expression with 116 molecules downregulated and 92 upregulated after treatment (Figure 3A). Of note, four of these molecules, NIBAN2, NHLRC2, PPIF, and JUN, were common in both gene and protein datasets. As shown in Figure 4, gene and protein expression levels of NIBAN2, PPIF and JUN were significantly downregulated in PBMCs by cladribine in vitro. In contrast, gene and protein expression levels for NHLRC2 were significantly upregulated by the effect of cladribine (Figure 4). These four candidates were selected as potential biomarkers for ex vivo validation.

To select individual miRNAs involved in the regulation of the genes differentially expressed associated with cladribine exposure, we quantified the miRNAs using microarrays and further performed the differential expression analysis. The final miRNA dataset consisted of 6,350 probesets that included mature miRNAs, miRNA precursor stem-loops, snoRNA, HAcaBox, CDBox, scaRNA, and probesets targeting 5.8s rRNA. Fifty-five out of the 6,350 miRNAs were differentially expressed (p-values < 0.05 and absolute Log2Fold-Change > 0.5) between the cladribine treated and untreated conditions (Figure 3B). Compared to the untreated condition, 37 miRNAs were downregulated and 18 miRNAs upregulated after in vitro treatment with cladribine (Figure 3B). Next, analysis of the validated miRNA-target gene interactions led to the identification of 13,751 validated targets for the 55 dysregulated miRNAs (Figure 3C). Twenty-two out of these 55 miRNAs were known to regulate the four molecules that were differentially expressed both at the gene and protein levels between the cladribine treated and untreated conditions (Table 2). Five of the most dysregulated miRNAs were selected for ex vivo validation based on their biological and functional relevance to MS, as reported in the literature (Table 2).

As a final step, we approached the validation of the in vitro findings in PBMCs collected from MS patients treated with cladribine. For this ex vivo validation, the expression levels of the four selected molecules (NIBAN2, NHLRC2, PPIF, JUN) and the five selected miRNAs targeting these molecules (miR-484, miR-30b-5p, miR-30e-5p, miR-21-5p, miR-766-3p) were determined by qPCR at baseline and after 3 and 12 months of treatment. As depicted in Figure 5A, gene expression levels of NHLRC2 and PPIF increased and decreased respectively over time after cladribine treatment and differences reached statistical significance at 12 months compared to the baseline condition. In contrast, no significant differences were observed for NIBAN2 and JUN between the baseline and the different treated time points (Supplementary Figure S2). Interestingly, miR-30b-5p, a miRNA known to target NHLRC2, was significantly downregulated by the effect of treatment after 3 and 12 months (Figure 5B). Cladribine also significantly downregulated miR-30e-5p and miR-21-5p expression at 3 months, although expression levels for these two miRNAs were later increased at 12 months by the effect of treatment (Figure 5B; Supplementary Figure S3). No significant differences associated with cladribine treatment were observed for miR-766-3p and miR-484 (Supplementary Figure S2). To notice, expression levels for these selected candidates did not differ following stratification of patients into naïve and previously treated patients (data not shown).