Dysregulation of MS risk genes and pathways at distinct stages of disease.

OBJECTIVE: To perform systematic transcriptomic analysis of multiple sclerosis (MS) risk genes in peripheral blood mononuclear cells (PBMCs) of subjects with distinct MS stages and describe the pathways characterized by dysregulated gene expressions.
METHODS: We monitored gene expression levels in PBMCs from 3 independent cohorts for a total of 297 cases (including clinically isolated syndromes (CIS), relapsing-remitting MS, primary and secondary progressive MS) and 96 healthy controls by distinct microarray platforms and quantitative PCR. Differential expression and pathway analyses for distinct MS stages were defined and validated by literature mining.
RESULTS: Genes located in the vicinity of MS risk variants displayed altered expression in peripheral blood at distinct stages of MS compared with the healthy population. The frequency of dysregulation was significantly higher than expected in CIS and progressive forms of MS. Pathway analysis for each MS stage-specific gene list showed that dysregulated genes contributed to pathogenic processes with scientific evidence in MS.
CONCLUSIONS: Systematic gene expression analysis in PBMCs highlighted selective dysregulation of MS susceptibility genes playing a role in novel and well-known pathogenic pathways.

MS is an inflammatory demyelinating disease of the CNS with evidence of immune dysfunction.[1] The first clinical episode with features suggestive of MS is classified as clinically isolated syndrome (CIS), unless lesion dissemination in time and space ratifying MS diagnosis is evidenced.[2] Approximately, 85% of individuals develop the relapsing-remitting form of MS (RR-MS), whereas 15% of MS individuals experience a progressive course (primary progressive MS, PP-MS). After a variable time, most RR-MS subjects advance to a secondary progressive (SP) phase (SP-MS), where neurologic worsening takes place without periods of remission.[3] Several MS risk variants have been uncovered by genome-wide association studies (GWAS), and many of them are located close to immunologically relevant genes,[4] suggesting that immune dysfunction may be partly genetically determined. Although effective in highlighting MS risk alleles, GWAS have failed in dissecting the genetic components of the susceptibility to distinct MS clinical forms. In fact, MS is a multifactorial disorder determined by the complex interaction between genetic and environmental factors, whose integration occurs at the epigenetic level and determines gene expression. Our group has recently shown the importance of blood transcriptomics in uncovering gene expression changes and transcriptional regulators in MS.[5] In this study, we have (1) examined the expression levels of known MS susceptibility genes in peripheral blood mononuclear cells (PBMCs) of CIS, RR-MS, PP-MS, SP-MS, and control individuals, (2) identified a panel of blood transcriptional signatures for distinct MS forms, and (3) explored the pathways contributed by the dysregulated susceptibility genes.

METHODS

Human subjects and blood sampling.

Investigations were conducted according to the principles expressed in the Declaration of Helsinki, and peripheral blood was drawn after signing of the institutional informed consent. We recruited 142 patients with MS (46 CIS, 52 RR-MS, 23 PP-MS, and 21 SP-MS) and 40 healthy controls (HCs) of Italian origin for the generation of the main transcriptomic data set by Illumina platform. Demographic and clinical parameters are shown in table 1. A second distinct cohort comprising 21 RR-MS, 15 PP-MS, 13 SP-MS, and 27 HCs was enrolled according to the same inclusion criteria and used for validation with a distinct array platform (table 1). Finally, a third independent cohort comprising 31 CIS, 30 RR-MS, 24 PP-MS, 21 SP-MS, and 29 HCs was included and used for validation experiments by quantitative PCR (q-PCR) (table 1). Patients with MS were diagnosed according to McDonald criteria[6] and were not suffering from any other acute or chronic inflammatory diseases or other autoimmune disorders. Furthermore, they had not started any immunomodulatory therapy for MS yet, as recruitment was performed over the last 15 years, in a period when decision to treat was not established and widespread as nowadays. Blood sampling was performed between 30 and 90 days after the first clinical attack in patients with CIS, and at least 4 weeks after the last clinical attack or steroid treatment for RR-MS subjects. All participants had peripheral blood counts within the reference range. All blood samplings were performed between 9 and 12 am

Table 1

