Skip to main content

ORIGINAL RESEARCH article

Front. Pharmacol., 07 August 2023
Sec. Pharmacology of Anti-Cancer Drugs

Comprehensive analysis of ferroptosis-related genes in immune infiltration and prognosis in multiple myeloma

Quanqiang Wang&#x;Quanqiang Wang1Misheng Zhao&#x;Misheng Zhao2Tianyu ZhangTianyu Zhang1Bingxin ZhangBingxin Zhang1Ziwei ZhengZiwei Zheng1Zhili LinZhili Lin1Shujuan ZhouShujuan Zhou1Dong ZhengDong Zheng1Zixing ChenZixing Chen3Sisi ZhengSisi Zheng1Yu ZhangYu Zhang1Xuanru LinXuanru Lin1Rujiao DongRujiao Dong1Jingjing ChenJingjing Chen1Honglan QianHonglan Qian1Xudong HuXudong Hu1Yan ZhuangYan Zhuang1Qianying ZhangQianying Zhang1Songfu Jiang
Songfu Jiang1*Yongyong Ma,,
Yongyong Ma1,4,5*
  • 1Department of Hematology, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, Zhejiang, China
  • 2Department of Clinical Laboratory, Wenzhou People’s Hospital, Wenzhou, China
  • 3Department of Hepatobiliary Surgery, The Second Affiliated Hospital and Yuying Children’s Hospital of Wenzhou Medical University, Wenzhou, Zhejiang, China
  • 4Key Laboratory of Intelligent Treatment and Life Support for Critical Diseases of Zhejiang Province, Wenzhou, Zhejiang, China
  • 5Zhejiang Engineering Research Center for Hospital Emergency and Process Digitization, Wenzhou, Zhejiang, China

Background: One particular type of cellular death that is known as ferroptosis is caused by the excessive lipid peroxidation. It is a regulated form of cell death that can affect the response of the tumor cells. Currently, it is not known if the presence of this condition can affect the prognosis of patients with multiple myeloma (MM).

Methods: In this study, we studied the expression differences and prognostic value of ferroptosis-related genes (FRGs) in MM, and established a ferroptosis risk scoring model. In order to improve the prediction accuracy and clinical applicability, a nomogram was also established. Through gene enrichment analysis, pathways closely related to high-risk groups were identified. We then explored the differences in risk stratification in drug sensitivity and immune patterns, and evaluated their value in prognostic prediction and treatment response. Lastly, we gathered MM cell lines and samples from patients to confirm the expression of marker FRGs using quantitative real-time PCR (qRT-PCR).

Results: The ability to predict the survival of MM patients is a challenging issue. Through the use of a risk model derived from ferroptosis, we were able to develop a more accurate prediction of the disease’s prognosis. They were then validated by a statistical analysis, which showed that the model is an independent factor in the prognosis of MM. Patients of high ferroptosis risk scores had a much worse chance of survival than those in the low-risk groups. The calibration and power of the nomogram were also strong. We noted that the link between the ferroptosis risk score and the clinical treatment was suggested by the FRG’s significant correlation with the immune checkpoint genes and the medication sensitivity. We validated the predictive model using qRT-PCR.

Conclusion: We demonstrated the association between FRGs and MM, and developed a new risk model for prognosis in MM patients. Our study sheds light on the potential clinical relevance of ferroptosis in MM and highlights its potential as a therapeutic target for patients with this disease.

1 Introduction

Uncontrolled plasma cell proliferation is the main cause of MM myeloma, which is a type of cancer that affects the bone marrow. Despite the advancements that have been made in treating this condition, the median overall survival rate for patients with this disease is still poor at only 5 years (Akhmetzyanova et al., 2020). Due to the existence of genetic abnormalities in MM patients, the importance of accurate risk stratification has been acknowledged. These abnormalities have been well-studied and are known to predict the prognosis of patients with this disease. (Abdallah et al., 2020) T (Palumbo et al., 2015; Bordini et al., 2020), T (Geeleher et al., 2014; Bordini et al., 2020), T (Bordini et al., 2020; Gonsalves et al., 2020), T (Rajkumar et al., 2014; Bordini et al., 2020), del (17/17p) and Nonhyperdiploid karyotype, karyotype del (Bordini et al., 2017) and gain (1q) were identified as important factors associated with poor prognosis. They can be used to form a predictive or prognostic score based on various clinical factors. Some of these include age, albumin, gender, and β2-microglobulin levels (Sonneveld et al., 2016). The Revised International Staging System (R-ISS) is also widely used to predict the prognosis of MM patients (Palumbo et al., 2015). The prognosis may vary significantly within each risk group indicated by the R-ISS, though. One study revealed that the R-ISS is more accurate than the International Staging System (ISS) when it comes to predicting a patient’s survival. However, it did not perform well when it came to assigning grades to the individuals (Cho et al., 2017). One study revealed that the R-ISS is more accurate than the International Staging System (ISS) when it comes to predicting a patient’s survival. However, it did not perform well when it came to assigning grades to the individuals. Several methodologies that could potentially enhance the precision of the R-ISS include analyzing blood clone plasma cells (Gonsalves et al., 2020), and imaging techniques such as positron emission tomography (Jung et al., 2019), and next-generation sequencing (Bolli et al., 2020). Unfortunately, these did not significantly increase R-ISS prediction accuracy.

