Transcriptome analysis of normal-appearing white matter reveals cortisol- and disease-associated gene expression profiles in multiple sclerosis

Brain tissue

Quantification of CRH-producing neurons and cortisol levels

Numbers of corticotropin-releasing hormone (CRH)-expressing neurons in the paraventricular nucleus (PVN) were quantified in fixed tissue as described previously [35, 58]. In short, serial 6-μm frontal sections were cut on a microtome. Delineation of the PVN was determined in thionine-stained sections. Each 100th section through the PVN was stained for CRH. Neurons that showed a nucleolus and expressed CRH were counted blinded. The total number of CRH-expressing neurons in the PVN was calculated on the basis of cell counts and the distance between the sections. Cortisol was measured by radioimmunoassay (Diagnostic Products Corporation, Los Angeles, CA, USA), using a radioactively labeled antibody that enables highly sensitive detection of cortisol levels in various types of fluid, including serum, CSF, and saliva.

Tissue processing and RNA isolation

Series of 10 cryostat sections (20 μm each) of subcortical NAWM were homogenized in Trizol (Invitrogen, Carlsbad, CA, USA). Sections preceding and following these series were stained by immunohistochemistry for proteolipid protein (PLP; Serotec, Oxford, UK) and HLA-DP, −DQ, −DR (DakoCytomation, Glostrup, Denmark) to confirm the absence of MS-lesion pathology, respectively, by ruling out demyelination and microglia/macrophage activation. RNA isolation and assessment of its quality by RNA integrity number (RIN) was performed as described previously [40]. Along with this, RNA was extracted from snap-frozen tissue dissected from various anatomical regions of control and MS brains, including MS lesions and NAWM, as well from tonsil. This RNA which was pooled to create common reference complementary RNA (cRNA), which was co-hybridized to every microarray slide to enable accurate comparison of expression levels across different cDNA microarray experiments.

Microarray hybridization

Labeling of isolated RNA was done using the Low Input Quick Amp Labeling kit (Agilent Technologies, Palo Alto, CA, USA), according to the manufacturer’s instructions. For whole-genome expression analysis, samples were hybridized to Agilent 4x44K v2 Whole Human Genome arrays (G4845A; Agilent Technologies), covering 27,958 genes. In brief, equal amounts of total RNA (50 ng) were amplified and labeled with either Cy3-CTP (experimental samples) or Cy5-CTP (reference material, obtained as described above) using the Low Input Quick Amp Labeling kit (Agilent Technologies). For hybridization, equal amounts (825 ng) of labeled samples were fragmented in Fragmentation Buffer (Agilent Technologies) for 30 min at 60 °C. Labeled and fragmented complementary RNA (cRNA) was hybridized to the array and incubated in a rotating hybridization chamber for 17 h at 60 °C. After hybridization, the array was washed subsequently for 5 min in 6 x saline sodium phosphate-EDTA (SSPE)/0.005% N-lauroylsarcosine, 1 min in 0.006 x SSPE/0.005% N-lauroylsarcosine, and 30 s in acetonitrile and dried quickly in nitrogen flow. The arrays were scanned at a resolution of 5 mm and at 5 and 100% photomultiplier tube settings using the Agilent DNA Microarray Scanner (Agilent Technologies). Scan data were extracted using Agilent Feature Extraction software (version 8.5.1; Agilent Technologies).

Normalization of gene expression and gene extraction

