Hypoxia–Immune-Related Gene SLC19A1 Serves as a Potential Biomarker for Prognosis in Multiple Myeloma

Background: Multiple myeloma (MM) remains an incurable malignant tumor of plasma cells. Increasing evidence has reported that hypoxia and immune status contribute to the progression of MM. In this research, the prognostic value of the hypoxia–immune-related gene SLC19A1 in MM was evaluated by bioinformatics analysis.

Method: RNA-sequencing (RNA-seq) data along with clinical information on MM were downloaded from the Gene Expression Omnibus (GEO) database. Consistent clustering analysis and ESTIMATE algorithms were performed to establish the MM sample subgroups related to hypoxia and immune status, respectively, based on the GSE24080 dataset. The differentially expressed analysis was performed to identify the hypoxia–immune-related genes. Subsequently, a hypoxia–immune-gene risk signature for MM patients was constructed by univariate and multivariate Cox regression analyses, which was also verified in the GSE4581 dataset. Furthermore, the mRNA expression of SLC19A1 was determined using qRT-PCR in 19 MM patients, and the correlations between the genetic expression of SLC19A1 and clinical features were further analyzed.

Result: A total of 47 genes were identified as hypoxia–immune-related genes for MM. Among these genes, SLC19A1 was screened to construct a risk score model that had better predictive power for MM. The constructed prognostic signature based on SLC19A1 was verified in the GSE4581 dataset. All independent prognostic factors (age, β₂-microglobulin, LDH, albumin, MRI, and gene risk score) were used to develop a nomogram that showed a better performance for predicting the survival probability of MM patients for 1–5 years. Furthermore, SLC19A1 was highly expressed in newly diagnosed and relapsed MM patients, and high expression of SLC19A1 is correlated with higher bone marrow aspiration plasma cells and β₂-microglobulin levels in MM patients.

Conclusion: In conclusion, our results suggest that SLC19A1 is aberrantly expressed in MM and highly expressed SLC19A1 might be a biomarker correlated with inferior prognosis. More importantly, we identified SLC19A1 as a hypoxia–immune-related gene in MM. Future functional and mechanistic studies will further clarify the roles of SLC19A1 in MM.

Introduction

Multiple myeloma (MM) is a plasma cell malignancy and remains incurable (1, 2). The main clinical manifestations of MM are renal insufficiency, hypercalcemia, bone marrow hematopoietic suppression, osteolytic bone lesions, and bone pain (3, 4). Diagnostic and prognostic biomarkers are used for predicting clinical outcomes and helping therapy selection and optimization so that MM patients can benefit from specific therapies (5, 6). Precise medical treatment of MM will require appropriate supporting diagnostic tests to detect the targetable lesions (5, 7); however, assay cost and availability remain major obstacles. This emphasizes the need to identify new molecules to provide prognostic biomarkers and/or therapeutic targets.

Hypoxia is one of the hallmarks of cancer, which is significantly associated with the occurrence and development of cancer (8). Hypoxia/oxygen sensing signal plays a vital role in the regulation of tumor progression (9–11). Previous studies have proven that the demethylases, KDM6A, JMJD3, JMJD6, and LSD1, have an affinity for oxygen, indicating their important significance in oxygen sensing and hypoxic reprogramming (12–14). In a variety of cancer models, the adaptation of tumor cells to the imbalance of oxygen demand and supply is related to inferior clinical outcomes (9, 15, 16). Through regulating gene and protein expression, hypoxia mediates genetic instability, malignant progression, and resistance to standard therapies (17).

In physiological immunological niches, such as the bone marrow, physiological hypoxia controls adaptive and innate immunity through transcriptional changes driven by the hypoxia-inducible factor (HIF) (18–20). On the contrary, in pathological immunological niches, such as tumors, pathological hypoxia can drive disease development and tissue dysfunction, via leading to immune cell imbalance (18). Through regulating the immunosuppressive effect, hypoxia also plays a vital role in anticancer immune responses (17). HIF signaling in liver cancer cells and innate immune cells together is conducive to the recruitment and maintenance of protumorigenic immune cells, inhibiting antitumorigenic immune cells and promoting immune evasion (21).

The reduced folate carrier (RFC, SLC19A1), thiamine transporter-1 (THTR1, SLC19A2), and thiamine transporter-2 (THTR2, SLC19A3) evolved from the same solute carrier family. SLC19A1 transports folates and SLC19A3 mediates intestinal thiamine absorption. SLC19A1 and SLC19A2 deliver their substrates to systemic tissues. SLC19A2 and SLC19A3 transport thiamine (22). THTR2 (SLC19A3) mRNA levels are downregulated in breast cancer, and the sensitivity to adriamycin was increased in the breast cancer cell line when transfected with SLC19A3 (23, 24). The SLC19A3-transfected cells were more sensitive to doxorubicin but not to paclitaxel or methotrexate (23). SLC19A1 is the primary means of cellular uptake for antifolate chemotherapeutic drugs;. therefore, the antifolate membrane transport of SLC19A1 is considered to be an important factor affecting antitumor activity (25). The reduced SLC19A1 plays a vital role in the transport of 5-methyltetrahydrofolate into the cells, and its lower expression has been related to the drug resistance of folate antagonists (26). In addition, SLC19A1 is the major transporter of cyclic dinucleotides (CDNs) into cells and has implications for the immunotherapeutic treatment of cancer (27–29). However, SLC19A1 has been identified but not yet studied for its impact on prognostic prediction or function in MM. In the present research, we tried to identify a group of prognostic genes related to hypoxia and immunity and construct a risk score model. Furthermore, the hypoxia–immune-related risk groups and several vital clinical parameters were included in a nomogram, which contributes to risk stratification and prognosis. The overall design of this research is shown in Figure 1.