Recently, people are starting to notice a new type of cell death known as ferroptosis, which is associated with various types of cancer. This mode of cell death is triggered by the imbalance between reactive oxygen species (ROS) and intracellular lipid peroxide. GPX4 and SLC7A11 are important regulators of ferroptosis and are highly expressed in MM cells. It was found that the immunosuppressant fingolimod (FTY720) could promote ferroptosis by reducing the mRNA and protein levels of GPX4 and SLC7A11 in U266 cells (Zhong et al., 2020). The process of lipid peroxide accumulation in cells requires the participation of iron ions, so ferroptosis is iron-dependent (Yan et al., 2021). In the study of MM model, plasma cells are very sensitive to excess iron. By producing antibodies, a large number of by-products such as H2O2 are synthesized, and finally ROS are stimulated by Fenton reaction (Tu and Weissman, 2004; Shimizu and Hendershot, 2009; Bordini et al., 2017). As a result, inducing iron overload may impede malignant plasma cell proliferation and improve the abilities bortezomib and carfilzomib to control disease development. Different MM cell lines were treated with high doses of ferrous ammonium citrate (FeAC), and untreated and other cell lines were used as controls for in vitro experiments. Iron caused cell death in MM cell lines, but not in control cells, and decreased proliferation was detected in all cell lines. Multiple myeloma cells exposed to iron undergo lipid oxidation and inhibit proteasome function, which means that the role of proteasome inhibitors will be enhanced (Campanella et al., 2013; Bordini et al., 2020). Our predictive model was developed using the FRGs public dataset and has demonstrated independent precision, and the R-ISS’ prognostic evaluation accuracy was further improved by incorporating it.

2 Materials and methods

2.1 Data acquisition

The gene expression profiles and clinicopathological information of three MM datasets GSE136337, GSE4204 and GSE24080 were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), and the gene expression profiles were normalized between different arrays. Subsequently, we obtained 252 ferroptosis-related genes (FRGs) using the FerrDb database (www.zhounan.org/ferrdb).

2.2 Construction and validation of a ferroptosis risk score

The GSE136337 data set is moderate in size and has sufficient clinical data, including but not limited to R-ISS staging, ISS staging, genetic mutation information, etc., so we decided to use the GSE136337 data set as the final training data set for the generated ferroptosis risk score. In order to discover which FRGs were substantially linked with prognosis (p < 0.001), we utilized the GSE136337 data set and carried out a univariate Cox regression analysis. For the purpose of integrating the various data related to the study, such as the survival status, gene expression, and survival time, we made use of the “glmnet” R software package. Following that, we carried out a regression analysis utilizing the Lasso-Cox method in order to obtain the best possible model. Finally, the lambda value was determined to be 0.00387132108721305, and five FRGs were selected. We use GSE4204 and GSE24080 datasets to verify the model. The subjects were divided into two risk groups based on the median risk score of each data set. To assess the risk of ferroptosis, we used univariate and multivariate Cox regression analysis on the training set (GSE136337) and the validation set (GSE24080).

2.3 Constructing a predictive nomogram

The ROC curves were used to compare the prediction value of the R-ISS and the joint model with respect to the 3- and 5-year survival of patients. The AUC, which is a measure of the models’ accuracy, was used to analyze the differences between them.

When compared to the R-ISS, we discovered that the combined model had a significantly higher probability of correctly predicting the survival rate of multiple myeloma patients over the next three to five years. This suggests that the use of both the ferroptosis risk score and the R-ISS stage parameter can improve the accuracy of the joint model when it comes to identifying patients with this condition.

2.4 Gene set enrichment analysis

The Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis is a widely used approach for pathway enrichment analysis. It helps to understand the biological function and interaction of genes and proteins by categorizing them into known pathways. Gene Set Enrichment Analysis (GSEA) is another widely used method for identifying differentially expressed genes and enriched pathways. It involves comparing the distribution of a pre-defined set of genes (e.g., a KEGG pathway) in two different groups of samples (high-risk vs. low-risk group) to determine whether the pathway is significantly enriched in one group compared to the other. In this study, the KEGG pathway was used to reveal the potential causes of ferroptosis risk score. GSEA was used to test the enrichment pathway in different ferroptosis risk score groups (GSEA v4.2.3 software, http://software.broadinstitute.org/gsea/login.jsp). Statistical significance was defined as p < 0.05, false discovery rate q < 0.25.

2.5 Drug sensitivity prediction

The “pRRophetic” (Geeleher et al., 2014) R package is a tool for predicting the half-maximal inhibitory concentration (IC50) of chemotherapeutic medicines using gene expression data. It uses a linear regression model that incorporates gene expression data from a training set of cell lines with known IC50 values for a given drug. The resulting model is then used to predict the IC50 values for new samples based on their gene expression profiles. This approach allows for the prediction of drug sensitivity or resistance in cancer patients and may be used to guide treatment decisions. We used the “pRRophetic” R package to test the drug sensitivity of the training set GSE136337 ferroptosis risk score high and low groups, and explored which drugs were different in the ferroptosis risk groups <0.05 was considered statistically significant.

2.6 Immune cell infiltration and immune microenvironment evaluation

The study utilized different algorithms for its analysis, such as the XCELL, TIDE, and CIBERSORT. The use of the XCELL algorithm can help determine the presence of various immune cells in the tumor microenvironment. This information can then be utilized to analyze the patterns of immune cell infiltration in the patients with MM of different risk group (Aran et al., 2017). The use of the CIBERSORT algorithm can also help identify the subsets of immune cells that are differentially expressed in patients with different risk levels. This method could provide a deeper analysis of the immune system’s response in patients (Chen et al., 2018). This can provide a deeper understanding of the immune response in different patient groups and can help identify potential targets for immunotherapy. The use of the TIDE algorithm can also help predict the response of the immune system to certain drugs in different patient groups. This method is useful in guiding the development of personalized treatment plans for patients with different types of cancer (Jiang et al., 2018). We also performed IPS scores on high- and low-risk groups to compare the different immune responses in each group.

2.7 Cell culture and reagents

The MM cell lines, namely H929, I9.2, U266, and RPMI 8226, were procured from Hunan Fenghui Biotechnology Co., Ltd. Gibco provided the cell culture and reagents. Short tandem repeat analysis was used to identify the cell lines, and they were routinely screened for mycoplasma contamination with Vazyme’s Myco-Blue Mycoplasma Detection Kit. The cells were grown in an incubator at 37 degrees Celsius and 5% carbon dioxide in a medium called RPMI 1640. The medium contained 10% fetal bovine serum, 100 international units per milliliter of penicillin, and 100 mg per milliliter of streptomycin.

2.8 Patient

31 MM patients were involved in the study at the Department of Clinical Hematology of the First Affiliated Hospital of Wenzhou Medical University, between January 2022 and November 2022. The International Myeloma Working Group (IMWG) criteria from 2014 were used to make the diagnoses (Rajkumar et al., 2014). Meanwhile, normal bone marrow samples from 14 healthy donors were collected as controls. Supplementary Table S1 depicts the distribution of clinical indicators and clinicopathological characteristics among the patients. The acquired samples were collected with the subjects’ informed consent. The study was authorized by the Ethics Committee of Wenzhou Medical University’s First Affiliated Hospital, and all operations were performed in accordance with the Helsinki Declaration.

2.9 Real-time fluorescence quantitative PCR

We validated the expression of five ferroptosis genes in MM using four human MM cell lines. Both the patients and the people in the control group had an aspirate of 5 mL taken from their bone marrow. The technique of density gradient centrifugation was utilized in order to isolate bone marrow mononuclear cells (BMMNC). RNA was extracted from the subjects ' BMMNC and myeloma cell lines using an RNA extraction kit (Yuanqi Bio, Shanghai, China) and reverse transcribed using HiScript ® III RT SuperMix (Vazyme, Nanjing, China). The primers were purchased from Yimai Biotechnology Co. Ltd. β-actin was used as a relatively quantitative endogenous control. The relative mRNA expression of various molecules was calculated using the 2−δCt method. The Ct value of the molecules was then normalized to the β-Actin value. Each experiment was repeated three times. Then different test methods are used to compare the experimental results. The forward and reverse primers of each molecule are placed in Supplementary Table S2.

2.10 Statistical analyses

R 4.1.2 was used for the analysis of the data. Normalization and lognormality tests were also performed using the software known as Prism 9.0. Mann-Whitney and the t-test were used to compare two groups, while Kruskal-Wallis and ANOVA were utilized for the comparison of multiple groups.

3 Results

3.1 Baseline covariates and subject selection

The analysis included survival data from 1,512 subjects from three datasets. There were sufficient clinical covariates in GSE136337 (N = 417) and GSE24080 (N = 557) for univariate and multivariate Cox regression analysis, but not in GSE4204 (N = 538). The clinical covariates of all three datasets were converted into categorical variables, as shown in Table 1. 252 FRGs were extracted from the FerrDb database to determine the risk score of ferroptosis in the GSE136337 cohort. Finally, the GSE4204 and GSE24080 datasets are used to verify the model.

TABLE 1
www.frontiersin.org

TABLE 1. Clinical messages of the training and validation cohorts.

3.2 The establishment of a prognostic ferroptosis risk score

We discovered five FRGs that were significant predictors of survival (p < 0.001) using GSE136337 as the training data set for univariate Cox regression analysis (Supplementary Figure S1). The ferroptosis risk score was then generated by incorporating various data elements, such as survival status, gene expression, and the time of the survival. The Lambda value was 0.00387132108721305, and five genes were chosen (Figures 1A, B). Ferroptosis risk score = (0.611358×YY1AP1 expression) + (0.26426×AURKA expression) - (0.36283×CDKN1A expression) + (0.07901×RRM2 expression) + (0.14758×STEAP3 expression). The algorithm was then used to estimate each individual’s risk score, and based on the median value of the data set, it was used to divide each individual into low-risk and high-risk categories.

FIGURE 1
www.frontiersin.org

FIGURE 1. Construction of ferroptosis model. (A) Lasso Cox regression analysis was used for 10-fold cross-validation of variable selection. (B) LASSO coefficient of ferroptosis genes. Each curve represents a ferroptosis gene.

3.3 Validation of prognostic ferroptosis risk score in each cohort

In order to evaluate the levels of specificity and sensitivity exhibited by the ferroptosis risk scores, a time-ROC analysis was carried out. The AUCs of the 1, 3, and 5-year prognostic survival models based on FRGs expression to predict the survival of plasma cell myeloma in the GSE136337 training dataset were 0.61 (95% CI: 0.50, 0.72), 0.69 (0.62, 0.79), and 0.70 (0.65, 0.76), respectively (Figure 2A). Kaplan-Meier curves were used to compare the survival of high-risk and low-risk cohorts in the training set (Figure 2D) and validation sets GSE4204 and GSE24080 (Figures 2E, F). Compared with the low-risk cohort (61% [54, 68%] vs.84% [78, 88%]); 71% [64, 79%] vs. 87% [82, 92%]; 58% [52, 65%] vs. 75% [69, 81%]; subjects in the high-risk group had worse survival than those in the low-risk group. A dot plot was created to compare the survival rates of individuals in high-risk and low-risk populations, and it was found that in each data set (Figures 2G, I), low-risk populations had better survival than high-risk populations. We also compared the expression of these five FRGs using heat maps. Although the expression is slightly different in each data set, it is still relatively stable (Figures 2G, I).

FIGURE 2
www.frontiersin.org

FIGURE 2. Validation of the ferroptosis risk-scoring model. (A–C) Sensitivity and specificity of the ferroptosis risk score model were assessed in each dataset by time-dependent ROC analysis. (D–F) Survival differences between high- and low-risk cohorts in each dataset. (G–I) Dot plots comparing outcomes of subjects in the high- and low-risk cohorts. The heat map displays results for the seven genes in both the training and validation cohorts. (A), (D) and (G) display GSE136337; (B), (E) and (H) display GSE4204; (C), (F) and (I) display GSE24080.

3.4 Uni-variable and multi-variable analyses

We analyzed the gender, age, albumin, β2-microglobulin, LDH, t [4; 14], t [14; 16], del [17p], ISS and R-ISS by univariate Cox regression analysis. In the following multivariate COX regression analysis, albumin, β2-microglobulin and LDH were excluded due to their collinearity with ISS and R-ISS. The same analysis was also performed on the clinical covariates in the validation dataset (GSE24080) (Table 2). According to multivariate analysis, ferroptosis risk scores were independently associated with survival. The hazard ratio (HR) in the GSE136337 was 2.896 (95% CI: 2.250–3.728, p < 0.001), while the HR in the GSE24080 was 2.086 (95% CI: 1.623–2.681, p < 0.001).

TABLE 2
www.frontiersin.org

TABLE 2. Univariate analysis and multi-variate regression analysis of overall survival in the training and validation cohorts by Cox regression analysis.

3.5 Creating a combined nomogram

To predict survival, a joint nomogram model was constructed using R-ISS stage and disease-related ferroptosis risk score (Figure 3A). The calibration plot revealed that the nomogram was a great predictor of 1-, 3-, and 5-year survival (Figure 3B) The graph shows an increase in the prediction accuracy of 3-year survival rate and 5-year survival rate in the training cohort increased from 0.55 (0.49, 0.62) and 0.60 (0.56, 0.65) of R-ISS to 0.68 (0.61, 0.76) of AUC. p = 0.001), 0.70 (0.65, 0.76); p = 0.0016) (Figure 3C; Supplementary Figure S2). However, the validation data for the study did not have R-ISS. This means that the risk score could not be combined with the R-ISS. The 3- and 5-year AUCs calculated from the ISS for the GSE24080 dataset increased from 0.60 to 0.67 and 0.58 to 0.61, respectively (Supplementary Figure S3). In both the training and validation sets, the pooled risk score outperformed other factors like sex, age, albumin, and LDH.