To allow for accurate comparison of expression levels across different cDNA microarray experiments, common reference cRNA was co-hybridized to every microarray. By using the reference cRNA, a ratio between the experimental and reference material could be calculated for every spot, and expression levels across different hybridizations were compared. These data were used for analysis by weighted gene co-expression network analysis (WGCNA) and by group-wise comparison between subgroups, which is further described below. For an overview of the experimental setup, see Fig. 1. Raw data from the extraction software was imported to the R statistical processing environment using linear models for microarray data (LIMMA)/Bioconductor (version 3.12.3) package. Feature and background non-uniformity outliers, as determined by the Agilent Feature Extraction software using default settings, were removed. The data were normalized using “between array normalization” with the “Gquantile” algorithm (http://www.bioconductor.org).

cDNA synthesis and quantitative real-time PCR

Using the same RNA from which cRNA was generated for hybridization, cDNA was synthesized with the Quantitect Reverse Transcription Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. Quantitative real-time polymerase chain reaction (PCR) was performed and analyzed as described elsewhere, with minor adaptations [40]. The amount of cDNA used per reaction was based on an input of 5 ng original RNA in a final volume of 20 μl.

Statistical analysis

To obtain lists of genes associated with either disease severity or HPA-axis activity, a WGCNA was performed on the ²log intensities to find gene modules [44, 61]. WCGNA defines, in an unbiased way, clusters of genes that show a co-expression pattern amongst samples, creating (a network with) modules of genes that are co-regulated. A soft power threshold of 9 was determined in the recommended way [44], a signed network was used, and dendrogram cluster detection was performed using a minimum cluster size of 27 and a deepSplit value of 1.

In addition to the WCGNA, linear regression analysis was performed for group-wise comparisons, using the LIMMA/bioconductor package in R (version 3.12.3). This was carried out using the limma function “lmscFit” for single channel analysis on two-color arrays with correction for intraspot correlation. A false discovery rate of 0.05 was used by multiple testing correction with the Benjamini-Hochberg method. All models were pairwise comparisons of two conditions. These group-wise comparisons were made for gene expression between all MS patients and control subjects. Furthermore, the group of MS patients was separated into a high and a low cortisol subgroup by the median for this parameter for objective group-wise comparison to identify cortisol-associated gene expression. To perform unbiased group-wise comparisons to find gene expression profiles associated with severity of MS, the patient population was divided into two subgroups by the median disease duration. Correlation coefficients were interpreted as follows: r < 0.1: no correlation; 0.1 < r < 0.3: weak correlation; 0.3 < r < 0.5: moderate correlation; r > 0.5: strong correlation. The significance level was set to 0.05. For all analyses that did not include the microarray data, we used SPSS software version 24 (IBM, Armonk, NY, USA).

All modules identified by WGCNA were subjected to the online tool DAVID, version 6.7 (http://david.abcc.ncifcrf.gov) for functional gene ontology (GO) analysis, using the whole human genome as background. For enrichment of top GO classes, the significance level was set to 0.05 (unadjusted). The minimum number of genes in each GO term was 3.

Weighted gene co-expression network analysis (WGCNA)

A first analysis on the microarray data was done by WGCNA, which in an unbiased way identifies clusters of genes that display a similar expression pattern across the different samples in the whole data set; clusters of genes that show the same expression pattern are referred to as modules. WGCNA can also be used to investigate whether such modules are related to traits of particular interest. Each module has an eigengene (the first principal component) that represents the average expression pattern of all genes within the module. Moreover, for each gene within a module, a correlation value to the module eigengene is calculated, which indicates to what extent the expression profile of an individual gene resembles that of the whole module [44, 61].

Two approaches were taken for WGCNA. The main analysis was performed on a dataset that included the samples from both control subjects and MS patients. Another analysis was performed on data from MS patients only, which was done to confirm findings made in the first approach and to identify additional gene expression profiles associated with HPA-axis activity and disease severity in MS. From both datasets, we selected out of a total of 27,958 genes covered by the microarray the 15,000 genes that showed the biggest variance across all samples for application of WGCNA. By hierarchical clustering using these data we found two outliers, one in the control group and one in the MS patient group (donor 01–069 and 98–158 respectively), which were excluded from further analyses. All data described below in the text, tables, and figures were generated by our main WGCNA analysis that included both control subjects and MS patients, unless stated otherwise. To find gene clusters associated with disease severity and HPA-axis activity in MS, we studied the correlation between gene modules identified by WCGNA and several clinical traits: disease duration and time to EDSS6 as parameters of disease severity, CSF cortisol levels and numbers of CRH-producing neurons as indicators of HPA-axis activity.

Modules generated by WCGNA correlate to HPA-axis activity and MS severity

Using the WGCNA method, we identified a total of 19 modules, of which the biggest module contained 2495 genes and the smallest 65 genes (Fig. 2). Of the identified modules, 9 showed an expression pattern that correlated with one or more of the traits of interest: cortisol levels in CSF, number of CRH-positive neurons, disease duration, and time to EDSS6. Notably, all modules that were correlated to CSF cortisol levels also showed a correlation in the same direction to numbers of CRH-expressing neurons, substantiating their association with activity of the HPA axis. In total, eight modules correlated only to indicators of HPA-axis activity, i.e. cortisol levels in CSF and numbers of CRH-expressing neurons in the PVN: the black module (r = − 0.47, p = 0.020; r = − 0.31, p = 0.100, respectively), yellow module (r = − 0.63, p < 0.001; r = − 0.57, p = 0.003, respectively), red module (r = 0.43, p = 0.030; r = 0.46, p = 0.020, respectively), salmon module (r = 0.40, p = 0.050; r = 0.52, p = 0.008, respectively), cyan module (r = 0.48, p = 0.010; r = 0.65, p < 0.001, respectively), pink module (r = 0.56, p = 0.003; r = 0.77, p < 0.001, respectively), tan module (r = − 0.63, p < 0.001; r = − 0.61, p = 0.001, respectively), and turquoise module (r = 0.25, p = 0.200; r = 0.46, p = 0.020, respectively).

A module strongly correlated to both severity of MS and HPA-axis activity

Among the gene clusters identified by the WCGNA, as depicted in Fig. 2, the lightgreen module had an eigengene that was positively correlated to cortisol levels (r = 0.56, p = 0.003), disease duration (r = 0.50, p = 0.010), and time to EDSS6 (r = 0.54, p = 0.006), indicating its positive association with both HPA-axis activity and disease severity. The module contained a total of 98 genes, encoding molecules, such as the S100 proteins S100A8, S100A9, and S100A12, which play a variety of roles in inflammation [3, 13, 21]. Functional annotation clustering indicated that the lightgreen module was enriched for genes involved in classes labeled by GO as ‘immune response’ (p = 6.0E-09), ‘negative regulation of cytokine biosynthetic process’ (p = 4.3E-03), and ‘regulation of IL-1β production’ (p = 4.3E-03), amongst others (Table 4). These GO classes contained several genes that are directly or indirectly implicated in MS, such as cystatin F (CST7), ghrelin (GHRL), IL-18 receptor accessory protein (IL18RAP), myeloperoxidase (MPO), matrix-metalloproteinase 9 (MMP9), and FMS-like tyrosine-3 (FLT3) [4, 11, 17, 24, 27, 29, 38, 43, 55, 70, 79]. The top 10 genes most strongly connected to the lightgreen module are involved in various molecular mechanisms, such as myelination (CST7), neurogenesis (ASPM), and inflammation (DEFA4, VSTM1, BPI, S100P, WISP3, NLRP12). Overall, the lightgreen module was characterized by molecules that actively regulate inflammation, though it also indicated changes in molecular pathways that affect multiple aspects of MS disease activity to slow down clinical progression.

A module positively correlated to HPA-axis activity, independent of disease severity

The pink module was positively correlated to cortisol levels (r = 0.56, p = 0.003) and numbers of CRH neurons (r = 0.77, p < 0.001). The module contained 310 genes and was enriched for molecules with the GO label ‘regulation of caspase activity’ (p = 9.4E-04), ‘heat shock protein binding’ (p = 8.0E-04), ‘regulation of lymphocyte activation’ (p = 2.0E-02), and ‘oxidation reduction’ (p = 1.1E-02), amongst others (Table 4). Moreover, the module was enriched for molecules involved in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway ‘antigen processing and binding’ (p = 1.7E-03). The molecules involved in regulation of caspase activity may be related to apoptosis of neurons and oligodendrocytes, but also play a role in microglia activation and neurotoxicity [7].

Interestingly, the pink module contained the interleukin-7 receptor (IL7R, CD127) and tumor necrosis factor ligand superfamily member 14 (TNFSF14), two molecules that are highly associated with genetic susceptibility to MS [28, 36]. In addition, the pink module contained various heat-shock proteins (HSP), which have been associated with MS pathogenesis, including HSP90, which functions as a chaperone protein for the glucocorticoid receptor (GR) [14, 66]. Also, sirtuin 7 (SRT7) was present in the GO class ‘cellular response to oxidative stress’, which has been described as a protective molecule in the brain and has been proposed as a therapeutic agent [72].

Similar to the observations regarding the lightgreen module, the pink module is also characterized by a combination of genes involved in active regulation of immunity and those involved in more neuroprotective mechanisms.

Modules negatively correlated to HPA-axis activity, independent of disease severity

The tan module had a module eigengene that was negatively correlated to cortisol levels (r = − 0.63, p = 7.0E-04) and numbers of CRH neurons (r = − 0.61, p = 1.0E-03). The module contained 253 genes and was enriched for molecules in the GO class ‘generation of precursor metabolites and energy’ (p = 7.6E-10) and ‘mitochondrion’ (p = 1.3E-10), which are likely related to changes in mitochondrial functioning and energy metabolism known to occur in MS (Table 4) .

The yellow module had a module eigengene that was strongly correlated negatively to cortisol levels (r = − 0.63, p = 7.0E-04) and numbers of CRH neurons (r = − 0.57, p = 3.0E-03). The module contained 922 genes and was enriched for molecules belonging to the GO class ‘lipid biosynthetic process’ (p = 1.7E-05) and ‘steroid metabolic process’ (p = 2.4E-02) (Table 4). One particularly interesting member of both enriched GO classes in the tan module was aconitase (ACO2). ACO2 is activated by iron and is involved in glutamate production [57], while both iron accumulation and glutamate excitotoxicity are implicated as pathogenic mechanisms in MS [30, 68].

Genes associated with severity of MS and HPA-axis activity, independent of WGCNA modules

Several molecules among the genes positively correlated to cortisol levels are well-known for their prominent role in shaping adaptive immune responses. For example, retinoic acid receptor alpha (RXRA) showed the strongest positive correlation with cortisol levels. Interestingly, RXRA-mediated signaling enhances differentiation of functionally competent CD4⁺ regulatory T cells and potently inhibits the formation of CD4⁺ T helper 17 cells, the latter of which have been strongly implicated in various autoimmune pathologies, including MS [6, 15, 18, 19, 22, 39, 62, 78]. Indeed, retinoic acid and other vitamin-A derivatives were shown to be protective in animal models of autoimmune disease [56]. Notably, IL7R was also present among the 10 genes most positively correlated to cortisol levels. In combination with its presence in the pink module, which showed the strongest positive correlation with cortisol levels in the CSF and numbers of CRH-expressing neurons on the PVN, this further confirms the positive association of IL7R expression with HPA-axis activity in MS. Another interesting members of the 10 genes most positively correlated to cortisol levels is gap junction alpha-4 protein (GJA4), as it is strongly expressed by endothelial cells and may affect leukocyte trafficking [71]. Expression of ankyrin repeat domain-containing protein 16 (ANKRD16) showed the strongest negative correlation to levels of cortisol in CSF. However, not much is known about this protein. The gene most strongly correlated to disease duration is the transcription factor NF-E2 45 kDa subunit (NFE2), which has been shown to play a role in differentiation and maturation of erythroid cells [49]. The strongest negative association to duration was seen for the gene coding for methionine sulphoxide reductase A (MSRA), which is involved in protection against oxidative stress [85].

Group-wise comparisons using linear regression (LIMMA)

The median cortisol level (220.5 nmol/l) was used to objectively divide the MS-patient population into two subgroups, allowing comparison of groups for gene expression profiles that might reveal associations with HPA-axis activity. Similarly, the median duration of MS was used to create an objective subdivision between patients with relatively severe or mild disease (i.e., having a disease course shorter or longer than 24 years, respectively), to compare these groups for identification of differences in gene expression profiles related to disease severity. Two patients in the study population had a disease duration of exactly 24 years. These were included in the group of patients with severe MS, as their disease duration was shorter than the median disease duration of 24.5 years in the entire MS-donor population of the Netherlands Brain Bank. The results of the different comparisons are shown in Fig. 3 and Additional file 2, with the 10 most strongly regulated genes in either direction highlighted. The GO classes enriched in the comparisons according to functional annotation clustering analysis are depicted in Table 6.

Patients – Controls
	GO class	Genes present
Up	Positive regulation of apoptosis (GO:0043065)	SIVA1, HTATIP2, ZAK, PML, TNFSF14, RPS27L, TLR4, ITSN1, CTNNBL1, RPS3, CASP3, HTRA2, CD44, CDKN2C, RPS3A, SOS1, MTCH1, TGM2, AATF, RUNX3, RPS27A, DEDD2, CEBPG, FADD, BAD, SOD1, TXNDC12, EI24, RNF7, NAIF1, HSPD1
	Activation of caspase activity (GO:0006919)	SIVA1, MTCH1, PML, HSPE1, HSPD1, RPS3
	Regulation of T cell activation (GO:0050863)	CD47, CASP3, IL6ST, NCK1, TGFBR2, TNFSF14, BAD, HSPD1, SOD1
	Cytokine binding (GO:0019955)	TNFRSF1A, ACVRL1, IL10RB, IL6ST, LEPR, TGFBR2, ENG, IFNAR1, ACVR1
Down	Neuron differentiation (GO:003082)	EGR2, GNAO1, ATL1, TBCE, NTNG2, DSCAML1, APP, RASGRF1, GHRL, MAPK8IP3, SEMA3B, NTM, C17ORF28
	Cell projection (GO:0042995)	RTN4, MYO5A, BBS5, KIAA1598, SSH1, ATL1, NEDD9, GIPC1, DNAH2, CPEB1, APP, CTTN, DNAI1, ARHGEF4, DBNL, ARHGEF7, FSCN1, LDB3, DNAI2, PCM1, CAMK2N1, RASGRF1, IFT172, MAPK8IP3, GHRL, SPEF1
	Regulation of lipid metabolic process (GO:0019216)	TNF, DHCR7, SF1, ACACB, PPARGC1A
MS High Cortisol – Low Cortisol
	GO class	Genes present
Up	Immune response (GO:0006955)	LAIR1, CEBPB, IL1RL1, IFITM2, TLR2, TNFSF14, CALCOCO2, SLC11A1, C1QB, UNC13D, APOL1, XBP1, LILRB3, IL4R, BCL3, HSPD1, VSIG4
	Regulation of cytokine production (GO:0001817)	INHBA, SLC11A1, CEBPB, TLR2, BCL3, NFKB1, BCL6, HSPD1, VSIG4, SRGN
	Myeloid cell differentiation (GO:0030099)	INHBA, RPS19, JMJD6, RPS14, BCL6, ZBTB16, RUNX1, CBFB, TIMP1
	Negative regulation of myeloid cell differentiation (GO:0045638)	INHBA, NFKBIA, ZBTB16, RUNX1
Down	Apoptosis (GO:0006915)	RTN4, CKAP2, POLR2G, DNM1L, TM2D1, EGLN3, RRAGA, PIGT, BAG1, NGFRAP1, EIF2AK2, MAGEH1, NDUFS1, PUF60, ZIM2
	Negative regulation of neurogenesis (GO:0050768)	RTN4, NOG, NF1, OMG
MS Mild Disease – Severe Disease
	GO class	Genes present
Up	Neuron projection development (GO:0031175)	APP, GNAO1, EGR2, ATL1, RASGRF1, TBCE, MAPK8IP3, NTNG2, DSCAML1, GHRL, SEMA3B
	Neuron differentiation (GO:0030182)	EGR2, GNAO1, ATL1, TBCE, NTNG2, DSCAML1, APP, RASGRF1, GHRL, MAPK8IP3, SEMA3B, NTM, C17ORF28
	Inflammatory response (GO:0006954)	HDAC5, YWHAZ, TNF, NDST1, TOLLIP, CCL3L3, ITIH4, CCL4L1, CCL4
Down	Lysosome (GO:0005764)	HGSNAT, SGSH, NAGLU, MFSD8, LIPA, GM2A, PPT1, CD63, ASAH1
	Steroid metabolic process (GO:0008202)	HSD17B11, MBTPS2, LIPA, CYB5R2, INSIG2, PRKAA1, CLN6

Comparison of MS patients and controls

Compared to control subjects, a total of 778 genes was significantly upregulated in MS patients, whereas 544 genes were downregulated. The four most strongly upregulated genes in NAWM of MS patients code for antibody subunits, which is possibly related to the synthesis of auto-antibodies in MS (Fig. 3a) [20]. Furthermore, NAWM of MS patients showed an increased expression of heat-shock proteins, HSPA1A and HSPA6. Among the genes upregulated in NAWM of MS patients compared to that of control subjects, there was an enrichment for molecules in GO classes associated with induction of apoptosis, activation of caspases, regulation of T-cell activation, and cytokines binding (Table 6).

Compared to control NAWM, the gene showing the most strongly decreased expression in MS NAWM was transthyretin (TTR), which was downregulated 4-fold (Fig. 3a). TTR is an important carrier in serum and CSF for the thyroid hormone thyroxin (T4) and for retinol, a form of vitamin A. Importantly, both TTR and T4 have been implicated in MS. Oxidative modifications of TTR protein and decreased levels of T4 were present in the CSF, and not in the serum, of MS patients and were correlated with disease duration [67]. Moreover, T4 was shown to play an important role in activating oligodendrocyte precursor activation and instill myelination [9].

GO classes overrepresented among the genes downregulated in NAWM of MS relative to control NAWM were related to neuron differentiation, cell projection, and regulation of lipid metabolism (Table 6).

Comparison of MS patients with high and low cortisol

In MS patients with high cortisol a total of 270 genes was upregulated, whereas 472 genes were downregulated compared to patients with low cortisol. Of note, the third most-strongly upregulated gene in MS patients with high cortisol by more than 4-fold was CD163 (Fig. 3b), which is a glucocorticoid-responsive gene that can be induced in myeloid immune cells, such as macrophages and microglia [59, 80].

Neuronal pentraxin-2 (NPTX2) was the gene most strongly increased in NAWM of MS patients with high cortisol, showing an almost 6-fold higher expression in comparison to NAWM of patients with low cortisol (Fig. 3b). NPTX2 is an immune-related molecule with structural similarities to several acute phase proteins that is thought to be essential in the compensatory synaptic response that occurs during prolonged neuronal inactivity [69]. An interesting finding regarding the 10 most strongly downregulated genes in NAWM of MS patients with high cortisol was the presence of the purinergic receptor P2RY12 (Fig. 3b), a microglia signature gene shown to play a major role in microglial activation, synaptic plasticity, and closure of the injured blood-brain barrier [8, 52, 76].

Interestingly, analysis by functional annotation clustering on the genes that are upregulated in MS patients with high cortisol levels, in comparison to patients with low cortisol levels, are enriched for several GO classes associated with (regulation of) inflammation (Table 6). In contrast, there is an overrepresentation of GO classes associated with apoptosis and negative regulation of neurogenesis among the genes showing decreased expression in MS patients with high cortisol.

Comparison of MS patients with severe and mild disease

In total 202 genes were upregulated in patients with mild MS, whereas 154 genes were downregulated compared to those with severe MS. There is a clear overlap between genes upregulated in MS patients with mild disease and those with high cortisol, as for example S100A8 and solute carrier organic anion transporter family member 4A1 (SLCO4A1) are present in the top 10 most strongly upregulated genes in both groups. IL-1 receptor like 1 (IL1R1) is among the 10 most strongly upregulated genes in mild MS (Fig. 3c), and is known to be a receptor for IL-33 and has been shown to be involved in induction of Th2 responses during allergic inflammation [73]. The oxidized low-density lipoprotein receptor 1 (OLR1) gene, also known as LOX-1, was the most downregulated gene in NAWM of MS patients mild disease (Fig. 3c). Interestingly, this gene was found to be strongly related to the extent of demyelination in white matter MS lesions [33]. The upregulated genes in NAWM of MS patients with mild disease were enriched for GO classes involved in neuron development and differentiation as well as the inflammatory response.

Distinct expression profile of inflammatory and GC-related genes in patients with high cortisol or mild MS

For the next analysis, we selected genes included in GO classes related to either inflammation or glucocorticoid signaling, including: ‘regulation of acute inflammatory response’ (GO:0002673), ‘chronic inflammatory response’ (GO:0002544), ‘macrophage differentiation’ (GO:0030225), ‘microglial cell activation’ (GO:0001774), ‘glucocorticoid-receptor signaling pathway’ (GO:0042921), and ‘cellular response to glucocorticoids’ (GO:0071385) and ‘glucocorticoid biosynthetic process’ (GO:0006704). These GO classes were selected to specifically study to what extent HPA-axis activity impacts on molecular mechanisms in NAWM of MS and how this affects neuroinflammation and disease severity. To this end, the selected genes were compared for their expression between five groups: control subjects, MS patients with high and low cortisol levels, and MS patients with severe and mild disease. This was done by cluster analysis, to find out which groups most strongly resemble each other in expression profiles for the selected genes. Interestingly, a distinct expression profile was present for inflammatory and glucocorticoid-associated genes in MS patients with high cortisol and those with mild MS (Fig. 4a and b). These subgroups differed in this respect from MS patients with low cortisol or severe MS, and control subjects. Remarkably, the subgroups of MS patients with high cortisol and those with mild MS were especially similar with respect to higher average expression values for inflammatory genes, when compared to the other three subgroups.

Molecules selected for validation

Expression of in total 14 genes strongly associated with HPA-axis activity and/or disease duration was validated by qPCR. For this purpose, the same tissue was used as for generating the original microarray data. For several of the validated targets no alterations have been described in the context of MS before, such as leucine-rich repeat-containing protein 32 (LRRC32) and nidogen-1 (NID1). Others have been implicated in MS pathogenesis, but were never studied in human subjects for their association with disease severity or HPA-axis activity, such as IL7R, tissue transglutaminase (TGM2), MAC-inhibitory protein (CD59), and gap junction alpha-4 (GJA4) [10, 41, 77, 83]. For 11 out of 14 genes the expression profile was indeed confirmed by qPCR (Fig. 5). Genes for which the qPCR data did not confirm the expression profile detected by microarray were SLC4, CD44, and CST7. For CST7, the expression detected by the microarray as well as the qPCR was very low. Therefore, the expression differences for CST7 observed in the microarray data may have been a technical artifact caused by background variation. For SLC4 and CD44, quite some variation in expression between subgroups was detected by qPCR, but apparently those differences were not large enough to be significant.