FIGURE 1

Figure 1 Overall design to screen and validate the hypoxia–immune-related gene SLC19A1 for multiple myeloma.

Methods and Materials

Data Collection

In our study, the gene expression profiles and corresponding clinical information of MM were retrieved from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database (GSE24080 and GSE4581) (30–33). The GSE24080 dataset was composed of 554 MM patients. The GSE4581 dataset was utilized as a validation set that included 414 MM samples to verify the prognostic risk signature.

Acquisition of Hypoxia-Related Differentially Expressed Genes

The 198 hypoxia-related genes were downloaded from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/) database. To determine the hypoxia status (hypoxia^high and hypoxia^low) of MM samples from the GSE24080 dataset, a consistent clustering analysis was conducted based on the 198 hypoxia-related genes using the R package ConsensusClusterPlus. The similarity distance between MM samples was computed by the Euclidean distance, and the optimal number of the clusters was determined by the cumulative distribution function (CDF). The different overall survival (OS) among the identified different clusters was detected by the K-M curve analysis using the survival package in R.

Differentially Expressed Analysis

p < 0.05 and |log₂ Fold Change| ≥0.5 were set as the cutoffs to identify hypoxia-related differentially expressed genes (DEGs) between the hypoxia^high and hypoxia^low group using the limma package in R. The same selection standard was utilized to identify the immune-related DEGs and hypoxia–immune-related DEGs.

Functional Enrichment Analysis of Hypoxia-Related DEGs

For the hypoxia-associated DEGs, the functions and pathway enrichment were investigated by the gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis utilizing the R package clusterProfiler. The cutoff criterion was the false-discovery rate (FDR) of <0.05.

Selection of Immune-Related DEGs

The immune scores of MM samples were assessed by the ESTIMATE algorithm using the R package estimate. The optimal cutoff of the immune score calculated by the surv cutpoint algorithm of the R package survival was used to classify MM samples into two groups: immune^high group and immnue^low group.

Identification of Hypoxia–Immune-Related Genes

To explore whether there is an internal relationship between immunity and hypoxia in the pathogenesis and development of MM disease, hypoxia and immune statuses were further combined into a two-dimension indicator that separated MM samples into low hypoxia and high immune groups (hypoxia^low/immune^high), high hypoxia and low immune groups (hypoxia^high/immune^low), and a mixed group including high hypoxia and the high immune group as well as low hypoxia and low immune group (hypoxia^high/immune^high + hypoxia^low/immune^low). Survival analysis of these three groups was executed by K-M plots, and the p < 0.05 was considered statistically significant.

Moreover, the DEGs between hypoxia^low/immune^high and hypoxia^high/immune^low groups were screened using the “limma” R package, with p-value <0.05 and |log2FC|≥0.5. Also, the volcano plot and heatmap were applied to the visualization of DEGs’ expression. The key hypoxia–immune-related genes were obtained by overlapping the DEGs (hypoxia^low/immune^high VS hypoxia^high/immune^low) and hypoxia/immune-related DEGs.

Constructed Prognostic Risk Signature Associated With Hypoxia and Immune Statuses

Firstly, the MM samples in the GSE24080 dataset were randomly stratified into an internal training set (n = 388) and an internal testing set (n = 166) at a ratio of 7:3. In the training set, K-M curve and univariate Cox regression analyses were operated to screen genes related to the OS of MM (p < 0.05), which were regarded as prognosis-based genes. Multivariate COX regression analysis was then used to identify the most predictive prognostic genes. The risk score was calculated according to the gene expression and coefficient value of samples. The formula was as follows:

$$\text{Risk score} = e^{\text{sum }\left(\text{each gene's expression levels} \times \text{corresponding coefficient}\right)} / e^{\text{sum }\left(\text{each gene's mean expression levels }\times \text{\hspace{0.17em}corresponding coefficient}\right)}$$

Furthermore, the risk score was determined by the predicted Cox ph algorithm of the R survival package, and the patients were classified into a high- and low-risk group based on the median risk score. A comparison of the survival between the two risk groups was tested by the K-M curve analysis. The receiver operating characteristic (ROC) analysis was employed to measure the predictive power of the risk model using the survival-ROC package in R. In addition, the areas under the curve (AUC) of each of the resulting ROC curves were used to assess the model’s discrimination.

Independent Prognostic Analysis and Construction of a Nomogram

Stepwise, univariate and multivariate COX regression analyses were performed to assess whether clinicopathological features (PROT, AGE, SEX, PACE, B2M, CPR, CREAT, LDH, ALB, HGB, ASPC, BMPC, MRI, CPS1, CPR1, FLC, lgA, lgD, lgG, nonsecretory, NSE) and gene risk score were independent prognostic factors for MM patients (p < 0.05). Hazard ratios (HRs) estimated from Cox models were reported as relative risks with corresponding 95% confidence intervals (CIs). HRs >1 were considered a risk factor, and HRs <1 were considered a protective factor. All identified independent prognostic factors were then included to construct a nomogram using the R package rms to predict the survival probability. Concordance indexes (C-index) and calibration curves were used to assess the discrimination and calibration of the nomogram.