FIGURE 3
www.frontiersin.org

FIGURE 3. Building the combined nomogram to predict the overall survival (OS) of patients with multiple myeloma (MM). (A) The nomogram plot was built based on R-ISS staging, ferroptosis risk score and total points in the training cohort. (B) Timedependent receiver operating characteristic (ROC) curves of nomograms were compared based on co-variates 5-year survival. (C) Five-year ROC curves of the merged risk score compared with clinical covariates in the training dataset.

3.6 Analyses of gene set enrichment

We performed GSEA in all datasets, exploring ferroptosis-associated pathways and KEGG pathways associated with other ferroptosis covariates. Significantly enriched pathways, such as base excision repair and DNA replication, cell cycle, and mismatch repair, were concentrated in high-risk groups (Figure 4). Other ferroptosis-related pathways enriched in the high-risk group of the training and validation cohorts included mismatch repair, p53 signaling pathway, proteasome nucleotide excision repair, oocyte meiosis, pathogenic Escherichia coli infection, and spliceosome.

FIGURE 4
www.frontiersin.org

FIGURE 4. Pathways with significant enrichment in each dataset. The top 10 pathways were enriched in the training cohort (A) and validation datasets (B–C).

3.7 Drug sensitivity drug testing

Furthermore, we compared the drug sensitivity of the high-risk with that of low-risk groups, and p < 0.05 was deemed statistically significant. The half inhibitory concentration (IC50), which symbolizes the concentration of the inhibitor needed to inhibit 50% of the drug, is used to measure the effectiveness of the drug. We found that all-trans-retinoic acid (ATRA), cisplatin, cyclopamine, bexarotene, gefitinib and imatinib had lower IC50 in high-risk groups, that is, higher sensitivity (Figure 5).

FIGURE 5
www.frontiersin.org

FIGURE 5. Prediction of drug sensitivity in high-risk group and low-risk group. Comparison of drug sensitivity between high-risk and low-risk groups was conducted. Statistical significance was determined at p < 0.05. The effectiveness of the drugs was measured using IC50 for the following drugs: ATRA, cisplatin, cyclopamine, bexarotene, gefitinib, and imatinib.

3.8 Immune-Related Analysis of MM Patients Using the Ferroptosis Score

To investigate the association between FRGs and antitumor immunity in MM patients, we utilized the CIBERSORT algorithm to evaluate immune cell infiltration in the GSE136337 dataset. The results revealed the proportional percentages of several immune cell types (Figure 6A). Stromal, immunological, and microenvironment scores were calculated to compare immune cell infiltration between individuals of different risk groups. Our analysis indicated that the scores of stromal, immune, and microenvironment were significantly higher in the high-risk group (p < 0.05) (Figure 6B). Furthermore, significant differences in activated CD4 memory T cells and Macrophage M0 were observed when comparing the proportions of different immune cell types in the two groups (Figure 6C). According to Pearson correlation analysis, the expression levels of FRGs showed a significant correlation with immune cell populations (p < 0.05). Within the risk model, three FRGs were found to be associated with T cell CD4 memory activation. Specifically, CDKN1A, the only protective gene in the model, exhibited a negative correlation with T cells CD4 memory activated, while AURKA and RRM2 showed a positive correlation (Figure 6D). Additionally, we investigated the relationship between FRGs and the tumor immune dysfunction and exclusion (TIDE) score. Our findings revealed significant correlations between AURKA, CDKN1A, and STEAP3 and tumor exclusion (Figure 6E). We also evaluated the immune prognostic signature (IPS) scores in two groups and found higher IPS scores in the high-risk group, indicating greater responsiveness to immunotherapy (Figure 6F).

FIGURE 6
www.frontiersin.org

FIGURE 6. Expression of ferroptosis-related genes in normal human and MM cell lines. The relative mRNA expression levels of YY1AP1, AURKA, RRM2, STEAP3, and CDKN1A were analyzed using qRT-PCR in RPMI8226, H929, U266, and I9.2 cell lines.

3.9 Relative mRNA expression levels of cell lines and patients