Subjects demographics and clinical information

PBMC isolation and RNA extraction.

PBMCs were isolated using a discontinuous density gradient (Lymphoprep; Nycomed, Oslo, Norway). Viable cells were counted by Trypan blue (Sigma-Aldrich, Milan, Italy) exclusion. Total RNA was extracted using Tri Reagent (Ambion; Applied Biosystems, Monza, Italy) and stored at −80°C.

Generation of gene expression data sets.

RNA quality was checked using Bioanalyzer 2100 (Agilent, Milan, Italy). Complementary RNA (cRNA) synthesis was performed using the Illumina TotalPrep RNA Amplification Kit (Ambion) according to the manufacturer's protocol. Hybridization of the cRNA relative to the first casistics was performed on Illumina Human Ref-8 v2 arrays (Illumina, Son, Netherlands). GenomeStudio GX Software (Illumina, San Diego, CA) was used to extract the array raw data. All the raw data were background subtracted by NEC method and normalized using cubic spline normalization as implemented in the software. We performed batch correction using COMBAT to remove batch effects in the PBMC expression data.[7] The second independent transcriptomic data set was generated by Affymetrix GeneChip Human Genome U133 Plus 2.0 Array. All the necessary preprocessing steps were performed using robust multiarray average algorithm present in the Affy package in Bioconductor and keeping genes with multiple probes as individual transcripts. The MIAMI compliant microarray data have been deposited in the EBI ArrayExpress database (Accession numbers E-MTAB-4890 and E-MTAB-5151). Probes with a mean intensity value lower than 100 in all experimental groups were filtered out. Furthermore, we applied Pearson correlation to identify and remove probes correlating with age in the healthy population of the Illumina study, where age ranged from 22 to 57 years. To identify differentially expressed genes (DEGs) between each clinical form and the healthy population, we performed a 2-sample t test with a p value threshold of 0.05. A detailed flow diagram of this study is depicted in figure 1. To verify whether stage-specific transcriptomes were enriched in dysregulated risk genes, we measured the frequency of differential expression of GWAS genes and compared it to that of DEG in the global transcriptome by χ2 test with Yates' correction in GraphPad.

Figure 1

Flow chart of the study

CIS = clinically isolated syndrome; DEG = differentially expressed gene; FDR = false discovery rate; GWAS = genome-wide association studies; HC = healthy control; MS = multiple sclerosis; PP-MS = primary progressive MS; RR-MS = relapsing-remitting MS; RT-PCR = real-time PCR; SP-MS = secondary progressive MS.

Collection of GWAS genes and their mapping to gene expression data.

Genes located in the vicinity of MS susceptibility loci were collected from the supplementary tables 3 and 8 published by the International Multiple Sclerosis Genetics Consortium[4] and reporting 48 gene variants from previous GWAS studies[8-10] and 110 non–major histocompatibility complex (MHC) susceptibility variants described in 2013.[4] Figure e-1 at Neurology.org/nn summarizes the complete list of GWAS genes considered in our study. To map MS susceptibility genes to the relative Illumina probe identifiers, all the filtered probes in the Illumina data set were submitted to DAVID Gene Accession Conversion Tool,[11] using Illumina HumanRef-8 v2 as a background reference so to retrieve Entrez and gene symbol identifiers for each probe. The same probe list was submitted to BioMart ID,[12] conversion package in R Bioconductor,[13] as a further control. In this way, gene expression levels were retrieved for 182 Illumina probes corresponding to 146 GWAS genes (figure 1). Similarly, GWAS genes were annotated to the corresponding Affymetrix probes using BioMart ID,[12] conversion package in R Bioconductor,[13] and gene expression levels were retrieved for 344 Affymetrix probes relative to 149 MS susceptibility genes.

cDNA synthesis and real-time PCR.

