Identification of CD4+ T cell biomarkers for predicting the response of patients with relapsing‑remitting multiple sclerosis to natalizumab treatment.

Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system of autoimmune etiopathogenesis, and is characterized by various neurological symptoms. Glatiramer acetate and interferon‑β are administered as first‑line treatments for this disease. In non‑responsive patients, several second‑line therapies are available, including natalizumab; however, a percentage of MS patients does not respond, or respond partially. Therefore, it is of the utmost importance to develop a diagnostic test for the prediction of drug response in patients suffering from complex diseases, such as MS, where several therapeutic options are already available. By a machine learning approach, the UnCorrelated Shrunken Centroid algorithm was applied to identify a subset of genes of CD4+ T cells that may predict the pharmacological response of relapsing‑remitting MS patients to natalizumab, before the initiation of therapy. The results from the present study may provide a basis for the design of personalized therapeutic strategies for patients with MS.

Introduction

Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system of autoimmune etiopathogenesis, and is characterized by neurological symptoms (1). A total of ~90% of patients with MS are diagnosed with relapsing-remitting disease, that involve acute periods of neurological dysfunctions followed by a period of recovery (2). Glatiramer acetate and interferon-β (IFN-β) are administered as first-line therapies. In patients non-responsive to treatment, several second-line therapies are available, including natalizumab and fingolimod (3). Natalizumab is a monoclonal antibody against the α4 subunit (CD49d) of α4 integrins [α4β1 (VLA-4) and α4β7), that prevents the transmigration of leukocytes across the blood-brain barrier by inhibiting interactions between α4β1 integrin/vascular cell adhesion molecule-1 and mucosal vascular addressin cell adhesion molecule-1 (4). In 2012, a longitudinal study assessing the effects of natalizumab on 333 patients with MS revealed that 69–88% of patients exhibited a positive outcome in all Patient-Reported Outcomes measures assessed (5). In the Natalizumab Safety and Efficacy in Relapsing Remitting Multiple Sclerosis (AFFIRM) trial (ClinicalTrials.gov Identifier: NCT00027300), of 942 patients, 627 were randomly selected for treatment with natalizumab and 315 were administered a placebo. The results revealed that natalizumab reduced the risk of sustained disability progression by 42% in a 2-year time-frame and decreased the rate of clinical relapse in 1 year by 68%, leading to an 83% reduction in the accumulation of new or enlarging T-2 hyperintense lesions (6). Additionally, the results of a SENTINEL trial indicated that 67% of patients receiving natalizumab plus IFN-β-1a remained free of new or enlarging T2-lesions compared with 30% of patients receiving IFN-β-1a alone (7). These findings indicate that despite the high efficacy, a percentage of patients with MS do not respond, or respond partially to natalizumab. Therefore, it is of the utmost importance to develop a diagnostic test to predict drug response in patients suffering from complex diseases, such as MS, in which several therapeutic options are readily available. This could lead to a double-fold advantage: Patients would benefit by avoiding ineffective therapies and healthcare costs would be notably reduced. In the present study, a machine learning approach was utilized to identify a subset of genes that may predict the response of patients with MS to natalizumab prior to the initiation of therapy.

Materials and methods

Molecular patterns of pharmacological resistance to natalizumab

