- 1The Second Clinical Medical College, Zhejiang Chinese Medical University, Hangzhou, China
- 2General Surgery, Cancer Center, Department of Breast Surgery, Zhejiang Provincial People’s Hospital (Affiliated People’s Hospital, Hangzhou Medical College), Hangzhou, China
- 3Laboratory Medicine Center, Department of Transfusion Medicine, Zhejiang Provincial People’s Hospital (Affiliated People’s Hospital, Hangzhou Medical College), Hangzhou, China
- 4Hangzhou Medical College, Hangzhou, China
- 5College of Pharmacy, Zhejiang University of Technology, Hangzhou, China
- 6Laboratory of Tumor Molecular Diagnosis and Individualized Medicine of Zhejiang Province, Zhejiang Provincial People’s Hospital (Affiliated People’s Hospital, Hangzhou Medical College), Hangzhou, China
Purpose: To identify molecular clusters associated with ferroptosis and to develop a ferroptosis-related signature for providing novel potential targets for the recurrence-free survival and treatment of breast cancer.
Methods: Ferroptosis-related gene (FRG) signature was constructed by univariate and multivariate Cox regression and least absolute shrinkage and selection operator (LASSO). Receiver operating characteristic curves, Kaplan–Meier survival analysis, principal component analysis, and univariate and multivariate Cox regression analyses in the training and test cohorts were used to evaluate the application of this signature. Quantitative reverse transcriptase–PCR (qRT-PCR) was employed to detect the expression of FRGs in the model. Furthermore, the correlations between the signature and immune microenvironment, somatic mutation, and chemotherapeutic drugs sensitivity were explored.
Results: Internal and external validations affirmed that relapse-free survival differed significantly between the high-risk and low-risk groups. Univariate and multivariate Cox regression analyses indicated that the riskScore was an independent prognostic factor for BRCA. The areas under the curve (AUCs) for predicting 1-, 2-, and 3-year survival in the training and test cohorts were satisfactory. Significant differences were also found in the immune microenvironment and IC50 of chemotherapeutic drugs between different risk groups. Furthermore, we divided patients into three clusters based on 18 FRGs to ameliorate the situation of immunotherapy failure in BRCA.
Conclusions: The FRG signature functions as a robust prognostic predictor of the immune microenvironment and therapeutic response, with great potential to guide individualized treatment strategies in the future.
Introduction
Breast cancer has surpassed lung cancer as being the most commonly diagnosed cancer with approximately 2.3 million new cases in 2020, accounting for 11.7% of all new cancer cases (1). Another scary truth is the drop in average onset age (2). Because breast cancer is a highly heterogeneous systemic disease, advancements in therapy are particularly crucial (3).
Neoadjuvant chemotherapy (NAC) is seen as the standard and first-line treatment for locally advanced breast cancer (4, 5), which not only is beneficial to breast-conserving surgery but also can detect tumor sensitivity to anticancer therapy for locally advanced breast cancer (6), and it could also be employed as a bridge to other therapies (7, 8). Anthracyclines and taxanes serve as the backbone of NAC regimens and are widely used clinically (9).
Ferroptosis is an emerging form of programmed cell death featured by the iron-dependent accumulation of lipid reactive oxygen species (ROS) of metabolic dysfunctions, iron accumulation, and antioxidant vulnerability (10–12). Accumulating evidence showed that the role of ferroptosis in carcinogenesis, progression, and chemoresistance had made progress. Fascin regulates SLC7A11 stability to induce ferroptosis (13). Renovation of SLC7A11 rescues miR-5096-mediated ferroptosis and antitumor effects of breast cancer (14). Ferroptosis-related gene (FRG) GPX4 promotes chemoresistance in nasopharyngeal carcinoma (15). Bufotalin induces ferroptosis by facilitating the ubiquitination and degradation of GPX4 in non-small cell lung cancer cells (16). Via ferroptosis, ETS1/miR-23a-3p/ACSL4 axis stimulates sorafenib resistance in HCC (17).
Ferroptosis has the characteristics of inhibiting chemoresistance and enhancing antitumor immunity (18), which may be a potential strategy to overcome the drug resistance mechanism of traditional cancer treatments (12). Previous studies prove the feasibility of ferroptosis-related prognostic markers to predict overall survival and immune characteristics. FRG signatures were constructed to predict overall survival in lung adenocarcinoma (19), colorectal cancer (20), and pancreatic adenocarcinoma (21). However, as far as we are aware, studies focusing on the correlation of ferroptosis with biochemical recurrence and antitumor immunology of BRCA were rather limited. Thus, it is an urgent need to discover a robust biomarker to predict relapse-free survival (RFS) in BRCA.
In this study, we constructed an FRG prognostic signature and identified three ferrClusters in predicting the RFS internally and externally, exploring the status of immune infiltrates and drug sensitivity of BRCA patients receiving NAC for guiding clinical practice. This signature may also serve as a novel and robust prediction tool for evaluating whether BRCA patients can benefit from immunotherapy.
Methods
Data Acquisition and Processing
Open expression matrix of mRNA (FPKM values) and clinical files of BRCA samples were downloaded from The Cancer Genome Atlas (TCGA) database. Datasets GSE25055 in the Gene Expression Omnibus (GEO) database were used to acquire RNA-sequencing (RNA-Seq) and clinical data of BRCA patients receiving NAC as a training cohort and GSE16446 and GSE25065 as test cohorts. Gene expression file of GSE25055 and GSE25065 was collected using platform GPL96 [HG-U133A] Affymetrix Human Genome U133A Array, and GSE16446 using platform GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. Batch effects and other unwanted variations in high-throughput experiments were eliminated using the “combat” function in the “sva” package (22) in R 4.1.1. Copy number variation (CNV) data were collected from the University of California, Santa Cruz (UCSC) website.
Construction of the Ferroptosis-Related Signature for Predicting Recurrence-Free Survival
FRGs including 150 drivers, 109 suppressors, and 123 markers were collected from FerrDb (19, 23, 24). A univariate Cox proportional hazards regression analysis was conducted to filtrate prognostic FRGs in the GSE25055 cohort with p < 0.05 considered to be statistically significant using the “coxph” function. Subsequently, with the help of the “cv.glmnet” function, the least absolute shrinkage and selection operator (LASSO) was performed for the dimension reduction and K-fold cross-validation, which was multiplied by ten, and the optimal parameter was the λ value that corresponded to the lowest deviation. The optimal penalty parameter was defined as the value within one SD of the minimum cross-validated partial likelihood deviance to obtain the best model. The proteins with non-zero regression coefficients were chosen for subsequent multivariate Cox regression analyses. The LASSO regression model was as follows:
Validation of the Prognostic Signature
First, in the GSE25055 dataset, the Kaplan–Meier (K-M) survival analysis using the “Surv” function in the “survival” package and univariate and multivariate Cox regression analyses between gene expression and clinical characters using the “coxph” function in the “survival” package were performed to certify that riskScore served as an independent predictor in predicting recurrence-free survival (RFS). Principal component analysis (PCA) using the “prcomp” function was used to visualize sample distribution. Receiver operating characteristic (ROC) using the “timeROC” package was done, and area under the curve (AUC) plots were generated for the 1-year, 2-year, and 3-year survival rates to assess the sensitivity and specificity of the prognostic model. Then, the prognostic signature was validated in the GSE25065 and GSE16446 datasets via the above methods.
Cell Culture
Normal breast epithelial cell line MCF-10A and the epithelial BRCA cell lines MCF-7, T47D, MDA-MB-231, MDA-MB-468, and BT-549 were acquired from the American Type and Culture Collection (ATCC; Manassas, VA, USA). MDA-MB-231 and BT-549 cells were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM) (ATCC; Manassas, VA, USA) supplemented with 10% fetal bovine serum (HyClone, Logan, UT, USA) and 1% antibiotic (100 IU/ml of penicillin and 100 µg/ml of streptomycin; HyClone, Logan, UT, USA). MCF-10A cells were cultured in DMEM/F12 medium supplemented with 20 ng/μl of epidermal growth factor, insulin, hydrocortisone, non-essential amino acid (NEAA), 5% horse serum (HS), and 1% penicillin/streptomycin (P/S) solution (Procell, Wuhan, China). MCF-7 and MDA-MB-468 were cultured in Minimum Essential Medium (MEM) (Gibco BRL, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (HyClone, Logan, UT, USA) and 1% antibiotic (100 IU/ml of penicillin and 100 µg/ml of streptomycin; HyClone, Logan, UT, USA). T-47D cells were cultured in Roswell Park Memorial Institute (RPMI) 1640 (HyClone, Logan, UT, USA) with 10% fetal bovine serum (HyClone, Logan, UT, USA). All the cell lines were incubated at 37°C, with a humidified atmosphere of 5% CO2.
Quantitative Reverse Transcriptase–PCR
Total RNAs were isolated from cells using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). PrimeScript™ RT reagent Kit (Takara, Maebashi, Japan) was employed to reverse transcribe into cDNA following the manufacturer’s protocol. Then SYBR Green PCRMaster Mix (Applied TaKaRa, Otsu, Japan) was used to conduct Real-time PCR on Applied Biosystems 7500 Fast Real-Time RCR System (Applied Biosystems, Foster City, CA, USA). The primers of FRGs for qRT-PCR utilized in this research were as follows:
Drug Sensitive Analysis
With the use of the “pRRophetic” package, the half-maximal inhibitory concentration (IC50) of BRCA patients was calculated on Genomics of Drug Sensitivity in Cancer (GDSC) (25) (https://www.cancerrxgene.org/) based on the given gene expression profiles in these datasets to evaluate the drug sensitivities (26–28).
Immune Infiltration Analyses
The CIBERSORT algorithm was used to explore the proportion of different types of immune cells in BRCA patients using CIBERSORT R script v1.04 (29–31). Based on the expression level of immune cell-related genes, the ESTIMATE algorithm was conducted to calculate the stromal score (SS), estimate score (ES), and immune score (IS) (the SS represents the level of stroma content in a tumor; the IS reflects the infiltration of immune cells in a tumor; the estimated score infers tumor purity) among the high- and low-risk groups using the “estimate” package (32). Single-sample gene set enrichment analysis (ssGSEA) was performed to calculate scores for antitumor immunity and protumor suppression for each sample (33) using “GSEABase” and “GSVA” packages.
Consensus Clustering Analyses for Identifying BRCA Subtypes
Consensus clustering based on Euclidean distance and Ward’s linkage was performed for hierarchical clustering to identify different subtypes using the “ConsensusClusterPlus” package and repeated the procedures 1,000 times to guarantee the stability of the classification (34). In consideration of a high consistency of clusters, a low coefficient of variation, and no significant increase in the CDF curve, the optimum cluster number could be determined (35).
Statistical Analysis
Correlation coefficients were calculated by Spearman’s and distance correlation analyses. For comparison of more than two groups, the Kruskal–Wallis and one-way ANOVAs were chosen as non-parametric and parametric methods, while Wilcoxon’s t-test was used for two groups. Student’s t-test was used to explore the statistical significance of quantitative data. The K-M and log-rank tests were employed to confirm the significance of prognostic differences (22). R 4.1.1 software was the main tool to conduct the statistical analysis. For all statistical results, a p-value of <0.05 was considered to be statistically significant.
Results
Construction of the Ferroptosis-Related Signature Associated With Recurrence-Free Survival
GSE25055 dataset was used as a training cohort; meanwhile, GSE25065 and GSE16446 datasets were used as test cohorts. Batch effects were removed for further study (Figures 1A, B). First, we performed a univariate Cox regression analysis in GSE25055. Among 382 FRGs retrieved from the FerrDB database, 76 FRGs were identified to be associated with RFS, with the standard of p < 0.01 (Figure 1C). Pearson’s correlation analysis revealed a correlation among these genes (Figure 1D). Then, LASSO regression analysis was used to establish the FRGs prognostic signature (Figures 1E, F):
Figure 1 RNA-sequencing (RNA-Seq) data of training and test cohorts before (A) and after (B) removing batch effects. (C) The hazard ratio (HR) and p-value of selected ferroptosis-related genes (FRGs) using the univariable Cox HR regression (criteria: p-value <0.01). (D) Expression interaction of the 76 FRGs in BRCA. The lines connecting the FRGs show how they are correlated with each other, with positive associations in red and negative associations in green. (E) The least absolute shrinkage and selection operator (LASSO) Cox analysis identified 18 FRGs most related to prognostics. (F) The 10-round cross-validation determined the optimal values of the penalty parameter. (G) In all, 40 of 983 (4.07%) BRCA patients experienced 18 FRG genetic alterations. (H) Copy number variation (CNV) mutation frequency of the 18 FRGs. This column represents the frequency of change. Deletion frequency is represented by green dots, while amplification frequency is represented by pink dots. (I) The location of the 18 FRGs in chromosomes. Blue point represents the genes that mainly had CNV deletion; red point represents the genes that mainly had CNV amplification. (J) Expression of the 18 FRGs in normal tissues and BRCA tissues. Genes with red color represent the differentially expressed genes. (K) The value of logFC of the 18 FRG genes. (L–N) qRT-PCR results showed the expression value of the three FRGs in the normal breast and five breast cancer cell lines. *, P < 0.05; **, P < 0.01; ***, P < 0.001.
Landscape of Gene Mutations and Expression in Ferroptosis-Related Genes in the Model in BRCA
Genomic mutations were common in these genes with 40 (4.07%) of 983 patients having experienced genetic changes, and a mutation frequency of 1% was observed in ENPP2 and NEDD4L (Figure 1G). We also found that CNV is prevalent among the 18 FRGs. ENPP2, CHMP6, MAFG, VEGFA, VDAC2, LPCAT3, CYP4F8, SLC1A4, XBP1, PIR, and ACO1 showed copy number amplification, while deletion happened in the other FRGs (Figure 1H). The location of the 18 FRGs in human chromosomes could be seen in Figure 1I. The result of differential analysis in normal breast tissue and tumor tissue showed that ACO1, CHMP6, ENPP2, MAFG, NEDD4L, PIR, SLC1A4, SLC7A5, and VEGFA had significant differential expression in breast cancer with p-value <0.001; ACADSB and NOX3 with p-value <0.01; and PEX12 with p-value <0.05 (Figure 1J). SLC7A5 was seen as a significantly upregulated gene, while ENPP2 and ACO1 were seen as significantly downregulated genes with |log FC| > 1 (Figure 1K). The result of RT-PCR provided strong support for our conclusion (Figures 1L–N). As described above, FRGs had significant heterogeneity of genetic variation and transcriptomic alteration landscape in BRCA patients, which played an important part in regulating the happening, aggravation, and prognosis of BRCA.
External Validation of the Ferroptosis-Related Gene Model
After the riskScore of each patient based on the risk model was calculated, we divided patients into the high- and low-risk groups with the standard median score in GSE25055 (Figure 2A). With the use of the median of GSE25055, patients in GSE25065 (Figure 2H) and GSE16446 (Figures 2H, O) were separated into the high-risk and low-risk groups in the same manner. The result of PCA showed significant heterogeneity between high-risk and low-risk patients in GSE25055 (Figures 2B, C), GSE25065 (Figures 2I, J), and GSE16446 (Figures 2P, Q), which certified the superior discrimination of the FRG model. For the purpose of exploring whether the signature could represent its prognostic value independently of other clinical factors, we conducted univariate and multivariate Cox regression analyses in the training and test cohorts. In univariate analyses, this risk score was able to independently predict survival outcomes in GEO cohorts (GSE25055, hazard ratio (HR) = 4.690, p < 0.001; GSE25065, HR = 6.350, p < 0.001; GSE16446, HR = 7.648, p < 0.001) (Figures 2D, K, R). The same conclusion could be drawn in multivariate analyses (Figures 2E, L, S). The results revealed that riskScore and pathologic response served as independent factors affecting receiving NAC BRCA patients’ prognosis. The AUCs of the time-dependent ROC curves at 1, 2, and 3 years were 0.818, 0.824, and 0.783 in GSE25055 (Figure 2F); 0.812, 0.824, and 0.783 in GSE25065 (Figure 2M); and 0.715, 0.725, and 0.723 in GSE16446 (Figure 2T). The AUCs in different years and cohorts were relatively high compared with those of other published literature, which suggested high sensitivity and specificity of the signature for predicting RFS. The K-M survival curve showed that patients in the high-risk group had a higher recurrence rate than those in the low-risk group using log-rank tests with p < 0.001 (Figure 2G), p < 0.001(Figure 2N), p = 0.029 (Figure 2U). Ferroptosis is a recently recognized form of regulated cell death that is characterized by lipid peroxidation, which mediates cell death in breast cancer. Among genes in our signature, ferroptosis driver genes such as NOX3 and PEBP1 had negative coefficients, while ferroptosis suppressor genes such as PIR and VDAC2 had positive coefficients. Therefore, high riskScore indicated that ferroptosis was suppressed in breast cancer, which might imply a worse prognosis.
Figure 2 Distribution of riskScore, scatterplot, and heatmap in the high-risk group and the low-risk group in GSE25055 (A), GSE25065 (H), and GSE16446 (O). PCA (principal component analysis) for BRCA based on the riskScore in GSE25055 (B, C), GSE25065 (I, J), and GSE16446 (P, Q). (D) Univariate and (E) multivariate Cox regression analyses of age, estrogen receptor (ER), progesterone receptor (PR), T, N, stage, grade, pathologic response, pam50 classification, and riskScore in GSE25055. (K) Univariate and (L) multivariate Cox regression analyses of age, ER, PR, T, N, stage, grade, pathologic response, pam50 classification, and riskScore in GSE25065. (R) Univariate and (S) multivariate Cox regression analyses of age, T, N, grade, her2-fish, top2atri, erbb2, and riskScore in GSE16446. Time-dependent receiver operating characteristic (ROC) curves for predicting 1-, 2-, and 3-year RFS in GSE25055 (F), GSE25065 (M), and GSE16446 (T). Kaplan–Meier curves of the high- and low-risk subgroup patients in GSE25055 (G), GSE25065 (N), and GSE16446 (U).
Clinicopathological Parameter Relevance Analysis
We further anatomized the association between riskScore and clinical parameters of BRCA patients. The detailed results depicted that the riskScore had a positive correlation with T stage, N stage, American Joint Committee on Cancer (AJCC) stage, and grade (Figure 3). BRCA patients with higher T, N, AJCC stage, and grade, combined with lower age, and negative status of progesterone receptor and estrogen receptor seemed to have higher riskScore, indicating a higher incidence rate of relapse, which was consistent with the conclusions of current accumulated literature. In other words, the results implied that the riskScore had a correlation with clinicopathological parameters.
Figure 3 RiskScore is correlated with clinicopathological features of BRCA. T stage (A), N stage (B), AJCC-stage (C), Grade (D), age (E), ER status (F), PR status (G), PAM50 subtypes (H). NA, P>0.05.
Chemotherapeutic Response Analysis
In order to improve the therapeutic benefit of BRCA patients from neoadjuvant therapy, we further explored whether FRG signature could predict the sensitivity to several chemotherapy drugs widely used in BRCA between two groups. According to the results calculated based on the GDSC database, IC50 values of chemotherapy drugs covering axitinib, bicalutamide, bleomycin, bortezomib, dasatinib, doxorubicin, gefitinib, lapatinib, and paclitaxel were evaluated. Compared with the low-risk group, IC50 values of paclitaxel, gefitinib, doxorubicin, bleomycin, and bortezomib were lower in the high-risk groups, which indicated that high-risk patients were more sensitive to these drugs (Figures 4A–I). The above results demonstrated that the riskScore had potential predictive value for chemotherapy and targeted therapy in breast cancer.
Figure 4 Correlation between ferroptosis-related gene (FRG) signature and drug sensitivity. Box plots for estimated IC50 of drugs between high- and low-risk BRCA patients. Paclitaxel (A), Lapatinib (B), Gefitinib (C), Dasatinib (D), Doxorubicin (E), Bleomycin (F), Bicalutamide (G), Bortezomib (H), Axitinib (I).
Comprehensive Analysis Between Ferroptosis-Related Gene Signature and Immune Microenvironment
We calculated the constitution of tumor-infiltrating immune cells in BRCA through the CIBERSORT algorithm (Figure 5A). Compared with the low-risk groups, the proportion of resting mast cells was lower in the high-risk groups (Figures 5B, C).
Figure 5 (A) Relative percent of different immune cells in each sample. (B) Different immune cell contents in low-risk and high-risk patients. (C) Correlation between immune cells and riskScore. (D–F) Kaplan–Meier curves of the high and low stromal score (SS), immune score (IS), and estimate score (ES) group patients. (G–L) Correlation between riskScore and SS, IS, and ES. (M, N) Comparisons of the expression levels of immune checkpoints between two groups. ns, P>0.05. *, P < 0.05; **, P < 0.01; ***, P < 0.001.
Then, the IS, SS, and ES of patients were evaluated using the ESTIMATE algorithm. Based on the optimum cutoff value of ISs or SSs respectively, BRCA patients were divided into the high and low IS/SS/ES groups. The K-M curves showed that patients with high IS/SS/ES exhibited significantly worse RFS as compared to the ones with low IS/SS/ES (Figures 5D–F). We further explored the relationships between the IS/SS/ES and riskScore. The result of Wilcoxon’s rank-sum test displayed that there is no significant difference between the high-risk and low-risk groups in SS (p = 0.53, Figure 5G) but significant in IS (p = 5.7e−09, Figure 5H) and ES (p = 1.9e−05, Figure 5I). Pearson’s correlation analysis showed that riskScore was positively associated with IS (R = 0.26, p = 1.6e−10, Figure 5K) and ES (R = 0.2, p = 5.1e−07, Figure 5L). However, the riskScore was not significantly correlated to the SS (P = 0.14, Figure 5J)
Furthermore, we dissected the role of riskScore in immune checkpoint blockade (ICB) treatment. We noticed that the expression levels of all immune checkpoints were significantly higher in the high riskScore group (Figures 5M, N). Taken together, the prognostic signature could predict the potential response to immunotherapy in BRCA patients, which provided guidance on whether or what to use for immunotherapy in clinical practice.
Identification of Three Consensus Clustering Subtypes
On the basis of the expression of 18 FRGs in the signature, we employed the “Partition Around Medoids” algorithm, along with Pearson’s distance to estimate similarity among patients to identify three clusters. We noticed that K = 3 seemed to be an optimal selection by clustering variable (k) increasing from 2 to 9, in which the intergroup correlations were the lowest and the intragroup correlations were the greatest (Figure 6B), indicating the optimal clustering stability of the three molecular phenotypes. The consensus cumulative distribution function (CDF) diagram showed that when k = 3, distribution reached an approximate maximum (Figure 6C), implying robust clustering for all samples (Figure 6A). The delta area plot depicts the relative change compared to k − 1 showing that the delta area was optimum when k = 3 (Figure 6D). Prognostic analysis of the three clusters revealed that patients in ferrCluster A were the least likely to relapse, while in ferrCluster B, they were the most likely to relapse (Figure 6E).
Figure 6 Consensus clustering of 18 ferroptosis-related genes (FRGs) identified three clusters of patients. (A) The tracking plot for k = 2 to k = 9. (B) The heatmap for K = 3. (C) Consensus clustering cumulative distribution function (CDF) with k = 2 to k = 9. (D) Relative change in area under CDF curve for k = 2–9. (E) Kaplan–Meier (K-M) curve of the survival difference among clusters 1–3. (F) Single-sample gene set enrichment analysis of immune status among three ferrClusters. ns, P>0.05. *, P < 0.05; **, P < 0.01; ***, P < 0.001.
We then performed an ssGSEA to quantify the scores of various immune cell subpopulations to further compare the differences in the number of immune cells among the three types of ferrClusters. The results indicated that the contents of monocyte cells were not significantly different. The proportion of immune cells was significantly different among the three clusters. Contents of nearly all types of immune cells in ferrCluster A seemed to be the poorest. The levels of activated CD4+ T cells, CD8+ T cells, dendritic cells, CD56 bright and dim NK cells, γδ-T cells, Tregs, and T helper cells were relatively the highest in ferrCluster B. Hence, we could draw the conclusion that ferrCluster A was a type of immune failure, ferrCluster B was a type of immune-activated characterized by T-cell subset enrichment, and ferrCluster C was a type of immune-activated characterized by B-cell subset enrichment (Figure 6F). These results indicated that the FRGs play key roles in immune cell infiltration and characteristic tumor immune microenvironment (TME) formation and affect the prognosis of BRCA patients.
Development of ferrScore to Quantify Individual Ferroptosis Pattern
With a view to the individual heterogeneity and complexity of BRCA patients, we calculated ferrScore to assess the ferroptosis pattern of each patient based on the PCA on the 18 FRGs in the model. The scoring framework was defined as ferrScore = PC1 + PC2 to quantify individual ferroptosis patterns of BRCA patients (36), further facilitating precise treatment. As indicated from the K-M curve, patients with lower ferrScore had a lower probability of relapse (Figure 7A). The ferrScore was closely related to immune cells (Figure 7B). We also observed that ferrScores of patients in ferrCluster A were significantly lower than those in ferrCluster B and C, while there was no significant difference between ferrClusters B and C (Figure 7D). The Sankey diagram shows the attribute changes in riskScore, ferrCluster, ferrScore, and recurrence status, indicating that the higher the riskScore and ferrScore, the higher the risk of relapse after receiving NAC (Figure 7C). The above results enriched treatment strategies for BRCA patients not only in targeted therapy and chemotherapy but also in immunotherapy. At last, the expression of CTLA4 was examined to elucidate a potential response to immunotherapy, and the high ferrScore group showed relatively high levels of expression (p = 8.6e−11, Figure 7B).
Figure 7 (A) Kaplan–Meier (K-M) curve of the survival difference between high and low ferrScore groups. (B) Correlation between immune cells and ferrScore. (C) Alluvial diagram of riskScore group, ferrCluster group, ferrScore group, and relapse-free status. (D) Correlation between ferrCluster and ferrScore *, P < 0.05.
Discussion
Ferroptosis is a newfound programmed cell death pattern distinguished from traditional cell death such as apoptosis, necrosis, and autophagy (37). Accumulating evidence demonstrated that dysregulated expression and genetic variations of FRGs were closely related to cell death, tumor carcinogenesis, and progression (22, 38).
TME is a cradle for tumorigenesis and cancer progression, in which immune infiltrating cells affect therapeutic outcomes (39). The relationships between TME infiltration immune cells and ferroptosis modifications have become a hotspot in the mechanism of tumorigenesis and development (40, 41). MIF secreted by nasopharyngeal carcinoma could suppress ferroptosis of macrophages and then increase the rate of metastasis (42). BEBT-908 induces immunogenic ferroptosis to potentiate cancer immune checkpoint therapy (43). SCD1 and FABP4 could drive ferroptosis, thereby leading to tumor resistance (44). Ferroptotic cells could also release chemotaxis to interact with immune cells, such as CD + T cells, and then modulate the anticancer immunity (45).
High-throughput genomic studies provided cutting-edge sights into the molecular mechanisms and identified new potential targets of breast cancer. Our research developed and verified a stepwise multivariate Cox regression model including 18 FRGs using LASSO and multivariate Cox regression for removing redundant factors to forecast the RFS of individual patients in GSE25055. The expression of FRGs in the signature was higher in BRCA tissues than in adjacent normal tissues, which was verified in several breast cancer cell lines using real-time PCR. Meanwhile, CNVs and mutation frequencies of FRGs were prevalent. Internal and external validations exhibited an excellent ability to predict the prognosis of BRCA patients. Specifically, a higher riskScore indicated a higher rate of recurrence. Moreover, riskScore was associated closely with clinicopathological features.
With a view to the significance of the immune system in antiviral and antitumor responses, we calculated the proportion of different types of tumor-infiltrating immune cells in BRCA using CIBERSORT and used ESTIMATE to explore IS, SS, and tumor purity. Higher SSs and ISs were observed in high-risk patients, leading to an unfavorable prognosis, which was consistent with a line of evidence from previous research (46, 47).
Extensive interest in cancer immunotherapy is reported according to the clinical importance of CTLA-4 and PD-1/PD-L1 in immune checkpoint therapies (48). The main immune checkpoints for breast cancer include CTLA-4, PD-1/PD-L1, lymphocyte activation gene 3 (LAG-3), T-cell immunoglobulin domain and mucin 3 (TIM-3), and other molecules (49). Clinical trials like SOLTI-1503 PROMETEO TRIAL (50), KEYNOTE-086 (51), NIMBUS (52), KEYNOTE-173 (53), and KEYNOTE-522 (54) showed that immunological checkpoint inhibitors have made significant progress in breast cancer immunotherapy, which is expected to become a new treatment for breast cancer.
Furthermore, for the purpose of exploring the response to chemotherapy sensitivity of patients, we calculated the IC50 value. The sensitivities of chemotherapeutic drugs widely used in BRCA showed a significant difference between the two groups.
In accordance with the expression matrix of the 18 FRGs in the signature, we identified three ferroptosis-related molecular clusters via consensus clustering analysis. The rate of relapse was significantly different among the three clusters. ssGSEA identified that the three ferrClusters as three immune types of immune failure, immune-activated characterized by T-cell subset enrichment, and immune-activated characterized by B-cell subset enrichment.
Inevitably, numerous limitations of our study should be included in the consideration. First, although our conclusion came through internal and external validation in TCGA, GSE25055, GSE25065, and GSE16446 cohorts, when it comes to its clinical application, caution is advised. Multicenter large-scale prospective clinical studies were needed rather than only retrospective data from public open databases to verify the signature. Second, the expression matrix of patients in GSE25055 and GSE25065 was extracted via platform GPL96 [HG-U133A] Affymetrix Human Genome U133A Array in 2010, which only included 12,549 genes, while GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array for GSE16446 contained 21,655 genes. Due to the relatively small number of detectable genes, bias may be amplified. Finally, detailed molecular mechanisms in the BRCA of the FRGs in the signature had not been fully revealed. Further in-depth studies were required to confirm relationships between FRGs and tumor microenvironment, and between ferroptosis and chemoresistance.
Conclusion
In brief, we constructed a novel FRG signature and identified three molecular subtypes for predicting the RFS of BRCA patients, which could predict the immune status of the tumor microenvironment and RFS of patients. It is worth noting that our conclusions provided more clues for the rational choices of chemotherapeutic drugs for patients with BRCA, provided a new immunological perspective and a new basis for immunotherapy of BRCA in the clinic, and had the potential possibility to coach and guide individualized healthcare decisions.
Data Availability Statement
The expression and clinical datasets presented in this study are available in databases of TCGA and Gene Expression Omnibus (GSE25055, GSE25065 and GSE16446), and the copy number variation data are presented in the University of California, Santa Cruz (UCSC) website https://xenabrowser.net/datapages/?dataset=TCGA-BRCA.gistic.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443.
Author Contributions
XM and QX substantially contributed to the conception of the work. YX, QZ, and TZ contributed to the data collection. YX, YD, and QX wrote the manuscript. QZ and YD helped to perform the enrichment and network analysis. YX, YD, and XM drafted and revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This research was supported by the Zhejiang Provincial Natural Science Foundation of China (Grant No. LQ21H200007), Zhejiang Provincial Ministry Medical and Health Co-construction Major Project (Grant No. 20214355173), and National Natural Science Foundation of China (Grant No. 81973861).
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.895110/full#supplementary-material
Supplementary Table 1 | Ferroptosis suppressor, driver and marker genes from FerrDB database.
Abbreviations
GEO, Gene Expression Omnibus; TCGA, The Cancer Genome Atlas; NAC, neoadjuvant chemotherapy; CNV, copy number variation; UCSC, University of California, Santa Cruz; ROC, receiver operating characteristic; GDSC, Genomics of Drug Sensitivity in Cancer; IC50, half-maximal inhibitory concentration; AUC, area under the curve; RFS, recurrence-free survival; FRG, ferroptosis-related gene; IS, immune score; SS, stromal score; ES, estimate score; ssGSEA, single-sample gene set enrichment analysis; TIME, tumor immune microenvironment; PCA, principal component analysis; ICI, immune checkpoint inhibitor.
References
1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660
2. Isik A, Firat D. Bilateral Intra-Areolar Polythelia. Breast J (2018) 24(1):89–90. doi: 10.1111/tbj.12838
3. Adams S, Gatti-Mays ME, Kalinsky K, Korde LA, Sharon E, Amiri-Kordestani L, et al. Current Landscape of Immunotherapy in Breast Cancer: A Review. JAMA Oncol (2019) 5(8):1205–14. doi: 10.1001/jamaoncol.2018.7147
4. Early Breast Cancer Trialists' Collaborative G. Long-Term Outcomes for Neoadjuvant Versus Adjuvant Chemotherapy in Early Breast Cancer: Meta-Analysis of Individual Patient Data From Ten Randomised Trials. Lancet Oncol (2018) 19(1):27–39. doi: 10.1016/S1470-2045(17)30777-5
5. Cortazar P, Zhang L, Untch M, Mehta K, Costantino JP, Wolmark N, et al. Pathological Complete Response and Long-Term Clinical Benefit in Breast Cancer: The CTNeoBC Pooled Analysis. Lancet (2014) 384(9938):164–72. doi: 10.1016/S0140-6736(13)62422-8
6. Hamy AS, Pierga JY, Sabaila A, Laas E, Bonsang-Kitzis H, Laurent C, et al. Stromal Lymphocyte Infiltration After Neoadjuvant Chemotherapy Is Associated With Aggressive Residual Disease and Lower Disease-Free Survival in HER2-Positive Breast Cancer. Ann Oncol (2017) 28(9):2233–40. doi: 10.1093/annonc/mdx309
7. Kuijer A, Straver M, den Dekker B, van Bommel ACM, Elias SG, Smorenburg CH, et al. Impact of 70-Gene Signature Use on Adjuvant Chemotherapy Decisions in Patients With Estrogen Receptor-Positive Early Breast Cancer: Results of a Prospective Cohort Study. J Clin Oncol (2017) 35(24):2814–9. doi: 10.1200/JCO.2016.70.3959
8. Harbeck N, Gnant M. Breast Cancer. Lancet (2017) 389(10074):1134–50. doi: 10.1016/S0140-6736(16)31891-8
9. Chaudhary LN. Early Stage Triple Negative Breast Cancer: Management and Future Directions. Semin Oncol (2020) 47(4):201–8. doi: 10.1053/j.seminoncol.2020.05.006
10. Xie Y, Hou W, Song X, Yu Y, Huang J, Sun X, et al. Ferroptosis: Process and Function. Cell Death Differ (2016) 23(3):369–79. doi: 10.1038/cdd.2015.158
11. Stockwell BR, Friedmann Angeli JP, Bayir H, Bush AI, Conrad M, Dixon SJ, et al. Ferroptosis: A Regulated Cell Death Nexus Linking Metabolism, Redox Biology, and Disease. Cell (2017) 171(2):273–85. doi: 10.1016/j.cell.2017.09.021
12. Sun LL, Linghu DL, Hung MC. Ferroptosis: A Promising Target for Cancer Immunotherapy. Am J Cancer Res (2021) 11(12):5856–63.
13. Chen C, Xie B, Li Z, Chen L, Chen Y, Zhou J, et al. Fascin Enhances the Vulnerability of Breast Cancer to Erastin-Induced Ferroptosis. Cell Death Dis (2022) 13(2):150. doi: 10.1038/s41419-022-04579-1
14. Yadav P, Sharma P, Sundaram S, Venkatraman G, Bera AK, Karunagaran D. SLC7A11/ xCT is a Target of miR-5096 and its Restoration Partially Rescues miR-5096-Mediated Ferroptosis and Anti-Tumor Effects in Human Breast Cancer Cells. Cancer Lett (2021) 522:211–24. doi: 10.1016/j.canlet.2021.09.033
15. Yuan L, Li S, Chen Q, Xia T, Luo D, Li L, et al. EBV Infection-Induced GPX4 Promotes Chemoresistance and Tumor Progression in Nasopharyngeal Carcinoma. Cell Death Differ (2022). doi: 10.1038/s41418-022-00939-8
16. Zhang W, Jiang B, Liu Y, Xu L, Wan M. Bufotalin Induces Ferroptosis in non-Small Cell Lung Cancer Cells by Facilitating the Ubiquitination and Degradation of GPX4. Free Radic Biol Med (2022) 180:75–84. doi: 10.1016/j.freeradbiomed.2022.01.009
17. Lu Y, Chan YT, Tan HY, Zhang C, Guo W, Xu Y, et al. Epigenetic Regulation of Ferroptosis via ETS1/miR-23a-3p/ACSL4 Axis Mediates Sorafenib Resistance in Human Hepatocellular Carcinoma. J Exp Clin Cancer Res (2022) 41(1):3. doi: 10.1186/s13046-021-02208-x
18. Lin HY, Ho HW, Chang YH, Wei CJ, Chu PY. The Evolving Role of Ferroptosis in Breast Cancer: Translational Implications Present and Future. Cancers (Basel) (2021) 13(18):4576. doi: 10.3390/cancers13184576
19. Wang S, Xie Z, Wu Z. Establishment and Validation of a Ferroptosis-Related Gene Signature to Predict Overall Survival in Lung Adenocarcinoma. Front Genet (2021) 12:793636:793636. doi: 10.3389/fgene.2021.793636
20. Du S, Zeng F, Sun H, Liu Y, Han P, Zhang B, et al. Prognostic and Therapeutic Significance of a Novel Ferroptosis Related Signature in Colorectal Cancer Patients. Bioengineered (2022) 13(2):2498–512. doi: 10.1080/21655979.2021.2017627
21. Yang J, Wei X, Hu F, Dong W, Sun L. Development and Validation of a Novel 3-Gene Prognostic Model for Pancreatic Adenocarcinoma Based on Ferroptosis-Related Genes. Cancer Cell Int (2022) 22(1):21. doi: 10.1186/s12935-021-02431-8
22. Liu DL, Wu MY, Zhang TN, Wang CG. Ferroptosis Regulator Modification Patterns and Tumor Microenvironment Immune Infiltration Characterization in Hepatocellular Carcinoma. Front Mol Biosci (2022) 9:807502:807502. doi: 10.3389/fmolb.2022.807502
23. Ke ZB, You Q, Sun JB, Zhu JM, Li XD, Chen DN, et al. A Novel Ferroptosis-Based Molecular Signature Associated With Biochemical Recurrence-Free Survival and Tumor Immune Microenvironment of Prostate Cancer. Front Cell Dev Biol (2021) 9:774625:774625. doi: 10.3389/fcell.2021.774625
24. Li X, Huang J, Chen J, Zhan Y, Zhang R, Lu E, et al. A Novel Prognostic Signature Based on Ferroptosis-Related Genes Predicts the Prognosis of Patients With Advanced Bladder Urothelial Carcinoma. Front Oncol (2021) 11:726486:726486. doi: 10.3389/fonc.2021.726486
25. Chen Z, Zou Y, Zhang Y, Chen Z, Wu F, Jin H, et al. A Pyroptosis-Based Prognostic Model for Immune Microenvironment Estimation of Hepatocellular Carcinoma. Dis Markers (2022) 2022:8109771. doi: 10.1155/2022/8109771
26. Qin M, Ma Y, Wang Z, Fang D, Wei J. Using Immune-Related lncRNAs to Construct Novel Biomarkers and Investigate the Immune Landscape of Breast Cancer. Transl Cancer Res (2021) 10(6):2991–3003. doi: 10.21037/tcr-21-783
27. Zhao Z, Liu H, Zhou X, Fang D, Ou X, Ye J, et al. Necroptosis-Related lncRNAs: Predicting Prognosis and the Distinction Between the Cold and Hot Tumors in Gastric Cancer. J Oncol (2021) 2021:6718443. doi: 10.1155/2021/6718443
28. Jia X, Chen B, Li Z, Huang S, Chen S, Zhou R, et al. Identification of a Four-Gene-Based SERM Signature for Prognostic and Drug Sensitivity Prediction in Gastric Cancer. Front Oncol (2021) 11:799223:799223. doi: 10.3389/fonc.2021.799223
29. Zhang S, Zhang W, Zhang J. Comprehensive Analysis of Immune Cell Infiltration and Significant Genes in Head and Neck Squamous Cell Carcinoma. Oral Oncol (2022) 126:105755. doi: 10.1016/j.oraloncology.2022.105755
30. Zhang H, Wang X, Hou C, Yang Z. Identification of Driver Genes and Interaction Networks Related to Brain Metastasis in Breast Cancer Patients. Dis Markers (2022) 2022:7631456. doi: 10.1155/2022/7631456
31. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells With CIBERSORT. Methods Mol Biol (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12
32. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
33. Kong G, Wang Y, Huang Y, Shi Z. Identification and Verification of Tumor Immune Microenvironment-Related Prognostic Genes in Kidney Renal Clear Cell Carcinoma. BioMed Res Int (2022) 2022:5563668. doi: 10.1155/2022/5563668
34. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A Class Discovery Tool With Confidence Assessments and Item Tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
35. Quan Y, Zhang X, Ping H. Construction of a Risk Prediction Model Using M6a RNA Methylation Regulators in Prostate Cancer: Comprehensive Bioinformatic Analysis and Histological Validation. Cancer Cell Int (2022) 22(1):33. doi: 10.1186/s12935-021-02438-1
36. Zhang X, Wei X, Wang Y, Wang S, Ji C, Yao L, et al. Pyroptosis Regulators and Tumor Microenvironment Infiltration Characterization in Clear Cell Renal Cell Carcinoma. Front Oncol (2021) 11:774279:774279. doi: 10.3389/fonc.2021.774279
37. Dixon SJ, Lemberg KM, Lamprecht MR, Skouta R, Zaitsev EM, Gleason CE, et al. Ferroptosis: An Iron-Dependent Form of Nonapoptotic Cell Death. Cell (2012) 149(5):1060–72. doi: 10.1016/j.cell.2012.03.042
38. Sun X, Ou Z, Chen R, Niu X, Chen D, Kang R, et al. Activation of the P62-Keap1-NRF2 Pathway Protects Against Ferroptosis in Hepatocellular Carcinoma Cells. Hepatology (2016) 63(1):173–84. doi: 10.1002/hep.28251
39. Hassannia B, Vandenabeele P, Vanden Berghe T. Targeting Ferroptosis to Iron Out Cancer. Cancer Cell (2019) 35(6):830–49. doi: 10.1016/j.ccell.2019.04.002
40. Xu H, Ye D, Ren M, Zhang H, Bi F. Ferroptosis in the Tumor Microenvironment: Perspectives for Immunotherapy. Trends Mol Med (2021) 27(9):856–67. doi: 10.1016/j.molmed.2021.06.014
41. Shi ZZ, Fan ZW, Chen YX, Xie XF, Jiang W, Wang WJ, et al. Ferroptosis in Carcinoma: Regulatory Mechanisms and New Method for Cancer Therapy. Onco Targets Ther (2019) 12:11291–304. doi: 10.2147/OTT.S232852
42. Chen W, Zuo F, Zhang K, Xia T, Lei W, Zhang Z, et al. Exosomal MIF Derived From Nasopharyngeal Carcinoma Promotes Metastasis by Repressing Ferroptosis of Macrophages. Front Cell Dev Biol (2021) 9:791187:791187. doi: 10.3389/fcell.2021.791187
43. Fan F, Liu P, Bao R, Chen J, Zhou M, Mo Z, et al. A Dual PI3K/HDAC Inhibitor Induces Immunogenic Ferroptosis to Potentiate Cancer Immune Checkpoint Therapy. Cancer Res (2021) 81(24):6233–45. doi: 10.1158/0008-5472.CAN-21-1547
44. Luis G, Godfroid A, Nishiumi S, Cimino J, Blacher S, Maquoi E, et al. Tumor Resistance to Ferroptosis Driven by Stearoyl-CoA Desaturase-1 (SCD1) in Cancer Cells and Fatty Acid Biding Protein-4 (FABP4) in Tumor Microenvironment Promote Tumor Recurrence. Redox Biol (2021) 43:102006. doi: 10.1016/j.redox.2021.102006
45. Wang W, Green M, Choi JE, Gijon M, Kennedy PD, Johnson JK, et al. CD8(+) T Cells Regulate Tumour Ferroptosis During Cancer Immunotherapy. Nature (2019) 569(7755):270–4. doi: 10.1038/s41586-019-1170-y
46. Turley SJ, Cremasco V, Astarita JL. Immunological Hallmarks of Stromal Cells in the Tumour Microenvironment. Nat Rev Immunol (2015) 15(11):669–82. doi: 10.1038/nri3902
47. Huang Y, Chen S, Xiao L, Qin W, Li L, Wang Y, et al. A Novel Prognostic Signature for Survival Prediction and Immune Implication Based on SARS-CoV-2-Related Genes in Kidney Renal Clear Cell Carcinoma. Front Bioeng Biotechnol (2021) 9:744659:744659. doi: 10.3389/fbioe.2021.744659
48. Son HY, Jeong HK. Immune Evasion Mechanism and AXL. Front Oncol (2021) 11:756225:756225. doi: 10.3389/fonc.2021.756225
49. Zhang W, Kong X, Ai B, Wang Z, Wang X, Wang N, et al. Research Progresses in Immunological Checkpoint Inhibitors for Breast Cancer Immunotherapy. Front Oncol (2021) 11:582664:582664. doi: 10.3389/fonc.2021.582664
50. Pascual T, Cejalvo JM, Oliveira M, Vidal M, Vega E, Ganau S, et al. SOLTI-1503 PROMETEO TRIAL: Combination of Talimogene Laherparepvec With Atezolizumab in Early Breast Cancer. Future Oncol (2020) 16(24):1801–13. doi: 10.2217/fon-2020-0246
51. Adams S, Loi S, Toppmeyer D, Cescon DW, De Laurentiis M, Nanda R, et al. Pembrolizumab Monotherapy for Previously Untreated, PD-L1-Positive, Metastatic Triple-Negative Breast Cancer: Cohort B of the Phase II KEYNOTE-086 Study. Ann Oncol (2019) 30(3):405–11. doi: 10.1093/annonc/mdy518
52. Harris KM, Smilek DE, Byron M, Lim N, Barry WT, McNamara J, et al. Effect of Costimulatory Blockade With Abatacept After Ustekinumab Withdrawal in Patients With Moderate to Severe Plaque Psoriasis: The PAUSE Randomized Clinical Trial. JAMA Dermatol (2021) 157(11):1306–15. doi: 10.1001/jamadermatol.2021.3492
53. Schmid P, Salgado R, Park YH, Munoz-Couselo E, Kim SB, Sohn J, et al. Pembrolizumab Plus Chemotherapy as Neoadjuvant Treatment of High-Risk, Early-Stage Triple-Negative Breast Cancer: Results From the Phase 1b Open-Label, Multicohort KEYNOTE-173 Study. Ann Oncol (2020) 31(5):569–81. doi: 10.1016/j.annonc.2020.01.072
Keywords: breast cancer, ferroptosis, relapse-free survival, neoadjuvant chemotherapy, immune microenvironment
Citation: Xu Y, Du Y, Zheng Q, Zhou T, Ye B, Wu Y, Xu Q and Meng X (2022) Identification of Ferroptosis-Related Prognostic Signature and Subtypes Related to the Immune Microenvironment for Breast Cancer Patients Receiving Neoadjuvant Chemotherapy. Front. Immunol. 13:895110. doi: 10.3389/fimmu.2022.895110
Received: 13 March 2022; Accepted: 06 April 2022;
Published: 04 May 2022.
Edited by:
Qihui Shi, Fudan University, ChinaReviewed by:
Yu Dong, Shanghai Jiao Tong University, ChinaMeiyi Li, The Chinese University of Hong Kong, China
Copyright © 2022 Xu, Du, Zheng, Zhou, Ye, Wu, Xu and Meng. 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: Xuli Meng, bXhsbWFpbEAxMjYuY29t; Qiuran Xu, d2luZHdheTYyNkBzaW5hLmNvbQ==
†These authors have contributed equally to this work