The purity and the quality of RNA were analyzed using NanoDrop8000 (Thermo-Scientific, Waltham, MA). Two-microgram aliquot of total RNA was reverse transcribed to cDNA using random primers and Superscript III reverse transcriptase (all from Invitrogen, Carlsbad, CA) following the manufacturer's protocol. Real-time (RT-PCR) was performed using TAQMan Universal Master Mix (Applied BioSystems). PSMB1 was used as housekeeping gene based on its low coefficient of variation in the Illumina transcriptomic data set. Transcript levels of the target gene CD86 were measured by the TAQMan assay corresponding to the gene sequence detected by the Illumina probe and expressed as percentage of the housekeeping gene. Normality of data distribution was assessed by Kolmogorov-Smirnov statistics, and nonparametric Mann-Whitney test was used to compare mean values. All p values were 2-sided and subjected to a significance level of 0.05.

Pathway analysis and literature mining.

Pathway enrichment analysis was performed starting from the distinct dysregulated expression profiles isolated in each clinical form of disease. MS specialty modules in MetaCore (Thomson Reuters, New York, NY) are manually curated pathways associated with MS pathology, which are then ranked based on statistical significance. We selected pathways passing the false discovery rate (FDR)-corrected p value threshold of 0.05 and containing 3 DEGs and at least 1 GWAS DEG. For literature mining, PubMed database was used to collect scientific evidence about the role of the identified biological pathways in MS pathogenesis by searching the name of the pathway and “Multiple Sclerosis” or “EAE.”

RESULTS

Transcriptomic analysis in MS.

Genes located in the vicinity of MS susceptibility loci (hereon called as GWAS genes) were collected from the supplementary tables 3 and 8 published by the International Multiple Sclerosis Genetics Consortium[4] (see also figure e-1). To examine the expression levels of the 162 GWAS genes, we generated high-throughput gene expression data from PBMC of CIS, RR-MS, PP-MS, SP-MS individuals, and HCs (table 1) by Illumina arrays, recovered the gene symbols for GWAS genes, and mapped them to the relative Illumina probes. As shown in figure 1, there were 72 probes for 68 expressed GWAS genes passing signal intensity filter. Next, we tested whether these probes were differentially expressed in any disease group compared with healthy controls. We observed 29 GWAS DEGs in CIS, 19 in RR-MS, 35 in PP-MS, and 35 in SP-MS (table 2). Notably, 6 GWAS DEGs were concordantly dysregulated in all disease groups compared with HCs, and several other GWAS DEGs were common to at least 2 MS forms (figure 2A). Among them, NDFIP1, AHI1, FCRL1, and STAT4 were concordantly regulated in CIS and RR-MS but not in progressive forms of disease, whereas 11 MS risk genes, including TYK2, IL7R, and TNFRSF1A, were differentially expressed in PP-MS and SP-MS but not in CIS and RR-MS subjects (figure 2A). Finally, some GWAS transcripts were altered at single stages of disease, suggesting distinct pathogenic alterations in distinct disease courses (figure 2A). To verify whether MS transcriptome was enriched in dysregulated GWAS genes, we measured the frequency of differential gene expression in the list of expressed GWAS genes and in the global transcriptome for each disease subtype (figure 1, table 2) and verified whether the frequency of GWAS DEG was significantly higher in the GWAS gene list than the expected frequency of dysregulation in a random selection of transcripts using χ2 statistics. As shown in table 2, we found a clear and significant enrichment of GWAS DEG in CIS, PP-MS, SP-MS, and not in RR-MS, demonstrating that the frequency of dysregulation in genes located in the vicinity of MS susceptibility variants is higher than expected in at least 3 distinct stages of disease. Then, we checked the expression of the identified GWAS DEG in PBMC of an independent cohort of subjects including 21 RR-MS, 15 PP-MS, 13 SP-MS, and 27 HCs, for which gene expression was generated using the Affymetrix array platform (table 1). After applying the same statistical threshold of the Illumina data set, several dysregulated GWAS DEGs were identified for each MS form in the Affymetrix data set. In fact, a total of 26 GWAS DEGs were concordantly upregulated (or downregulated) compared with the healthy population in 2 independent transcriptomic data sets (figure 2B). Finally, we used a third case-control cohort to validate the expression of CD86, one of the GWAS genes upregulated in all forms of disease in the Illumina and Affymetrix data sets, by q-PCR (figures 2C and e-2).

Table 2

Probability of GWAS DEG enrichment at distinct stages of MS

Figure 2

Dysregulated MS risk genes and pathways at distinct stages of disease