Real-Time Quantitative Polymerase Chain Reaction

Bone marrow (BM) samples derived from healthy candidate donors (n = 5), newly diagnosed MM patients (n = 13), and relapsed MM patients (n = 6) were obtained with consent at Pingxiang People’s Hospital and the Third Xiangya Hospital of Central South University (CSU) during the period 2018–2020. The experimental protocols were approved by the ethical committee of CSU. BM mononuclear cells were isolated with density gradient separation using Ficoll-Hypaque (Sigma-Aldrich, St Louis, MO). CD138⁺ cells were enriched using the Miltenyi microbead separation system (Miltenyi BioTech, Auburn, CA). Bone marrow processing was done within 2 h after collection. The clinical information of the specimens is listed in Supplementary Table S1. RNA was isolated from 1 × 10⁶ cells using TRIzol reagent (Vazyme, Nanjing, China). Then, 2 μg total RNA was reverse transcribed into cDNA with the HiScript III RT SuperMix (Vazyme). cDNA was amplified in a total volume of 20 μl using the SYBR qPCR Master Mix (Vazyme) and the amplification was performed by LightCycler 480 (Roche, Burgess Hill, UK) with specific primers (Table 1). The Gene expression levels were calculated using the 2^−ΔΔCt method (35).

TABLE 1

Table 1 Characteristics of gene-specific real-time PCR assays.

Statistical Analysis

Statistical analysis was conducted using the R software. Comparisons in proportions of variables between two groups were analyzed using the χ² test. An unpaired t-test or the Mann–Whitney U test was performed for the differences between continuous variables. Unless otherwise stipulated, the p < 0.05 was considered statistical significance.

Results

Identification of Hypoxia-Related Genes for MM

Hypoxia is well known as a critical hallmark of the tumor microenvironment (36). With the expression matrix of 198 hypoxia-related genes from the MSigDB database, the Euclidean distance was calculated between two cases in the GSE24080 dataset to identify the hypoxia status in the tumor environment of MM patients. As shown in Figure 2A, k = 2 was found to be more optimal and classified MM samples into two groups with less correlation: clusters 1 and 2. The result of the survival comparison analysis suggested that cluster 2 patients presented worse OS than cluster 1 patients (p = 0.0011, Figure 2B). The target gene expression changes of the HIF-1 signaling pathway were analyzed to identify the hypoxia status between the two clusters. Nine of 15 genes (60%) associated with the process of increasing oxygen delivery were overexpressed in cluster 1, while 2 of 12 genes (16.7%) associated with the process of reducing oxygen consumption were increased in cluster 2 (Figure 2C). Taken together, we uncovered that patients in clusters 1 and 2 were classified as hypoxia^low and a hypoxia^high, respectively. Moreover, we identified 639 DEGs between clusters 2 and 1 that were considered hypoxia-related DEGs for MM, containing 71 upregulated and 568 downregulated genes (Figures 2D, E; Supplementary Table S2). The results of functional enrichment analysis demonstrated that those DEGs were significantly enriched in the immune-related terms and pathways, including humoral immune response, acute inflammatory response, and myeloid leukocyte migration (Figure 2F; Supplementary Figure S1).

FIGURE 2

Figure 2 Differentially expressed analysis of hypoxia-associated 639 genes for multiple myeloma (MM). (A) The Euclidean distance was used to calculate the similarity distance between MM samples, and the optimal number of the clusters was determined by the cumulative distribution function (CDF). (B) Kaplan–Meier analysis of the different clusters to predict the survival of the patient. Upper: Kaplan–Meier curve of the overall survival between clusters 1 and 2 groups. Lower: the number of patients at risk in clusters 1 and 2 groups at different time points. (C) Identified the hypoxia status between the two clusters by analyzing the target gene expression changes of the HIF-1 signaling pathway. (D, E) In total, 639 hypoxia-related differentially expressed genes (DEGs) were identified with hierarchical cluster analysis for MM. (F) Functional enrichment analysis of the 639 DEGs.

Identification of Immune-Related Genes for MM

Considering the importance of the immune response in the MM tumor environment, the immune status was identified by ESTIMATE. The MM samples were clustered into high and low immunity groups based on the calculated median immune score (Figure 3A), which showed a significant difference in survival between the two groups (Figure 3B). Patients with high immune scores represented a higher survival rate than those with low immune scores (p = 0.0082, Figure 3B). Moreover, 443 genes were overexpressed in the immune^high group and 27 genes were overexpressed in the immune^low group, which were considered immune-related DEGs (Figures 3C, D; Supplementary Table S3).

FIGURE 3

Figure 3 Differentially expressed analysis of immune-related 470 genes for MM. (A) The distributions of immune scores across multiple myeloma samples in the GSE24080 datasets. (B) Kaplan–Meier analysis of the different immune scores to predict patient survival. (C, D) In total, 470 immune-related differentially expressed genes were identified with hierarchical cluster analysis for MM.

Identification of Hypoxia–Immune-Related Genes for MM

