Novel peripheral blood mononuclear cell mRNA signature for IFN-beta therapy responsiveness prediction in multiple sclerosis

Abstract

Interferon-beta (IFN-$\beta$) is one of the classical drugs for immunomodulatory therapy in relapsing-remitting multiple sclerosis (RRMS) patients, but the drug responsiveness of different patients varies. Currently, there is no valid model to predict IFN-$\beta$ responsiveness. This research attempted to develop an IFN-$\beta$ responsiveness prediction model based on mRNA expression in RRMS patient peripheral blood mononuclear cells. Peripheral blood mononuclear cell mRNA expression datasets including 50 RRMS patients receiving IFN-$\beta$ treatment were obtained from GEO. Among the datasets, 24 cases from GSE24427 were included in a training set, and 18 and 9 cases from GSE19285 and GSE33464, respectively, were adopted as two independent test sets. In the training set, blood samples were collected immediately before first, second, month 1, 12, and 24 IFN-$\beta$ injection, and the mRNA expression data at four time points, namely, two days, one month, one year and two years after the onset of IFN-$\beta$ treatment, were compared with pre-treatment data to identify IFN-stimulated genes (ISGs). The ISGs at the one-month time point were used to construct the drug responsiveness prediction model. Next, the drug responsiveness model was verified in the two independent test sets to examine the performance of the model in predicting drug responsiveness. Finally, we used CIBERSORTx to estimate the content of cell subtypes in samples and evaluated whether differences in the proportions of cell subtypes were related to differences in IFN-$\beta$ responsiveness. Among the four time points, one month was the time point when the training set GSE24427 and test set GSE33464 had the highest number of ISGs. Functional analysis showed that these one-month ISGs were enriched in biological functions such as the innate immune response, type-I interferon signalling pathway, and other IFN-$\beta$-associated functions. Based on these ISGs, we obtained a four-factor prediction model for IFN-$\beta$ responsiveness including MX1, MX2, XAF1, and LAMP3. In addition, the model demonstrated favourable predictive performance within the training set and two external test sets. A higher proportion of activated NK cells and lower naive CD4/total CD4 ratio might indicate better drug responsiveness. This research developed a polygene-based biomarker model that could predict RRMS patient IFN-$\beta$ responsiveness in the early treatment period. This model could probably help doctors screen out patients who would not benefit from IFN-$\beta$ treatment early and determine whether a current treatment plan should be continued.

Graphical Abstract

Keywords:

• Multiple sclerosis
• IFN-$\beta$ therapy
• drug responsiveness
• prognostic model construction
• survival analysis

1. Introduction

Relapsing-remitting multiple sclerosis (RRMS) is a common autoimmune disease of the central nervous system, and its clinical manifestations alternate between relapse and remittance [1]. Since the neural functional disorder caused by disease relapse each time is partially irreversible, patients might experience serious neural functional disorders after suffering from relapse many times [2]. Therefore, a main focus of MS treatment is reducing the number of relapses and extending the relapse interval in patients. As an immunomodulatory therapy used to prevent relapse in RRMS, interferon (IFN)-$\beta$ has a history of more than 30 years, and was once used as a first-line treatment option, but is now not the most widely used drug because it is not as effective, on average, as other disease-modifying treatments, such as newer forms of IFN, monoclonal antibodies, and oral medications [3,4]. Responsiveness to this drug varies significantly among patients, and it has been reported that approximately 30% to 50% of patients may not respond to the drug [5,6]. The mechanism underlying poor drug responsiveness has not yet been clarified but is probably related to neutralising antibodies (NAbs) [7,8]. However, depending on factors such as the specific MS subtype, severity of disease, side effect profile, and therapeutic goals, IFN-$\beta$ may still be preferred in some cases, especially for women considering pregnancy, breastfeeding, and older patients [9–11]. Currently, patients with a poor treatment response cannot be screened out before therapy or in the early therapy phase. Therefore, for a long time, researchers have been attempting to identify biomarkers that can predict the responsiveness of RRMS patients to IFN-$\beta$ treatment and guide individual patients’ adoption of IFN-$\beta$ treatment [12–15].