(A) MS susceptibility genes dysregulated at distinct MS stages compared with as measured by Illumina microarrays. (B) Dysregulated MS susceptibility genes validated in 2 distinct array platforms. For A and B, fold-change heatmaps represent upregulated (red), downregulated (green), and unchanged (gray) expressions in the disease group compared with the healthy population. (C) Real-time PCR validation of CD86 transcript in peripheral blood mononuclear cells of a third case-control cohort. Refer also to figure e-2 for data elaboration after removal of a few outliers (samples with expression levels above 500). *p < 0.05, ***p < 0.005. (D) Venn diagram representing number of common and unique enriched pathways for each disease stage. CIS = clinically isolated syndrome; HC = healthy control; MS = multiple sclerosis; PP-MS = primary progressive MS; RR-MS = relapsing-remitting MS; SP-MS = secondary progressive MS.

Biological pathways contributed by MS GWAS genes.

To verify whether the dysregulated GWAS genes were involved in specific biological pathways, we applied the MetaCore pathway analysis program which uses a curated annotation database of biological functions and identifies the significantly enriched pathways for a given DEG list. Critically considering the limited number of dysregulated GWAS genes for a pathway analysis, for each disease group, we enriched the list of GWAS DEG with those dysregulated genes which were not located in the vicinity of MS susceptibility variants, performed pathway analysis using the global DEG list as input, and selected those pathways passing FDR-corrected p value cut-off of 0.05 and containing 3 global DEGs and at least 1 GWAS DEG for each disease form (figure 1). We found 165, 95, 184, and 223 significant pathways for CIS, RR-MS, PP-MS, and SP-MS respectively, with 65 pathways shared by all the disease courses and others common to at least 2 or 3 clinical forms (figure 2D, table e-1). Of interest, the 65 shared pathways belonged to major biological themes relative to transcription and translation, cell signaling and proliferation, cancer, cytoskeleton remodeling and cell migration, apoptosis, immunity, and nervous system (table 3). Functions related to immunity and cell signaling were predominant and included pathways regulating innate and adaptive immunity, and triggered by lipids, cytokines, and growth factors. An interactive link in the table e-1 opens the graphical maps representing the overall genes and the dysregulated MS transcripts involved in each pathway. It is important that when searching for scientific evidence about a role of the pathways in MS or its animal model by literature mining, we found it in most cases (table 3). Thus, bioinformatical annotation of gene expression data led to the identification of a core of pathogenetic functions common to all MS stages and characterized by dysregulation in genes including MS susceptibility genes.

Table 3

Shared pathways across MS stages

DISCUSSION

In this study, we provide evidence that genes located in the vicinity of MS risk variants display dysregulated expression in peripheral blood of subjects with distinct MS stages compared with healthy individuals, that the frequency of dysregulation is significantly higher than expected, and that those genes contribute to pathogenic pathways in MS.