A previous study has reported the interaction between hypoxia and immune status in the MM tumor microenvironment (36). The above hypoxia and immune statuses were combined as a two-dimension indicator that split MM patients into three groups: hypoxia^high/immune^low group, hypoxia^low/immune^high group, and mixed group (hypoxia^high/immune^high + hypoxia^low/immune^low). K-M curve analysis demonstrated a remarkably different survival rate among the three groups. As shown in Figure 4A, the patients in the hypoxia^low/immune^high group exhibited the best survival rate, while patients in the hypoxia^high/immune^low group exhibited the worst survival rate, manifesting the mutual effects of hypoxia and immune response on the development of MM (p < 0.0001). A total of 1,909 DEGs were selected between the hypoxia^high/immune^low group and the hypoxia^low/immune^high group, in which 177 genes were upregulated in the hypoxia^high/immune^low group and 1,765 genes were upregulated in the hypoxia^low/immune^high group (Figures 4B, C). To obtain the hypoxia–immune-related genes, the above selected 1,909 hypoxia–immune-related DEGs, 639 hypoxia-related DEGs, and 470 immune-related DEGs were intersected to obtain 47 genes for subsequent study (Figure 4D; Supplementary Table S4).

FIGURE 4

Figure 4 Differentially expressed analysis of hypoxia–immune-related 47 genes for MM. (A) Kaplan–Meier analysis of the different groups based on the hypoxia and immune status to predict patient survival. (B, C) In total, 1,909 hypoxia–immune-related differentially expressed genes (DEGs) were identified with hierarchical cluster analysis for MM. (D) Forty-seven hypoxia–immune-related genes were obtained from intersecting with 639 hypoxia-related DEGs, 470 immune-related DEGs, and 1,909 hypoxia–immune-related DEGs.

Establishment of Risk Signature Related to the Hypoxia and Immune Status for the Prognosis of MM

To explore the prognostic value of the hypoxia–immune-related genes in MM, we first performed the K-M curve and univariate Cox regression analyses in the internal training set. IGHM, LYZ, LST1, and SLC19A1 were identified as significantly associated with the overall survival of MM (Figure 5A; Table 2). Afterward, LST1 (p = 0.016, HR = 0.8556) and SLC19A1 (p = 1e−04, HR = 1.2922) were further identified by the multiple stepwise Cox regression analysis to generate a prognostic risk signature for MM (Figure 5B; Table 3). The individual-level gene risk score for a patient was calculated, and the cutoff risk score divided patients into a high- and low-risk group in the internal training set. Survival analysis indicated that patients in the high-risk group showed a poor prognosis compared to those in the low-risk group (Figure 5C). As shown in Figure 5D, the patients died with an increasing gene risk score. The AUC for survival rate was 0.5754 at 1 year, 0.6078 at 3 years, and 0.6498 at 5 years, implying a better prognostic ability for MM (Figure 5E). Moreover, SLC19A1 was overexpressed and LST1 was downregulated in the high-risk group (Figure 5C).

FIGURE 5

Figure 5 Univariate and multiple stepwise Cox regression analyses of the hypoxia–immune-related genes in MM. (A, B) Multivariate Cox regression analysis of the genes IGHM, LYZ, LST1, and SLC19A1. (C) Kaplan–Meier analysis of the risk score model based on SLC19A1 and LST1 to predict patient survival in the internal training set. Upper: Kaplan–Meier curve. Lower: the number of patients at different time points. (D) The risk score distribution and the survival status. According to the cutoff value, patients with multiple myeloma were divided into high- and low-risk groups. (E) Time-dependent ROC curves for the two-gene model to predict patient survival.

TABLE 2

Table 2 Hypoxia-immune-related genes in K-M curve and univariate Cox regression analysis.

TABLE 3

Table 3 LST1 and SLC19A1 in the multiple stepwise Cox regression analysis.

Genetic Overexpression of SLC19A1 Is Associated With Poor Prognosis in MM

The constructed hypoxia–immune-gene prognosis signature based on SLC19A1 upregulation and LST1 downregulation was further validated in the internal testing set and GSE4581 dataset. Consistent with the findings in the internal training set, patients in the high-risk group showed shorter survival times than those in the low-risk group in both testing sets (Figures 6A, D). Additionally, ROC analysis verified the results that the risk signature had good performance for prognosing MM patients (Figures 6B, E). The detailed information on the distribution of risk score, survival status, and the two prognostic genes SLC19A1 and LST1 expression profiles in the two sets was displayed in Figures 6C, F. In conclusion, genetic overexpression of SLC19A1 is correlated with inferior prognosis in MM.

FIGURE 6

Figure 6 Genetic overexpression of SLC19A1 is correlated with poor prognosis in MM. (A, D) Kaplan–Meier analysis of the risk score model based on SLC19A1 and LST1 to predict patient survival in the internal testing set (A) and GSE4581 dataset (D). (B, E) Time-dependent ROC curves for the two-gene model to predict patient survival in the internal testing set (B) and GSE4581 dataset (E). (C, F) The risk score distribution and the survival status of the high- and low-risk groups in the internal testing set (C) and GSE4581 dataset (F).

Independent Prediction Capacity of the Risk Score Based on SLC19A1