The qRT-PCR results showed that the relative mRNA expression for YY1AP1, AURKA, RRM2, and STEAP3 in RPMI8226, H929, U266 and I9.2 was statistically significantly higher than the control, while CDKN1A was decreased (Figure 7). In addition, compared with the control, the relative mRNA expression of YY1AP1, AURKA, RRM2, STEA3 were significantly increased in BMMNC of 31 MM patients while the expression of CDKN1A was decreased, which is consistent with our previous analysis results (Figure 8).

FIGURE 7
www.frontiersin.org

FIGURE 7. Expression of ferroptosis-related genes in normal human and MM patients. The relative mRNA expression levels of YY1AP1, AURKA, RRM2, STEAP3, and CDKN1A were analyzed using qRT-PCR in BMMNC derived from 31 MM patients.

FIGURE 8
www.frontiersin.org

FIGURE 8. Immune-Related Analysis of MM Patients Using the Ferroptosis Score. (A) CIBERSORT algorithm determined the immune cell infiltration of all MM patients with GSE136337, showing the proportion of each typical immune cell. (B) The XCELL algorithm compared the matrix score, immune score and microenvironment score between the high-risk group and the low-risk group. (C) Compare the proportion of each immune cell between high-risk and low-risk groups. (D) The correlation between ferroptosis-related genes and immune cells was analyzed by Spearman statistical method. (E) Correlation analysis between model ferroptosis genes and TIDE score. (F) IPS model scores of high-risk group and low-risk group.

4 Discussion

We developed a new prognostic risk score for ferroptosis based on 5-FRGs, using the GSE136337 dataset as the training dataset. We validated this score in two independent datasets, and in multivariate analysis with R-ISS, the ferroptosis risk score was confirmed as an independent prognostic factor.

Our ferroptosis risk score identified CDKN1A as a positive prognostic factor and YY1AP1, AURKA, RRM2, and STEAP3 as negative prognostic factors. AURKA, a serine/threonine kinase that regulates cell division during mitosis, was overexpressed in several tumor types compared to normal tissues according to TCGA database. It has been shown to promote the development of hematological malignancies and solid tumors (Du et al., 2021).

RRM2, a stress response factor, has been associated with the initiation and advancement of various cancers. Analysis of the Oncomine database showed elevated expression of RRM2 in MM patients in comparison to healthy controls. In addition, downregulation of RRM2 by targeting the Wnt/-catenin signaling pathway may promote MM cell death and activate DNA damage (Liu et al., 2019).

CDKN1A is responsible for encoding P21, which is a founding member of the cyclin-dependent kinase inhibitor family (CKI) (Xiong et al., 1993). CKI is a crucial regulator of the cell cycle that maintains genomic stability, and is frequently dysregulated in human cancers (Roninson, 2002; Abbas and Dutta, 2009; Kreis et al., 2015). Studies have shown that downregulation of p21 expression can promote the proliferation of MM cells (Xiang et al., 2021). Highly expressed CDKN1A may play a tumor suppressor role by encoding p21.

Prostate family member 3 (STEAP3) six transmembrane epithelial antigen was originally discovered on mouse M1 bone marrow cells. According to research, STEAP3 plays a critical role in the regulation of ferroptosis via modulating iron metabolism (Mou et al., 2019). The proliferation of cancer cells, including, pancreatic cancer (Yu et al., 2020), colorectal cancer (Barresi et al., 2016) and bladder cancer (Kim et al., 2016) may be facilitated by the overexpression of STEAP3, which promotes iron uptake and storage.

Based on the results of GSEA, the high-risk score group demonstrated a significant enrichment of ferroptosis-related pathways. Among them, the base excision repair (BER) pathway was notably enriched, which may suggest that cancer cells depend on heightened BER activity to withstand oxidative stress. The proteasome pathway was also enriched, which is a myeloma cell-dependent proteasome complex that manages protein overload by collecting misfolded proteins in the endoplasmic reticulum. Additionally, DNA replication-enriched pathways were found to promote cell biosynthesis.

We performed a drug sensitivity test with high and low risk groups. ATRA is a vitamin A derivative. The anti-IL-6 drug ATRA can inhibit the growth of human MM cell lines. It has also been shown that the combination of this drug and the other medication, daratumumab, can effectively treat patients with this type of cancer (Frerichs et al., 2021). The ATRA works by increasing the CD38 expression in cancer cells and boosting the complement-dependent and antibody-mediated cytotoxic properties of daratumumab. It also increased the activity of the drug in a mouse model of MM (Nijhof et al., 2015). After large-scale data analysis of the gene spectrum of MM patients, it can be confirmed that ATRA can make human myeloma cells sensitive to carfilzomib, a proteasome inhibitor, by activating the retinoic acid receptor (RAR) γ and interferon-β response pathways, thereby effectively improving their survival (Wang et al., 2022). Cyclopamine has the ability to make myeloma cells more susceptible to apoptosis by circulating replacement TRAIL (Wang et al., 2021). Imatinib mesylate (STI571) has the ability to limit the proliferation of multiple myeloma (MM) cells and boost the effectiveness of anti-myeloma medicines like dexamethasone, which are routinely utilized (Pandiella et al., 2003). Understanding the drug sensitivity genomics of MM can better arrange targeted anti-myeloma treatment, which is also of great significance for myeloma drugs such as bortezomib (Ghasemi et al., 2016).