The most recent Immunochip study has highlighted more than 100 novel non-MHC risk variants in MS and has underlined that several genes close to the identified polymorphisms play a key role in immunity, thus suggesting genetic susceptibility to a dysimmune state.[4] On the other hand, genetics contributes only to part of the risk of developing MS and has been unsuccessful in determining predisposition to distinct clinical courses of MS, indicating an important role for environmental factors. Thus, the potential effect of the identified GWAS genes during distinct stages of MS remains elusive and needs clarification. As a first step into this direction, we have verified whether the GWAS genes postulated by the Immunochip analysis have altered transcript expression in peripheral mononuclear cells in MS subjects starting from the early phase of disease (CIS) to the distinct clinical courses most of the patients experience. Here, we show that although a core of GWAS genes is concordantly dysregulated in all MS stages when compared with the healthy population, each disease group is characterized by a distinct set of dysregulated GWAS genes, clearly indicating that alteration in the expression of MS risk genes is not a stable feature reproduced throughout the distinct phases and forms of disease, but is a selective event for each disease stage. Notably, 26 MS susceptibility genes were concordantly regulated in 2 independent cohorts measured with distinct array platforms which use different probes to detect gene expression, emphasizing the robustness of several dysregulations. A GWAS DEG in our study was CD86 (B7-2), a costimulatory molecule expressed by a wide range of leukocytes[14] and critically involved in T-lymphocyte activation and proliferation.[15,16] Monocytes of patients with RR-MS display higher expression of CD86 compared with healthy controls and interferon β treatment triggers reversal regulation of CD86.[17] Here, we demonstrate the overexpression of CD86 transcripts in PBMC at all the stages of MS compared with HCs in 2 independent microarray data sets and in a third independent cohort using q-PCR. Further, we report that the expression of JDP2, MAF, MAPK3, and RGS1 is upregulated in MS. Of interest, RGS1 is also upregulated by pathogenic Th17 T cells during experimental autoimmune encephalomyelitis (EAE).[18] Conversely, BACH2, a critical transcriptional factor regulating the balance between tolerance and immunity by repressing CD4+ T-cell differentiation,[19] is downregulated in pathogenic T cells in EAE mice[18] and in the blood of RR-MS subjects.[20] Our study confirms this observation and extends it to the other MS stages. Very limited information is present in the literature about the expression of other GWAS genes in MS blood cells. A recent report has highlighted the downregulation of the transcription factor EOMES in whole blood of CIS and a pooled group of MS subjects with various disease courses compared with HCs by both RNA-seq and RT-PCR analysis.[21] We confirm the low transcript levels in CIS but not in the distinct clinical forms when analyzed separately. Furthermore, our study dissects the previously reported differential regulation of IKZF3, FOXP1, ZNF438 in CIS, and all MS clinical forms[21] and reproduces IL7R downregulation in PP-MS.[22] The regulation of the other GWAS genes during MS finds the first description in our study. Concerning the frequency of dysregulation of MS susceptibility genes, a recent work shows that CD4-positive pathogenic T cells during EAE are significantly enriched in dysregulated GWAS genes.[18] It is important that this observation is conserved in human PBMC, as we demonstrate a significantly high frequency of dysregulated GWAS genes in CIS and progressive forms of MS but not in RR-MS. As the potential to find transcriptional dysregulation in RR-MS data, however, was comparable with that in CIS (25.68% and 23.67% global dysregulation in CIS and RR-MS, respectively), the lack of enrichment in GWAS DEG in RR-MS has most probably a biological reason, whose explanation remains obscure.

Thus, to our knowledge, we report the first systematic expression analysis for MS GWAS genes using large cohorts of untreated subjects with distinct MS stages and HCs, and describe several novel GWAS dysregulations, whose single functional effects on MS pathogenesis deserve further investigations in appropriate animal models. In addition, genome-wide expression profiling has provided relevant information about the overall quality and frequency of dysregulation in PBMC transcriptomes at distinct MS stages. Previous studies have searched for pathways in MS using the full list of GWAS genes under the assumption that all GWAS genes are dysregulated in MS in all tissues and at all disease stages.[23,24] Our study demonstrates that this hypothesis may be wrong, as MS GWAS genes are not all dysregulated at the transcriptional level in PBMC and that there may be differences among disease courses. To investigate the combined effects of multiple dysregulations, we have searched for biological processes significantly enriched in DEG, including MS risk genes, as we have considered that MS is not a pure genetic disorder and pathogenic pathways may be contributed by transcriptional dysregulation determined by environmental factors as well. We report here the complete list of pathways obtained for each disease stage. Notably, several of them were detected in single forms of disease, whereas others were common to more than 1 MS stage with 65 pathways appearing in all the distinct MS subtypes. The shared pathways involve well-known functions in MS pathogenesis as cell signaling and migration, immunity, and apoptosis. Scientific evidence about a pathogenic role for most of these shared pathways is already available in the literature; however the contribution of the single pathways to the distinct MS stages still needs to be clarified. Furthermore, we describe a series of novel candidate pathogenic pathways in MS, which deserve investigation and validation.

Overall, we offer the scientific community a large transcriptomic data set for expression studies in MS, the evidence about dysregulation of several known MS risk genes in blood at distinct MS stages, and the complete list of significant pathways contributed by GWAS genes for each disease form. This information may serve as the basis for novel translational medicine studies investigating the molecular framework underlying MS pathogenesis.