Next, to determine whether the risk score based on SLC19A1 was an independent prognostic factor for MM patients, the univariate and multivariate Cox regression analyses were conducted in the GSE24080 and GSE4581 datasets. In the GSE24080 dataset, age-identified clinical prognostic factors (β₂-microglobulin, C-reactive protein, creatinine, LDH, albumin, hemoglobin, ASPC, BMPC, and MRI) and gene risk score were obviously connected with OS by univariate Cox regression analysis (Figure 7A) (38). The result of multivariate Cox regression analysis suggested that β₂-microglobulin, LDH, albumin, and gene risk score were independent prognostic factors for MM patients (Figure 7B). The prognostic signature contained all harvested independent factors that had better predictive power than alternative options (Figure 8A). In addition, age, β₂-microglobulin, LDH, albumin, MRI, and gene risk score were used to construct a nomogram to predict the 1-, 2-, 3-, 4, and 5-year overall survival of MM patients. The calibration plot manifested that the nomogram showed a good agreement with the ideal model, suggesting the constructed nomogram had great performance (Figures 8B, C).

FIGURE 7

Figure 7 Multivariate Cox regression analysis of the SLC19A1 risk score and clinical risk parameters. (A, B) B2M, LDH, ALB, and SLC19A1 risk scores were markedly correlated with the overall survival of patients with multiple myeloma.

FIGURE 8

Figure 8 The SLC19A1 risk score and clinical risk parameters with multivariate Cox regression analysis. (A) ROC curves are utilized to compare the predictive performance for AGE, B2M, LDH, ALB, MRI, riskScore, and ALL to predict patient survival. (B) The calibration curves for predicting 1-, 2-, 3-, 4-, and 5-year overall survival. (C) Establishment of a nomogram for prediction of 1-, 2-, 3-, 4-, and 5-year survival of MM patients.

High SLC19A1 Expression Is Correlated With Risk Stratification, Bone Marrow Aspiration Plasma Cell Numbers, and β₂-Microglobulin in MM Patients

The results indicated the expression levels of SLC19A1 in MM, including newly diagnosed and relapsed MM, were higher than those in normal control (Figures 9A, B). Moreover, we found that patients with high-risk cytogenetics (39) had much higher SLC19A1 expression compared with standard-risk patients (p = 0.017, Figure 9C). To determine the correlation of SLC19A1 expression with clinical features in MM patients, we classified the patients into a low SLC19A1 expression group (SLC19A1^low, n = 9) and a high SLC19A1 expression group (SLC19A1^high, n = 10) based on the median SLC19A1 expression level. As shown in Table 4, the high expression of SLC19A1 was found to be correlated with a higher number of bone marrow aspirated plasma cells (p = 0.037). Moreover, marked differences were observed in the group of β₂-microglobulin ≥5.5 mg/L (p = 0.029). However, no significant differences were observed in sex, age, albumin, C-reactive protein, LDH, creatinine, hemoglobin, cytogenetic abnormalities, and relapse MM patients between SLC19A1^high and SLC19A1^low patients (Table 4).

FIGURE 9

Figure 9 Expression of SLC19A1 in primary MM samples. (A) qRT-PCR analysis of SLC19A1 mRNA expression in MM patient samples (n = 19) and normal controls (n = 5). (B) SLC19A1 mRNA expression in newly diagnosed MM patients (n = 13), relapsed MM patients (n = 6), and normal controls (n = 5). (C) SLC19A1 mRNA expression in different risk stratifications. High risk defined as: del 17p, t(14;16), t(14;20), +1q, and del 1p. Intermediate risk defined as: t(4;14), del 13q, and hyperdiploidy. Standard-risk defined as: t(11;14), t(6;14), del 14q (IgH), and others.

TABLE 4

Table 4 Correlation of SLC19A1 expression with clinical and laboratorial parameters in MM patients.

Discussion

Hypoxia and HIFs are known to affect cancer cells by regulating the expression of many genes that control tumor cell metabolism, angiogenesis, tumor invasion, and metastasis (40, 41). Hypoxia in the tumor microenvironment regulates the function of immune cell effectors and plays a positive role in tumor development (18). Although there is a clear link between hypoxia and immunity (42, 43), fewer studies pay attention to the comprehensive molecular mechanism.

In the present study, we performed the differentially expressed analysis and Cox regression analysis to screen the potential hypoxia–immune-related genes and establish a gene-related risk score model to predict the overall survival (OS) of MM patients. Moreover, we included gene risk groups and clinical parameters in our nomogram to precisely predict patient survival. A total of 47 hypoxia–immune-related genes were analyzed by univariate and multiple stepwise Cox regression in our study, and the overexpressed gene SLC19A1 and the downregulated gene LST1 were identified as markedly correlated with the OS of MM. Additionally, the risk score model was established based on the two prognostic genes SLC19A1 and LST1. Through further univariate and multivariate Cox regression analyses, we found that the risk score model was well validated in the GSE24080 and GSE4581 datasets. These two genes showed good predictive performance, indicating that SLC19A1 and LST1 may act as independent prognostic factors for MM patients. Overall, our data suggested that genetic overexpression of SLC19A1 is correlated with inferior prognosis in MM.