For the identification of the molecular patterns underlying the pharmacological resistance to natalizumab in MS, we selected the GSE44964 microarray dataset, available from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). The dataset comprised whole-genome expression data from CD4+ T cells isolated from patients with MS and stimulated in vitro with precoated anti-CD3/-CD28 monoclonal antibodies for 48 h. The Agilent Sureprint G3 Human Gene Expression 8×60k platform was used to generate the dataset (Agilent Technologies, Inc., Santa Clara, CA, USA). The raw data were quantile normalized and batch effect-corrected using ComBat v2 (8). The patients were of Swedish origin and had relapsing-remitting disease. The GSE44964 dataset comprises data generated from two different microarray platforms. To avoid obtaining biased results, the data of two platforms were not combined; thus, analysis was conducted using the largest set of samples only. Patients were diagnosed with MS according to the McDonald criteria (9) and prospectively classified as low responders (LRs, n=6), if at least one period of relapse occurred during the follow-up period (3 years) and as high responders (HRs, n=6), providing no relapse was observed. Other parameters, such as magnetic resonance imaging could be used for the classification of LRs and HRs; however, the relapse rate is considered as a primary endpoint of several phase 2/3 clinical trials (6), and classifying patients as responsive and non-responsive on the basis of whether relapse had occurred or not is appropriate for a preliminary transcriptomic analysis. All samples were collected for gene expression analysis prior to the initiation of natalizumab treatment. The LR and HR groups were matched for sex, age, Expanded Disability Status Scale score (10) and disease duration (11). A total of 5/6 patients in the groups were males; the age of patients was 36±6.3 and 33.7±7.1 years old for LRs and HRs, respectively. Statistical differences between HR and LR patients were assessed using LIMMA version 3.26.8 (Linear models for microarray data) in R version 3.2.3 (12). P<0.01 was considered to indicate a statistically significant difference. Statistical analysis and principal component analysis (PCA) were performed using MultiExperiment Viewer software (http://mev.tm4.org/). Gene Ontology (GO) analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 web-based tool (13,14). Functional annotation provided by DAVID comprises >40 annotation categories, including GO terms, protein-protein interactions, protein domains, disease associations, pathways, homology, gene function, gene tissue expression and literature. Network analysis was performed using the GeneMania utility (15).

Identification of biomarkers for natalizumab responsiveness

In order to identify a specific gene expression signature for predicting patient responsiveness to natalizumab treatment, the UnCorrelated Shrunken Centroid (UCSC; http://home.cc.umanitoba.ca/~psgendb/birchhomedir/BIRCHDEV/doc/MeV/manual/usc.html) algorithm was used. UCSC analysis was performed with the probes that were determined to be significantly modulated in HRs compared with the LR group. Cross-validation was conducted using the following parameters: 5-fold and 10-fold cross-validation; each cross-validation run was divided five-fold and therefore, a total 10 cross-validation runs were performed. Δ-(shrinkage threshold) and ρ-(correlation threshold) values were empirically selected so that the smallest number of classification errors were obtained using the fewest genes. Subsequently, PCA and Hierarchical Clustering (HCL) was performed using only the set of the identified predictors. For HCL, Euclidean distance and average linkage degree were used.

Results

HR and LR patients have different transcriptomic patterns

Statistical analysis of the transcriptomic differences between CD4+ T cells from patients of the HR and LR groups revealed 45 significant probes (Table I). PCA (Fig. 1A) produced two main clusters that contained HRs and LRs, respectively. The results indicated that the mRNA expression levels of several genes notably differed between natalizumab-responsive and non-responsive patients, and that a distinct pattern of gene expression could be associated with natalizumab resistance. Functional annotation revealed that the most enriched categories and their associated genes were: ‘Lipid-binding’ [oxysterol binding protein like 6 (OSBPL6), fatty acid binding protein 3 (FABP3), estrogen receptor 1 (ESR1) and sorting nexin 10 (SNX10)], ‘estrogen-responsive protein Efp controls cell cycle and breast tumors growth’ [ESR1 and stratifin (SFN)], ‘cytoplasm’ [OSBPL6, ESR1, SFN, amyloid β precursor-like protein 1 (APLP1), inorganic pyrophosphatase PPA1), tropomyosin 3 (TPM3), 20S proteasome subunit β-2 (PSMB7), RAB28, member RAS oncogene family (RAB28), regulator of G-protein signaling 5, FABP3, TBC1 domain family member 32 (TBC1D32), SNX10 and coiled-coil domain-containing protein 8] and ‘protein localization to cilium’ (SNX10 and TBC1D32) (Fig. 1B). A regulatory network comprising the significant genes and the top 20 related genes, was presented as Fig. 1C. The computational gene network prediction tool GeneMania identified the DNA topoisomerase II α gene to interact with PPA1, TNFAIP3 interacting protein 3, PSMB3, RAB28, APLP1, ESR1 and ring finger protein 113A. Other nodes were represented by PSMB7, alanyl-tRNA synthetase, APLP1 and ESR1 (Fig. 1C).

Table I.

List of genes significantly modulated between LRs and HRs CD4+ T cells.

Probe	Gene accession no.	Gene symbol	Gene name	P-value	Log fold change
A_19_P00317412		XLOC_000787		0.0048	−0.4088
A_19_P00317731		XLOC_002473		0.0078	0.6845
A_19_P00319311				0.0046	0.6016
A_19_P00321466				0.0079	0.5516
A_19_P00805263		XLOC_001851		0.0088	0.3698
A_19_P00807336				0.0030	−0.5004
A_23_P108823	NM_032523	OSBPL6	Oxysterol binding protein-like 6	0.0054	−0.3303
A_23_P112187	NM_032843	FIBCD1	Fibrinogen C domain containing 1	0.0059	−1.1805
A_23_P164958	NM_032040	CCDC8	Coiled-coil domain containing 8	0.0044	0.6849
A_23_P27983	NM_005166	APLP1	Amyloid β (A4) precursor-like protein 1	0.0077	−0.7417
A_23_P381017	NM_152559	WBSCR27	Williams Beuren syndrome chromosome region 27	0.0051	0.8273
A_23_P386478	NM_024873	TNIP3	TNFAIP3 interacting protein 3	0.0065	−0.9948
A_23_P398172	NM_020819	FAM135A	Family with sequence similarity 135, member A	0.0099	0.7579
A_23_P401547	NM_015480	PVRL3	Poliovirus receptor-related 3	0.0076	0.4645
A_23_P54205	NM_017926	C14orf118	Chromosome 14 open reading frame 118	0.0067	0.6402
A_23_P77135	NM_080650	ATPBD4	ATP binding domain 4	0.0084	1.0763
A_23_P90523	NM_024578	OCEL1	Occludin/ELL domain containing 1	0.0025	−0.3799
A_24_P62783	NM_004102	FABP3	Fatty acid binding protein 3, muscle and heart (mammary-derived growth inhibitor)	0.0016	−0.6532
A_24_P920048	AK092807	LOC100127972	Uncharacterized LOC100127972	0.0043	0.9088
A_32_P204376	NM_001012421	ANKRD20A2	Ankyrin repeat domain 20 family, member A2	0.0084	0.3258
A_32_P7204	NM_004249	RAB28	RAB28, member RAS oncogene family	0.0057	0.8683
A_32_P797019	XM_003403482	NPEPL1	Aminopeptidase-like 1	0.0060	0.6953
A_32_P80245	NM_001109809	ZFP57	Zinc finger protein 57 homolog (mouse)	0.0042	1.0614
A_33_P3221925	NM_152730	C6orf170	Chromosome 6 open reading frame 170	0.0018	0.6414
A_33_P3227269	NR_024256	FLJ45983	Uncharacterized LOC399717	0.0047	0.9719
A_33_P3231297	NM_003851	CREG1	Cellular repressor of E1A-stimulated genes 1	0.0073	−0.3014
A_33_P3233273	NM_001142928	LRRC61	Leucine rich repeat containing 61	0.0067	0.4877
A_33_P3243093	NM_003617	RGS5	Regulator of G-protein signaling 5	0.0042	−0.7003
A_33_P3271241	NM_021129	PPA1	Pyrophosphatase (inorganic) 1	0.0079	−0.3602
A_33_P3273552	NM_002282	KRT83	Keratin 83	0.0004	−0.9603
A_33_P3289536	NM_001199835	SNX10	Sorting nexin 10	0.0082	0.9177
A_33_P3292724	BC017576			0.0042	0.6890
A_33_P3333600		LOC400950	Uncharacterized LOC400950	0.0078	−0.6855
A_33_P3338674	NR_033298	CCDC163P	Coiled-coil domain containing 163, pseudogene	0.0006	−1.1794
A_33_P3341474	NM_001080412	ZBTB38	Zinc finger and BTB domain containing 38	0.0067	0.7430
A_33_P3346193	NM_001043351	TPM3	Tropomyosin 3	0.0022	0.8160
A_33_P3354514	BC047507	SLC2A13	Solute carrier family 2 (facilitated glucose transporter), member 13	0.0004	0.4126
A_33_P3379356	NM_001122742	ESR1	Estrogen receptor 1	0.0062	0.5516
A_33_P3382489				0.0056	0.6209
A_33_P3387110	NR_015419	LOC145783	Uncharacterized LOC145783	0.0081	−0.5672
A_33_P3389286	NM_006142	SFN	Stratifin	0.0059	0.7813
A_33_P3400292				0.0094	−0.5263
A_33_P3407429		PSMB7	Proteasome (prosome, macropain) subunit, β type, 7	0.0052	−0.4364
A_33_P3541279	AF116649			0.0021	0.8729
A_33_P3876414	AJ272176	C17orf6	Chromosome 17 open reading frame 6	0.0090	0.5759

Figure 1.

Molecular patterns of pharmacological resistance to natalizumab. For the identification of the molecular patterns underlying the pharmacological resistance to natalizumab in MS, we selected the GSE44964 microarray dataset that included expression data from CD4+ T cells isolated from patients with MS, which were stimulated in vitro with anti-CD3/-CD28 monoclonal antibodies for 48 h. Statistical differences between HR (n=6) and LR (n=6) patients were assessed using LIMMA. (A) Principal component analysis was performed using MultiExperiment Viewer software and the genes significantly modulated between HR and LR patients. The distribution of the samples of the first two components was presented. (B) Gene Ontology analysis for the genes significantly modulated between HR and LR patients was performed using the Database for Annotation, Visualization and Integrated Discovery v6.8 web-based tool. (C) Network analysis for genes significantly modulated between HR and LR patients was performed using the GeneMania. HR, high responder; LR, low responder; NF-κB, nuclear factor-κB.

Machine learning-identified genes for predicting natalizumab responsiveness

To identify a specific gene signature for natalizumab responsiveness, the UCSC algorithm was applied to the significant probes identified. The following parameters were selected: Δ=1 and ρ=1 for UCSC analysis. A total of 17 predictors of the 45 probes (Fig. 2A) were identified from UCSC analysis that were able to classify HR and LR samples with 89.2% agreement with the clinical data (Fig. 2B). Consistent with these findings, HCL and PCA based on the 17 markers were able to accurately separate HR and LR patients (Fig. 2C and D). The 17 identified predictors were presented in Table II.

Figure 2.

Prediction of pharmacological response to natalizumab. In order to determine a specific gene expression signature for predicting responsiveness to natalizumab treatment, the UCSC algorithm was applied to the genes significantly modulated between HR and LR patients. (A) The experimental procedure. (B) The number of classification errors obtained from UCSC analysis using different Δ-values (C) Hierarchical clustering of the 17 predictor genes obtained via UCSC analysis. (D) Principal component analysis using the 17 predictor genes obtained via UCSC analysis. The distribution of the samples on the first three components. HR, high responder; LR, low responder; UCSC, UnCorrelated Shrunken Centroid.

Table II.

Predictors of natalizumab responsiveness.

Probe	Gene accession no.	Gene symbol	Description
A_33_P3273552	NM_002282	KRT83	Keratin 83
A_33_P3338674	NR_033298	CCDC163P	Coiled-coil domain containing 163, pseudogene
A_24_P62783	NM_004102	FABP3	Fatty acid binding protein 3, muscle and heart (mammary-derived growth inhibitor)
A_33_P3221925	NM_152730	C6orf170	Chromosome 6 open reading frame 170
A_33_P3541279	AF116649
A_33_P3346193	NM_001043351	TPM3	Tropomyosin 3
A_33_P3243093	NM_003617	RGS5	Regulator of G-protein signaling 5
A_32_P80245	NM_001109809	ZFP57	Zinc finger protein 57 homolog (mouse)
A_24_P920048	AK092807	LOC100127972	Uncharacterized LOC100127972
A_33_P3227269	NR_024256	FLJ45983	Uncharacterized LOC399717
A_23_P381017	NM_152559	WBSCR27	Williams Beuren syndrome chromosome region 27
A_32_P7204	NM_004249	RAB28	RAB28, member RAS oncogene family
A_23_P112187	NM_032843	FIBCD1	Fibrinogen C domain containing 1
A_33_P3389286	NM_006142	SFN	Stratifin
A_23_P386478	NM_024873	TNIP3	TNFAIP3 interacting protein 3
A_33_P3289536	NM_001199835	SNX10	Sorting nexin 10
A_23_P77135	NM_080650	ATPBD4	ATP binding domain 4

Discussion

Several studies have investigated expression array phenotyping as a means for predicting drug response and clinical prognosis (16–18), and for the classification of diseases (18–25)the most common lymphoid malignancy in adults, is curable in less than 50% of patients. Prognostic models based on pre-treatment characteristics, such as the International Prognostic Index (IPI. Predicting diagnostic classes based on a sample using its expression profile is known as supervised learning or classification. The use of microarray data, although practical, poses the problem of predicting diagnostic classes using a number of genes that is notably higher than the number of sample types available. Therefore, it is necessary to select subsets of genes that are relevant for the characterization of the different diagnostic classes. In addition, the identification of specific subsets of genes may improve the classification accuracy, allow the development of cost-effective diagnostic tests and may provide novel biological insight into certain diseases. Classification can be defined as a supervised learning approach, in which the classes of a series of samples are inputted to an algorithm. This is distinct from unsupervised clustering, in which no prior knowledge of the samples is available. The aim of classification is to identify the smallest possible subgroup of genes highly associated with the known sample classes. The UCSC algorithm is based on the ‘Shrunken Centroid’ algorithm reported by Tibshirani et al (25). Briefly, genes are considered one at a time and the difference between the class centroid (the mean expression in a class) of a gene and the overall centroid (the mean expression level across all classes) of a gene is compared with the within-class standard deviation plus a Δ-value, which is determined by cross-validation, in order to minimize classification errors. In the present study, the UCSC algorithm was applied for the identification of a subset of genes that could predict the pharmacological response to natalizumab treatment among patients with relapsing-remitting MS. Natalizumab is a disease-modifying drug that can effectively reduce the frequency of relapse and short-term disability progression in relapsing-remitting MS, and it is often used as second-line treatment in patients exhibiting active disease, despite treatment with glatiramer acetate or IFN-β (26).

To the best of our knowledge, the present study is the first to identify, at the whole-genome level, the genes that were significantly modulated in HR patients compared with the LR group. Our findings suggest that a specific gene expression profile of CD4+ T cells may characterize the pharmacological responsiveness to natalizumab in patients with MS. Interestingly, no significant differences in the transcription levels of CD49d and CD29, which encode the target of natalizumab comprising the α4 and β1 subunits of VLA-4, were observed between the HR and LR groups of patients (data not shown). Additionally, we applied machine learning to select the minimum number of genes able to predict the response to natalizumab. The results indicate the genes that may be relevant for P4 medicine, which constitutes predictive, preventative, personalized and participatory medicine (27). At present, the mechanisms of resistance to natalizumab remains largely unknown. Recently, Cavaliere et al (28) applied molecular dynamics simulation to determine whether a polymorphism could induce conformational changes in VLA-4, affecting the binding affinity with natalizumab; expression profiling of circulating blood cells should be conducted for the identification of biomarkers of natalizumab resistance. The role of the genes identified in our study requires further investigation; however, certain genes may be involved in immunity and the pathology of MS. In particular, the Ras-related protein, Rab-28, has been detected in the serum of Alzheimer's disease and MS patients (29). TPM3 was determined to be phosphorylated following T cell costimulation, resulting in downregulated interleukin-2-stimulated T cells (30). The locus 3 kb upstream of the zinc finger 57, that encodes a protein likely to act as a transcriptional repressor, has been reported to be hypomethylated in CD4+ T cells from patients with MS compared with healthy controls (31). Increasing efforts are required to validate our findings determine of the role of the genes involved.

The use of biomarkers to predict natalizumab resistance in MS would lead to notable therapeutic and economic benefits; however, our study has several limitations. The number of patients is limited and no external validation could be performed. In addition, it has not been disclosed by the original authors of the microarray datasets whether natalizumab treatment was administered as first-line therapy or after failure with other medications, such as glatiramer acetate or IFN-β. Furthermore, a comparison with other second-line drugs, such as fingolimod, should be conducted. Finally, the expression profiles of unsorted and unstimulated circulating cell populations, such as whole blood cells or peripheral blood mononuclear cells should be determined for the development of simple and economically viable diagnostic tests. Despite these limitations of the present study, findings may serve as a basis for the design of personalized therapeutic options for patients with MS.