Using the XCELL algorithm, the data collected by the study revealed that the immune infiltration levels of high-risk individuals were more significant than those of the lower-risk group. It suggested that treating them with immunotherapies could improve their cancer responses. Using the XCELL algorithm, the data collected by the study revealed that the immune infiltration levels of high-risk individuals were more significant than those of the lower-risk group. It suggested that treating them with immunotherapies could improve their cancer responses (MacLeod et al., 2010). These mechanisms can be utilized by tumor-specific CD4 + T cells to deliver effective anti-cancer immune responses. Some of these include the licensing of CD8 + T cells (Schietinger et al., 2010), cytotoxic killing of tumor cells expressing MHC-II (Lundin et al., 2003; Quezada et al., 2010), macrophage (Corthay et al., 2005) and natural killer (NK) cell activation (Corthay et al., 2005) and cytokine-mediated effects on the tumor vascular system (Qin and Blankenstein, 2000). And the correlation analysis of FRGs with TIDE score and T cell exclusion and dysfunction score was significant. It is possible that these factors are responsible for the higher immunological score exhibited by the high-risk group. And the IPS model score also fully proves this finding. The MM ferroptosis risk score can help predict the disease’s prognosis and play a crucial role in the development of immunotherapies for this condition.

Our study is not without limitations. Firstly, some potential prognostic factors could not be accounted for in all three datasets due to missing information, with validation cohorts only including ISS data. Secondly, treatment effects were not adjusted for in the data, making our scores more prognostic than predictive. Thirdly, it is not known how much each gene contributes individually to the ferroptosis score, which is something that has to be investigated further. Lastly, the specific role of YY1AP1 in MM remains unclear. Future studies should focus on adjusting for treatment effects to enhance the predictive capabilities of the scores and improve patient outcome assessment. Additionally, investigating the function and significance of YY1AP1 in multiple myeloma, including its interactions with other genes and pathways, can provide valuable insights for disease progression and treatment strategies. Moreover, integrating phenomics into the research can also deepen our understanding of how genes, epigenetics, symbiotic microorganisms, diet, and environmental exposure contribute to multiple myeloma, potentially leading to the discovery of novel biomarkers and therapeutic targets (Liu et al., 2022; Ying, 2023).

In conclusion, we have created and validated a prognostic model for MM that is based on the expression of FRGs in MM cells. Our ferroptosis risk score proved to be a distinct predictor of survival after subjecting it to multivariate analysis. Individuals with high-risk scores can also have better options in the selection of chemotherapy drugs. And ferroptosis risk score may play an important role in the guidance of immunotherapy in clinical MM patients. In vitro experiments confirmed the expression of FRGs in MM patients. The addition of our ferroptosis risk score to R-ISS improved the accuracy of survival prediction and could potentially be used in combination with R-ISS for prognostic purposes.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving human participants were reviewed and approved by the First Affiliated Hospital of Wenzhou Medical University’s ethics committee. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

QW designed the study, conducted data analysis and interpretation, carried out experiments, and wrote the manuscript. As a newly added author, MZ participated in the revision of the manuscript and completed a supplementary experiment. The other authors provided input on data analysis. YM contributed vital comments and performed critical revisions on the article. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by grants from The National Natural Science Foundation (grant no. 82270212), the Natural Science Foundation of Zhejiang province (grant no. LY20H080003) and the Wenzhou Municipal Science and Technology Bureau (grant no. Y20220716).

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/fphar.2023.1203125/full#supplementary-material

References

Abbas, T., and Dutta, A. (2009). p21 in cancer: intricate networks and multiple activities. Nat. Rev. Cancer 9 (6), 400–414. doi:10.1038/nrc2657

PubMed Abstract | CrossRef Full Text | Google Scholar

Abdallah, N., Rajkumar, S. V., Greipp, P., Kapoor, P., Gertz, M. A., Dispenzieri, A., et al. (2020). Cytogenetic abnormalities in multiple myeloma: Association with disease characteristics and treatment response. Blood Cancer J. 10 (8), 82. doi:10.1038/s41408-020-00348-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Akhmetzyanova, I., McCarron, M. J., Parekh, S., Chesi, M., Bergsagel, P. L., and Fooksman, D. R. (2020). Dynamic CD138 surface expression regulates switch between myeloma growth and dissemination. Leukemia 34 (1), 245–256. doi:10.1038/s41375-019-0519-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Aran, D., Hu, Z., and Butte, A. J. (2017). xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18 (1), 220. doi:10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Barresi, V., Trovato-Salinaro, A., Spampinato, G., Musso, N., Castorina, S., Rizzarelli, E., et al. (2016). Transcriptome analysis of copper homeostasis genes reveals coordinated upregulation of SLC31A1,SCO1, and COX11 in colorectal cancer. FEBS Open Bio 6 (8), 794–806. doi:10.1002/2211-5463.12060

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolli, N., Genuardi, E., Ziccheddu, B., Martello, M., Oliva, S., and Terragna, C. (2020). Next-generation sequencing for clinical management of multiple myeloma: Ready for prime time? Front. Oncol. 10, 189. doi:10.3389/fonc.2020.00189

PubMed Abstract | CrossRef Full Text | Google Scholar

Bordini, J., Galvan, S., Ponzoni, M., Bertilaccio, M. T. S., Chesi, M., Bergsagel, P. L., et al. (2017). Induction of iron excess restricts malignant plasma cells expansion and potentiates bortezomib effect in models of multiple myeloma. Leukemia 31 (4), 967–970. doi:10.1038/leu.2016.346

PubMed Abstract | CrossRef Full Text | Google Scholar

Bordini, J., Morisi, F., Cerruti, F., Cascio, P., Camaschella, C., Ghia, P., et al. (2020). Iron causes lipid oxidation and inhibits proteasome function in multiple myeloma cells: A proof of concept for novel combination therapies. Cancers (Basel). 12 (4), 970. doi:10.3390/cancers12040970

PubMed Abstract | CrossRef Full Text | Google Scholar

Campanella, A., Santambrogio, P., Fontana, F., Frenquelli, M., Cenci, S., Marcatti, M., et al. (2013). Iron increases the susceptibility of multiple myeloma cells to bortezomib. Haematologica 98 (6), 971–979. doi:10.3324/haematol.2012.074872

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol. Biol. 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, H., Yoon, D. H., Lee, J. B., Kim, S-Y., Moon, J. H., Do, Y. R., et al. (2017). Comprehensive evaluation of the revised international staging system in multiple myeloma patients treated with novel agents as a primary therapy. Am. J. Hematol. 92 (12), 1280–1286. doi:10.1002/ajh.24891

PubMed Abstract | CrossRef Full Text | Google Scholar

Corthay, A., Skovseth, D. K., Lundin, K. U., Røsjø, E., Omholt, H., Hofgaard, P. O., et al. (2005). Primary antitumor immune response mediated by CD4+ T cells. Immunity 22 (3), 371–383. doi:10.1016/j.immuni.2005.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, R., Huang, C., Liu, K., Li, X., and Dong, Z. (2021). Targeting AURKA in cancer: Molecular mechanisms and opportunities for cancer therapy. Mol. Cancer 20 (1), 15. doi:10.1186/s12943-020-01305-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Frerichs, K. A., Minnema, M. C., Levin, M-D., Broijl, A., Bos, G. M. J., Kersten, M. J., et al. (2021). Efficacy and safety of daratumumab combined with all-trans retinoic acid in relapsed/refractory multiple myeloma. Blood Adv. 5 (23), 5128–5139. doi:10.1182/bloodadvances.2021005220

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghasemi, M., Alpsoy, S., Türk, S., Üy, M., Atakan, Ş., Haznedaroğlu, İ. C., et al. (2016). Expression profiles of the individual genes corresponding to the genes generated by cytotoxicity experiments with bortezomib in multiple myeloma. Turk J. Haematol. 33 (4), 286–292. doi:10.4274/tjh.2015.0145

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonsalves, W. I., Jevremovic, D., Nandakumar, B., Dispenzieri, A., Buadi, F. K., Dingli, D., et al. (2020). Enhancing the R-ISS classification of newly diagnosed multiple myeloma by quantifying circulating clonal plasma cells. Am. J. Hematol. 95 (3), 310–315. doi:10.1002/ajh.25709

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jung, S-H., Kwon, S. Y., Min, J-J., Bom, H-S., Ahn, S-Y., Jung, S-Y., et al. (2019). 18F-FDG PET/CT is useful for determining survival outcomes of patients with multiple myeloma classified as stage II and III with the Revised International Staging System. Eur. J. Nucl. Med. Mol. Imaging 46 (1), 107–115. doi:10.1007/s00259-018-4114-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. H., Ho, J-N., Jin, H., Lee, S. C., Lee, S. E., Hong, S-K., et al. (2016). Upregulated expression of BCL2, MCM7, and CCNE1 indicate cisplatin-resistance in the set of two human bladder cancer cell lines: T24 cisplatin sensitive and T24R2 cisplatin resistant bladder cancer cell lines. Investig. Clin. Urol. 57 (1), 63–72. doi:10.4111/icu.2016.57.1.63

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreis, N. N., Louwen, F., and Yuan, J. (2015). Less understood issues: p21(Cip1) in mitosis and its therapeutic potential. Oncogene 34 (14), 1758–1767. doi:10.1038/onc.2014.133

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Peng, J., Zhou, Y., Xie, B., and Wang, J. (2019). Initiation and persistence with antiplatelet agents among the patients with acute coronary syndromes: A retrospective, observational database study in China. Mol. Med. Rep. 20 (3), 2159–2169. doi:10.2147/PPA.S228065

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Zhao, H., Fu, B., Jiang, S., Wang, J., and Wan, Y. (2022). Mapping cell phenomics with multiparametric flow cytometry assays. Phenomics 2 (4), 272–281. doi:10.1007/s43657-021-00031-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Lundin, K. U., Hofgaard, P. O., Omholt, H., Munthe, L. A., Corthay, A., and Bogen, B. (2003). Therapeutic effect of idiotype-specific CD4+ T cells against B-cell lymphoma in the absence of anti-idiotypic antibodies. Blood 102 (2), 605–612. doi:10.1182/blood-2002-11-3381

PubMed Abstract | CrossRef Full Text | Google Scholar

MacLeod, M. K. L., Kappler, J. W., and Marrack, P. (2010). Memory CD4 T cells: Generation, reactivation and re-assignment. Immunology 130 (1), 10–15. doi:10.1111/j.1365-2567.2010.03260.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mou, Y., Wang, J., Wu, J., He, D., Zhang, C., Duan, C., et al. (2019). Ferroptosis, a new form of cell death: Opportunities and challenges in cancer. J. Hematol. Oncol. 12 (1), 34. doi:10.1186/s13045-019-0720-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Nijhof, I. S., Groen, R. W. J., Lokhorst, H. M., van Kessel, B., Bloem, A. C., van Velzen, J., et al. (2015). Upregulation of CD38 expression on multiple myeloma cells by all-trans retinoic acid improves the efficacy of daratumumab. Leukemia 29 (10), 2039–2049. doi:10.1038/leu.2015.123

PubMed Abstract | CrossRef Full Text | Google Scholar