Previous studies revealed that solute carrier (SLC) transporters play an important role in MM drug resistance, including resistance to novel therapies (44). SLC19A1 belongs to the solute carrier (SLC) group of transporters and is the main mechanism by which folates and antifolate drugs are delivered to mammalian cells and tissues (45). High levels of SLC19A1 transcripts are detected in the liver and placenta, with appreciable levels in other tissues, including bone marrow (46). SLC19A1 was recognized to drive uptake of the new generation antifolate pemetrexed (PMX) into tumor cells (47). Recent studies showed that the mRNA expression levels of SLC19A1 were markedly increased in bladder tumor specimens (25). SLC19A1 was significantly highly expressed in osteosarcoma cells (48). Maggini et al. reported that the SLC19A1 T-233T genotype is associated with an improved response after treatment with melphalan and ASCT in MM. In this study, we identified SLC19A1 as a hypoxia–immune-related gene in MM, which has been found to be the major importer of 2′3′-cyclic-GMP-AMP (cGAMP) and other CDNs (29). cGAMP is an immunotransmitter secreted by cancer cells that, when taken up by host cells, can trigger an antitumoral immune response (29, 49). Due to the antitumoral effects of cGAMP, SLC19A1’s role in MM is worthy of further study. Furthermore, our data suggested that the expression of SLC19A1 in newly diagnosed or relapsed MM patients was significantly increased. Moreover, the upregulation of SLC19A1 has been found to be associated with higher risk stratification, bone marrow aspiration plasma cells, and β₂-microglobulin level in MM patients.

In conclusion, we established a hypoxia–immune-associated gene risk score model for predicting the survival of MM patients by differential expression analysis and univariate and multivariate Cox regression analyses. The overexpressed gene SLC19A1 may be a prognostic biomarker for the inferior survival of MM patients, which is supported by previous studies (34). In addition, SLC19A1 is a potential hypoxia–immune-related gene for MM, which is really worthy of further study about the function of SLC19A1 in multidrug resistance and the molecular mechanism in hypoxic immune niches.

Data Availability Statement

The data is publicly accessible in the Gene Expression Omnibus (GEO) and can be found at the website: https://www.ncbi.nlm.nih.gov/geo/, accession numbers (GSE24080 and GSE4581).

Ethics Statement

The studies involving human participants were reviewed and approved by the Cancer Research Institute Review Board of Central South University (CSU). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

XW and DJ conceived, designed and supported the study. WJL and PY performed experiments and wrote the manuscript. WJL and QM conducted most of the data analysis. WQL, LX, and CX participated in collecting data and helped in data analysis. SX and YH revised the paper. All authors reviewed and approved the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (82060044) and the Natural Science Foundation of Guangxi Province (2020GXNSFAA159018).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.843369/full#supplementary-material

References

1. Bazarbachi AH, Al Hamed R, Malard F, Malard F, Harousseau JL, Mohty M. Relapsed Refractory Multiple Myeloma: A Comprehensive Overview. Leukemia (2019) 33(10):2343–57. doi: 10.1038/s41375-019-0561-2

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Pawlyn C, Davies FE. Toward Personalized Treatment in Multiple Myeloma Based on Molecular Characteristics. Blood (2019) 133(7):660–75. doi: 10.1182/blood-2018-09-825331

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Brigle K, Rogers B. Pathobiology and Diagnosis of Multiple Myeloma. Semin Oncol Nurs (2017) 33(3):225–36. doi: 10.1016/j.soncn.2017.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Bhutani M, Foureau DM, Atrash S, Voorhees PM, Usmani SZ. Extramedullary Multiple Myeloma. Leukemia (2020) 34(1):1–20. doi: 10.1038/s41375-019-0660-0

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Wallington-Beddoe CT, Mynott RL. Prognostic and Predictive Biomarker Developments in Multiple Myeloma. J Hematol Oncol (2021) 14(1):151. doi: 10.1186/s13045-021-01162-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Levin A, Hari P, Dhakal B. Novel Biomarkers in Multiple Myeloma. Transl Res (2018) 201:49–59. doi: 10.1016/j.trsl.2018.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Durer C, Durer S, Lee S, Chakraborty R, Malik MN, Rafae A, et al. Treatment of Relapsed Multiple Myeloma: Evidence-Based Recommendations. Blood Rev (2020) 39:100616. doi: 10.1016/j.blre.2019.100616

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Hanahan D, Weinberg RA. Hallmarks of Cancer: The Next Generation. Cell (2011) 144(5):646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yang G, Shi R, Zhang Q. Hypoxia and Oxygen-Sensing Signaling in Gene Regulation and Cancer Progression. Int J Mol Sci (2020) 21(21):8162. doi: 10.3390/ijms21218162

CrossRef Full Text | Google Scholar

10. Jing X, Yang F, Shao C, Wei K, Xie M, Shen H, et al. Role of Hypoxia in Cancer Therapy by Regulating the Tumor Microenvironment. Mol Canc (2019) 18(1):157. doi: 10.1186/s12943-019-1089-9

CrossRef Full Text | Google Scholar

