Discovery of grey matter lesion-related immune genes for diagnostic prediction in multiple sclerosis

Data collection

In the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), the keywords “multiple sclerosis”, “grey matter lesion”, “expression profiling by array” and “Homo sapiens” were used as screening criteria for the initial search. Datasets were then further filtered by applying the criteria below: (1) contained MS and normal group; (2) the samples were derived from the grey matter of the brain; and (3) the sample size was not less than 10. According to the above criteria, the GSE135511 dataset annotated by GPL6883 (Illumina HumanRef-8 v3.0 expression beadchip) and GSE131282 annotated by GPL10558 (HumanHT-12 V4.0 expression beadchip) were selected and downloaded. The GSE135511 dataset contained several subgroups including control GM, normal appearing GM, and GM lesion, etc. A total of ten samples of control GM tissues (GSM4013310 to GSM4013319) and ten samples of GM lesions (GSM4013300 to GSM4013309) were finally selected for subsequent study and analysis. GSE131282 was used as a validation dataset.

Differentially expressed gene identification

The “R” software (R v4.1.3) was adopted for the analysis. By using the “limma” package (Ritchie et al., 2015), a filter of differentially expressed genes (DEGs) in GSE135511 was performed with cutoff values of adjusted p-value < 0.05 and |log2 fold change (FC)| ≥ 1. Volcano plots were created to show the differences in gene expression levels for the two groups. The clustering results are displayed by using the “pheatmap” package, which shows the top 20 up-regulated genes and the top 20 down-regulated genes.

Functional enrichment analysis of DEGs

To reveal the biological functions of DEGs in GSE135511, enrichment analyses of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were carried out using the R package “clusterProfiler” (Yu et al., 2012). By employing the “enrichplot” package, the results were visualized as bubble plots. The critical values of GO and KEGG were determined as p < 0.05.

Weighted gene co-expression network analysis

WGCNA was performed on the gene expression profiles and clinical data in GSE135511 by using the “WGCNA” package to determine potential functioning modules (Langfelder & Horvath, 2008). All samples were clustered by average linkage and Pearson correlation value. A suitable soft threshold β was defined for the establishment of the scale-free network. The adjacency matrix was transformed into the topological overlap matrix (TOM). A hierarchical clustering tree was structured, and the branches of the clustering tree represented gene modules (minimum gene number of gene modules is 30). The correlation between the obtained modules and the clinical trait was analyzed, and the modules significantly correlated with the clinical trait were selected. A heatmap was also constructed to illustrate the correlation among different modules. The scatter plots were constructed to show the correlation between gene significance and module membership.

Screening of immune-related genes