In recent years, an increasing number of post-injection predictive biomarkers in predicting the IFN treatment response in the plasma or cerebrospinal fluid have been discovered and reported, such as NAbs, myxovirus resistance protein (Mx) A or Mx1, and NLRP3 [16–20]. However, the specificity of these biomarkers in predicting the treatment response of individual patients is limited. The construction of prognostic prediction models by integrating multiple biomarkers has been shown to be effective in many diseases [21–23]. At the same time, the timing of biomarker assessment is also a key but unclear factor [24,25]. Although some studies found that the peak mRNA levels of MxA and other IFN-$\beta$ therapy-related genes in the plasma occurred within 6-12 h after IFN-$\beta$ injection, the classification of MS patients based on the MxA expression level at this time point could not effectively predict the risk of recurrence after IFN-$\beta$ therapy [26,27]. We hypothesised that IFN-$\beta$ may have some long-term pharmacological effects on the body in addition to this short-term response, so more time points need to be tested. This research evaluated mononuclear cell gene expression in peripheral blood before and after IFN-$\beta$ treatment of RRMS patients. By comparing the mRNA expression data at different time points following IFN-$\beta$ treatment with data collected before treatment, we identified differentially genes, namely, IFN- stimulated genes (ISGs), at different time points after treatment administration to define the optimal detection window after a drug response. Then, the ISGs in the optimal time window were used to build a drug response prediction model, which was verified with two mutually independent and external datasets. The flow chart is shown in Figure 1.

Figure 1. Flowchart of the development of an IFN-$\beta$ responsiveness prediction model based on the mRNA expression level of RRMS patients’peripheral blood mononuclear cells. DEGs, differential expressed genes.

2. Materials and methods

2.1. Selection and collection of peripheral blood mononuclear cell expression datasets from RRMS patients

