Article Text
PDF
Abstract
Background Pancreatic cancer (PAC) is one of the most malignant cancer types and immunotherapy has emerged as a promising treatment option. PAC cells undergo metabolic reprogramming, which is thought to modulate the tumor microenvironment (TME) and affect immunotherapy outcomes. However, the metabolic landscape of PAC and its association with the TME remains largely unexplored.
Methods We characterized the metabolic landscape of PAC based on 112 metabolic pathways and constructed a novel metabolism-related signature (MBS) using data from 1,188 patients with PAC. We evaluated the predictive performance of MBS for immunotherapy outcomes in 11 immunotherapy cohorts from both bulk-RNA and single-cell perspectives. We validated our results using immunohistochemistry, western blotting, colony-formation assays, and an in-house cohort.
Results MBS was found to be negatively associated with antitumor immunity, while positively correlated with cancer stemness, intratumoral heterogeneity, and immune resistant pathways. Notably, MBS outperformed other acknowledged signatures for predicting immunotherapy response in multiple immunotherapy cohorts. Additionally, MBS was a powerful and robust biomarker for predicting prognosis compared with 66 published signatures. Further, we identified dasatinib and epothilone B as potential therapeutic options for MBS-high patients, which were validated through experiments.
Conclusions Our study provides insights into the mechanisms of immunotherapy resistance in PAC and introduces MBS as a robust metabolism-based indicator for predicting response to immunotherapy and prognosis in patients with PAC. These findings have significant implications for the development of personalized treatment strategies in patients with PAC and highlight the importance of considering metabolic pathways and immune infiltration in TME regulation.
- metabolic landscape
- pancreatic cancer
- immunotherapy
- immune infiltration
- prognosis
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ON THIS TOPIC
The existing biomarkers presently do not provide sufficient accuracy in predicting the prognosis and effectiveness of immunotherapy in pancreatic cancer (PAC). Consequently, there is an urgent necessity to discover reliable and precise clinical biomarkers that can predict the prognosis and responsiveness to immunotherapy in patients with PAC through minimally invasive methods.
WHAT THIS STUDY ADDS
Our study has created a machine learning-based platform called the metabolism-related signature (MBS), which enables the prediction of prognosis and immunotherapy outcomes for patients with PAC. Additionally, we have identified alternative therapeutic agents for patients with PAC who are deemed unsuitable for immunotherapy.
HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY
The study conducted a systematic investigation on the significance of metabolic patterns in the anticancer immune response of PAC. The MBS demonstrated exceptional performance in predicting the response to immunotherapy, surpassing established signatures. This signifies the potential of MBS to serve as a valuable biomarker for guiding personalized immunotherapy strategies for patients with PAC. By leveraging MBS, clinicians can make more informed decisions regarding immunotherapy treatment and improve patient outcomes.
Key points
We characterized the metabolic landscape of pancreatic cancer based on 112 metabolic pathways and constructed a novel metabolism-related signature (MBS) using machine learning methods.
MBS was a powerful and robust biomarker for predicting prognosis of pancreatic cancer compared with 66 published signatures.
MBS outperformed other acknowledged signatures for predicting immunotherapy response in multiple immunotherapy cohorts and in-house cohorts.
We identified dasatinib and epothilone B as potential therapeutic options for MBS-high patients, which were validated through experiments.
Background
Pancreatic cancer (PAC) is one of the most aggressive and lethal cancer types worldwide. Although most patients have already developed metastasis at initial diagnosis, surgical resection is the primary treatment for PAC.1 While chemotherapy and radiotherapy represent alternative treatment options, little progress has been made in improving patient outcome in recent decades. Immunotherapy has ushered in a novel approach, with immune checkpoint inhibition relieving tumor-infiltrating lymphocytes (TILs) suppression, thus leading to the enhanced activation of TILs and subsequent tumor cell clearance.2 3 Nevertheless, there are less than 20% of patients with cancer benefit from immunotherapy as single agents.4 PAC are among the immunogenically “cold” tumors that does not respond to immunotherapy as single agents because they lack of natural infiltration of antitumor effector T cells and fibrotic stroma in PAC tumor microenvironment (TME).5 6 According to the previously published papers, increasing antigenicity, enhancing effector T cells function, and overcoming immunosuppressive factors in TME are crucial strategies to convert PAC into “hot” tumors.7 8
Besides, the analysis of TILs alone is not sufficient to fully characterize the complex tumor immune microenvironment. Furthermore, immunotherapy resistance can still occur in patients with high levels of TILs. As a result, various indicators have been developed for predicting response to immunotherapy, including tumor mutational burden (TMB), programmed death-ligand 1 (PD-L1) expression, and microsatellite instability (MSI).9–11 Therefore, identifying powerful and robust biomarkers, which aid clinicians in identifying patients suitable for immunotherapy.
To proliferate rapidly and cope with stress, tumor cells undergo metabolic reprogramming, opting for aerobic glycolysis and lactic acid production (Warburg effect).12 TILs are usually subjected to metabolic stress due to tumor cell metabolic dysregulation, resulting in an impaired antitumor immune response. Further tumor cells suppress TILs function by expressing checkpoint molecules.13 Accumulating evidence suggests that immune cell metabolisms within the TME greatly influences immunotherapy response. TME conditions regulate immune cell energy consumption and metabolic reprogramming, often inducing, which compromises cancer cell clearings.14 Using metabolic regulator drugs or antibodies against specific immune receptors to enhance nutrient usage represents a strategy to promote the energetic rewiring of immune cells, enhancing antitumor efficacy.15 Taken together, metabolic changes influence immune function and undermine immunotherapy.16 Targeting cancer cell metabolism shows promise in overcoming immunotherapy resistance and identifying potential markers for predicting immunotherapy response. Therefore, targeting metabolism holds promise for overcoming immunotherapy resistance and identifying potential indicators of immunotherapy response. Systematic understanding of anticancer immune response has emphasized the importance based on immunometabolism. Herein, we characterized metabolic patterns and explored their relationship with the TME in 1,188 patients with PAC, constructing a metabolism-related score (MBS) based on these metabolic patterns for predicting both prognosis and immunotherapy outcomes.
Methods
Collection of PAC data sets and preprocessing
Publicly available gene expression data and clinical annotations of data sets obtained from the Gene Expression Omnibus, The Cancer Genome Atlas (TCGA), International Cancer Genome Consortium (ICGC), and ArrayExpress were used as previously described.17 Patients without survival information were removed from further analysis. For the TCGA-pancreatic adenocarcinoma (PAAD) cohort, RNA sequencing data (fragments per kilobase of transcripts per million mapped reads (FPKM) values) and clinical information were downloaded via the TCGAbiolinks18 package of R software. FPKM values were then transformed into transcripts per kilobase million values similar to microarray results. Batch effects were corrected using the ComBat method from the “SVA” package.19 Further, batch effects were confirmed via principal component analysis (PCA). In total, we analyzed GSE28735, GSE57495, GSE62452, MTAB-6134, and TCGA-PAAD data sets, including 635 patients as the training cohort. We included 295 patients from ICGC-AU and ICGC-CA data sets as the ICGA-validation cohort. Besides, clinical proteomic tumor analysis consortium (CPTAC)-PAAD (n=135), GSE131050 (n=44), and GSE85916 (n=79) were included in the independent validation cohort. In total, 1,188 patients with PAC with survival information were included in this study (online supplemental figure S1A,B and online supplemental table S5).
Supplemental material
Supplemental material
Statistical analysis
The detailed methods and statistics were described in online supplemental methods. All data processing, plotting, and statistical analyses were performed using R V.4.0.4. The Kaplan-Meier method was used to analyze the survival probability, and the log-rank test was used to calculate the significant differences. Continuous variables were assessed for normality first. A Student’s t-test was used to analyze the difference between the two groups for normally distributed data, and the Wilcoxon matched-pairs signed-rank test was used for non-normally distributed data. Two-sided p values<0.05 were considered statistically significant.
Supplemental material
Results
Metabolism and the immune microenvironment are heavily involved in PAC progression
Metabolic changes in tumor cells can affect TME to limit the immune response and current obstacles to cancer treatment. Improving our understanding of these changes can reveal the opportunity to advance the translation of metabolic pathway and immunity. Moreover, it has been established that there exists a complex interplay between these two categories.20 Online supplemental figure S1 illustrates the workflow of the present study. To investigate the significance of metabolism and TME in PAC, the authors analyzed the enrichment score of 112 metabolism-related pathways in normal and cancer tissues using the GSE71729 data set. The results revealed differential regulation of most metabolic pathways between normal and tumor tissues, and most metabolic pathways were closely related to the survival of patients with PAC (online supplemental figure S2C,D and online supplemental table S22). We also analyzed infiltrating immune cells in seven independent cohorts with a total of 930 patients with PAC, and found that patients with high immune cell infiltration had the longest overall survival(OS) (online supplemental figure S3A,B). The infiltrating immune cell types were related to the clinical characteristics of patients with PAC (online supplemental figure 1). These results confirmed the importance of both metabolism and the immune microenvironment in PAC progression, but the relationship between these two therapeutic targets remains largely unknown. Therefore, a systematic and global analysis is needed to reveal the complex relationship between metabolism and TME in PAC.
The metabolic landscape of PAC
First, we constructed metabolic clusters (MB clusters) based on the enrichment score of metabolic pathways in the combined data set with 930 patients to classify these patients into distinct subtypes (online supplemental figure S4A–K). A value of k=3 was identified to have optimal clustering stability based on the similarity displayed by the pathway scores and the proportion of ambiguous clustering measures (figure 1A and online supplemental figure S4L). Principal coordinates analysis (PCoA) revealed a clear distinction among these three cluster subgroups (p<0.001, figure 1B). Cluster 3 had better OS than the other two clusters (figure 1C). Most metabolic pathways were differentially regulated among these three subgroups (figure 1D). In particular, fatty acid degradation, pyruvate metabolism, tyrosine metabolism, and tryptophan metabolism were mostly enriched in cluster 3, whereas glycogen degradation and glycogen biosynthesis were mostly enriched in cluster 1. Furthermore, cluster 3 exhibited the highest levels of infiltration by cytotoxic cells, CD8 T cells, T cells, and plasmacytoid dendritic cells (pDCs) (figure 1E). These cell types play a crucial role in tumor cell clearance and are generally associated with a more favorable prognosis (online supplemental figure 3C). In the present study, CD274, CD276, CTLA4, HAVCR2, IDO1, LAG3, PDCD1, and PDCD1LG2 were selected as immune checkpoint-related factors, while CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF CD4, CD8A, CXCL10, GZMB, PRF1, PTPRC, STAT1, STAT2, STAT3, STAT4, and STAT5A were used as immune activity-related factors.21–23 Cluster 3 had the lowest CD274, IDO1 and CTLA-4 expression but the highest CD8A and LAG3 expression (figure 1F and online supplemental figure S5A–F).
Construction of MB clusters and the landscape of metabolism and immune-cell infiltration in PAC. (A) Consensus clustering matrix of all 930 patients with PAC for k=3. (B) The principle coordinate analysis confirmed the three MB clusters. The circles and error bars indicate the mean and SEs of the mean. (PERM-analysis of variance test with 10,000 permutations). (C) Kaplan-Meier survival analysis of patients with PAC by metabolic subgroup. (D) Metabolic pathways dysregulated among three cluster subgroups. (E) The comparison of the fraction of tumor-infiltrating immune cells in three MB clusters. (F) Immune-activation-relevant genes expression in three MB subgroups. (G) TME-related pathways among three MB cluster subtypes. *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001, ns represents p>0.05. CAFs, cancer-associated fibroblasts; DC, dendritic cell; ICB, immune checkpoint blockade; MB, metabolic clusters; MDSCs, myeloid-derived suppressor cells; mRNA, messenger RNA; NK, natural killer; PAC, pancreatic cancer; PCoA, principal coordinate analysis; pDCs, plasmacytoid dendritic cells; TAMs, tumor-associated macrophages; TLSs, tertiary lymphoid structures; EMT, epithelial mesenchymal transition; TReg, regulatory T cell.
As expected, the panfibroblast TGFβ response characteristics (Pan-F-TBRs), cancer-associated fibroblasts (CAFs), myeloid-derived suppressor cells (MDSCs), tumor-associated macrophages (TAMs), epithelial–mesenchymal transition (EMT), and immune checkpoint blockade (ICB)_resistance signatures, which are the main signatures contributing to immune therapy resistance, were generally higher in patients from cluster 1 compared with cluster 2 and three patients, whereas tertiary lymphoid structures (TLSs) were the highest in cluster 3 (figure 1G). Additionally, cluster 3 had the lowest Tumor Immune Dysfunction and Exclusion (TIDE) score, which was consistent with the survival data (online supplemental figure S5G,H). Furthermore, cluster 3 had the lowest level of T-cell dysfunction and the highest level of CD8 and MSI scores, which contribute to anticancer immunity (online supplemental figure S5I–K). Moreover, the study used the cancer immunity cycle to uncover which step is most involved in killing tumor cells. The results showed that cluster 3 had high levels of dendritic cells (DCs), macrophage, and natural killer (NK) cell recruitment (online supplemental figure S6A), in addition to the highest interferon (IFN)-gamma expression as well as the lowest cell cycle and DNA replication rates (online supplemental figure S6B). These results indicated that patients in cluster 3 had the best survival, with enhanced immune cell infiltration potentially underlying their prognosis and indicating immunotherapy sensitivity.
Construction of metabolic gene clusters and association with the tumor immune microenvironment
To further characterize distinct MB clusters, we identified DEGs between normal and PAC tumor tissues from gene expression profiling interactive analysis database (GEPIA). Then, we obtained 3,876 differentially expressed genes among the three MB clusters by using limma packages in R. Lastly, we took intersections from DEGs and 3,876 genes, and we harvested 2,831 differentially expressed metabolism-related genes (DEMBGs, online supplemental figure S7A). We performed unsupervised clustering of DEMBGs in the training cohort and divided the latter into gene clusters A, B, and C (online supplemental figure S7B–L). The 305 genes positively correlated with the gene cluster were named as MB gene signature A, while the remaining DEMBGs were termed gene signature B (online supplemental figure S7M), (online supplemental table S2). We performed functional enrichment analysis of gene signatures A and B using over representation analysis (online supplemental figure S7M). In order to reduce the noise or redundant genes, we used the Boruta algorithm to perform dimension reduction in the gene signatures A and B. Online supplemental figure S8A presents the expression of the 308 most abundant DEMBGs identified across gene clusters A, B, and C (online supplemental table S3).
PCoA confirmed these three gene clusters (p<0.001, online supplemental figure S8B), and survival analysis indicated that patients in gene cluster C had the best OS, whereas those in A and B had worse prognosis (p<0.001, online supplemental figure S8C). We then explored the relationship between gene clusters and infiltrating immune cells. Gene clusters B and C had higher levels of B cells, DC cells, T cells, and Th17 cells, which play an antitumor role (online supplemental figure S8D). Moreover, gene cluster C had the lowest expression of immunoinhibitory factors IDO1, CD274, HAVCR2 and CTLA-4 (online supplemental figure S8E). The cancer immunity cycle results implied that gene cluster C had the highest level of basophil recruitment and the lowest level of MDSC recruitment (online supplemental figure S8f), in addition to the lowest TIDE score (online supplemental figure S9A). Immunotherapy-related pathways in gene clusters were also analyzed and consistent with above-presented results (online supplemental figure S9B). The consistency among prognostic performance, immune profile, and immunotherapeutic prediction in the three MB-related gene clusters validated our classification.
Development of MBS and performance in prognostic prediction
To quantify the MB landscape, we used a PCA method to construct MBS based on MB gene signatures A and B (figure 2A). We found that MBS was lowest in gene cluster C (online supplemental figure S10A). Using the surv_cutpoint function in the survminer package in R, we stratified patients in the training cohort into two subgroups, MBS-high and MBS-low, and this was further confirmed via PCoA (figure 2B, p<0.001). In the training cohort, the MBS-low group, which greatly overlaps with MB cluster 3, had better overall survival than the MBS-high group (figure 2C, p<0.001), as well as in the CPTAC-PAAD cohort (figure 2D, p<0.01), GSE131050 cohort (figure 2E, p<0.01), SE85916 cohort (figure 2F, p<0.05), and the other two validation cohorts (online supplemental figure S10B,C, p<0.05).
Construction and validation of MBS. (A) The association landscape among MBS and molecular characteristics (MB clusters, gene clusters, immune infiltration subtypes, survival status, and cohorts). Columns showed PAC samples sorted by MBS from low to high. (B) The principal coordinate analysis confirmed the high and low MBS groups. The circles and error bars indicate the mean and SEs of the mean. (PERM-analysis of variance test with 10,000 permutations). (C–F) Kaplan-Meier curves of OS for patients with PAC based on the MBS in the training cohort, CPTAC-mRNA cohort, GSE131050_Yeh cohort, and GSE85916 cohorts. (G) Differences of MBS between two Moffitt subtypes in the GSE131050 data set. (H–J) Comparison of C-index among MBS and 66 published signatures in TCGA-PAAD, ICGC-AU, and meta-cohorts. ICGC, International Cancer Genome Consortium; MB, metabolic clusters; MBS, metabolism-related score; mRNA, messenger RNA; PAC, pancreatic cancer; PCoA, principal coordinates analysis; TCGA, The Cancer Genome Atlas; CPTAC, Clinical Proteomic Tumor Analysis Consortium; OS, overall survival.
In addition, we investigated the relationship between MBS and molecular subtypes of PAC. Previous studies have shown that basal-like subtype patients have worse outcomes and are molecularly similar to basal tumors in bladder and breast cancer.24 Quasi-mesenchymal subtype patients exhibit high expression of mesenchymal-associated genes and also have poor prognosis.25 Interestingly, consistent with previous findings, MBS was significantly upregulated in basal-like and quasi-mesenchymal subtypes of patients with PAC (figure 2G, S10D, p<0.01).
To compare the performance of MBS with other signatures, we comprehensively collected 66 published signatures, including long noncoding RNA (lncRNA) and messenger RNA, associated with various biological processes such as immunotherapy response, autophagy, hypoxia, epigenetic modification, RNA-binding, stemness, immune cell phenotypes, KRAS and TP53 mutated signatures, glycolysis, and drug resistance. These signatures have been reported to be promising prognostic biomarkers for PAC. In this study, we compared the C-index of MBS with these published signatures in TCGA-PAAD, ICGC-AU, and combined cohorts. Notably, MBS exhibited outstanding performance in all three cohorts, which demonstrated its stability and robustness (figure 2H–J). Some models showed weak performance across data sets, possibly due to overfitting during model development. However, our MBS model was dimensionally reduced by multiple algorithms and therefore had better extrapolation potential.
Immune significance and validation in a clinical in-house cohort
Immune cells and immune-related genes play a critical role in cancer immunotherapy, with numerous immune agonists and antagonists being evaluated in clinical oncology. Consistent with this, a low MBS was associated with greater infiltration of B cells, eosinophils, NK cells, T cells, and Th17 cells, and lower regulatory T cell (Treg) infiltration (online supplemental figure S10E). To further advance this research, it is important to understand the expression of immune-related genes and modes of control in different states of the TME. Therefore, we examined the expression of these genes, as well as somatic copy-number alterations (SCNAs) and epigenetic mechanisms. Our results showed that the expression of most immune-related genes varied across MBS subtypes, potentially indicating their role in shaping the TME based on the TCGA-PAAD cohort (figure 3A). To further confirm the clinical applicability of MBS, we employed IPS, ESTIMATE, MCPcounter, xCell, and TIMER algorithms to quantify the immune cell infiltration landscape in our internal cohort. Consistent with the previous results, the low MBS group exhibited higher levels of CD4+T cells, CD8+T cells, B cells, and immunophenotype score (IPS) scores (figure 3B). Notably, the protein expression of PD-1, CD8A, and CD4 were significantly higher in the low-MBS group than in the high-MBS group (figure 3C,D). Furthermore, we found that there was a strong correlation between cytotoxic cells and T cells in both MBS-high and MBS-low subtypes (online supplemental figure S10F). Most metabolic pathways were differentially regulated between the two MBS subgroups (online supplemental figure S10H). Fatty acid degradation, pyruvate metabolism, tyrosine metabolism, and tryptophan metabolism were mostly enriched in the MBS-low group, whereas glycogen degradation and glycogen biosynthesis were most enriched in the MBS-high group.
The landscape of tumor microenvironment analysis between high and low MBS groups. (A) mRNA expression of IM genes (median normalized expression levels); expression versus methylation (gene expression correlation with DNA-methylation beta-value); amplification frequency (the difference between the fraction of samples in which an IM is amplified in a particular subtype and the amplification fraction in all samples); and the deletion frequency (as amplifications) for 59 IM genes by MBS subtypes in The Cancer Genome Atlas-PAAD cohort. (B) Heatmap for immune landscape in Tangdu in-house cohort based on IPS, ESTIMATE, MCPcounter, xCell, and TIMER algorithms between high and low MBS group. (C) Boxplot displays the IHC scores of PD-1, CD8A, and CD4 between two MBS groups based on internal data from Tangdu cohort (n=49) according to staining results. *p<0.05. (D) Representative IHC staining images of PD-1, CD8A, and CD4 between two MBS groups (n=49). Scale bars=100 µm. DC, dendritic cell; IPS, immunophenotype score; MB, metabolic clusters; MBS, metabolism-related score; mRNA, messenger RNA; NK, natural killer; PD-1, programmed cell death-1.
Next, we analyzed the correlation between the MBS and anticancer-related signatures. As expected, Pan-F-TBRs, CAFs, MDSCs, EMT, ICB_resistance, T-cell exclusion, and TIDE score were obviously higher in MBS-high patients than MBS-low patients, whereas the opposite was noted for TLSs, MSI_expr_Sig, and IPS, which contribute to anticancer immunity (figure 4A–J). MBS most significantly correlated with the cancer immunity cycle and immunotherapy-related pathways (figure 4K,L). The alluvial diagram indicated that MBS-low patients overlapped with those in MB cluster 3, whereas MBS-high overlapped with those in MB cluster 1 (online supplemental figure S10E). The findings described above demonstrate the consistent prognostic performance and immune profiles in MB clusters and MBS subgroups, suggesting that MBS may serve as a robust indicator for predicting immunotherapy response in PAC. Moreover, these results suggest that patients with low MBS may have more alternative resources for immune checkpoint inhibitors (ICI) treatment.
The association between MBS and tumor microenvironment. (A–J) Differences of immunotherapy-related score, including Pan_F_TBRs, MDSCs, ICB_resistance, CAFs, EMT, exclusion, TIDE, MSI_Expr, TLS, and IPS signature scores, respectively, in two MBS groups. (L) Correlation between MBS and the steps of the cancer immunity cycle. (K) Correlation between MBS and the enrichment score of immunotherapy-predicted pathways. (M) Raincloud plots depicted the differences of MBS in non-responders (NR) and treatment-naive patients (TN) groups in the SKCM-GSE115978 immunotherapy cohort. (N) Raincloud plots depicted the differences of MBS in NR and responders (R) groups in the BCC-GSE123813 immunotherapy cohort. *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001. BCC, basal cell carcinoma; CAFs, cancer-associated fibroblasts; ICB, immune checkpoint blockade, IFN, intrerferon; IPS, immunophenotype score; MB, metabolic clusters; MBS, metabolism-related score; MDSCs, myeloid-derived suppressor cells; MSI, microsatellite instability; TAMs, tumor-associated macrophages; TIDE, Tumor Immune Dysfunction and Exclusion; TLSs, tertiary lymphoid structures; EMT, epithelial–mesenchymal transition; SKCM, Skin Cutaneous Melanoma.
Immunotherapy outcome prediction by MBS
We sought to evaluate the predictive performance of MBS for immunotherapy outcome in 15 cohorts from bulk-RNA and scRNA-seq aspects. First, we analyzed the two scRNA-seq ICI cohorts, including melanoma (SKCM-GSE115978) and basal cell carcinoma (BCC-GSE123813). As the previous paper described,26 we used 24 patients from the melanoma cohort, including 11 non-responders (NR) and 13 treatment-naive patients (TN), a total of 2,142 malignant cells were obtained after removal of non-malignant cells. Ideally, it would be better to compare the MBS between responders (R) and NR. However, data on R were not available in this cohort, and since untreated patients may include both potential R and NR, MBS comparisons were made between NR and TN. As the results showed that MBS in the NR group was significantly higher than that in the TN group (p<0.001, figure 4M). In the BCC cohort, a total of four NR and six R were included, and 1,826 malignant cells were obtained after removing non-malignant cells. We also found that MBS in the NR group was significantly higher than that in the R subgroup in the BCC cohort (p<0.001, figure 4N). These results further indicated that MBS was robust in predicting the immunotherapy responsiveness of patients with tumor.
Next, we systematically analyzed thirteen bulk-RNA immunotherapy cohorts. MBS-low group patients had a better prognosis and immunotherapy response in the IMvigor210 cohort (figure 5A–C). Moreover, the MBS was negatively related to FMOne mutation and neoantigen burden (figure 5D,E). Patients with nivolumab-treated advanced melanoma cohort were divided based on MBS, and it was found that MBS was significantly higher in the stable disease (SD)/progressive disease (PD) group than in the complete response (CR)/partial response (PR) groups (figure 5G,H). Additionally, MBS-low patients had a greater overall survival than MBS-high patients (figure 5I). In the GSE78220 cohort of patients with anti-PD-1-treated melanoma, MBS was negatively correlated with the cytolytic score and was significantly higher in the PD group than in the PR/CR groups, with MBS-low patients achieving better overall survival (figure 5F and figure 5J,L). The predictive role of MBS was validated in other cohorts, such as mice receiving anti-CTLA-4 and anti-PD-L1 (GSE117358, online supplemental figure S11A,B), patients with immunotherapy-treated chronic lymphocytic leukemia (GSE148476, online supplemental figure S11C,D), and anti-CTLA-4-treated mices (GSE63557, online supplemental figure S11E–G), as well as GSE173839, GSE165252, and GSE168204 (online supplemental figure S11H–J). Besides, we verified the performance of MBS in four immunotherapy PAC cohorts. There was significant difference of MBS between NR compared with R in both CRC and PDAC patients, and higher MBS indicated worse OS and PFS (online supplemental figure S12A–F) in Parikh cohort. In terms of PDAC, NR had higher MBS compared with R, however, there was no statistical difference between NR and R (online supplemental figure S12G,H). The main reason for this problem, we suspect, is that the sample size of patients is too small. Our result showed that short-term survivors of PDAC has higher MBS compared with those long-term survivors of PDAC in Balachandran cohort (online supplemental figure S12I,J). Vaccinated patients showed lower MBS compared with non-vaccinated patients (online supplemental figure S12K,L) in Lutz cohort and CD11b cells had highest MBS, and CD4+T cell and CD8+T cell had lower MBS in Li cohort (online supplemental figure S12M,N). Additionally, we further predicted the immunotherapeutic responsiveness using the TIDE algorithm in nine PAC data sets. A lower TIDE score indicates a better response to immunotherapy. Consistent with previous research findings, we observed a significant positive correlation between MBS and TIDE scores. Furthermore, the MBS in the R group was significantly lower than that in the NR group (online supplemental figure S13A–E).
The performance of MBS in the prediction of immunotherapeutic benefits. (A) Violin plots depicted the differences of MBS in CR, PR, SD, and PD groups in the IMvigor210 cohort with metastatic urothelial cancer. (B) Rate of clinical response (CR/PR and SD/PD) to anti-programmed death-ligand 1 immunotherapy in high and low MBS subgroups in the IMvigor210 cohort. (C) Kaplan-Meier curves for patients with high and low MBS in the IMvigor210 cohort. (D–E) The correlation between MBS and FMOne/neoantigen TMB in the IMvigor210 cohort. (F) The correlation between MBS and cytolytic score in the GSE91061 cohort with melanoma. (G) Violin plots depicted the differences of MBS in CR/PR and SD/PD groups in the GSE91061 cohort. (H) Rate of clinical response (CR/PR and SD/PD) to anti-PD-1 immunotherapy in high and low MBS subgroups in the GSE91061 cohort. (I) Kaplan-Meier curves for patients with high and low MBS in the GSE91061 cohort. (J) Violin plots depicted the differences of MBS in CR/PR and PD groups in the GSE78220 cohort with melanoma. (K) Rate of clinical response (CR/PR and PD) to anti-PD-1 immunotherapy in high and low MBS subgroups in the GSE78220 cohort. (L) Kaplan-Meier curves for patients with high and low MBS in the GSE78220 cohort. (M) Circus plot depicting the performance of MBS in predicting immunotherapy outcome in four data sets compared with already acknowledged signatures. The vertical axis indicated AUC values. (N) The heatmap depicts the performance of MBS in predicting immunotherapy outcome with already acknowledged signatures in four data sets. CR, complete response; GEP, T-cell-inflamed gene expression profile; IMPRES, immune-predictive score; MB, metabolic clusters; MBS, metabolism-related signature; MIAS, MHC 1 association immunoscore; PD, progressive disease; PR, partial response; SD, stable disease; PD-1, programmed cell death-1; AUC, area under curve; CTLA-4, cytotoxic T lymphocyte-associated antigen-4; IDO1, indoleamine2,3-dioxygenase1.
Finally, we compared the performance of MBS with other previously recognized signatures for predicting immunotherapy outcomes in four different data sets. The results demonstrated that MBS outperformed other signatures for predicting immunotherapy response in all four cohorts, indicating its stable and robust predictive performance (figure 5M,N). These findings suggest that MBS has the potential to serve as a valuable biomarker for guiding personalized immunotherapy strategies for patients with PAC.
Potential therapeutics for MBS-high patients
As mentioned above, patients with low MBS are highly sensitive to immunotherapy and have a better prognosis. However, we also need to identify alternative therapeutic targets for patients with high MBS. To address this, we applied three different approaches—Cancer Therapeutics Response Portal (CTRP), profiling relative inhibition simultaneously in mixtures (PRISM), and Genomics of Drug Sensitivity in Cancer (GDSC) to identify candidate drugs for patients with PAC with high MBS. Differential compound sensitivity analysis between MBS-high (top decile) and MBS-low (bottom decile) groups was first performed to identify drugs with lower estimated AUC values in the former (log2 FC>0.1). Then we performed Spearman’s correlation analysis between MBS and AUC value to identify drugs with negative correlation coefficient (r<−0.2 for CTRP and PRISM, p<0.05). We identified six drugs via CTRP (dasatinib, PD318088, selumetinib, paclitaxel, SB-743921, and triazolothiadiazine) and five via PRISM (LY2606368, cobimetinib, epothilone B, dasatinib, and trametinib), with estimated AUC values negatively correlated with MBS (figure 6A,B). Furthermore, we analyzed the maximal inhibitory concentration (IC50) of these compounds based on the GDSC database. Three drugs (dasatinib, paclitaxel, and epothilone B) had a lower estimated IC50 in the MBS-high group (figure 6C–E), suggesting that they may hold promise for the treatment of patients with PAC with high MBS. Based on our analysis of the combat data set derived from GSE14701 and GSE45757, we observed that the MIAPaCa-2 cell line exhibited the highest MBS compared with PATu-8988S and PATu-8988T (online supplemental figure S4A–C), therefore, we generated colony formation using these three PAC cell lines to validate the effect of dasatinib and epothilone B with five different concentrations. Our results demonstrated the rationality of our findings (figure 6G–I).
Identification of potential therapeutic drugs for patients with PAC with high MBS. (A) The results of Spearman’s correlation and differential drugs response analysis of six CTRP-obtained drugs. (B) The results of Spearman’s correlation and differential drugs response analysis of five PRISM-obtained drugs. (C) Comparison of estimated IC50 levels of dasatinib between high and low MBS groups. (D) Comparison of estimated IC50 levels of paclitaxel between high and low MBS groups. (E) Comparison of estimated IC50 levels of epothilone B between high and low MBS groups. (F) Colony-formation assays were conducted based on two PAC cell lines (MIAPaCa-2 and PATu-8988T) treated with dasatinib using gradient concentrations (0 nM, 6.25 nM, 12.5 nM, 25 nM, and 50 nM). (G) Colony-formation assays were conducted based on two PAC cell lines (MIAPaCa-2 and PATu-8988T) treated with epothilone B using gradient concentrations (0 nM, 0.1 nM, 0.2 nM, 0.4 nM, and 0.8 nM). Note: Lower estimated AUC and IC50 values indicate greater drug sensitivity. *p<0.05; **p<0.01; ***p<0.001. (H) Relative colony-formation rate of three PAC cell lines treated with dasatinib using gradient concentrations (0 nM, 6.25 nM, 12.5 nM, 25 nM, and 50 nM). (I) Relative colony-formation rate of three PAC cell lines treated with epothilone B using gradient concentrations (0 nM, 0.1 nM, 0.2 nM, 0.4 nM, and 0.8 nM). CTRP, Cancer Therapeutics Response Portal; PAC, pancreatic cancer; MBS, metabolism-related score; AUC, area under curve.
To explore the underlying mechanisms, we performed high-throughput sequencing on the MIAPaCa-2 cell line, MIAPaCa-2 cells treated with dasatinib, and cells treated with epothilone B. Through high-throughput sequencing, we discovered that dasatinib primarily acts through inhibiting signaling pathways such as tumor necrosis factor (TNF)-α, hypoxia, DNA repair, and E2F_target signaling pathways in PAC cell lines. On the other hand, we found that epothilone B mainly affects signaling pathways including MTORC1, TNF-α, and unfolded_protein_response (online supplemental figure S15E,F). These findings are consistent with the activation of these signaling pathways observed in patients with PAC with high MBS (online supplemental figure S19A–I). Furthermore, dasatinib is an oral, once-daily SRC kinase inhibitor commonly used in the treatment of chronic myeloid leukemia and Philadelphia chromosome-positive acute lymphoblastic leukemia. Our analysis revealed that SRC family genes, including SRC, HCK, FGR, and FYN are significantly upregulated in patients with PAC with high MBS (online supplemental figure S15A–D). This further explains why dasatinib is better efficient in patients with PAC with high MBS. In conclusion, our study suggests that dasatinib and epothilone B may be promising candidates for treating patients with PAC with high MBS, and we provide potential mechanistic explanations for the actions of these drugs. However, further research and clinical trials are needed to validate these results and evaluate the efficacy and safety of these drugs in patients with PAC.
Clinical performances of the MBS
We sought to further validate the clinical performance of MBS in the TCGA-PAAD cohort, which had the most extensive clinical data. Lower MBS was related to better OS, disease-free interval, disease-specific survival, and progression-free interval (figure 7A–D). In patients undergoing their first course of treatment, a low MBS was associated with better clinical responses (figure 7E,F). Next, we analyzed the correlation between the MBS and clinical characteristics of patients with PAC. A low MBS was significantly related to age, survival status, MSI status, TP53 mutation, and KRAS mutation status (figure 7G and online supplemental figure S16A–C). Furthermore, we also analyzed the relationship between MBS and response to radiation therapy and chemotherapy. TCGA samples were classified into radiation-sensitive and radiation-resistant classes based on their reported sensitivity to radiation therapy using the response evaluation criteria in solid tumours (RECIST) classification method and data obtained from published papers.27 The results showed that patients with PAC who were sensitive to both chemotherapy and radiotherapy had lower MBS values compared with those who were insensitive to both therapies (online supplemental figure S16D–F). This finding further emphasizes the importance of MBS as a predictor of treatment response in PAC and suggests that incorporating MBS into personalized treatment strategies may be beneficial.
The clinical performance of MBS for patients with PAC. (A–D) Kaplan-Meier curves of OS, DFI, DSS, and PFI for patients with PAAD based on the MBS groups in the TCGA-PAAD cohort. (E) Rate of high and low MBS subgroups in the clinical response (complete response (CR)/partial response (PR) and stable disease (SD)/progressive disease (PD) according to treatment outcome after first course treatment in the TCGA-PAAD cohort). (F) The MBS difference in the two groups according to the treatment outcome after first course treatment. (G) Correlation between MBS and clinical characteristics of patients with PAAD. DFI, disease-free interval; DSS, disease-specific survival; MBS, metabolism-related score; MSI, microsatellite instability; PFI, progression-free interval; TCGA, The Cancer Genome Atlas; OS, overall survival; MSI, microsatellite instability; PAAD, pancreatic adenocarcinoma.
MBS-associated mechanisms in pancreatic cancer
To further explore the mechanisms underlying MBS, the study compared PAC-specific cancer driver gene expression, genetic alterations, and epigenetically driven transcriptional networks between MBS-high and MBS-low subgroups. Differentially regulated cancer driver genes were observed, such as ACVR2A, BAP1, BRCA2, and TP53 (online supplemental figure S17A). Strong correlations were observed between CD8A and GZMA expression, as well as CXCL9 and IDO1 expression (online supplemental figure S17B). Additionally, somatic mutations were analyzed in the TCGA-PAAD data set, and it was found that more somatic mutations occurred in the MBS-high group, including synonymous and non-synonymous mutations (figure 8A–C and online supplemental figure 18A,B). Fourteen genes were more frequently mutated in MBS-high patients, including KRAS, TP53, and ATOB (figure 8D). Significant co-occurrences were observed between KRAS and TP53 mutations, as well as FBN3 and FLG mutations (figure 8F,G and online supplemental figure S18C,D). Finally, we analyzed the prognosis of patients with PAC with and without mutations in these four genes using the cBioPortal database and found that those with mutations in these genes had a poor prognosis (online supplemental figure S18E–G). We also explored the regulon activity of 20 regulator profiles between high and low MBS groups,28 as well as the correlation between MBS and the expression of 28 known epigenetic regulators29 in patients with PAC. Our results suggested that epigenetically driven transcriptional networks might be important factors for MBS subtypes (figure 8H,I). This further supports the idea that alterations in gene regulation and epigenetic modifications may contribute to the differences in clinical outcomes observed between MBS-high and MBS-low subtypes.
The epigenetically driven transcriptional networks and tumor mutation status in high and low MBS subgroups. (A–C) The correlation between MBS and the synonymous mutations, non-synonymous mutations, and all mutations counts and their distribution in high and low MBS subtypes. (D) Forest plot of differences in gene mutations between low MBS and high MBS groups. (E) Interaction of differentially mutated genes in patients with low and high MBS groups. (F–G) Lollipop diagram showing mutation sites of TP53 and KRAS proteins, respectively. (H) Regulon activity of 20 regulators profiles between high and low MBS groups. (I) The correlation between MBS and the expression of 28-known epigenetic regulators *p<0.05; **p<0.01; ***p<0.001. MBS, metabolism-related score.
The study also revealed that several pathways, including TGF-β, TNF-α, hypoxia, inflammatory response, interleukin-6-JAK-STAT3, the unfolded protein response (UPR), and EMT pathways, were significantly activated in the MBS-high group, while the pancreas β cell pathway was activated in the MBS-low group (online supplemental figure S19A–H). Additionally, the expression of UPR and EMT regulators was upregulated in PAC samples (online supplemental figure S19I). To determine the relative importance of the 308 MBS genes, the researchers used the randomForestSRC and randomSurvivalForest algorithms for feature selection and ranked their importance. The results were verified, as shown in online supplemental figure S20A–D. These findings suggest that dysregulation of multiple signaling pathways and immune evasion may contribute to the differences in clinical outcomes observed between MBS-high and MBS-low groups. Furthermore, these pathways may represent potential targets for personalized treatment strategies.
Potential links between MBS and immune resistance using pan-cancer cohort
The study also investigated the relationship between MBS and immunosuppression using TCGA pan-cancer data sets. First, the researchers analyzed the correlation between MBS and the expression level of immunosuppressive genes. As expected, there was a significantly positive association across multiple cancer types (online supplemental figure S21A). Second, we explored the association between MBS and hallmark pathway enrichment to determine whether immunosuppressive biological functions were enriched in the high MBS group. The results showed that TGF-β, EMT, and Notch signaling pathways, which have been reported to contribute to immunosuppression, were significantly upregulated in the high MBS group (online supplemental figure S21B). Additionally, cancer stemness and intratumoral heterogeneity have also been reported to confer immunosuppressive properties. The study found a positive association between MBS and Intra-tumor heterogeneity (ITH) and five stemness signatures across 30 cancers (online supplemental figure S21C,D), suggesting that these factors may contribute to the observed associations between MBS and immunosuppression.30 In conclusion, the study found that MBS is negatively correlated with anticancer immunity. The researchers also conducted a systematic exploration of potential therapeutic targets in synergy with MBS using 17-CRISPR data sets. We ranked 22,505 genes according to their average Z-score in the 17-CRISPR data sets and found that 23 MBS genes were among the top 10% ranked genes across multiple different CRISPR data sets (online supplemental figure S21E). The top-ranked genes (with negative Z-scores) are immune resistance genes that may promote antitumor immunity when knocked out, while immunosensitive genes are located at the bottom of the list. These findings suggest that targeting these MBS-related genes may represent potential therapeutic strategies for enhancing anticancer immunity in PAC.
Discussion
Clinical trials of immunotherapy for PAC confirmed its critical role in eradicating tumors and improving quality of life of patients.31 However, only a minority of patients are sensitive to immunotherapy owing to mechanisms of immune evasion and suppression as well as competition for basic nutrients and the suppression of immune cell metabolisms.32 Increasing antigenicity, enhance effector T cells function, and overcome immunosuppressive factors in TME are crucial strategies to convert PAC into “hot” tumors. To satisfy the enormous energy demands, cancer cells preferentially use glycolysis rather than oxidative phosphorylation.33 34 The former metabolic pathways provide energy in a rapid manner, although glucose is not fully used. Lactate is one of the main products of glycolysis, and the lactate concentration in tumors is 20−30 times greater than that in normal tissues. This acidic microenvironment suppresses infiltrating immune cells, thus compromising the efficacy of immunotherapy.35 Cancer cells and infiltrating immune cells reprogram their metabolism to adapt to the specific TME, with amino acid metabolism in latter having a similar suppressive effect on immune cell function. Altogether, metabolic pathways modulate the tumor immune microenvironment, suggesting the potential for improving immunotherapy outcomes by targeting metabolic pathways or specific metabolites.3 36 Glycolysis–cholesterol synthesis axis was reported to be related to PAC prognosis and prognostic subtype classifier gene expressions.37 However, a comprehensive characterization of the metabolic landscape in relation to immunotherapy response in PAC is lacking.
In this study, we analyzed 112 metabolism-related pathways to characterize the metabolic landscape of PAC. We classified patients into three metabolic cluster subtypes, with cluster 3 exhibiting the best overall survival and enriched metabolic pathways, including fatty acid degradation, pyruvate metabolism, tyrosine metabolism, and tryptophan metabolism. It is currently believed that different immune cell subsets play distinct roles in antitumor immunity.38 Importantly, cluster 3 also had the highest infiltration levels of cytotoxic cells, CD8+T cells, and T cells, indicating a potent antitumor immune response and potential sensitivity to immunotherapy. In contrast, cluster 1 had the highest expression of PD-L1 and lowest CD8 T cells, suggesting a potential for immune escape and resistance to immunotherapy. The Pan-F-TBRs signature is positively correlated with poor immunotherapy responses and unfavorable prognosis.39 Immunosuppressive cell types, including CAFs, MDSCs, and TAMs, suppress T cells in tumors and promote immune escape.40 In contrast, TLSs improve antigen presentation and are correlated with immunotherapy response.41 42 Pan-F-TBRs, CAFs, MDSCs, TAMs, EMT, TIDE, T-cell dysfunction, and ICB_resistance were mostly lower in cluster 3 than clusters 2 and 1, whereas TLSs were the highest in cluster 3. The TIDE score, which integrates T-cell dysfunction and elimination characteristics, is a poor indicator when predicting immunotherapy responses and prognosis. PD-L1 played the structural carcinogenic roles in “cold tumors”, known as innate immune drug-resistant tumors, including those that are PD-L1 positive in the absence of CD8+T cells.43 Patients with this type of tumor emphasize the importance of considering the presence of TILs in TME in conjunction with the state of PD-L1 in order to predict the immunotherapy response. While recently clinical trials of anti-PD-L1 antibodies as monotherapy have not shown clinical benefit in the majority of PAC.44 As supported by previous literature, metabolic activity has been shown to significantly impact the differentiation and fate of effector T cells.16 Therefore, it is crucial to identify strategies that can modulate metabolic patterns in order to enhance the immune response against tumors.
To evaluate quantitative indicators of PAC metabolic landscape, we developed the MBS using machine learning methods. Patients in the MBS-low subgroup, which highly overlapped with those of metabolic cluster 3, had favorable prognosis. To demonstrate the prognostic prediction performance of MBS, we conducted a comparison with other published signatures using the C-index as a measure of performance. The results showed that MBS had the most robust potency to predict prognosis, outperforming the other signatures. This suggests that MBS may be a valuable tool for predicting prognosis in the future. With regard to the TME, lower MBS was related to greater infiltration of antitumor immune cells, while immunosuppressive indicators were significantly decreased. Further, dysregulated metabolism, TIDE and IPS have good performance in predicting immunotherapy response.45 46 Liu et al reported29 demethylase fat mass and obesity-associated protein (FTO) as an essential regulator used by tumors to escape immune surveillance through regulation of glycolytic metabolism. Motivated by above results, we hypothesized that MBS is a promising predictor of immunotherapy response, with a lower MBS indicating greater sensitivity. Therefore, we validated the predictive performance of MBS for immunotherapy outcome in 15 cohorts from bulk-RNA and single-cell aspects. The results confirmed its value in the clinical decision-making process and low MBS was correlated with a better immunotherapy response. We compared the performance of MBS and other already acknowledged scores for predicting immunotherapy response in kinds of cancer data sets using receiver operating characteristic (ROC) curves. And the results showed that MBS displayed outstanding performance for predicting immunotherapy response in multiple immunotherapy-cohorts, which further demonstrated the stability of predictive performance of MBS. However, despite the robust predictive ability of MBS for immunotherapy response that has been evaluated in 15 immunotherapy cohorts, including four PAC immunotherapy cohorts, the limited sample size in PAC immunotherapy cohorts necessitates further validation of the immunotherapeutic predictive value and prognostic significance of MBS in larger cohorts of patients with PAC who receive immunotherapy. Moreover, it is important to consider the combination or comparison of MBS with other potential biomarkers, such as TMB and tumor-infiltrating immune cells, to enhance its predictive capabilities and explore their complementary roles in PAC immunotherapy. Promising biomarkers would improve therapeutic selection for patients. Therefore, we identified dasatinib and epothilone B were identified as promising therapeutic compounds in these patients and we further validated using experiments.
In addition to glucose and glutamine, solid tumors often consume tryptophan. This amino acid is metabolized to kynurenine by the IDO. Kynurenine promotes the survival ability of tumor cells and supports the generation of Treg cells, which inhibit the antitumor response. Therefore, the increased expression of IDO is associated with tumor progression and poor prognosis.47 Conversely, blocking IDO1 can inhibit the generation of Treg cells or convert these cells into non-inhibitory Th17 cells, thus weakening the immunosuppressive TME.48 49 In the present study, IDO1 expression was significantly upregulated in patients from cluster 1 and the MBS-high subgroup, who also had higher Treg cell infiltration. Chang et al50 reported that excessive glucose consumption by tumor cells can restrict the anticancer capacity of T cells, which is related to a decrease in mammalian target of rapamycin (mTOR) activity, glycolysis, and the production of IFN-γ. Antibodies against immune checkpoint-relevant molecules can restore glucose level in TME and enhance TIL glycolysis and IFN-γ production. Further, Treg cells, which are affected by glucose availability, are associated with unfavorable prognosis in multiple cancer types, with their depletion enhancing effective anticancer immunity.51 52
We further explored the potential links between MBS and immune resistance in patients with PAC. We found that TGF-β, EMT, and Notch signaling pathways, which pathways contributed to the immunosuppression, were significantly upregulated in high MBS tumors. EMT could promote cancer progression and metastasis by modulating the immune escape in the local TME.53 Besides, cancer stemness and intratumoral heterogeneity were reported to confer immunosuppressive properties.30 Remarkably, positive association between MBS and ITH and stemness signatures were found across 30 cancers. Taken together, all these data indicated that MBS was negatively correlated with anticancer immunity, which corroborates the predictive value of MBS.
In summary, our study offers novel insights for the stratification of patients with PAC for immunotherapy based on a metabolism-related indicator. The limitation of the current study was that the selected immunotherapy cohort included melanoma, patients with bladder cancer, and only in a small group of patients with PAC. In the follow-up study, we will further validate the MBS in a larger series of patients with PAC treated by immunotherapy and will explore associated mechanisms through in vivo and in vitro experiments.
Ethics approval
Ethics approval was granted by the Ethics Committee of the Fourth Military Medical University with number 2021-03-082. Participants gave informed consent to participate in the study before taking part.
Acknowledgments
Authors thank National Natural Science Foundation of China for giving financial support for this study.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
YG, RW and JS contributed equally.
Contributors HS, JL, and YS designed and supervised research and obtained research funding, and acts as guarantor. YG performed research, analyzed data and prepared a first draft of manuscript, YG, RW and JS drafted the manuscript. PM and CY collected the data. JM and TZ analyzed the data and edited the manuscript. All authors reviewed and approved the final manuscript.
Funding This study was supported by a grant from the National Natural Science Foundation of China to HS (No. 81372255) and to RW (No. 81902523).
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.