A list of immune-related genes was obtained from two databases: ImmPort (https://www.immport.org/) and InnateDB (https://www.innatedb.ca/). Modules most related to MS obtained by WGCNA were selected. The SangerBox (http://vip.sangerbox.com/) was used to draw a Venn diagram and show the intersection of immune-related genes and gene modules. The intersection genes were defined as immune-related mRNAs.

Least absolute shrinkage and selection operator model construction

The least absolute shrinkage and selection operator (LASSO) method could obtain a refined model by employing a penalty function to reduce variable numbers. Applying the “glmnet” R package, LASSO regression analysis was performed, and penalty parameters were tuned by 10-fold cross-validation to filter immune-related candidate genes.

Protein-protein interaction network analysis

The online STRING database (https://string-db.org) was used to identify potential interaction among the genes screened by LASSO. Then the result of STRING was uploaded to Cytoscape software. The cytoHubba plug-in was used to explore the hub genes in the PPI network. Top ten genes were generated using six common topological analysis methods in cytoHubba plug-in. The common hub genes were selected to perform subsequent analyses.

Immune cell infiltration analysis

CIBERSORT is a widely used algorithm that uses linear support vector regression modeling to deconstruct the expression matrix of immune cell subtypes. To identify the immune characteristics of MS samples in GSE135511, the CIBERSORT was applied for assessing the proportions of 22 immune cell types (Newman et al., 2015). The relationship between immune cell subtypes and candidate genes were conducted by Pearson correlation analysis with “psych” and “corrplot” packages. The absolute value of correlation coefficient greater than 0.5 and p-value less than 0.05 were set as the cut-off.

Nomogram model construction

In order to predict the onset of MS, a diagnostic nomogram model was developed by applying the “regplot” package. The accuracy and consistency of the nomogram model were assessed using calibration curves and decision curve analysis (DCA) curve. Receiver operating characteristic (ROC) curve analysis was performed with the “pROC” package. Diagnostic power of the hub gene in the GSE131282 dataset was checked with area under the curve (AUC).

Gene set enrichment analysis

GSEA was carried out using the “clusterProfiler” package in order to identify biological functions of the hub genes in GSE135511. We chose “c2.cp.kegg.v7.5.1.symbols.gmt” as the reference gene list, |NES| > 1 and FDR < 0.25 indicated significant enrichment.

Cuprizone-induced MS model

Male C57BL/6 mice (six-to-eight-week-old) were purchased from Zhejiang Vital River Laboratory Animal Technology Co., Ltd. (Zhejiang, China). Twenty mice were randomly divided into two groups (control and CPZ) of ten animals and each five mice were kept in one cage. All mice were housed in the specific pathogen-free facility under a 12-h light/dark cycle. The ambient temperature was 20–25 °C and the humidity was 50% ± 10%. For the induction of MS model, the CPZ mice were fed with 0.2% cuprizone chow for 5 weeks. The control mice were kept on a normal diet. Water was available ad libitum. Upon withdrawal of the cuprizone diet, the mice were euthanized by carbon dioxide asphyxiation, then the brain samples were collected. All experiments were approved by Experimental Animal Ethics Committee of Henan University of Chinese Medicine (Code: DWLL202103124). No excess experimental animals were used in this experiment.

Luxol fast blue staining and immunohistochemical staining

Paraffin brain sections were deparaffinized, and LFB staining was applied to analyze demyelination according to standard protocol. For immunohistochemical (IHC) staining, sections were boiled in EDTA buffer for 20 min for antigen repair. The blocking solution was a PBS solution containing 5% BSA. Sections were then incubated with primary antibodies for MBP (1:50; Millipore), PDGFRB (1:1,000; Proteintech), TLR9 (1:500; Servicebio), and CCL5 (1:50; Huabio) overnight at 4 °C. After incubating the corresponding secondary antibody, images were obtained using Olympus BX53 microscopy. Three animals for each group were tested.

Western-blotting analyses

Total proteins of cerebral cortex tissues of three animals for each group were extracted with RIPA buffer containing 1% PMSF and separated by 10% SDS-PAGE. Then, proteins were transferred to PVDF membranes. After blocking with 5% delipidated lipids for 1 h, the membranes were incubated overnight at 4 °C with the following primary antibodies: MBP (1:3,000; Millipore), PDGFRB (1:2,000; Proteintech), TLR9 (1:1,000; Servicebio), CCL5 (1:50; Huabio) and GAPDH (1:5,000; Proteintech). After incubation with HRP-conjugated secondary antibody, bands were visualized using ECL reagent on a Chemidoc detection system (Bio-Rad, Hercules, CA, USA). Three animals for each group were tested.

Quantitative RT-PCR

Total RNA of cerebral cortex tissue was extracted using Trizol reagent (Ambion, Austin, TX, USA), followed by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) for quantitative analysis. Reverse transcription was performed using the Fastking RT Kit (Tiangen, Beijing, China). KAPA SYBR FAST (KAPA Biosystems, Wilmington, MA, USA) was used to detect the mRNA expression. The following PCR conditions were used: initial denaturation at 95 °C for 20 s, followed by 40 cycles of 95 °C for 3 s and 60 °C for 20 s. The primer sequences were as follows: TLR9 (forward: 5′-TATCCACCACCTGCACAACT-3′, reverse: 5′-TTCAGCTCCTCCAGTGTACG-3′), CCL5 (forward: 5′-TGCCCACGTCAAGGAGTATTTC-3′, reverse: 5′-AACCCACTTCTTCTCTGGGTTG-3′), PDGFRB (forward: 5′-AGCTACATGGCCCCTTATGA-3′, reverse: 5′-GGATCCCAAAAGACCAGACA-3′), GAPDH (forward: 5′-AGGTCGGTGTGAACGGATTTG-3′, reverse: 5′-GGGGTCGTTGATGGCAACA-3′). The results were gathered and analyzed by Quant Studio 7 Flex (Applied Biosystems, Foster City, CA, USA). The relative levels of target genes were calculated by using the 2^−ΔΔCT method, and an internal reference, GAPDH, was utilized. Four animals for each group were tested.

Statistical analysis

Bioinformatic analyses were conducted using R software. Student’s t or Wilcoxon test was used to evaluate differences between two groups and was deemed significant at p < 0.05.

Identification and enrichment analysis of DEGs

The results of mRNA differential expression analysis showed a total of 502 DEGs, including 265 down-regulated genes and 237 up-regulated genes (Table S1). All these DEGs were displayed in the form of volcano plot (Fig. 1A). The heatmap depicts the top 20 up-and-down-regulated DEGs, with gene expression levels expressed in color, red for high expression, and blue for low expression (Fig. 1B). The results of the GO analysis showed that the enriched biological process (BP) terms were viral response, digestion, and synaptic vesicle cycles; the enriched cell composition (CC) terms were transport vesicles, synaptic membranes, and distal axon; the enriched molecular function (MF) terms were GTP enzyme binding, protein kinase regulator activity and sodium channel regulator activity (Fig. 1C). The top 10 enriched signaling pathways for the DEGs were determined by KEGG analysis (Fig. 1D), and DEGs identified were enriched in MAPK signaling pathways, neuroactive ligand-receptor interactions, and human tumor virus infection.

Screening modules associated with GM lesions by WGCNA

To uncover genes that significantly related to GM lesions, WGCNA was performed to identify potential functional modules. Genes with high expression correlations tend to be grouped into the same module, and each module is assigned a unique color, which implies that they are involved in a related biological pathway. Then, based on the scale-free R² = 0.9 and high average connectivity, we set the soft threshold to 6 (Fig. 2A). The identified modules were shown under a clustering tree (Fig. 2B), and the reliability of module delineation was elucidated by correlation analysis between the modules, which showed that there was no substantial association between them (Fig. 2C). Correlations between modules and clinical traits were analyzed. The results showed that turquoise (|R| = 0.84, p = 4e−06) and green modules (|R| = 0.78, p = 5e−05) showed a strong correlation with MS, and therefore, these two modules were considered as key modules for further analysis. Correlation analysis of gene significance vs. module membership showed that they were highly correlated in the turquoise (R = 0.81, p < 1e−200) and green modules (R = 0.87, p = 4.6e−135) (Fig. 2E).

Screening for immune-related genes

To explore immune-related gene modules, the PPI network was employed. We defined 2,660 immune-related genes from ImmPort and InnateDB (Table S2). In total 1,271 genes were identified after combining turquoise and green modules in WGCNA, and the intersection of these genes with immune genes was filtered (Fig. 3A).

We selected these 116 genes for LASSO analysis. With the increase of lambda value, the coefficient of some genes decreased to zero (Fig. 3B). Then tenfold cross-validation was carried out; when lambda = 0.018, the optimal model consisting of a total of 28 genes was developed (Fig. 3C).

PPI networks reflecting the association of the proteins encoded by LASSO screened genes were generated from the STRING database (Fig. 3D). Next, by taking the intersection of six algorithms of cytoHubba, a total of nine genes were selected, including CXCL8, CCL5, TLR9, MPO, KNG1, CRLF2, PDGFRB, SDC4 and GATA3 (Table 1). The significant module obtained from the PPI network using cytoHubba was shown in Fig. 3E.

Additionally, the expression values of these nine genes were verified using the data of 42 healthy control and 37 GM lesions tissue samples in the GSE131282 dataset, and eight of them (CXCL8, CCL5, TLR9, MPO, KNG1, CRLF2, PDGFRB and SDC4) were differentially expressed and exploited for further study (Fig. 3F).

Correlation analysis of immune-related genes and immune cells

Meanwhile, we calculated the infiltration of immune cells in MS samples, and the results showed that proportions of seven types of immune cells were significantly different between the two groups (Fig. 4A). There was an obvious increase in the proportion of naive CD4⁺ T cells, activated NK cells, M0 macrophages, M1 macrophages and eosinophils in MS group relative to that in normal group, while plasma cells and activated dendritic displayed the opposite proportion. The correlation of 22 types of immune cells was calculated (Fig. 4B). Naive B cells were significantly negatively correlated with CD4 memory resting T cells (R = −0.47, p = 0.040), but significantly positively correlated with resting NK cells (R = 0.64, p = 0.003), M1 macrophages (R = 0.51, p = 0.022) and resting dendritic cells (R = 0.62, p = 0.003). Plasma cells were significantly negatively correlated with activated NK cells (R = −0.56, p = 0.011), monocytes (R = −0.45, p = 0.044), M1 macrophages (R = −0.44, p = 0.050) and resting mast cells (R = −0.57, p = 0.008). CD8 T cells were significantly positively correlated with resting dendritic cells (R = 0.49, p = 0.029). CD4 memory resting T cells were significantly negatively correlated with CD4 memory activated T cells (R = −0.47, p = 0.038), regulatory Tregs T cells (R = −0.48, p = 0.032) and monocytes (R = −0.58, p = 0.007), but significantly positively correlated with M0 macrophages (R = 0.49, p = 0.027) and eosinophils (R = 0.48, p = 0.030). CD4 memory activated T cells were significantly negatively correlated with follicular helper T cells (R = −0.55, p = 0.011), but significantly positively correlated with resting NK cells (R = 0.64, p = 0.003), M1 macrophages (R = 0.51, p = 0.022) and resting dendritic cells (R = 0.62, p = 0.003). Follicular helper T cells were significantly negatively correlated with regulatory Tregs T cells (R = −0.57, p = 0.009) and resting NK cells (R = −0.54, p = 0.013). Regulatory Tregs T cells were significantly positively correlated with resting NK cells (R = 0.61, p = 0.004), M1 macrophages (R = 0.53, p = 0.017) and resting dendritic cells (R = 0.60, p = 0.005). Resting NK cells were significantly negatively correlated with activated NK cells (R = −0.44, p = 0.050), but significantly positively correlated with M1 macrophages (R = 0.53, p = 0.016) and resting dendritic cells (R = 0.48, p = 0.031). Activated NK cells were significantly positively correlated with resting mast cells (R = 0.51, p = 0.020). M0 macrophages were significantly negatively correlated with M2 macrophages (R = −0.46, p = 0.041) and activated dendritic cells (R = −0.61, p = 0.004). Resting mast cells were significantly positively correlated with eosinophils (R = 0.73, p < 0.001).

We further conducted correlation analysis to determine the association between the immune cells and the seven candidate genes (Fig. 4C). The results showed that CXCL8 displayed the most negative correlation with M1 macrophages (R = −0.67, p = 0.03) and the most positive correlation with activated dendritic cells (R = −0.62, p = 0.01), CCL5 displayed the most positive correlation with monocytes (R = 0.68, p = 0.03) and the most negative correlation with plasma cells (R = −.58, p = 0.03). TLR9 showed the most positive correlation with naive CD4⁺ T cells (R = 0.64, p = 0.01) and the most negative correlation with plasma cells (R = −0.51, p = 0.04), PDGFRB showed the most positive correlation with regulatory T cells (R = 0.74, p = 0.02) and the most negative correlation with plasma cells (R = −0.57, p = 0.04). Statistically significant correlations were not observed between the other four candidate genes (MPO, KNG1, CRLF2 and SDC4) and the 22 types of immune cells. Therefore, four genes (CXCL8, CCL5, TLR9 and PDGFRB) were finally determined to be immune-related hub genes.

Diagnostic model construction

Subsequently, we established a diagnosis nomogram model based on gene signatures (Fig. 5A). In order to detect the performance of the nomogram, we built a calibration curve for visualization (Fig. 5B); the nomogram prediction model displayed a C-index of 0.854 which revealed high consistency between the actual and predicted survival rates. According to the DCA curve (Fig. 5C), the nomogram prediction model has good performance in predicting the incidence of MS. Subsequently, we conducted ROC curve analysis in the third-party dataset GSE131282. The AUC values of hub genes were CCL5: 0.795, CXCL8: 0.614, PDGFRB: 0.815, TLR9: 0.759 (Fig. 5D). Hence, the nomogram based on a four-gene signature is a good model for MS diagnosis.

Analysis of GSEA pathways of immune hub genes in MS

To further explore the hub gene-related signaling pathway, the GSEA method was used. GSEA enrichment analysis showed a dominant enrichment of CCL5 in the pathways of adhesion molecules, regulation of actin cytoskeleton and cytokine-cytokine receptor interactions (Fig. 6A). CXCL8 was prominently enriched in the pathways of cell adhesion molecules, human T-cell leukemia virus 1 infection and Phagosome (Fig. 6B). PDGFRB exhibited a significant enrichment in the pathways of cytokine-cytokine receptor interaction, focal adhesion and JAK-STAT signaling pathway (Fig. 6C). TLR9 was significantly enriched in the pathways of Parkinson disease, ribosome, and spinocerebellar ataxia (Fig. 6D). Genes significantly enriched in hub gene-related pathways were shown in Fig. 7.

Verification of immune hub genes in MS mouse model

To analyze the myelin loss in the cerebral cortex of CPZ-fed mice, LFB staining for myelin sheath and IHC staining for MBP (a marker of mature myelin) were performed (Fig. 8A). We observed a strong myelin loss in mice cortical tissues after 5 weeks of treatment with CPZ. We further detected the protein level of MBP in cortical tissues using western blot, and in line with the IHC result, MBP level was significantly downregulated in the CPZ group, indicating that CPZ successfully induced cortical demyelination. As rodents lack the CXCL8 gene, we verified the protein and mRNA expressions of TLR9, PDGFRB and CCL5 in cuprizone-induced MS mice by IHC, western blot and RT-qPCR (Figs. 8B–8D). The results showed that compared with control group, the protein and mRNA expressions of TLR9 and PDGFRB were decreased, while the protein and mRNA expression of CCL5 was increased in CPZ group.