11. Liao C, Zhang Q. Understanding the Oxygen-Sensing Pathway and Its Therapeutic Implications in Diseases. Am J Pathol (2020) 190(8):1584–95. doi: 10.1016/j.ajpath.2020.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Chakraborty AA, Laukka T, Myllykoski M, Ringel AE, Booker MA, Tolstorukov MY, et al. Histone Demethylase KDM6A Directly Senses Oxygen to Control Chromatin and Cell Fate. Science (2019) 363(6432):1217–22. doi: 10.1126/science.aaw1026

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Batie M, Frost J, Frost M, Wilson JW, Schofield P, Rocha S. Hypoxia Induces Rapid Changes to Histone Methylation and Reprograms Chromatin. Science (2019) 363(6432):1222–6. doi: 10.1126/science.aau5870

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Qian X, Li X, Shi Z, Bai X, Xia Y, Zheng Y, et al. KDM3A Senses Oxygen Availability to Regulate PGC-1α-Mediated Mitochondrial Biogenesis. Mol Cell (2019) 76(6):885–895.e7. doi: 10.1016/j.molcel.2019.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lee P, Chandel NS, Simon MC. Cellular Adaptation to Hypoxia Through Hypoxia Inducible Factors and Beyond. Nat Rev Mol Cell Biol (2020) 21(5):268–83. doi: 10.1038/s41580-020-0227-y

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Xiong Q, Liu B, Ding M, Zhou J, Yang C, Chen Y. Hypoxia and Cancer Related Pathology. Cancer Lett (2020) 486:1–7. doi: 10.1016/j.canlet.2020.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Multhoff G, Vaupel P. Hypoxia Compromises Anti-Cancer Immune Responses. Adv Exp Med Biol (2020) 1232:131–43. doi: 10.1007/978-3-030-34461-0_18

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Taylor CT, Colgan SP. Regulation of Immunity and Inflammation by Hypoxia in Immunological Niches. Nat Rev Immunol (2017) 17(12):774–85. doi: 10.1038/nri.2017.103

PubMed Abstract | CrossRef Full Text | Google Scholar

19. McGettrick AF, O'Neill LAJ. The Role of HIF in Immunity and Inflammation. Cell Metab (2020) 32(4):524–36. doi: 10.1016/j.cmet.2020.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Colgan SP, Furuta GT, Taylor CT. Hypoxia and Innate Immunity: Keeping Up With the HIFsters. Annu Rev Immunol (2020) 38:341–63. doi: 10.1146/annurev-immunol-100819-121537

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Yuen VW, Wong CC. Hypoxia-Inducible Factors and Innate Immunity in Liver Cancer. J Clin Invest (2020) 130(10):5052–62. doi: 10.1172/JCI137553

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhao R, Goldman ID. Folate and Thiamine Transporters Mediated by Facilitative Carriers (SLC19A1-3 and SLC46A1) and Folate Receptors. Mol Aspect Med (2013) 34(2-3):373–85. doi: 10.1016/j.mam.2012.07.006

CrossRef Full Text | Google Scholar

23. Liu S, Huang H, Lu X, Golinski M, Comesse S, Watt D, et al. Down-Regulation of Thiamine Transporter THTR2 Gene Expression in Breast Cancer and its Association With Resistance to Apoptosis. Mol Cancer Res (2003) 1(9):665–73.

PubMed Abstract | Google Scholar

24. Liu S, Stromberg A, Tai HH, Moscow JA. Thiamine Transporter Gene Expression and Exogenous Thiamine Modulate the Expression of Genes Involved in Drug and Prostaglandin Metabolism in Breast Cancer Cells. Mol Cancer Res (2004) 2(8):477–87. doi: 10.1158/1541-7786.477.2.8

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Abdel-Haleem AM, El-Zeiry MI, Mahran LG, Abou-Aisha K, Rady MH, Rohde J, et al. Expression of RFC/SLC19A1 is Associated With Tumor Type in Bladder Cancer Patients. PloS One (2011) 6(7):e21820. doi: 10.1371/journal.pone.0021820

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Moscow JA. Methotrexate Transport and Resistance. Leuk Lymphoma (1998) 30(3-4):215–24. doi: 10.3109/10428199809057535

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Corrales L, Glickman LH, McWhirter SM, Kanne DB, Sivick KE, Katibah GE, et al. Direct Activation of STING in the Tumor Microenvironment Leads to Potent and Systemic Tumor Regression and Immunity. Cell Rep (2015) 11(7):1018–30. doi: 10.1016/j.celrep.2015.04.031

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Luteijn RD, Zaver SA, Gowen BG, Wyman SK, Garelis NE, Onia L, et al. SLC19A1 Transports Immunoreactive Cyclic Dinucleotides. Nature (2019) 573(7774):434–8. doi: 10.1038/s41586-019-1553-0

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Ritchie C, Cordova AF, Hess GT, Hess GT, Bassik MC, Li L. SLC19A1 Is an Importer of the Immunotransmitter cGAMP. Mol Cell (2019) 75(2):372–381.e5. doi: 10.1016/j.molcel.2019.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhu FX, Wang XT, Zeng HQ, Yin ZH, Ye ZZ. A Predicted Risk Score Based on the Expression of 16 Autophagy-Related Genes for Multiple Myeloma Survival. Oncol Lett (2019) 18(5):5310–24. doi: 10.3892/ol.2019.10881

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Ibata S, Kobune M, Kikuchi S, Yoshida M, Miura S, Horiguchi H, et al. High Expression of Nucleoporin 133 mRNA in Bone Marrow CD138⁺ Cells is a Poor Prognostic Factor in Multiple Myeloma. Oncotarget (2018) 9(38):25127–35. doi: 10.18632/oncotarget.25350

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yang J, Wang F, Chen B. HLA-DPA1 Gene Is a Potential Predictor With Prognostic Values in Multiple Myeloma. BMC Canc (2020) 20(1):915. doi: 10.1186/s12885-020-07393-0