Palumbo, A., Avet-Loiseau, H., Oliva, S., Lokhorst, H. M., Goldschmidt, H., Rosinol, L., et al. (2015). Revised international staging system for multiple myeloma: A report from international myeloma working group. J. Clin. Oncol. 33 (26), 2863–2869. doi:10.1200/JCO.2015.61.2267

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandiella, A., Carvajal-Vergara, X., Tabera, S., Mateo, G., Gutiérrez, N., and San Miguel, J. F. (2003). Imatinib mesylate (STI571) inhibits multiple myeloma cell proliferation and potentiates the effect of common antimyeloma agents. Br. J. Haematol. 123 (5), 858–868. doi:10.1046/j.1365-2141.2003.04706.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, Z., and Blankenstein, T. (2000). CD4+ T cell--mediated tumor rejection involves inhibition of angiogenesis that is dependent on IFN gamma receptor expression by nonhematopoietic cells. Immunity 12 (6), 677–686. doi:10.1016/s1074-7613(00)80218-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Quezada, S. A., Simpson, T. R., Peggs, K. S., Merghoub, T., Vider, J., Fan, X., et al. (2010). Tumor-reactive CD4(+) T cells develop cytotoxic activity and eradicate large established melanoma after transfer into lymphopenic hosts. J. Exp. Med. 207 (3), 637–650. doi:10.1084/jem.20091918

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajkumar, S. V., Dimopoulos, M. A., Palumbo, A., Blade, J., Merlini, G., Mateos, M. V., et al. (2014). International Myeloma Working Group updated criteria for the diagnosis of multiple myeloma. Lancet Oncol. 15 (12), e538–e548. doi:10.1016/S1470-2045(14)70442-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Roninson, I. B. (2002). Oncogenic functions of tumour suppressor p21(waf1/cip1/sdi1): Association with cell senescence and tumour-promoting activities of stromal fibroblasts. Cancer Lett. 179 (1), 1–14. doi:10.1016/s0304-3835(01)00847-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Schietinger, A., Philip, M., Liu, R. B., Schreiber, K., and Schreiber, H. (2010). Bystander killing of cancer requires the cooperation of CD4(+) and CD8(+) T cells during the effector phase. J. Exp. Med. 207 (11), 2469–2477. doi:10.1084/jem.20092450

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimizu, Y., and Hendershot, L. M. (2009). Oxidative folding: Cellular strategies for dealing with the resultant equimolar production of reactive oxygen species. Antioxid. Redox Signal 11 (9), 2317–2331. doi:10.1089/ars.2009.2501

PubMed Abstract | CrossRef Full Text | Google Scholar

Sonneveld, P., Avet-Loiseau, H., Lonial, S., Usmani, S., Siegel, D., Anderson, K. C., et al. (2016). Treatment of multiple myeloma with high-risk cytogenetics: A consensus of the international myeloma working group. Blood 127 (24), 2955–2962. doi:10.1182/blood-2016-01-631200

PubMed Abstract | CrossRef Full Text | Google Scholar

Tu, B. P., and Weissman, J. S. (2004). Oxidative protein folding in eukaryotes: Mechanisms and consequences. J. Cell Biol. 164 (3), 341–346. doi:10.1083/jcb.200311055

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Geng, C., Zhou, H., Zhang, Z., and Chen, W. (2021). Cyclopamine sensitizes multiple myeloma cells to circularly permuted TRAIL-induced apoptosis. Oncol. Lett. 21 (4), 295. doi:10.3892/ol.2021.12556

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Lin, Z., Wang, Z., Ye, L., Xian, M., Xiao, L., et al. (2022). RARγ activation sensitizes human myeloma cells to carfilzomib treatment through the OAS-RNase L innate immune pathway. Blood 139 (1), 59–72. doi:10.1182/blood.2020009856

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, P., Yeung, Y. T., Wang, J., Wu, Q., Du, R., Huang, C., et al. (2021). miR-17-3p promotes the proliferation of multiple myeloma cells by downregulating P21 expression through LMLN inhibition. Int. J. Cancer 148 (12), 3071–3085. doi:10.1002/ijc.33528

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, Y., Hannon, G. J., Zhang, H., Casso, D., Kobayashi, R., and Beach, D. (1993). p21 is a universal inhibitor of cyclin kinases. Nature 366 (6456), 701–704. doi:10.1038/366701a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, B., Ai, Y., Sun, Q., Ma, Y., Cao, Y., Wang, J., et al. (2021). Membrane damage during ferroptosis is caused by oxidation of phospholipids catalyzed by the oxidoreductases POR and CYB5R1. Mol. Cell 81 (2), 355–369.e10. doi:10.1016/j.molcel.2020.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Ying, W. (2023). Phenomic studies on diseases: Potential and challenges. Phenomics 3 (3), 285–299. doi:10.1007/s43657-022-00089-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Z., Wang, H., Fang, Y., Lu, L., Li, M., Yan, B., et al. (2020). Molecular chaperone HspB2 inhibited pancreatic cancer cell proliferation via activating p53 downstream gene RPRM, Bai1, and TSAP6. J. Cell Biochem. 121 (3), 2318–2329. doi:10.1002/jcb.29455

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, Y., Tian, F., Ma, H., Wang, H., Yang, W., Liu, Z., et al. (2020). FTY720 induces ferroptosis and autophagy via PP2A/AMPK pathway in multiple myeloma cells. Life Sci. 260, 118077. doi:10.1016/j.lfs.2020.118077

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: multiple myeloma, ferroptosis, prognostic marker genes, overall survival, gene set enrichment analysis, drug sensitivity, immune analysis

Citation: Wang Q, Zhao M, Zhang T, Zhang B, Zheng Z, Lin Z, Zhou S, Zheng D, Chen Z, Zheng S, Zhang Y, Lin X, Dong R, Chen J, Qian H, Hu X, Zhuang Y, Zhang Q, Jiang S and Ma Y (2023) Comprehensive analysis of ferroptosis-related genes in immune infiltration and prognosis in multiple myeloma. Front. Pharmacol. 14:1203125. doi: 10.3389/fphar.2023.1203125

Received: 10 April 2023; Accepted: 24 July 2023;
Published: 07 August 2023.

Edited by:

Chen Ling, Fudan University, China

Reviewed by:

Ibrahim C. Haznedaroglu, Hacettepe University Hospital, Türkiye
Xiuli Wu, Jinan University, China
Mei Luo, Huazhong University of Science and Technology, China

Copyright © 2023 Wang, Zhao, Zhang, Zhang, Zheng, Lin, Zhou, Zheng, Chen, Zheng, Zhang, Lin, Dong, Chen, Qian, Hu, Zhuang, Zhang, Jiang and Ma. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Songfu Jiang, jiangsongfu@189.cn; Yongyong Ma, mayy@wmu.edu.cn

These authors have contributed equally to this work

Disclaimer: 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.