- 1Department of Neurosurgery, West China Hospital of Sichuan University, Chengdu, China
- 2Department of Head and Neck Surgery, Sichuan Cancer Hospital and Institute, Sichuan Cancer Hospital, School of Medicine, University of Electronic Science and Technology of China, Chengdu, China
- 3Neuroscience & Metabolism Research, State Key Laboratory of Biotherapy, West China Hospital of Sichuan University, Chengdu, China
Glioma is one of the most malignant intracerebral tumors, whose treatment means was limited, and prognosis was unsatisfactory. Lactate metabolism patterns have been shown to be highly heterogenous among different tumors and produce diverse impact on the tumor microenvironment. To understand the characteristics and implications of lactate metabolism gene expression, we developed a lactate metabolism-related gene expression signature of gliomas based on RNA-sequencing data of a total of 965 patient samples from TCGA, CGGA, and our own glioma cohort. Sixty-three lactate metabolism-related genes (LMGs) were differentially expressed between glioma and normal brain tissue, and consensus clustering analysis identified two clusters distinct LMG expression patterns. The consensus clusters differed in prognosis, molecular characteristics and estimated immune microenvironment landscape involving immune checkpoint proteins, T cell dysfunction and exclusion, as well as tumor purity. Univariate Cox regression and Least Absolute Shrinkage and Selection Operator (LASSO) Cox hazard regression was applied in determining of prognosis-related lactate metabolism genes (PRLMGs), on which prognostic lactate metabolism risk score (PLMRS) was constructed. The high PLMRS group was associated with significantly poorer patient outcome. A nomogram containing PLMRS and other independent prognostic variables was established with remarkable predictive performance on patient survival. Exploration on the somatic mutations and copy number variations of the high- and low-PLMRS groups demonstrated their distinct genetic background. Together, our results indicated that the expression signature of LMG was associated with the prognosis of glioma patients and influenced the activity of immune cells in the tumor microenvironment, which may serve as a potential biomarker for predicting response of gliomas to immunotherapy.
Introduction
Glioma accounts for over 80% of all primary malignant intracerebral tumors in adult. Glioblastoma (GBM) accounts for 58.4% of newly diagnosed adult primary gliomas and is the most aggressive glioma subtype (1). Despite comprehensive anti-cancer treatment involving surgical resection, followed by radiotherapy and temozolomide (TMZ) chemotherapy, the 5-year-survival of GBM remained 6.8%-9.8%, while the relative survival for patients with astrocytoma and oligodendroglioma were 43.5% and 77.9%, respectively (1, 2). Warburg effect describes a metabolic characteristics of cancer cells, whose energy is mostly provided through high level of glycolysis, regardless of richness of oxygen supply, which is also named as “aerobic glycolysis” (3). Lactate, produced by active glycolysis, has received much attention in cancer research in the past decades, and lactate metabolism has been identified as a hallmark of cancer, which contributes to tumorigenesis (4). Clinical investigations have connected high level of lactate in the tumor microenvironment (TME) with poor prognosis in cervical cancer, breast cancer and head and neck squamous cell carcinomas (5–7). A magnetic resonance spectroscopy (MRS) study found that lactate concentration was significantly higher in isocitrate dehydrogenase (IDH) wild-type gliomas than the less aggressive IDH mutants (8).
Recent studies have revealed the important role of tumor lactate metabolism in the tumor immune microenvironment (9). For example, lactate accumulation in TME results in downregulation of CD4+ T cells, natural killer (NK) and natural killer T cells and upregulation of T-regular cells (Tregs) and myeloid-derived suppressor cells (MDSCs), leading to suppression of immune response (10–13). On the other hand, elevated serum lactate dehydrogenase (LDH) predicts shorter survival in melanoma, non-small cell lung cancer (NSCLC) and esophageal squamous cell carcinoma patients treated with immune checkpoint inhibitors (ICIs) (14–16). Additionally, extracellular acidic microenvironment induced by lactate is more toxic for normal cells compared with tumor cells (17, 18). Duan et al. reported lactate metabolic characteristics in glioblastoma and found hypoxia inducible factor 1α (HIF-1α), monocarboxylate transporters 1 (MCT1) and 4 (MCT4) indicating poor prognosis (19).
Lactate metabolism-related gene (LMG) expression signatures have been used in prognosis prediction and microenvironment of breast cancer, lung adenocarcinoma and hepatocellular carcinoma (20–22). However, such analysis hasn’t been applied in glioma. In the present study, we aimed to understand the relationship between lactate metabolism and glioma TME with bioinformatics methods and attempted to investigate its clinical implication by constructing a prognostic LMG expression signature.
Patients and methods
Data acquisition and definition of lactate metabolism-related genes
RNA-sequencing and clinical data of adult primary gliomas and normal control samples were collected from The Cancer Genome Atlas (TCGA)1 database (TCGA-GBM and TCGA-LGG), the Chinese Glioma Genome Atlas (CGGA)2 together with glioma patients from West China Hospital (WCH). Sample collection and sequencing process are available in Supplementary Materials Data Sheet 2. R package “TCGAbiolinks” was used to download data from TCGA (23, 24). Patients under 18-years-old at the time of diagnosis were excluded. The Molecular Signature Database (MSigDB) was used to search for genes involved in lactate metabolism with the key word “lactate” (25). To make the data comparable, genes unavailable or expressed at very low level (maximum fragment per kilobase million, FPKM under 0.1) in either of the dataset were excluded. R package “limma” was used to identify differentially expressed genes (DEG) comparing gliomas and normal brain tissues in TCGA database. Genes with adjusted P-values <0.05 were considered significant DEGs. Lactate metabolism-related genes (LMGs) referred to the overlap between genes included in the lactate metabolism pathways gene sets curated in the MSigDB and the differentially expressed genes between gliomas and normal brain samples in the TCGA dataset.
Unsupervised clustering analysis using LMGs
Consensus clustering analysis was applied to divide different lactate metabolism patterns based on expression of LMGs using R package “ConsensusClusterPlus” with 100 iterations (26). The optimal number of clusters was determined based on sample size of each cluster and the cumulative distribution function (CDF) curve of consensus index values in the consensus matrix (27). Briefly, we aimed to maximize sample size of each cluster while keeping a gradually increasing CDF. Principal component analysis (PCA) of all DEGs performed to visualize the difference in the transcriptome between clusters.
Clinical characteristics and biological function of consensus clusters
DEGs between the two clusters identified in consensus clustering with |log2FC| >0.5 and adjusted P-values <0.05 were selected for functional analysis. In the gene sets curated in MSigDB, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Hallmark gene sets (HALLMARK) were utilized to conduct Over-representation (ORA) and Gene set enrichment analysis (GSEA) using R package “clusterProfiler” (28). To analyze the differential expressed signaling pathways between the clusters, we first transformed the logFPKM matrix into pathway expression matrix using R package “GSVA” and identified differentially expressed pathways using “limma” (29).
Tumor immune microenvironment landscape
Cibersortx analysis3 was used to estimate the fraction of infiltrating immune cells of each tumor based on the RNA-seq data (30). The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) score was also utilized in analyzing the differences of tumor stromal and immune microenvironment (31).
For tumor purity estimates, we used the ESTIMATE tumor purity and consensus purity estimation (CPE) data published by Aran et al. (32). The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to evaluate T cell dysfunction score and T cell exclusion score (33).
Identification of prognosis-related lactate metabolism genes
To further investigate the role of LMGs in glioma prognosis and construct a prognostic model, prognosis-related lactate metabolism genes (PRLMGs) were identified. Data in TCGA database were randomly split into training set and validation set at a ratio of 6:4. Genes with adjusted P-value <0.05 in univariate Cox regression were selected in following analysis. Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression was performed in the training set to identify PRLMGs using R package “glmnet” (34). Genes whose coefficient was not 0 at lambda.min in 100 random repetition of LASSO Cox regression were defined as PRLMGs and selected for further analysis.
Prognostic model and clinical analysis based on PRLMGs
To better understand the prognostic role of PRLMGs, a lactate metabolism related prognostic signature was established with PRLMGs. The formula of prognostic lactate metabolism risk score (PLMRS) was:
Patients in training set were stratified into low-risk group and high-risk group accordingly. The optimal cut-off PLMRS was determined by function “surv_cutpoint” in R package “survminer” with a minimum group proportion of 0.3. Receiver operating characteristic (ROC) curves in validation sets of 1-year, 2-year and 3-year survival were illustrated and area under the ROC curve (AUC) were calculated by R package “timeROC”. Then, the prognostic model was tested in validation set.
Establishment and validation of nomogram
R package “rms” was used to construct and evaluate the nomogram. The nomogram was constructed by training a Cox-regression model using the cph function in the rms package with PLMRS and other potentially prognostic clinical factors to predict patient survival at 1, 2, and 3 years after diagnosis. Calibration curve was calculated by bootstrapping using calibrate function in the rms package to compare the actual survival rate and nomogram-predicted survival probability at these time points. Corrected concordance index (C-index) was calculated using validate function from the “rms” package.
Association between PLMRS and genetic alterations of gliomas
The somatic mutation and copy number variation (CNV) data of the glioma patients in the TCGA cohort was downloaded from cBioPortal4 (35). R package “maftools” were used to visualize most frequently mutated genes as well as PLMRGs (36). Gistic segment mean CNV values were used to generate the whole genome CNV profile. Amplification and homozygous deletion of genes were also reported in the cBioPortal data.
Statistical analysis
R software5 (Version 3.6.1) was used for data analysis and visualization. For continuous variables student t test was performed in two groups to compare differences, and one-way ANOVA was conducted to compare between three or more groups. Chi-square test was performed to detect the differences of categorical variables. Survival analysis was performed by R package “survminer”. The comparison of each Kaplan–Meier (KM) curve was accomplished using the log-rank test. Coefficients and their significance in the univariate and multivariate Cox regression analyses were computed using the coxph function. As for the statistical results, a two-sided P value <0.05 was considered a significant statistical difference. In the figures, a * indicated P < 0.05, while ** P < 0.01, *** P < 0.001, **** P < 0.0001.
Statement of ethics
Ethical review and approvement of this research were conducted by the Ethical Committee of Sichuan University. All principles in the Declaration of Helsinki (Ethic number: 2018.569) were strictly followed. Patients and their authorized trustees were informed and written informed consents were obtained before surgery.
Results
Identification of glioma LMGs
Flow chart of the current study was shown in Figure 1. RNA-seq data of 662 adult primary gliomas and five normal brain samples were downloaded from TCGA, 226 gliomas from CGGA, and generated from 77 glioma samples and 16 normal brain samples from our own patient cohort. The clinical information of these patients was recorded in Table 1.
Figure 1 Flow chart of study design. MSigDB, Molecular Signature Database; DEG, differentially expressed gene; LMG, lactate metabolism-related gene; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; TIDE, Tumor Immune Dysfunction and Exclusion; LASSO, Least Absolute Shrinkage and Selection Operator; PRLMG, prognosis-related lactate metabolism gene; PLMRS, prognostic lactate metabolism risk score.
Table 1 Clinicopathological characteristics of adult primary glioma patients in TCGA, CCGA, and WCH cohort.
A total of 223 genes related to lactate metabolism were found in MSigDB, of which 205 genes were included (Table S1). In the TCGA cohort, 7010 genes were differentially expressed in mRNA level between glioma samples and normal brain sample (Figure 2A; Table S2). The 63 genes overlapped between DEGs, and lactate metabolism gene sets were defined as LMGs (Figure 2B).
Figure 2 Consensus clustering of gliomas using LMG expression. (A) Volcano plot of gene fold change between adult primary glioma samples and normal brain in TCGA. (B) Overlap between tumor/normal DEGs and lactate metabolism related genes. (C) Consensus matrix of two LMGs clusters. (D, E) Concordance between two and three dimensional PCA of LMG expression and consensus clusters. LMG, lactate metabolism-related gene; TCGA, The Cancer Genome Atlas; DEG, differentially expressed gene; PCA, principal component analysis.
Consensus clustering of gliomas using LMG expression
Unsupervised consensus clustering was performed using RNA expression of LMGs in TGCA gliomas. Evaluation of the cluster sizes and CDFs found that the gliomas could be classified into two consensus clusters (Figure 2C; Figure S1). The two distinct LMG expression patterns of gliomas were also verified using two and three dimensional PCA (Figures 2D, E).
To investigate the stratification, we combined clinicopathologic and survival information with consensus clusters. We found that clustering of gliomas was significantly related to patients’ age at diagnosis, WHO grade, IDH mutational status, 1p19 codeletion status, Alpha Thalassemia Developmental Delay (ATRX) mutational status, O6-methylguanine-DNA methyltransferase (MGMT) promoter methylation status, and telomerase reverse transcriptase (TERT) promoter mutational status, but not gender (Figures 3A–H). K-M curve also revealed significant difference in patient outcome between the two clusters (Figure 3I).
Figure 3 Comparison of clinicopathological characteristics between two clusters. Differences between two clusters were compared in age at diagnosis (A), gender (B), tumor grade (C), IDH mutations (D), 1p19q codeletion (E), ATRX mutation (F), MGMT promoter methylation (G), TERT promoter mutation (H). (I) K-M survival curves of two clusters. TCGA, The Cancer Genome Atlas; WHO, World Health Organization; IDH, isocitrate dehydrogenase; ATRX, Alpha Thalassemia Developmental Delay; MGMT, O6-methylguanine-DNA methyltransferase; TERT, telomerase reverse transcriptase; K-M, Kaplan–Meier; WT, wild type. ****p < 0.0001. ns, not significant.
Functional analyses of DEGs between the two clusters found cell adhesion, immune response and epithelial mesenchymal transition (EMT) were over-represented in the KEGG and HALLMARK pathway annotation of these genes (Figures 4A, B). At the same time, in the GSEA, significant enrichment of the DEGs in cell cycle (normalized enrichment score = 2.564, adjusted-P value = 0.007) and EMT (normalized enrichment score = 4.413, adjusted-P value = 0.004) was also observed in the top 5 enriched gene sets of the KEGG and HALLMARK databases (Figures 4C, D). GSVA of two clusters showed enhanced expression of genes involved in the typical anabolic pathways of malignant tumors in the more aggressive cluster 1 (Figures 4E, F).
Figure 4 Biological function analysis of DEGs between the two consensus clusters. (A, B) Top five over-represented pathways in KEGG (A) and HALLMARK (B) annotation of the DEGs. (C, D) Top five pathways with the highest normalized enrichment score of GSEA in the KEGG (C) and the HALLMARK gene sets (D). (E, F) Ten pathways with smallest adjusted p-value of GSVA in KEGG (E) and HALLMARK pathways (F). TCGA, The Cancer Genome Atlas; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis; UP, upregulation; GSVA, gene set variation analysis.
Differential TME immune characteristics between consensus clusters
In silico analysis of immune cell infiltration using Cibersortx revealed a number of significant differences between the LMG expression clusters (Figure 5A). More specifically, in the more aggressive cluster 1 gliomas, the proportion of CD8+ T cells, which were the main target of ICIs, were higher, while that of Tregs, which have immunosuppressive properties, was also higher. In the meantime, higher fraction of resting NK cells was predicted to reside tumors in cluster 1. Additionally, higher stromal cell and immune cell contents were also found in the cluster 1 tumor by the ESTIMATE algorithm (Figure 5B). On the other hand, tumor purity, which was negatively correlated with the ESTIMATE scores was lower in cluster 1 (Figure 5C). In concordance with the TME immune cell content estimations, prediction of T cell exclusion and functional status using the TIDE algorithm revealed a higher T cell dysfunction and exclusion score in cluster 1, which suggests a potentially more complicated immune background in the cluster 1 gliomas and that these tumors were less likely to respond to immunotherapy in contrast to cluster 2 (Figure 5D). To shed light on the possible immune pathways that may interact with lactate metabolism, expression of 29 immune checkpoints (ICPs) mRNA was compared in TCGA database (Figure 5E) (37). The results indicated that most representative ICPs were significant higher in cluster 1, which further implied an immunosuppressive environment in cluster 1.
Figure 5 Immune microenvironment of two clusters. (A) Fraction of 22 types of immune infiltrating cells in Cibersortx in two clusters. (B) ESTIMATE score of two clusters. (C) Tumor purity estimation of two clusters using the ESTIMATE and the CPE algorithms (D) T cell dysfunction score and T cell exclusion score of two clusters. (E) mRNA expression of 29 ICPs in two clusters. TCGA, The Cancer Genome Atlas; NK, natural killer; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; CPE, consensus measurement of purity estimations. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.
Prognostic value of LMG expression signature
To understand if the LMG expression profile of gliomas could be utilized to evaluate prognosis of glioma patients, we attempted to construct an LMG-based prognostic measurement. We first split the TCGA glioma samples into training and testing datasets at a ratio of 6:4. In the training set, 63 LMGs were screened for their significant prognostic value using univariate Cox regression analysis. Then, 47 LMGs with adjusted P-value < 0.05 in univariate Cox regression were used to build multivariate LASSO Cox regression models for 100 random repetitions to select LMGs that showed unconfounded robust correlation with patient survival. In the end, 14 genes had non-zero coefficients in 100 LASSO Cox regression models and were qualified as prognostic LMGs (PRLMGs, Figures 6A, B).
Figure 6 LASSO Cox regression and identification of PRLMGs. (A) Average of means and 95% confidence intervals of the C-index of LASSO Cox regression model at each lambda value. (B) Average of coefficients of LMGs in the LASSO Cox regression model at each lambda value. (C) mRNA expression of each PRLMG in glioma and normal brain tissue. (D) Univariate Cox regression analysis of PRLMGs in glioma prognosis. (E) Alluvial plot showing relationships between LMGs clusters and PLMRS risk groups. LASSO, Least Absolute Shrinkage and Selection Operator; C-index, concordance index; LMG, lactate metabolism-related gene; PRLMG, prognosis-related lactate metabolism gene. *P < 0.05; **P < 0.01; ***P < 0.001.
Expression of PRLMGs was listed in Figure 6C. Among these genes, 13 of them had higher expression in tumor tissues, while MPL was expressed at a lower level in gliomas than the normal brain tissues (Figure 6C). In univariate cox regression analysis, 12 of PRLMGs were associated with worse prognosis, while the rest two were associated with prolonged OS (Figure 6D). All the PRLMGs were used to construct a final Prognostic LMG expression Risk Score (PLMRS) using multivariate Cox regression in the TCGA training set. The final calculation formula was:
Using optimal PLMRS cutoff determined by maximizing outcome difference in each cohort (1.76 for the TCGA cohort, 8.79 for the CGGA cohort, 0.042 for our own cohort), the gliomas was divided into low- risk and high-risk groups. An alluvial diagram demonstrated the gliomas in cluster 1 were more likely to be at the high-risk group (Figure 6E). ROC analysis found that the PLMRS showed high predictive performance for 1-, 2-, 3-year survival in not only the TCGA-validation group, but also the CGGA and our own patient cohort (Figures 7A–C). In each cohort, patients in the low-risk group had significantly better outcome compared to those in the high-risk group (Figures 7D–G).
Figure 7 Prognostic value of PLMRS. ROC curves and matched AUC of 1-, 2-, 3-year survival in TCGA validation group (A), CGGA (B) and WCH (C). K-M curves of high-risk and low-risk groups in TCGA training group (D), TCGA validation group (E), CGGA (F) and WCH (G). TCGA, The Cancer Genome Atlas; CGGA, Chinese Glioma Genome Atlas; WCH, West China Hospital; PLMRS, prognostic lactate metabolism risk score; ROC, receiver operating characteristic curve; AUC, area under the curve.
Development and validation of a PLMRS-based nomogram
To demonstrate the translational value of PLMRS in clinical settings, we attempted to construct a nomogram that involved the PLMRS and potential prognostic clinical information. Univariate Cox regression analysis of the TCGA cohort found that apart from the PLMRS, the age and Karnofsky Performance Score (KPS) of patients at diagnosis, the WHO grade, codeletion of chromosome arms 1p and 19q, mutational status of IDH genes, as well as chemotherapy and radiotherapy were significantly associated with patient survival (Figure 8A). In further multivariate analysis (Figure 8B) tumor grade, PLMRS, IDH mutational status and radiotherapy remained significant prognostic factors. As a result, we included these factors in the final nomogram construction (Figure 8C; Figure S2). The calibration plot of the nomogram for each cohort showed that the prediction of the nomogram was highly consistent with observed patient outcome in all three cohorts (Figures 8D–F). Compared with PLMRS alone, nomogram combining PLMRS, and other independent prognostic predictors showed superiority in their predictive performance of patient prognosis (C-index 0.837 for PLMRS only vs. 0.851 for combined nomogram, TCGA; 0.771 vs. 0.787, CGGA; 0.627 vs. 0.699, WCH).
Figure 8 Nomogram construction and validation. Univariate (A) and multivariate (B) Cox regression to identify independent prognostic factors. (C) Nomogram of 1-, 2-, 3-year survival of glioma established in TCGA. Calibration plots of nomogram developed by TCGA (D), CGGA (E) and WCH (F). TCGA, The Cancer Genome Atlas; CGGA, Chinese Glioma Genome Atlas; WCH, West China Hospital; WHO, World Health Organization; PLMRS, prognostic lactate metabolism risk score; RT, radiotherapy; CT, chemotherapy; KPS, Karnofsky Performance Score; IDH, isocitrate dehydrogenase; WT, wild type.
Associations between PLMRS and characteristics of gliomas
To understand how gliomas from two risk groups differed from each other, we analyzed the relationships between PLMRS and the clinicopathological features, as well as genetic background of gliomas. PLMRS was significantly associated with age at diagnosis, tumor grade, IDH mutational status, 1p19q codeletion status, ATRX mutational status, MGMT promoter methylation status, TERT promoter status and histology, but not gender (Figures 9A–H). Inspection at the expression of PRLMGs and the clinicopathological characteristics revealed apparent difference in these genes between different glioma subtypes (Figure 9I).
Figure 9 Comparison of clinicopathological characteristics in two PLMRS risk groups. Differences between the high and low risk groups were compared in age at diagnosis (A), gender (B), tumor grade (C), IDH mutations (D), 1p19q codeletion (E), ATRX mutation (F), MGMT promoter methylation (G), TERT promoter mutation (H). (I) Heatmap summarizing the clinicopathological feature and expressing PRLMGs in two risk groups. PLMRS, prognostic lactate metabolism risk score; TCGA, The Cancer Genome Atlas; WHO, World Health Organization; IDH, isocitrate dehydrogenase; ATRX, Alpha Thalassemia Developmental Delay; MGMT, O6-methylguanine-DNA methyltransferase; TERT, telomerase reverse transcriptase; WT, wild type. ****p < 0.0001. ns, not significant.
Next, we analyzed the differences in the genetic alterations between the two risk groups by comparing the top 20 somatic mutations of low-risk group and high-risk group (Figures 10A, B). While a high rate of TP53 mutations were present in both groups, IDH1 mutation occurred in 90% of gliomas in low-risk group, but less than 20% in high-risk group. EGFR was the second most common mutations in high-risk group. At larger genomic landscape, (Figure 10C), high-risk gliomas exhibited evident gain of chromosome 7 and coexisting loss of chromosome 10.
Figure 10 Somatic mutations and CNV of two risk groups. Top 20 somatic mutations in low-risk group (A) and high-risk group (B). (C) Landscape of autosomal CNV in two risk groups. CNV, copy number variation; TCGA, The Cancer Genome Atlas; PLMRS, prognostic lactate metabolism risk score; IDH, isocitrate dehydrogenase; TERT, telomerase reverse transcriptase; ATRX, Alpha Thalassemia Developmental Delay; MGMT, O6-methylguanine-DNA methyltransferase; GISTIC, Genomic Identification of Significant Targets in Cancer.
Discussion
Glioma is the most common primary malignant intracerebral tumor in adult (1). Despite the massive research investment and clinical experiments reported worldwide, the treatment of glioma has not yet to be revolutionized and the prognosis of glioma patients is not ideal, causing loss of quality of life (QOL) and financial stress for patients and their family (38, 39). Lactate is one of the components in the TME, whose concentration has been associated with the molecular features and prognosis of gliomas (8). Meanwhile, lactate metabolism has been considered to regulate tumor immunity through diverse mechanisms. With the assistance of public resources and our in-house dataset, we interrogated the relationships between the expression of lactate metabolism-related genes and clinicopathological characteristics of gliomas, as well as the potential impact of LMG expression on the landscape of glioma immune TME.
For tumor cells, aerobic glycolysis was considered more of a survival strategy than a compromise to insufficient oxygen supply. Nicotinamide adenine dinucleotide (NADH), the intermediate product of glycolysis, was used in synthesis of other organic molecules, while lactate was released, thus creating an acidic microenvironment (40). In gliomas, elevated lactate was associated with higher blood volume, and lower pH was connected with higher mitotic index (41, 42). Some researchers proposed a dynamic equilibrium of lactate in TME, where a group of tumor cells utilize lactate as energy source (43). In this study, we discovered that the unsupervised clustering with LMG expression signature could stratify glioma samples into two consensus clusters, in which one was associated with IDH wildtype gliomas and was more likely to have a more sophisticated and immunosuppressive microenvironment. This result was consistent with the findings reported by Wenger et al. that the higher level of glycolysis in the IDH wild type gliomas was associated with more lactate in TME, which could potentially disturb the distribution and functionality of immune cells (8).
Tumor immunotherapy has gained growing attention in the recent years. Current immunotherapies include cytokine therapy, ICIs, chimeric antigen receptor (CAR) T-cell therapy, oncolytic virus therapy, vaccine therapy and therapies targeting other immune related molecules or cells. However, none of the treatment showed promising clinical effect on gliomas so far (44). Among the numerous factors interfering anti-glioma immunity, lactate in the TME could play an important role in suppressing anti-tumor immune, by inducing Tregs and TAM M2, inhibiting NK and dendritic cells, activating MDSCs and decreasing cytotoxic T lymphocytes (CTLs) (45). Here, we investigated the association between lactate metabolism and the disturbed immune microenvironment of gliomas. In our result, both indolamine 2,3-dioxygenase 1(IDO1) and programmed cell death 1 ligand 1 (PD-L1 or CD 274) mRNA expression were higher in cluster 1. Overexpression of IDO1 resulted in recruitment of Tregs (46). This finding was consistent with the higher level of Tregs in cluster 1. PD-L1 expression was significant higher in TAM of GBM (47). Shan et al. reported that the expression of PD-L1 protein in TAM M2, which tended to be pro-tumor, was upregulated when treated with lactate (48). As expected, a higher fraction of TAM M2 score was observed in cluster 1. Programmed Cell Death 1 (PDCD-1 or PD-1) and Cytotoxic T-Lymphocyte Associated Protein 4 (CTLA4) were well-characterized ICPs. Inhibition of PD-1 leading to rescue of T lymphocyte and NK cell function, and CTLA-4 block resulted in activation of T cells (44). In our findings, both PD-1 and CTLA4 mRNA expression were elevated in cluster 1. In the meantime, NK cells were less activated in cluster 1. Although CD8+ T cells were higher in cluster 1, they could be inactivated due to the co-stimulation of PD-1 and CTLA4. In summary, our results suggest that the anti-tumor immunity is closely connected with the expression profile of LMGs.
To understand the clinical implications of LMG expression, we identified the LMGs exhibiting robust correlation with the clinical outcome of glioma patients, which we defined as PRLMGs. Among the PRLMGs, solute carrier family 7 member 7 (SLC7A7) and glycogen phosphorylase, liver form (PYGL) has been previously reported to be predictor of poor survival of GBM (49, 50). Solute carrier family 16 member 1 (SLC16A1), also known as monocarboxylate transporter 1 (MCT1), is responsible for the transmembrane delivery of lactate. Overexpression of SLC16A1 also predicted poor survival of high-grade gliomas (51, 52). The therapeutic and prognostic role of SLC16A1 has also been pointed out (53, 54). Tet methylcytosine dioxygenase 2 (TET2) catalyzed conversion of 5-methylcytosine (5mC) to 5-hydroxymethylcytocine (5hmC), which is important in homeostasis of hematopoietic cells (55). Of note, IDH mutation reversibly sabotaged the demethylation function of TET2 via oncometabolite 2-hydroxyglutarate (2HG) in acute myeloid leukemia (AML) (56). Similar mechanisms were reported in cholangiocarcinoma, mutation of IDH caused repression of TET2, which could induce immune suppression (57). In GBM, loss of TET2 was associated with poor survival, which could be mediated by sex determining region Y-box transcription factor 2 (Sox2) (58). Together, our results and previous findings indicate that, among the PRLMGs, SLC16A1 and TET2 are potentially important targets for intervention of glioma lactate metabolism that deserve priority in future investigation.
Despite the correlations between the expression of PRLMGs and aggressiveness of gliomas, a number of prognostic factors were also known to produce significant impact on the outcome of glioma patients. Besides patient comorbidity status and histological tumor grading, mutation of IDH genes and codeletion of chromosome arm 1p/19q have been well-recognized indicators for patient prognosis (59). Compared to the genetic background of gliomas, adjuvant treatments appeared to produce less impact on patient survival, but our results indicate that radiotherapy remained significant in prolonging patient survival. To translate this panel of clinical, pathological, and molecular information into succinct measurement, we combined the PLMRS with other independent prognostic factors to construct a nomogram that could accurately predict patient survival. Of note, in both the TCGA and CGGA cohort, PLMRS contributed most to the output of the nomogram, which emphasized the potential clinical value of the LMG expression signature in predicting outcome of gliomas.
In the present study, we comprehensively explore the connections of the expression of lactate metabolism-related genes with the clinical, pathological and immunological aspects of gliomas. However, there were still some limitations in this study. First, in spite of the usage of up to three independent data sets, the sequencing protocols and data preprocessing steps were different, therefore the PLMRS cutoffs and nomograms were established separately in each cohort. Nevertheless, the major trends identified in the study should not be influenced by the lack of consistency in these parameters. Secondly, in this study, the interactions between expression of LMGs and tumor immunity were derived from in silico analysis. The precise effects of different LMG expression profiles on the glioma immune TME and the mechanisms mediating these effects should be elaborated in the future experimental studies.
Conclusion
In conclusion, we found that the expression signature of LMGs could stratify gliomas into two consensus clusters which exhibited differential clinical, pathological, molecular characteristics, as well as the prognosis of the patients. Gliomas in two LMG expression clusters displayed distinct tumor immune microenvironment landscapes that could potentially imply their response to immunotherapies.
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. The sequencing data generated in this study are available at the Genome Sequence Archive for Humans (accession code HRA002839). The website of the Genome Sequence Archive for Humans is https://ttps://ngdc.cncb.ac.cn/gsa-human.
Ethics statement
The studies involving human participants were reviewed and approved by Ethical Committee of Sichuan University. The patients/participants provided their written informed consent to participate in this study.
Author contributions
Study design: ZW, SZ, YL, and MC; data retrieve: JL, SC, YY, WL, and ZW; statistical analysis: SZ and MZ; result interpretation: ZW, SZ, and WF; writing-original draft: all authors; writing-revise: YL, MC, and SZ. All authors contributed to the article and approved the submitted version.
Funding
This study is supported by Clinical Research Innovation Project, West China Hospital, Sichuan University (19HXCX009 to YL), Post-Doctor Research Project, West China Hospital, Sichuan University (20HXBH035 to SZ), Science and technology project, technology innovation research and development project, Chengdu (2022-YF05-01456-SN to MC).
Acknowledgments
Our sincere thanks to the developers and maintainers of TCGA, CGGA, the R software contributors, and all the related algorithms and websites.
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/fonc.2022.1008219/full#supplementary-material
Footnotes
- ^ https://portal.gdc.cancer.gov/
- ^ http://www.cgga.org.cn/index.jsp
- ^ https://cibersortx.stanford.edu/
- ^ http://cbioportal.org
- ^ http://www.R-project.org
References
1. Ostrom QT, Cioffi G, Waite K, Kruchko C, Barnholtz-Sloan JS. CBTRUS statistical report: Primary brain and other central nervous system tumors diagnosed in the united states in 2014-2018. Neuro Oncol (2021) 23(12 Suppl 2):iii1–iii105. doi: 10.1093/neuonc/noab200
2. Stupp R, Hegi ME, Mason WP, van den Bent MJ, Taphoorn MJ, Janzer RC, et al. Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial. Lancet Oncol (2009) 10(5):459–66. doi: 10.1016/S1470-2045(09)70025-7
3. Warburg O. On the origin of cancer cells. Science (1956) 123(3191):309–14. doi: 10.1126/science.123.3191.309
4. Hirschhaeuser F, Sattler UG, Mueller-Klieser W. Lactate: a metabolic key player in cancer. Cancer Res (2011) 71(22):6921–5. doi: 10.1158/0008-5472.CAN-11-1457
5. Walenta S, Wetterling M, Lehrke M, Schwickert G, Sundfør K, Rofstad EK, et al. High lactate levels predict likelihood of metastases, tumor recurrence, and restricted patient survival in human cervical cancers. Cancer Res (2000) 60(4):916–21.
6. Cheung SM, Husain E, Masannat Y, Miller ID, Wahle K, Heys SD, et al. Lactate concentration in breast cancer using advanced magnetic resonance spectroscopy. Br J Cancer (2020) 123(2):261–7. doi: 10.1038/s41416-020-0886-7
7. Ziebart T, Walenta S, Kunkel M, Reichert TE, Wagner W, Mueller-Klieser W. Metabolic and proteomic differentials in head and neck squamous cell carcinomas and normal gingival tissue. J Cancer Res Clin Oncol (2011) 137(2):193–9. doi: 10.1007/s00432-010-0875-y
8. Wenger KJ, Steinbach JP, Bähr O, Pilatus U, Hattingen E. Lower lactate levels and lower intracellular pH in patients with IDH-mutant versus wild-type gliomas. AJNR Am J Neuroradiol (2020) 41(8):1414–22. doi: 10.3174/ajnr.A6633
9. Certo M, Tsai CH, Pucino V, Ho PC, Mauro C. Lactate modulation of immune responses in inflammatory versus tumour microenvironments. Nat Rev Immunol (2021) 21(3):151–61. doi: 10.1038/s41577-020-0406-2
10. Husain Z, Huang Y, Seth P, Sukhatme VP. Tumor-derived lactate modifies antitumor immune response: effect on myeloid-derived suppressor cells and NK cells. J Immunol (2013) 191(3):1486–95. doi: 10.4049/jimmunol.1202702
11. Comito G, Iscaro A, Bacci M, Morandi A, Ippolito L, Parri M, et al. Lactate modulates CD4(+) T-cell polarization and induces an immunosuppressive environment, which sustains prostate carcinoma progression via TLR8/miR21 axis. Oncogene (2019) 38(19):3681–95. doi: 10.1038/s41388-019-0688-7
12. Xie D, Zhu S, Bai L. Lactic acid in tumor microenvironments causes dysfunction of NKT cells by interfering with mTOR signaling. Sci China Life Sci (2016) 59(12):1290–6. doi: 10.1007/s11427-016-0348-7
13. Angelin A, Gil-de-Gómez L, Dahiya S, Jiao J, Guo L, Levine MH, et al. Foxp3 reprograms T cell metabolism to function in low-glucose, high-lactate environments. Cell Metab (2017) 25(6):1282–93.e7. doi: 10.1016/j.cmet.2016.12.018
14. Schouwenburg MG, Suijkerbuijk KPM, Koornstra RHT, Jochems A, van Zeijl MCT, van den Eertwegh AJM, et al. Switching to immune checkpoint inhibitors upon response to targeted therapy; the road to long-term survival in advanced melanoma patients with highly elevated serum LDH? Cancers (Basel) (2019) 11(12):1940. doi: 10.3390/cancers11121940
15. Zhang Z, Li Y, Yan X, Song Q, Wang G, Hu Y, et al. Pretreatment lactate dehydrogenase may predict outcome of advanced non small-cell lung cancer patients treated with immune checkpoint inhibitors: A meta-analysis. Cancer Med (2019) 8(4):1467–73. doi: 10.1002/cam4.2024
16. Wang X, Zhang B, Chen X, Mo H, Wu D, Lan B, et al. Lactate dehydrogenase and baseline markers associated with clinical outcomes of advanced esophageal squamous cell carcinoma patients treated with camrelizumab (SHR-1210), a novel anti-PD-1 antibody. Thorac Cancer (2019) 10(6):1395–401. doi: 10.1111/1759-7714.13083
17. Zhang X, Lin Y, Gillies RJ. Tumor pH and its measurement. J Nucl Med (2010) 51(8):1167–70. doi: 10.2967/jnumed.109.068981
18. Webb BA, Chimenti M, Jacobson MP, Barber DL. Dysregulated pH: a perfect storm for cancer progression. Nat Rev Cancer (2011) 11(9):671–7. doi: 10.1038/nrc3110
19. Duan K, Liu ZJ, Hu SQ, Huo HY, Xu ZR, Ruan JF, et al. Lactic acid induces lactate transport and glycolysis/OXPHOS interconversion in glioblastoma. Biochem Biophys Res Commun (2018) 503(2):888–94. doi: 10.1016/j.bbrc.2018.06.092
20. Yang L, Tan P, Sun H, Zeng Z, Pan Y. Integrative dissection of novel lactate metabolism-related signature in the tumor immune microenvironment and prognostic prediction in breast cancer. Front Oncol (2022) 12:874731. doi: 10.3389/fonc.2022.874731
21. Wang Z, Embaye KS, Yang Q, Qin L, Zhang C, Liu L, et al. Establishment and validation of a prognostic signature for lung adenocarcinoma based on metabolism-related genes. Cancer Cell Int (2021) 21(1):219. doi: 10.1186/s12935-021-01915-x
22. Li Y, Mo H, Wu S, Liu X, Tu K. A novel lactate metabolism-related gene signature for predicting clinical outcome and tumor microenvironment in hepatocellular carcinoma. Front Cell Dev Biol (2021) 9:801959. doi: 10.3389/fcell.2021.801959
23. Zhao Z, Zhang KN, Wang Q, Li G, Zeng F, Zhang Y, et al. Chinese Glioma genome atlas (CGGA): A comprehensive resource with functional genomic data from Chinese glioma patients. Genomics Proteomics Bioinf (2021) 19(1):1–12. doi: 10.1016/j.gpb.2020.10.005
24. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res (2016) 44(8):e71. doi: 10.1093/nar/gkv1507
25. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics (2011) 27(12):1739–40. doi: 10.1093/bioinformatics/btr260
26. 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
27. Monti S, Tamayo P, Mesirov JP, Golub TR. Consensus clustering: A resampling-based method for class discovery and visualization of gene expression microarray data. Mach Learn (2003) 52(1-2):91–118. doi: 10.1023/A:1023949509487
28. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. Omics (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118
29. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7
30. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337
31. Yoshihara K, Shahmoradgoli M, Martínez 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
32. Aran D, Sirota M, Butte AJ. Systematic pan-cancer analysis of tumour purity. Nat Commun (2015) 6:8971. doi: 10.1038/ncomms9971
33. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1
34. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software (2010) 33(1):1–22. doi: 10.18637/jss.v033.i01
35. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discovery (2012) 2(5):401–4. doi: 10.1158/2159-8290.CD-12-0095
36. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118
37. Song W, Ren J, Xiang R, Kong C, Fu T. Identification of pyroptosis-related subtypes, the development of a prognosis model, and characterization of tumor microenvironment infiltration in colorectal cancer. Oncoimmunology (2021) 10(1):1987636. doi: 10.1080/2162402X.2021.1987636
38. Colen CB, Allcut E. Quality of life and outcomes in glioblastoma management. Neurosurg Clin N Am (2012) 23(3):507–13. doi: 10.1016/j.nec.2012.04.010
39. Ståhl P, Henoch I, Smits A, Rydenhag B, Ozanne A. Quality of life in patients with glioblastoma and their relatives. Acta Neurol Scand (2022) 146(1):82–91. doi: 10.1111/ane.13625
40. Vander Heiden MG, Cantley LC, Thompson CB. Understanding the warburg effect: the metabolic requirements of cell proliferation. Science (2009) 324(5930):1029–33. doi: 10.1126/science.1160809
41. Li X, Vigneron DB, Cha S, Graves EE, Crawford F, Chang SM, et al. Relationship of MR-derived lactate, mobile lipids, and relative blood volume for gliomas in vivo. AJNR Am J Neuroradiol (2005) 26(4):760–9.
42. Ferrauto G, Di Gregorio E, Auboiroux V, Petit M, Berger F, Aime S, et al. CEST-MRI for glioma pH quantification in mouse model: Validation by immunohistochemistry. NMR Biomed (2018) 31(11):e4005. doi: 10.1002/nbm.4005
43. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell (2011) 144(5):646–74. doi: 10.1016/j.cell.2011.02.013
44. Sener U, Ruff MW, Campian JL. Immunotherapy in glioblastoma: Current approaches and future perspectives. Int J Mol Sci (2022) 23(13):7046. doi: 10.3390/ijms23137046
45. Wang JX, Choi SYC, Niu X, Kang N, Xue H, Killam J, et al. Lactic acid and an acidic tumor microenvironment suppress anticancer immunity. Int J Mol Sci (2020) 21(21):8363. doi: 10.3390/ijms21218363
46. Wainwright DA, Balyasnikova IV, Chang AL, Ahmed AU, Moon KS, Auffinger B, et al. IDO expression in brain tumors increases the recruitment of regulatory T cells and negatively impacts survival. Clin Cancer Res (2012) 18(22):6110–21. doi: 10.1158/1078-0432.CCR-12-2130
47. Bloch O, Crane CA, Kaur R, Safaee M, Rutkowski MJ, Parsa AT. Gliomas promote immunosuppression through induction of B7-H1 expression in tumor-associated macrophages. Clin Cancer Res (2013) 19(12):3165–75. doi: 10.1158/1078-0432.CCR-12-3314
48. Shan T, Chen S, Chen X, Wu T, Yang Y, Li S, et al. M2−TAM subsets altered by lactic acid promote t−cell apoptosis through the PD−L1/PD−1 pathway. Oncol Rep (2020) 44(5):1885–94. doi: 10.3892/or.2020.7767
49. Zhao CY, Hua CH, Li CH, Zheng RZ, Li XY. High PYGL expression predicts poor prognosis in human gliomas. Front Neurol (2021) 12:652931. doi: 10.3389/fneur.2021.652931
50. Fan S, Meng D, Xu T, Chen Y, Wang J, Li X, et al. Overexpression of SLC7A7 predicts poor progression-free and overall survival in patients with glioblastoma. Med Oncol (2013) 30(1):384. doi: 10.1007/s12032-012-0384-8
51. Lin HH, Tsai WC, Tsai CK, Chen SH, Huang LC, Hueng DY, et al. Overexpression of cell-surface marker SLC16A1 shortened survival in human high-grade gliomas. J Mol Neurosci (2021) 71(8):1614–21. doi: 10.1007/s12031-021-01806-w
52. Miranda-Gonçalves V, Gonçalves CS, Granja S, Vieira de Castro J, Reis RM, Costa BM, et al. MCT1 is a new prognostic biomarker and its therapeutic inhibition boosts response to temozolomide in human glioblastoma. Cancers (Basel) (2021) 13(14):3468. doi: 10.3390/cancers13143468
53. Miranda-Gonçalves V, Honavar M, Pinheiro C, Martinho O, Pires MM, Pinheiro C, et al. Monocarboxylate transporters (MCTs) in gliomas: expression and exploitation as therapeutic targets. Neuro Oncol (2013) 15(2):172–88. doi: 10.1093/neuonc/nos298
54. Thakur A, Qiu G, Xu C, Han X, Yang T, Ng SP, et al. Label-free sensing of exosomal MCT1 and CD147 for tracking metabolic reprogramming and malignant progression in glioma. Sci Adv (2020) 6(26):eaaz6119. doi: 10.1126/sciadv.aaz6119
55. Li Z, Cai X, Cai CL, Wang J, Zhang W, Petersen BE, et al. Deletion of Tet2 in mice leads to dysregulated hematopoietic stem cells and subsequent development of myeloid malignancies. Blood (2011) 118(17):4509–18. doi: 10.1182/blood-2010-12-325241
56. Kernytsky A, Wang F, Hansen E, Schalm S, Straley K, Gliser C, et al. IDH2 mutation-induced histone and DNA hypermethylation is progressively reversed by small-molecule inhibition. Blood (2015) 125(2):296–303. doi: 10.1182/blood-2013-10-533604
57. Wu MJ, Shi L, Dubrot J, Merritt J, Vijay V, Wei TY, et al. Mutant IDH inhibits IFNγ-TET2 signaling to promote immunoevasion and tumor maintenance in cholangiocarcinoma. Cancer Discovery (2022) 12(3):812–35. doi: 10.1158/2159-8290.CD-21-1077
58. Lopez-Bertoni H, Johnson A, Rui Y, Lal B, Sall S, Malloy M, et al. Sox2 induces glioblastoma cell stemness and tumor propagation by repressing TET2 and deregulating 5hmC and 5mC DNA modifications. Signal Transduct Target Ther (2022) 7(1):37. doi: 10.1038/s41392-021-00857-0
Keywords: glioma, lactate metabolism, tumor microenvironment, immunity, prognosis
Citation: Wang Z, Zhang S, Li J, Yuan Y, Chen S, Zuo M, Li W, Feng W, Chen M and Liu Y (2022) Prognostic value of lactate metabolism-related gene expression signature in adult primary gliomas and its impact on the tumor immune microenvironment. Front. Oncol. 12:1008219. doi: 10.3389/fonc.2022.1008219
Received: 31 July 2022; Accepted: 26 August 2022;
Published: 20 September 2022.
Edited by:
Min Zhao, University of the Sunshine Coast, AustraliaReviewed by:
Ding Hu, Qingdao University, ChinaYufang Hou, Chinese Academy of Medical Sciences and Peking Union Medical College, China
Copyright © 2022 Wang, Zhang, Li, Yuan, Chen, Zuo, Li, Feng, Chen and Liu. 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: Yanhui Liu, bGl1eWhAc2N1LmVkdS5jbg==; Mina Chen, Y2hlbm1pbmEyMDEwQHNjdS5lZHUuY24=
†These authors have contributed equally to this work