The GEO database (https://www.ncbi.nlm.nih.gov/gds) was used for gene chip collection, and relevant chip datasets, including GSE24427 [28], GSE19285 [29], and GSE33464 [30], were downloaded under the keyword “relapsing-remitting multiple sclerosis”. All three GSE studies were from same lab. The peripheral blood mononuclear cell gene expression chip data for RRMS patients aged 23 years old to 63 years old who at least regularly received IFN-$\beta$ treatment for two years, or received IFN-$\beta$ treatment for one year and had not relapsed within the following year, and finished two years of follow-up visit were chosen from the aforementioned datasets. These datasets included expression data collected before and one month after the onset of IFN-$\beta$ treatment and contained clinical data for relapse during the two follow-up years and Expanded Disability Status Scale (EDSS) score changes.

We included 24 peripheral blood mononuclear cell specimens from RRMS patients receiving IFN-$\beta$ in GSE24427. Among these patients, 16 were females, and 8 were males; their average age was 40 years old. Gene expression data for five time points including before treatment and 2 days, 1 month, 1 year, and 2 years after the onset of IFN-$\beta$ treatment were chosen. We included 18 peripheral blood mononuclear cell specimens from RRMS patients receiving IFN-$\beta$ treatment in GSE19285. Among these patients, 13 were females, and 5 were males; their average age was 37 years old. The gene expression data collected before treatment and 1 month after the onset of treatment were chosen (no data for one year and two years after administration). We included 9 peripheral blood mononuclear cell specimens from RRMS patients receiving IFN-$\beta$ treatment in GSE33464. Among these patients, 5 were females, and 4 were males; their average age was 39 years old. The gene expression data for four time points including before treatment and 1 month, 1 year, and 2 years after the treatment were chosen. Detailed clinical information of the RRMS patients in the training and test sets is presented in Table S1.

2.2. Screening of ISGs and gene ontology (GO) functional enrichment and Kyoto encyclopedia of genes and gnomes (KEGG) signalling pathway analyses

GSE24427 was adopted as the training set. Using the limma package [31], we compared the gene expression data of patients at 2 days, 1 month, 1 year, and 2 years after IFN-$\beta$ treatment, respectively, with the pre-treatment data to obtain ISG lists at 4 different time points after treatment. The screening standard was p < 0.05 and fold change (FC) −2 to 2. STRING was used to establish the network of interactions among the ISG proteins. Enrichment of ISG gene clusters and signalling pathways was realised through DAVID (https://david.ncifcrf.gov/home.jsp) [32]. Specifically, under the background of anthropogenic genes, GO functional enrichment and KEGG signalling pathways were analysed.

2.3. Screening of key genes related to the response of RRMS patients to IFN-$\beta$ treatment

Twenty-four RRMS patients in GSE24427 were defined as the training set. Seventeen MS patients in GSE19285 and 9 MS patients in GSE33464 were defined as test set 1 and test set 2, respectively. The training set was used to construct the IFN-$\beta$ treatment response prediction model. The test set was used to verify the predictive performance of the response prediction model. The ISGs of RRMS patients were combined with the adverse response occurrence data and clinical rating changes of 24 patients. Patients who had adverse event, including relapse event, an EDSS score increase of $\ge 1$ (if baseline EDSS $\le 5.5$) or $\le 0.5$ (if baseline EDSS $>5.5$) [33] within two years of IFN-$\beta$ treatment were defined as poor response to treatment, otherwise defined as good response to treatment. The R software “survival” function package was adopted for single-factor Cox analysis and to screen out genes related to patients’ treatment responses. Next, the results obtained by single-factor Cox analysis were included in a LASSO regression analysis. LASSO regression, through construction of the regression coefficient of the penalty function towards independent variables, thus produced a stronger correlation with the RRMS patients’ IFN-$\beta$ treatment response data. Next, multifactor Cox analysis was conducted on the genes obtained through the LASSO regression analysis, and the key genes related to patients’ treatment responses were identified.

2.4. Construction of the IFN-$\beta$ response prediction model for RRMS patients in the training set

The regression coefficients of key factors related to treatment response were identified through LASSO and Cox regression analyses. In this way, the IFN-$\beta$ treatment response prediction model for the training set could be built, in which the response value equals $\beta$1Exp1 +$\beta$2Exp2+……$\beta$ nExpn. $$\begin{array}{l}\hspace{0.17em}\text{Model: Risk\hspace{0.33em}score}\hspace{0.17em}={\displaystyle \sum }_{i=1}^k\beta \text{i\hspace{0.33em}Expi}\hfill \end{array}$$

Exp denotes the expression of a gene, “$\beta$” denotes the regression coefficient in the multifactor Cox regression analysis of genes, and “n” denotes the gene dosage related to patient prognosis. The response value of each case in the training set could be computed with this formula. Then, the median of the response value was set as the cut-off value. According to the cut-off value, patients in the training set were divided into good and poor responsiveness groups. The “survival” package of R software was used to carry out exacerbation- and progression-free survival analysis of the good and poor responsiveness group. Additionally, the “Survival ROC” package was used to draw the receiver operating characteristic (ROC) curve and obtain an assessment of the predictive performance of the response prediction model.

2.5. Single-factor and multifactor prognostic analyses

Patient clinical information such as age, sex, baseline EDSS score before treatment, and relapse times before treatment were included. Then, the response rating was combined for single-factor and multi-factor Cox analyses to assess the impacts of clinical factors on the response rating.

2.6. Verification of the IFN-$\beta$ response prediction model for RRMS patients with the test sets

To test the accuracy of the treatment response prediction model, the risk values of 24 cases were calculated with the response prediction model built on the training set. The median of the risk value was adopted as the cut-off value, and the patients in the test set were divided into good and poor responsiveness groups. The R software “survival” package was adopted for Kaplan-Meier survival analysis of the high- and low-risk groups. The “survival ROC” package was used to draw the ROC curve to assess the performance of the response prediction model.

2.7. Subtype analysis of the proportions of mononuclear cells

To characterise the abundances of mononuclear cell subtypes based on RNA-seq data for blood, the CIBERSORTx (https://cibersortx.stanford.edu/) [34] web tool was applied. Using a deconvolution algorithm, CIBERSORTx computed the levels of different cell types. CIBERSORTx derived an empirical p-value for the deconvolution of each case using Monte Carlo sampling, and samples with p < 0.05 were adopted for analysis because of the high reliability of the inferred cell composition. Therefore, cases with a p-value $\ge 0.05$ were not retained for subsequent analysis. To validate the accuracy of the CIBERSORTx data, the Tumour Immune Estimation Resource (TIMER) database was also employed to calculate the abundances of immune cells. Subsequently, box plots were utilised to present the differences in infiltrating immune cells, T cell activating and inhibitory receptors, and macrophage-associated molecules between the high- and low-risk groups using the “ggplot2” package.

2.8. Statistical analysis

All statistical analyses were conducted using R (version 4.0.3). The Wilcoxon test was used to screen significantly differentially expressed genes and infiltrating immune cells. Kaplan-Meier curves and the log-rank test were used to evaluate the statistical significance of differences in survival rates between risk groups. The predictive accuracy of risk signatures was determined with ROC curves. The proportional hazards assumption was tested with Schoenfeld residuals. Then, univariate and multivariate Cox regression analyses were performed to evaluate significant prognostic factors. Finally, the results of multivariate Cox regression analyses were visualised with a nomogram. A p-value < 0.05 was considered statistically significant.

3. Results

3.1. The best time point for prediction of the IFN-$\beta$ treatment response in RRMS patients

RRMS patients in the training set, GSE24427, had the largest number of ISGs after being treated with IFN-$\beta$ for one month. ISGs were evaluated after treatment with IFN-$\beta$ at four time points, which identified 0 (2 days), 53 (1 month), 0 (1 year), and 0 (2 years) ISGs at the different time points (Figure 2). PCA showing the percentage variability captured by the first two principal components is displayed across PC1 and 2. The analysis plot also revealed a distinct separation before and after IFN-$\beta$ treatment for 1 month. In the test set GSE33464, the similar trend were detected (Figure S1). The GSE19285 dataset was not included in this analysis because it does not have data collected at two days, one year or two years after treatment. The above results indicated that one month after patients started receiving treatment seemed to be a good window to test the IFN-$\beta$ treatment response.

Figure 2. Gene expression at different time points in RRMS patients’ peripheral blood mononuclear cells following IFN-$\beta$ treatment in the training set. Volcano maps and PCA showed that several genes altered significantly in patients at one month of treatment.

3.2. Protein-protein interaction network of ISGs and gene function analysis

Fifty-three ISGs were detected at one month after the onset of IFN-$\beta$ therapy (Table S2). Functional analysis suggested that ISGs detected at one month after treatment were mainly enriched in the cellular response to type I interferon, response to virus, response to other organisms, response to biotic stimulus, etc. (Figure 3, Table S3).

Figure 3. Protein-protein interaction network and functional analysis of ISGs detected at one month after the onset of IFN-$\beta$ therapy. (a) The network of interactions among the ISG proteins. Each node represents a protein. Disconnected nodes are not shown. (b) GO functional (biological process and molecular function) and KEGG signalling pathway enrichment analysis of ISGs. The bubble size and colour correlate with the number of ISGs and the significance of the p-value, respectively.

3.3. Construction of the IFN-$\beta$ treatment response prediction model with the data for RRMS patients in the training set

Single-factor Cox regression analysis of the training set indicated that 13 ISGs induced by IFN-$\beta$ were related to prognosis after IFN-$\beta$ treatment. The 13 ISGs were CMPK2, HERC6, IFIT3, LAMP3, MX1, MX2, OASL, OTOF, RSAD2, SCO2, SERPING1, SIGLEC1, and XAF1 (Figure 4a). Then, LASSO and multifactor Cox regression analyses revealed that three genes, namely, MX1, MX2, and XAF1, were key genes related to the IFN-$\beta$ treatment response (p < 0.05) and that these three genes could be adopted as independent prognostic factors. Increasing expression of MX1 was an indicator of a poor response, whereas increasing expression of MX2 and XAF1 was an indicator of a good response. In addition, the multifactor analysis yielded a p-value of 0.066 for LAMP3, which we included in the rating formula to improve the applicability of the formula. Finally, we constructed a predictive model for IFN-$\beta$ treatment in RRMS patients consisting of four key prognosis-related genes using LASSO and Cox regression models. The risk score was written as “(15.50879$\times$ MX1) + (-8.94$\times$ MX2) + (-9.49369$\times$ LAMP3) + (-10.6426$\times$ XAF1)” (Figure 4b, Table S4).

Figure 4. The expression level of ISGs associated with the prognosis of RRMS patients. (a) Single-factor Cox regression analysis. (b) Multifactor Cox regression analyses. The four ISGs obtained by constructing the prognostic model shown in the forest map, red indicates high-risk gene and green denotes low-risk genes.

3.4. Verification of the RRMS patient IFN-$\beta$ treatment response prediction model with the training set and test sets

The risk score of each patient could be calculated with the formula, and the median of the risk score was adopted as the cut-off value, with the risk score based on the onset or progression of adverse events. The 24 cases in the training set, GSE24427, were divided into a high-risk group (12 cases) and low-risk group (12 cases). The clinical information of high-risk and low-risk group were shown in Table S5. The risk score and adverse events-free time corresponding to each sample and the expression value of each gene are shown in Figure 5a. Kaplan-Meier survival curve analysis demonstrated that the median time of adverse events in the high-risk group was significantly earlier than that in the low-risk group (p = 0.024) (Figure 5b). ROC curve analysis of the predictive performance for adverse events within 2 years of the responsiveness prediction model showed that the area under the curve for adverse event in MS patients within 2 years was 70.80% (Figure 5c), indicating that the model could be used to predict responsiveness to IFN-$\beta$ treatment.

Figure 5. Adverse events-free survival analysis and verification of IFN-$\beta$ treatment response prediction model. (a) The risk score and adverse events-free status distribution of training set (GSE24427) and test sets (GSE19285 and GSE33464). (b) Adverse events-free survival analysis curves of patients in the high- and low-risk group in training set and test sets. (c) ROC curves of the predictive performance for exacerbation and progression within 2 years of the responsiveness prediction model of training set and test sets.

Next, we verified the prediction model in the test sets. The 18 cases in test set 1, GSE19285, were divided into a high-risk group (8 cases) and low-risk group (10 cases). The 9 cases in test set 2, GSE33464, were divided into a high-risk group (5 cases) and low-risk group (4 cases). The median time of adverse events in the high-risk group was significantly earlier than that in the low-risk group in both sets (p = 0.002, p = 0.032) (Figure 5b). ROC curve analysis showed that the area under the curve for adverse events in MS patients within 2 years was 88.75% and 75%, respectively (Figure 5c). The model was proved to be effective in predicting the response to IFN-$\beta$ treatment.

3.5. Single-factor and multifactor prognostic analyses

Considering that age, gender, baseline EDSS score and relapse times before treatment might also affect patients’ response to IFN-$\beta$ treatment, this research included the above variables in the IFN-$\beta$ treatment response prediction model described above. The proposed response prediction model could still predict the response of patients to IFN-$\beta$ treatment. This finding was verified in single-factor and multifactor Cox models (p < 0.01) (Figure 6).

Figure 6. Single-factor (a) and multifactor (b) prognostic analyses of other demographic or clinical data of the patients.

3.6. Correlations between the IFN-$\beta$ response score and proportions of mononuclear cell subtypes

We combined the training set and test set data, used the classification of the high- and low-risk samples in the previous survival analysis and divided the 50 samples into a high-risk group containing 25 cases and a low-risk group containing 25 cases. As shown by the violin plot in Figure 7a, the abundances of mononuclear subtypes determined by using CIBERSORTx did not show much difference between the high-risk and low-risk group. There might have been a tendency for the proportion of activated NK cells in the high-risk group to be lower than that in the low-risk group (p = 0.05). However, it is worth noting that the proportion of naive CD4 T cells in the total CD4 T cell population was different between the high- and low-risk groups. Given the background that the overall number of CD4 T cells remained almost constant, the proportion of naive CD4 T cells in the total CD4 T cell population in the high-risk group was greater than that in the low-risk group (p = 0.018) (Figure 7).

Figure 7. Correlations between the IFN-$\beta$ response score and proportions of mononuclear cell subtypes. (a) Violin plot of the abundances of mononuclear subtypes between the high-risk (blue) and low-risk (red) groups. (b) The abundances of naive CD4 T cells in the total CD4 T cell population between the high- and low-risk groups.

4. Discussion

There are actually multiple disease-modifying therapies available today for the treatment of relapsing MS, however the development of a predictive model of IFN-$\beta$ reactivity could make the drug still an important first-line treatment option. In this study, we carried out a predictive study of the response to IFN-$\beta$ in a group of RRMS patients by comparing gene expression in peripheral blood mononuclear cells at different post-treatment times with respect to the time prior to the start of treatment. Several genes more altered in patients at one month of treatment, with no change at two days, 1 year, and 2 years post-treatment. Previous studies have described changes in gene expression in MS patients following interferon therapy. Serrano et al. identified differentially expressed genes at four time points during a 1-year period in MS patients receiving IFN-$\beta$ therapy [35]. Samples were collected from patients receiving IFN-$\beta$ treatment before, 2 days, 1 month, and 1 year after treatment. The 15 IFN response genes, including MX1 and XAF1, were differentially expressed at all time points, and the gene expression of most patients reached the highest at 1 month after treatment. Another study examined gene expression in peripheral blood mononuclear cells of seven MS patients at baseline as well as 3 months after IFN-$\beta$-1a treatment, the gene expression changes became less pronounced, suggesting a steady-state response to IFN-$\beta$ [36]. The above results are consistent with our findings. In addition, by studying the most altered genes, we made a predicitive analysis that allows us to identify patients with high and low risk of not being good responders to IFN-$\beta$. The scaling of high and low risk is not observed with other demographic or clinical data of the patients.

We found that one month after the onset of IFN-$\beta$ treatment seemed to be a good time point to identify early-stage prediction indexes for IFN-$\beta$ treatment responsiveness. The experiment revealed that the largest number of differentially expressed genes appeared one month after the onset of IFN-$\beta$ treatment instead of 2 days. The same phenomena could also be observed in the test set (GSE33464). These probably indicate that the immune processes involved in ISGs in peripheral blood mononuclear cells at one month of IFN-$\beta$ treatment are more complex. At one year and two years, we found that the number of ISGs was reduced again, which probably indicated that the ISGs originated from temporary physiological processes, or there is a feedback mechanism that ensures homeostasis. As a pleiotropic cell factor, IFN-$\beta$ can affect relevant molecular models of multiple cells in response to external pathogenic stimulus reactions, such as various Toll like receptors, and affect some biological processes, such as immune cell autophagy [37,38]. In this study, the ISGs identified at the one-month time point were enriched in functions including double-stranded RNA binding, the antiviral response, and 2′-5′-oligoadenylate synthetase activity.

The Mx protein is a specific protein generated during the IFN-$\beta$ antiviral reaction that can inhibit virus transcription and is often regarded as a reaction-based biological indicator to assess IFN-$\beta$ biological functions or patient responsiveness [39]. Both MX1 and MX2 are IFN-inducible proteins but have different antiviral mechanisms and subcellular localisations, especially in autoimmune diseases such as MS. MX1 is mainly localised in the nucleus and associated with the nuclear membrane, where it exerts its antiviral activity by blocking the entry of viral ribonucleoproteins into the nucleus. MX2 is mainly located in the cytoplasm and associated with the intracellular membrane, where it exerts its antiviral activity by preventing viral ribonucleoproteins from entering the nucleus [40]. In MS, immune processes occur mainly in the extra-nuclear, cytoplasmic and extra-cellular spaces within the central nervous system. The primary goal of MS treatment is not to directly target viruses but to modulate the immune system’s inflammatory response. Therefore, MX2, which is localised in the cytoplasm, may be more important in the treatment of MS.

It is worth noting that the biological effects of Mx can be blocked by anti-IFN-$\beta$ NAb, resulting in less effect [41]. NAbs are one of the most recognised biological markers for the therapeutic effect of IFN-$\beta$ treatment on RRMS and are related to the resistance of patients to IFN-$\beta$ treatment. However, there are some problems in the clinical application of the Mx protein and NAb. One of the problems is that the baseline titre of NAb varies greatly among patients [42]. In this way, the multiple-mRNA model we built in this study may help alleviate this problem. We comprehensively evaluated the mRNA transcripts of several genes (MX1, MX2, XAF1, and LAMP3) to improve the consistency of the predictive model when applied to different individuals. Furthermore, this study observed that low-level MX1 and high-level MX2 were positively correlated with patients’ drug responsiveness, which showed good agreement with previous studies on the correlations between MXs and IFN-$\beta$ treatment [43]. However, the predictive model of IFN-$\beta$ treatment responsiveness proposed in this study involves the cooperative interaction of several genes, but may not be fully applicable to accurately predict a single gene for treatment responsiveness.

X-chromosome-linked inhibitor of apoptosis proteins associated factor 1 (XAF1) is a kind of apoptosis-promoting tumour inhibitor that has been extensively investigated in tumour-related disease research [44]. Moreover, it has been reported that XAF1 can form a positive feedback link with IFN regulating factor-1 (IRF-1) to strengthen the role of IFN in vitro [45]. In this research, we found that patients with good IFN-$\beta$ responsiveness had a higher level of XAF1. Accordingly, this paper speculated that the stronger XAF1 expression and mutual positive feedback of IRF-1 might be reasons for good IFN-$\beta$ responsiveness. In addition, it is important to note that XAF-1 is an X chromosome-associated gene that may be copied more in females, and the role of XAF1 may be potentially related to the male/female ratio in each group. However, 70%-80% of X-linked genes are inactivated in females, maintaining similar levels of X-linked protein expression between the sexes [46,47]. Currently, we have found no relevant reports of gender differences in XAF1 [48], but studies with larger samples are still needed to rule out potential gender effects. Lysosomal membrane protein 3 (LAMP3) is one of the membrane proteins related to lysosomes. LAMP1, LAMP2, and LAMP3 together are known as the products of cell reactions towards IFN-$\beta$ [49]. LAMP3 can activate cell apoptosis in response to autoimmune diseases, such as sicca syndrome, and is related to its autoantigen [50]. However, the roles of LAMP3 in RRMS disease occurrence and development have not yet been reported. This research observed a positive correlation between an increase in LAMP3 levels and IFN-$\beta$ treatment responsiveness, which indicated that enhanced expression of LAMP3 could benefit RRMS patients’ IFN-$\beta$ treatment responsiveness.

We also tried to incorporate clinical factors into the IFN-$\beta$ treatment response prediction model, including gender, age, baseline disease severity degree, and relapse times before treatment. We found that none of these factors showed predictive power for IFN-$\beta$ responsiveness. It is widely recognised that gender plays a role in the incidence and course of MS, but there is a lack of reports on gender differences in responsiveness to IFN-$\beta$ therapy and Nab production in MS patients [51]. However, in this research, although males showed a trend towards a lower risk of adverse events, it was not statistically significant. We believe this trend may not be instructive and may be related to the small sample size. Future studies with larger samples may help to deepen the understanding of gender differences in IFN-$\beta$ therapy.

In addition, this study found that at one month after the initiation of IFN-$\beta$ therapy, there were differences in the proportions of mononuclear cell subtypes between the high and low IFN-$\beta$ responsiveness prediction score groups. A higher proportion of activated NK cells and a lower naive CD4/total CD4 ratio showed correlations with better drug responsiveness. Since the proportion of total CD4 T cells did not change significantly, we speculated that IFN-$\beta$ might promote the transformation of naive T cells into functional T cells without affecting the total number of CD4 T cells. The existence of this correlation indicated that the genes used in the responsiveness prediction model for IFN-β treatment might be involved in the changes in lymphocytes that occur after IFN-$\beta$ treatment. Additionally, LAMP3 is closely related to dendritic cells and antigen-presenting functions [52], which were not included in this study. More cell subtype studies are necessary in the future to predict IFN-$\beta$ responsiveness.

Although IFN-$\beta$ is no longer the most widely used treatment for MS because it is not as effective on average as some of the other disease-modifying drugs, a subset of patients identified through the current approach could be very high responders. This study shows that the levels of several genes induced by IFN-$\beta$ may predict future responses to IFN-$\beta$ therapy in RRMS patients and may enhance the selection of drug responders. This topic has been studied extensively before, and the main originality of our study lies in the use of other longer experimental times. All three GSE studies were from same lab, the arrays may be identical and thus more reproducible and compatible. However, this study has several limitations: (a) the sample size of this study is limited, especially the limited sample of test set 2, so the prognostic model needs more data to be validated; (b) the data of this study are from public databases, so further experiments are still needed to investigate the clinical value of the model genes.

Author contributions

YX designed the research and collected data. YX and PY led the statistical analysis and imaging data process, contributing equally to the work. All authors contributed to interpretation of the data, preparation of the manuscript, and read and approved the final manuscript.

Abbreviations
RRMS	=	Relapsing-remitting multiple sclerosis
IFN	=	interferon
NAbs	=	neutralising antibodies
Mx	=	myxovirus resistance protein
ISGs	=	IFN-stimulated genes
DEGs	=	differentially expressed genes
EDSS	=	Expanded Disability Status Scale
FC	=	fold change
GO	=	Gene Ontology
KEGG	=	Kyoto Encyclopaedia of Genes and Genomes
TIMER	=	the Tumour Immune Estimation Resource
ROC	=	receiver operating characteristic
LAMP3	=	Lysosomal membrane protein 3
XAF1	=	X-chromosome-linked inhibitor of apoptosis proteins associated factor 1
IRF-1	=	IFN regulatory factor-1

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Heilongjiang Postdoctoral Financial Assistance (No. LBH-Z21176) and Research and Innovation Fund of the First Affiliated Hospital of Harbin Medical University (2021203).