CrossRef Full Text | Google Scholar

33. Zhong Y, Liu Z, Li D, Liao Q, Li J. Identification and Validation of a Potential Prognostic 7-lncRNA Signature for Predicting Survival in Patients With Multiple Myeloma. BioMed Res Int (2020) 2020:3813546. doi: 10.1155/2020/3813546

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Jiang D, He Y, Mo Q, Liu E, Li X, Huang L, et al. PRICKLE1, a Wnt/PCP Signaling Component, Is Overexpressed and Associated With Inferior Prognosis in Acute Myeloid Leukemia. J Transl Med (2021) 19(1):211. doi: 10.1186/s12967-021-02873-8

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Sarasquete ME, Martínez-López J, Chillón MC, Alcoceba M, Corchete LA, Paiva B, et al. Evaluating Gene Expression Profiling by Quantitative Polymerase Chain Reaction to Develop a Clinically Feasible Test for Outcome Prediction in Multiple Myeloma. Br J Haematol 2013 163(2):223–34.

PubMed Abstract | Google Scholar

36. Semenza GL. Oxygen Sensing, Hypoxia-Inducible Factors, and Disease Pathophysiology. Annu Rev Pathol (2014) 9:47–71. doi: 10.1146/annurev-pathol-012513-104720

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Sarkar S, Germeraad WT, Rouschop KM, Steeghs EM, van Gelder M, Bos GM, et al. Hypoxia Induced Impairment of NK Cell Cytotoxicity Against Multiple Myeloma can be Overcome by IL-2 Activation of the NK Cells. PloS One (2013) 8(5):e64835. doi: 10.1371/journal.pone.0064835

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Qi T, Qu J, Tu C, Lu Q, Li G, Wang J, et al. Super-Enhancer Associated Five-Gene Risk Score Model Predicts Overall Survival in Multiple Myeloma Patients. Front Cell Dev Biol (2020) 8:596777. doi: 10.3389/fcell.2020.596777

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Mikhael JR, Dingli D, Roy V, Reeder CB, Buadi FK, Hayman SR, et al. Management of Newly Diagnosed Symptomatic Multiple Myeloma: Updated Mayo Stratification of Myeloma and Risk-Adapted Therapy (mSMART) Consensus Guidelines 2013. Mayo Clin Proc (2013) 88(4):360–76. doi: 10.1016/j.mayocp.2013.01.019

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Chappell JC, Payne LB, Rathmell WK. Hypoxia, Angiogenesis, and Metabolism in the Hereditary Kidney Cancers. J Clin Invest (2019) 129(2):442–51. doi: 10.1172/JCI120855

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Schito L, Semenza GL. Hypoxia-Inducible Factors: Master Regulators of Cancer Progression. Trends Canc (2016) 2(12):758–70. doi: 10.1016/j.trecan.2016.10.016

CrossRef Full Text | Google Scholar

42. Quäschling T, Friedrich D, Deepe GS Jr, Rupp J. Crosstalk Between Autophagy and Hypoxia-Inducible Factor-1α in Antifungal Immunity. Cells (2020) 9(10):2150. doi: 10.3390/cells9102150

CrossRef Full Text | Google Scholar

43. Krzywinska E, Stockmann C. Hypoxia, Metabolism and Immune Cell Function. Biomedicines (2018) 6(2):56. doi: 10.3390/biomedicines6020056

CrossRef Full Text | Google Scholar

44. Mynott RL, Wallington-Beddoe CT. Drug and Solute Transporters in Mediating Resistance to Novel Therapeutics in Multiple Myeloma. ACS Pharmacol Transl Sci (2021) 4(3):1050–65. doi: 10.1021/acsptsci.1c00074

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Matherly LH, Hou Z, Deng Y. Human Reduced Folate Carrier: Translation of Basic Biology to Cancer Etiology and Therapy. Cancer Metastasis Rev (2007) 26(1):111–28. doi: 10.1007/s10555-007-9046-2

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Whetstine JR, Flatley RM, Matherly LH. The Human Reduced Folate Carrier Gene is Ubiquitously and Differentially Expressed in Normal Human Tissues: Identification of Seven Non-Coding Exons and Characterization of a Novel Promoter. Biochem J (2002) 367(Pt 3):629–40. doi: 10.1042/bj20020512

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Qiu A, Jansen M, Sakaris A, Min SH, Chattopadhyay S, Tsai E, et al. Identification of an Intestinal Folate Transporter and the Molecular Basis for Hereditary Folate Malabsorption. Cell (2006) 127(5):917–28. doi: 10.1016/j.cell.2006.09.041

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Li C, Yuan B, Yu X, Ju L, Ma Y, Guo D, et al. SLC19A1 May Serve as a Potential Biomarker for Diagnosis and Prognosis in Osteosarcoma. Clin Lab (2020) 66(11). doi: 10.7754/Clin.Lab.2020.200246

CrossRef Full Text | Google Scholar

49. Marcus A, Mao AJ, Lensink-Vasan M, Wang L, Vance RE, Raulet DH. Tumor-Derived cGAMP Triggers a STING-Mediated Interferon Response in Non-Tumor Cells to Activate the NK Cell Response. Immunity (2018) 49(4):754–763.e4. doi: 10.1016/j.immuni.2018.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar