Zhenyu Wang1†Ping Tao2†Peidang Fan1†Jiongyuan Wang3†Tao Rong1Yingyong Hou4Yuhong Zhou5Weiqi Lu3Liang Hong6*Lijie Ma6,7*Yong Zhang3*
Hanxing Tong3*- 1First Affiliated Hospital, Anhui University of Science and Technology, Huainan, China
- 2Department of Laboratory Medicine, Shanghai Traditional Chinese Medicine-Integrated Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai, China
- 3Department of General Surgery, Zhongshan Hospital, Fudan University, Shanghai, China
- 4Department of Pathology, Zhongshan Hospital, Fudan University, Shanghai, China
- 5Department of Medical Oncology, Zhongshan Hospital, Fudan University, Shanghai, China
- 6Department of General Surgery, Shanghai Fifth People’s Hospital, Fudan University, Shanghai, China
- 7Department of Liver Surgery, Renji Hospital, School of Medicine, Shanghai Jiaotong University, Shanghai, China
Introduction: The exploration of lipid metabolism dysregulation may provide novel perspectives for retroperitoneal liposarcoma (RPLS). In our study, we aimed to investigate potential targets and facilitate further understanding of immune landscape in RPLS, through lipid metabolism-associated genes (LMAGs) based prognostic model.
Methods: Gene expression profiles and corresponding clinical information of 234 cases were enrolled from two public databases and the largest retroperitoneal tumor research center of East China, including cohort-TCGA (n=58), cohort-GSE30929 (n=92), cohort-FD (n=50), cohort-scRNA-seq (n=4) and cohort-validation (n=30). Consensus clustering analysis was performed to identify lipid metabolism-associated molecular subtypes (LMSs). A prognostic risk model containing 13 LMAGs was established using LASSO algorithm and multivariate Cox analysis in cohort-TCGA. ESTIMATE, CIBERSORT, XCELL and MCP analyses were performed to visualize the immune landscape. WGCNA was used to identify three hub genes among the 13 model LMAGs, and preliminarily validated in both cohort-GSE30929 and cohort-FD. Moreover, TIMER was used to visualize the correlation between antigen-presenting cells and potential targets. Finally, single-cell RNA-sequencing (scRNA-seq) analysis of four RPLS and multiplexed immunohistochemistry (mIHC) were performed in cohort-validation to validate the discoveries of bioinformatics analysis.
Results: LMS1 and LMS2 were characterized as immune-infiltrated and -excluded tumors, with significant differences in molecular features and clinical prognosis, respectively. Elongation of very long chain fatty acids protein 2 (ELOVL2), the enzyme that catalyzed the elongation of long chain fatty acids, involved in the maintenance of lipid metabolism and cellular homeostasis in normal cells, was identified and negatively correlated with antigen-presenting cells and identified as a potential target in RPLS. Furthermore, ELOVL2 was enriched in LMS2 with significantly lower immunoscore and unfavorable prognosis. Finally, a high-resolution dissection through scRNA-seq was performed in four RPLS, revealing the entire tumor ecosystem and validated previous findings.
Discussion: The LMS subgroups and risk model based on LMAGs proposed in our study were both promising prognostic classifications for RPLS. ELOVL2 is a potential target linking lipid metabolism to immune regulations against RPLS, specifically for patients with LMS2 tumors.
Introduction
Retroperitoneal liposarcoma (RPLS) is a rare type of mesenchymal tumor, but the most common subtype of retroperitoneal sarcoma (1). It is also characterized by accumulation of intracellular lipid, induction of adipocyte-specific genes (2), dismal immunotherapy response and poor prognosis (3), since the significant challenges of large tumor size and adjacent organ involvement (4). In addition, other therapeutic strategies, such as combination chemotherapy and molecular targeted drugs also have limited efficacy due to intrinsic chemo-resistance, even in histology-tailored neoadjuvant chemotherapy (5). Therefore, novel strategies are needed to improve the therapeutic condition of RPLS.
Aberrant lipid metabolism and lipid metabolism reprogramming are critically involved in drug resistance (6), the adaptation of immune microenvironment (7), energy supplement, cell signaling (8) and regarded as a new hallmark of tumor ecosystem (9). Emerging evidence indicated that targeting the lipid metabolism pathway was an attractive tumor treatment strategy (10). Previous studies have indicated that lipid metabolism-associated genes (LMAGs) exhibit potent capability in predicting the prognosis of various tumors (11–14). However, lipid metabolism dysregulation in patients with RPLS remains unknown.
Immunotherapy has been extensively studied as a promising treatment, but has had limited therapeutic benefit in RPLS, which is considered as a “non-immunogenic” and highly variable tumor (15). Increasing evidence suggest that the alterations in lipid metabolism, including metabolite abundance and accumulation of lipid biomolecules, lead to local immunosuppression in the tumor microenvironment (16). However, the association between abnormality of lipid metabolism and immune microenvironment remains obscure in RPLS.
Elongation of very long chain fatty acids protein 2 (ELOVL2), an enzyme that catalyzes the elongation of fatty acids with chain lengths greater than 18 carbons. Research has shown that ELOVL2 is involved in the maintenance of cellular homeostasis in normal cells (17). Specifically, ELOVL2 was implicated in the regulation of autophagy (18) and the activity of the mTOR signaling pathway (19), which played a key role in the regulation of cell growth and proliferation. It has been suggested that the decline in ELOVL2 expression with age may contribute to the aging process and age-related diseases. Furthermore, mutations in the ELOVL2 gene have been associated with intellectual disability and developmental delay (20). However, further research was still needed to fully elucidate the functions of ELOVL2 in normal cells and its potential implications for cancer treatment.
Therefore, in this study, we explored the role of lipid metabolism dysregulation in RPLS through LMAGs related immune landscape, using multiple bioinformatics methods. A novel LMAGs based prognostic risk model was established and validated in independent cohorts. To the best of our knowledge, this is the first study to promote the understanding and clinical applications about lipid metabolism dysregulation and serve as a reliable reference for further developing target in RPLS.
Materials and methods
Patient and clinical specimens
The cohort-FD consists of 50 RPLS patients (34% female and 66% male) with a mean age of 55 years. In total, 50 tumor samples were surgically resected and collected in Zhongshan Hospital, Fudan University between 2018 and 2020. For scRNA-seq, four fresh surgical specimens (four primary tumors and matched PBMC) were sequenced and incorporated in further analyses. All samples were confirmed by pathologists through both cytological detections during the surgery and the paraffin section after surgery. Clinical information, including demographics and tumor clinicopathologic characteristics of all cohorts were summarized in Supplementary Tables 1 and 2.
Database search
RNA transcriptome sequencing data of 150 RPLS patients and detailed characteristics were obtained from the TCGA (https://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases. We excluded RPLS samples with no complete expression profile data and unknown overall survival (OS) or living status, and included clinical features such as gender, pathological grade, tumor stage and survival status in this study. Additionally, all the obtained data, TPM (transcripts per kilobase of exon model per million mapped reads) values, were normalized using the log2 (TPM + 1) transformation.
Consensus clustering
Firstly, 135 genes were found to be associated with the prognosis of RPLS through the univariate Cox regression analysis. Consensus clustering was performed according to the expression matrix of the 135 genes using the R package “Consensus Cluster Plus”.
Construction, validation and evaluation of risk model based on LMAGs
Least absolute shrinkage and selection operator (LASSO) analysis was performed to downsize the OS and DFS related genes previously filtrated using “glmnet” R package. The minimum lambda value was defined as the optimal value. The genes applied for the establishment of risk model was enrolled by multivariate Cox regression analysis. OS related risk score of each patient in each cohort was calculated as: OS related risk score=1.866577428*ACOT7-0.040721477*ARSJ-0.451127087*ARSK+ 0.479584604*CPT1B-0.000114123*CYP21A2 + 0.424535615*ELOVL2-1.414705024*FDX2-0.32557149*GSTM4 + 0.692304756*HACD1-0.670006158*HSD17B14-0.763286635*MTMR8-0.734212004*ORMDL2 + 1.159963024*TNFRSF21. DFS related risk score of each patient in each cohort was calculated as: DFS related risk score=0.468299088*ACOT1-0.286593499*FABP6. Patients were divided into high and low risk groups according to the medium value. ROC and Martingale residuals method were used to evaluate the predictive efficiency of model.
Functional analysis
Differentially expressed genes (DEGs) between two LMSs were visualized using R package “Limma”. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed through “clusterProfifiler” R package and visualized by Metascape5 (21). Based on “GO biological process” gene dataset was downloaded from molecular signature database, Gene Set Enrichment Analysis (GSEA) was conducted to analyze the difference between subtypes. Compared with LMS2, the upregulated differential genes in LMS1 was visualized by PPI network through the Search Tool for the Retrieval of Interacting Genes (STRING) online tool and the minimum required interaction value was set as 0.7 (22).
cBioPortal analysis
cBioPortal for cancer genomics (23) (cBioPortal, http://www.cbioportal.org, version v3.2.11) is an open-access online tool integrating the raw data from large scale genomic projects. In this study, cBioPortal was used to visualize the gene alteration of potential antigens against tumors in cohort-TCGA, including the correlation between ELOVL2 gene expression and DNA methylation.
TIMER analysis
Tumor Immune Estimation Resource (24) (TIMER, https://cistrome.shinyapps.io/timer/) is a comprehensive resource for the systematical analysis of the immune infiltrates across diverse cancer types. In this study, TIMER was used to visualize the correlation between antigen-presenting cell (APC) infiltration and the expression of the identified potent antigens. The partial Spearman’s correlation was used to perform purity adjustment. Spearman correlation analysis was used to analyze the correlation between the abundance of APCs and the expression of the selected antigens. Statistical significance was set at P < 0.05.
Construction of immune landscape
The immune score and stromal score of each sample in cohort-TCGA was calculated by the “estimate” package in R. The proportion of the 22 types of immune cells in the tumor microenvironment (TIME) of each sample was evaluated via the CIBERSORT algorithm in R software (25).
Intra-cohort immune classifications
Unsupervised clustering of samples in each cohort was performed based on the metagene Z-score for the included populations of MCP-counter using R software, with the Euclidian distance and Ward’s linkage criterion, using the gplots package. Cohort-TCGA and Cohort-GSE30929 were further divided into 5 groups (SIC-A, B, C, D and E) (26).
Establishment of LMAGs based nomogram
Univariate Cox regression analysis was performed to evaluate the prognostic value of identified signatures and clinicopathological features. Multivariate Cox regression analysis was used to further determine the independent prognostic factors. Two nomograms were established by the “rms” package for predicting OS and DFS. The C-index and calibration plot were constructed to estimate the accuracy and consistency of the prognostic models.
Gene signatures for the functional orientation
The gene signatures used to determine the functional orientation were reported as previously described. Each signature was summarized as the following: immunosuppression (CXCL12, TGFB1, TGFB3 and LGALS1), T cell activation (CXCL9, CXCL10, CXCL16, IFNG and IL15), T cell survival (CD70 and CD27), regulatory T cells (FOXP3 and TNFRSF18), major histocompatibility complex class I (HLA-A, HLA-B, HLA-C, HLA-E, HLA-F, HLA-G and B2M), myeloid cell chemotaxis (CCL2), and tertiary lymphoid structures (CXCL13).
Estimation of TLS and immune cell enrichment
12 chemokines were highly expressed by TLS, including CCL2, CCL3, CCL4, CCL5, CCL8, CCL18, CCL19, CCL21, CXCL9, CXCL10, CXCL11 and CXCL13, was applied as the gene signature of TLS. The enrichment score of TLS was calculated by single-sample gene set enrichment analysis (ssGSEA) method as implemented by R-package (27).
Weighted gene co-expression network analysis
WGCNA algorithm was used to identify the hub genes among model genes (28). Gene co-expression modules were identified after a weighted gene co-expression network, and the association between gene network and clinical phenotype were also explored. The WGCNA-R package was applied to establish the co-expression network of all genes in the cohort-TCGA, and the genes with variance within the first 5000 were identified by the algorithm for subsequent analysis. The soft-threshold β was determined by the function “sft$powerEstimate”. The weighted adjacency matrix was transformed into a topological overlap matrix (TOM) to estimate the network connectivity, with hierarchical clustering being used to create the clustering tree structure of the TOM. Different branches of the clustering tree indicated different gene modules with different colors. Tens of thousands of genes were classified into modules based on having similar expression patterns (using their weighted correlation coefficients).
Single-cell RNA sequencing
The single cell suspensions were converted to barcoded scRNA-seq libraries using the Chromium Single Cell 30 Library, Gel Bead & Multiplex Kit, and Chip Kit (10x Genomics), aiming for 6,000 cells per library. Samples were processed using kits pertaining to V2 barcoding chemistry of 10x Genomics. Single samples were always processed in a single well of a PCR plate, allowing all cells from a sample to be treated with the same master mix and in the same reaction vessel. For each experiment, all samples were processed in parallel in the same thermal cycler. Libraries were sequenced on an Illumina HiSeq4000, and mapped to the human genome (buildGRCh38) or to the mouse genome (build mm10) using CellRanger software (10x Genomics, version 3.0.2).
Single-cell transcriptome data processing
The output of the cell-gene count matrix was processed with the Seurat (v 3.1.0) package of R software (version 3.6.1) for quality control and down-streaming analysis. Low-quality cells with < 200 genes or with > 40% mitochondrial genes were removed from the analysis. As the cells from tumor and adjacent normal tissues were loaded in batch for each patient, the data for each patient as individual Seurat Object. The Seurat object for each patient was integrated with the harmony algorithm (R package, Harmony, version 1.0). The top 50 principal components (PCAs) were used for graph-based clustering to identify a distinct group of cells at the indicated resolution. In the subgroup analysis, significant PCAs identified with the ElbowPlot() function were used for graph-based clustering for each cell cluster to identify subgroup cells based on the t-SNE analysis (29). The cell types of the identified cells were defined based on their expression of the canonical marker genes: T cells (CD3D, CD3E, CD4, and CD8A), B cells (MS4A1, CD79A, JCHAIN and CD19), NK cells (NCAM1, NKG7, KLRD1 and NCR1),Monocyte cells (CD14, FCN1, VCAN and CD300E), Neutrophil cells (FCGR3B, S100A9, S100A8 and CXCR1), Macrophages cells (CD163, CD68, C1QA and CSF1R), Dendritic cells (HLA-DQB2, XCR1, CD1C, and CLEC10A), Mast cells (TPSAB1, TPSB2, KIT and CPA3), Smooth muscle cells (ACAT2, TAGLN, MYL9 and MYH11), Endothelial cells (PECAM1, VWF, CLDN5 and PTPRB), Fibroblasts (DCN, COL1A1, PDPN and COLA2).
Immunohistochemistry
Serial FFPE sections (4 um) were deparaffinized in xylene and then rehy- drated in 100%, 90% and 70% alcohol successively. Antigen unmasking was performed with a preheated epitope retrieval solution (100X citrate buffer, pH 6.0), and endogenous peroxidase was inactivated by incubation in 3% H2O2 for 20 mins. Next, the sections were preincubated with 10% normal goat serum and then incubated overnight with the following primary antibodies: anti-ELOVL2 antibody (1:50, 20308-1-AP, Proteintech); anti-CD3 antibody (1:100, ab16669, Abcam); anti-CD8 antibody (1:100, ab17147, Abcam). Next, sections were incubated with the corresponding horseradish peroxidase (HRP)-conjugated secondary antibodies (ready-to-use, MP-7451 and MP-7452, VectorLab) for 30 mins at room temperature and and development with DAB substrate (Vector Laboratories). Sections were counterstained with hematoxylin. Slides were scanned using PANNORAMIC Digital Slide Scanners, and QuPath software was used to quantify positive staining cells.
Multiplexed immunohistochemistry
mIHC was performed in four DDLPS cases of cohort-sc-RNA seq. Briefly, Fresh tumor tissues were fixed in 4% paraformaldehyde solution and embedded in paraffin. FFPE slides were made using 4 μm sections of the tumor samples. Deparaffinization and rehydration were performed with xylene and ethanol respectively, followed by microwave antigen retrieval using heated citric acid buffer (pH 6.0) for 10 minutes and endogenous peroxidase blocking in 3% H2O2 for 20 minutes. Blocker/Diluent was used to block nonspecific binding sites. Afterward, relevant primary antibodies were incubated for 1 hour at room temperature, followed by the corresponding secondary antibodies for 30 minutes. Slides were then incubated with fluorescein TSA plus for 10 minutes (Akoya Boscience, NEL861001KT), after which microwave antigen retrieval was repeated with the above steps until the last antibody was added. After multiplexing, DAPI (Sigma, D9542) was used to stain the nuclei. Antibodies and fluorescent dyes used for multiplexing are listed in Supplementary Table 3. The slides were scanned by Vectra 3 automated high-throughput multiplexed biomarker imaging system (Akoya Phenoimager HT). Immune cells were classified into the following types: T cells (CD3+), B cells (CD20+), DC cells (CD11b+), NK cells (CD57+), Mac (CD68+) and ELOVL2+ cell.
Statistical analysis
SPSS 22.0 (SPSS Inc., Chicago, IL, USA) and R 4.0.4 (R Foundation for Statistical Computing, Vienna, Austria; http://www.r-project.org/) were used for all statistical analyses. Univariate and multivariate Cox regression analyses, ROC curve analysis and K-M survival analysis were performed by R software and the corresponding R packages. The continuous data are expressed as the mean ± standard deviation (SD). The Wilcoxon test was used for comparisons between the two groups, and the Kruskal-Wallis test was used for comparisons of prognosis between groups. Except for the special instructions, all statistical tests were two-tailed, and a P < 0.05 was considered to be statistically significant.
Results
Identification of multi-omics landscape and prognostic LMAGs in RPLS
The whole study process, including consensus clustering, immune landscape visualization, DEGs analysis, nomogram construction and WGCNA analysis were systematically evaluated and depicted in the workflow chart, as shown in Supplementary Table 4.
To systematic appraise the lipid metabolism dysregulation in RPLS, 741 LMAGs were obtained from the MSigDB database. Initially, we encapsulated the incidence of copy number variations and somatic mutations in cohort-TCGA, where 14862 amplified genes were screened to identify potent antigens (Figure 1A), while missense mutation was the most common variant classification (Figure 1B). Additionally, we investigated the incidence of CNV among these LMAGs. By evaluating the frequency of CNV status, existing CNV alterations of LMAGs were also clarified, and the top 35 genes in amplified or deleted CNV status were summarized (Supplementary Figure 1A).
Figure 1 Identification of multi-omics landscape and prognostic LMAGs in RPLS. (A) The chromosomal distribution of the aberrant copy number genes in RPLS of cohort-TCGA. (B) Genetic profile of cohort-TCGA. (C) Unsupervised clustering of cohort-TCGA based on the LMAGs. (D) Heatmaps integrating 135 LMAGs and associated clinicopathologic features in cohort-TCGA. (E) Kaplan-Meier survival curve of OS and DFS for LMSs in cohort-TCGA. (F) Circos plot depicting the location on chromosomes and expression level of 8 overlapping OS-related and DFS-related LMAGs. (G) Triangle heatmap showing the correlation features among 8 overlapping LMAGs in cohort-TCGA. (H) Correlation and prognosis analysis showing the inter-gene communications among eight LMAGs.
Subsequently, according to the transcriptomic data, consensus clustering was performed to cluster RPLS patients into two lipid metabolism subgroups (LMSs) (K = 2) (Figure 1C) based on 135 overall survival (OS) related genes by univariable Cox analysis (Supplementary Figure 1B). 31 and 27 patients were clustered into LMS1 and LMS2, respectively (Supplementary Figure 1C). Heatmap visualization also indicated that prognostic LMAGs profiles differed significantly between LMS1 and LMS2, while LMS2 was enriched with LMAGs (Figure 1D), but predicted poor prognosis (Figure 1E). Interestingly, two subtypes harbored heterogenous somatic mutations profiles, demonstrated ATRX and B4GALNT1 as the genes with the highest mutation frequency in LMS1, but MUC16 in LMS2 (Supplementary Figures 1D, E). However, not significant difference was found in genome altered fraction and mutation counts between LMS1 and LMS2 (Supplementary Figures 1F, G).
Moreover, we also investigated 26 disease free survival (DFS) related LMAGs through univariate Cox analysis (Supplementary Figure 1H). Intersecting the results of OS-related and DFS-related genes, 8 overlapping LMAGs (GSTM4, GRHL1, PI4K2B, GK3P, ARSK, ELOVL2, NR1H4, and ABHD4) were excavated and eligible for further screen prognostic relevant antigens (Supplementary Figure 1I). The location on chromosomes and expression levels were visualized in Figure 1F. The regulatory network described the comprehensive landscapes of the 8 LMAGs indicating their interactions, correlation feature (Figure 1G) and prognostic values (Figure 1H). These findings indicated that the LMAGs classified RPLS patients into two subtypes with different molecular features and prognosis.
Heterogeneous functional enrichment and immune landscape in LMAGs subtypes
To better understand the innate difference of survival and underlying signaling mechanisms between LMS1 and LMS2, DEGs and functional enrichment analyses were performed, respectively. A total of 4144 DEGs were identified, of which 3586 genes were downregulated and 558 genes were upregulated in LMS2, as compared with LMS1 (Figure 2A). GO enrichment analysis indicated that these up-regulated DEGs were involved in positive regulation of immune system process, immune response and regulation of immune system process (Figure 2B). Similarly, KEGG enrichment analysis also validated these pathways associated with cytokine-cytokine receptor interaction, chemokine signaling pathway and B cell receptor signaling pathway, of which are part of the immune response (Figure 2C). Meanwhile, PPI analysis further confirmed 15 submodels, all of which were closely associated with immune and metabolism (Supplementary Figure 2A). These results depict that the expression of LMAGs is closely related to immune-related biological processes, confirming that lipid metabolic reprogramming is significantly associated with tumor immune microenvironment (TIME) in RPLS.
Figure 2 Heterogeneous functional enrichment and immune landscape in LMAGs subtypes. (A) Volcano plot depicting the gene expression difference between LMS1 and LMS2. (B, C) Bubble diagram showing the biological processes enriched by GO analysis and the signaling pathways enriched by KEGG analysis. (D) The comparison of stromal score, immune score and ESTIMATE score calculated by ESTIMATE algorithm. (E) Abundance of immune-related molecular signatures evaluated by XCELL indicating significant differences between LMS1 and LMS2. (F) Estimation of TLS based on 12-chemokine signature between LMS1 and LMS2. (G) Heatmap showing the transcriptomic expression of 12-chemokine involving TLS according to LMSs. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.
To further evaluate the dysregulation of immune and metabolism remodeling involved in RPLS, a series of TIME profiles was conducted. Firstly, we illustrated the distribution feature of previously reported six pan-cancer immune subtypes (C1-C6) (30), of which LMS1 and LMS2 were mostly clustered into C1 (Wound Healing), C2 (IFN-γ Dominant), C3 (Inflammatory), C4 (Lymphocyte Depleted) and C6 (TGF-β Dominant) (Supplementary Figure 2B). Intriguingly, C3 and C6 presented a higher proportion in LMS1, while they were associated with better outcomes. Accordingly, immune score, stromal score and ESTIMATE score were also calculated using the ESTIMATE algorithm (Figure 2D). Remarkably, LMS1 demonstrated significantly higher TIME scores and better prognosis than those in LMS2, implying that favorable immune components and immune-related molecules were abundant in LMS1.
Next, the relationship between LMSs and 34 infiltrated immune cells was further explored. In concordance with previous TIME scores, LMS1 harbored more infiltrated DC, B cells, CD4+Tem, macrophages, monocytes and NKT than those in LMS2 (Figure 2E), indicating the significant impacts of these immune cells in the progression of RPLS. However, the complete immune response involves a close combination of multiple events, not only the infiltrated immune cells (4). Thereafter, we further calculated and compared the immune activity score of each step through TIP analysis. Similarly, the abundance of anti-tumor immune cells was significant higher in LMS1 than those in LMS2, as well as in Step1, Step3 and Step5 (Supplementary Figure 2C). In addition, we explored substantial differences in TLS-associated 12-chemokine signature between LMS1 and LMS2 (Figure 2F), while the expression of lymphoid-structures-associated B-cell-specific chemokine CXCL13 was notably higher in LMS1 (Figure 2G).
Accumulating evidence indicated that tumor with high TMB level predicted better efficacy of immunotherapy (31). We then estimated the value of TMB in both LMS1 and LMS2, but not significant difference was found (Supplementary Figure 2D). Intriguingly, patients in LMS1 with low TMB demonstrated a satisfactory survival benefit. Considering the importance of immune checkpoint inhibitors in the treatment of solid tumor (3), we further examined the differences in immune checkpoint profiles and found notably substantial differences in CD28, CD40, CD86, HAVCR2 and PD-1, between these two subtypes (Supplementary Figure 2E). Immunogenic cell death (ICD) has been classified as a form of regulated cell death (RCD) that is sufficient to activate an adaptive immune response (32). We next identified ICD-related genes and analyzed the expression patterns. Importantly, we discovered that significant higher expression of FPR1, TLR4 and CXCL10 were enriched in LMS1 (Supplementary Figure 2F). Taken together, these findings demonstrated the unique characteristics of TIME within two LMSs, offering a conducive complement to previous studies.
Identification of immune gene co-expression modules and immune hub genes of RPLS
The immune gene co-expression module was designed and applied to classify immune-related genes, whose expression may significantly influenced the effectiveness of potential targets (33). Therefore, we re-analyzed and enrolled immune-related genes to establish gene modules through WGCNA (Figure 3A). The soft threshold was set at four in the scale-free network (Supplementary Figure 3A). The representation matrix was converted to adjacency and next to a topological matrix. The average-linkage hierarchy clustering approach was applied with a minimum of 30 genes for each network according to the standard of a hybrid dynamic shear tree. Eigengenes of each module were calculated and the close modules were integrated (height = 0.25, deep split = 3 and min module size = 30). Notably, six gene modules were identified (Figures 3B, C), and correlation feature was also visualized (Supplementary Figure 3B). In addition, the module eigengenes in LMS1 were significantly higher in yellow and blue modules (Figure 3D). Moreover, the prognostic analysis indicated that eigengenes in the brown module was significantly associated with OS in RPLS (Figure 3E).
Figure 3 Identification of immune gene co-expression modules and immune hub genes of RPLS. (A, B) Gene co-expression network analysis based on the immune-related genes. (C) Dot plot of the co-expression gene modules. (D) The comparison of identified gene modules according to LMSs. (E) Forest maps depicting prognosis prediction value of six modules in RPLS. (F) Bubble plot showing the top 10 KEGG terms enrichment pathways from the brown module. (G) Kaplan-Meier survival curve of OS for two risk-score groups in cohort-TCGA. Distribution of survival status according to risk score in RPLS. Heatmap illustrating the expression panel of three immune hub genes according to the risk score. *P < 0.05, **P < 0.01, and ****P < 0.0001.
Further functional enrichment analysis illustrated that genes involved in cytokine-cytokine receptor interaction and T cell receptor signaling were enriched in the brown module (Figure 3F). The brown module (361 immune-related genes) was further selected, and three of them (PLCG1, ZC3HAV1L and NFAT5) were filtered to build the risk score through LASSO algorithm (Supplementary Figures 3C, D). Patients were classified into the high-risk and low-risk groups, while high-risk group predicted unfavorable OS (P = 0.02) (Figure 3G). The area under the receiver operating curve (AUC) was 0.75, indicating a good accuracy of the model (Supplementary Figure 3E). Taken together, this risk model may serve as a novel tool for prognostic predicting in RPLS, based on the immunogenic genes co-expression network.
Development of survival and relapse risk models and nomograms based on LMAGs
Given the significant biological roles of LMAGs in lipid metabolic reprogramming, the association between LMAGs-related risk score and the prognosis needed thoroughly study. Thus, two prognostic models were conducted for OS and DFS, respectively.
24 LMAGs were found to be considerably linked to the OS of patients through LASSO regression analysis (Figures 4A, B), 13 of which were tested and selected for the risk score model from multivariate Cox analysis (Supplementary Figure 4A). We also investigated the relationship between risk score and survival status, and the low-risk subgroup harbored significant more alive statuses (Figure 4C) and better OS than those in high-risk subgroup (Figure 4D). Specifically, this OS related model indicated a great accuracy with AUC values of 0.94 in 1 year, 0.97 in 3 years and 0.97 in 5 years (Figure 4E).
Figure 4 Development of survival and relapse risk models and nomograms based on LMAGs. (A, B) LASSO regression analyses for screening optimal OS-related LMAGs. (C) Distribution of survival status according to risk score in RPLS. Heatmap illustrating the expression panel of 13 OS-related LMAGs according to the risk score. (D) Kaplan-Meier survival curve of OS for two risk-score groups in cohort-TCGA. (E) Time-dependent ROC curve analyses of the LMAGs risk signature in cohort-TCGA. (F) The comparison of stromal score, immune score and ESTIMATE score calculated by ESTIMATE algorithm in two risk-score groups. (G) Nomogram predicting OS for RPLS. *P < 0.05.
To evaluate and validate the universality of this LMAGs-related prognostic model from cohort-TCGA, an independent dataset (cohort-FD) was performed as a validation cohort. The risk score of each case in cohort-FD was calculated using the same formula as that for the cohort-TCGA. Similarly, patients in high-risk group suffered unfavorable OS than those in low-risk group (Supplementary Figure 4B). In addition, the AUC values of this model according to ROC analysis were 0.7 in 1 year, 0.8 in 3 years, and 0.85 in 5 years (Supplementary Figure 4C).
Since this model predicted great potency in clinical prognosis based on LMAGs, we further investigated the correlation between risk score and TIME in cohort-TCGA. In concordance with previous findings, the high-risk subgroup was characterized with a significantly lower immune score (Figure 4F), and infiltrated less CD8+ T cells and plasma cells, but more resting memory CD4+ T cells and Tregs than those in low-risk subgroup (Supplementary Figure 4D), implying that high-risk subgroup was considered with the characteristics of immunosuppressive status.
Meanwhile, we also explored the value of model between the two subgroups stratified by different clinical features. Univariate Cox analysis indicated that patients with dismal OS were authenticated with larger tumor size, chemotherapeutic efficacy and high risk score (Table 1). Multivariate Cox analysis further confirmed that all of them were independent risk factors (Table 1). Subsequently, we developed a nomogram for OS prediction using these two clinical parameters and the LMAGs-based risk scores. A calibration plot for internal validation of the nomogram presented excellent consistency between the nomogram-predicted probability and actual observations of the 1-, 3-, and 5-year OS (Figure 4G).
Table 1 Univariate and multivariate analysis of OS in cohort-TCGA (n=58).
Given the significant high local recurrence rate in clinical treatment of RPLS (34), 21 LMAGs were also found to be considerably linked to the DFS of patients through LASSO regression analysis (Supplementary Figure 4E), 2 of which were tested and selected for the prediction model in the multivariate Cox analysis (Supplementary Figure 4F). The association between risk score and recurrence status was next evaluated, and the low-risk subgroup harbored significant more alive statuses (Supplementary Figure 4G) and better DFS than those in high-risk subgroup (Supplementary Figure 4H). Similarly, this DFS related model also presented a great accuracy with AUC values of 0.72 in 1 year, 0.89 in 3 years and 0.79 in 5 years (Supplementary Figure 4I). In addition, univariate Cox analysis indicated that patients with worse DFS were authenticated with tumor residue, poor chemotherapeutic efficacy and high risk score (Table 2). Multivariate Cox analysis further confirmed that dismal chemotherapeutic efficacy and high risk score were independent risk factors (Table 2). Furthermore, a nomogram for DFS prediction was conducted with excellent consistency (Supplementary Figure 4J).
Table 2 Univariate and multivariate analysis of DFS in cohort-TCGA (n=58).
Identification of lipid metabolism-associated targets
To explore key genes that functioned as potential candidates for RPLS, we further systematically screened and identified two candidates (NR1H4 and ELOVL2) with both gene amplification and mutation, which were also associated both OS and DFS from the 135 LMAGs (Figure 5A). Given the essential role of antigen-presenting cells (APCs) in the function of immunological reaction, we also evaluated the relationship of these two genes with APCs using TIMER analysis (35). Intriguingly, ELOVL2, but not NR1H4, was identified with closely related to APCs (Pearson correlation coefficient > 0.3; Figure 5B; Supplementary Figure 5A), which could serve as a potential target and triggered strong immune response. In addition, similar results were also found in TCGA-SARC (Supplementary Figure 5B). Notably, survival analysis demonstrated that high mRNA expression of ELOVL2 was associated with unfavorable OS and DFS (Figures 5C, D), suggesting ELOVL2 was of importance in RPLS development and progression. In concordance with cohort-TCGA, the mRNA expression of ELOVL2 validated similar prognostic efficiency in both cohort-FD and cohort-GSE30929 (Supplementary Figures 5C, D).
Figure 5 Identification of lipid metabolism-associated targets. (A) The intersection of amplification, mutation, OS-related and DFS-related genes were identified as potential lipid metabolism-associated targets. (B) Identification of targets associated with APCs. Correlation between ELOVL2 expression and infiltration of macrophages, dendritic cells and B cells in RPLS. (C, D) Kaplan-Meier survival curves of OS and DFS for ELOVL2 expression in cohort-TCGA. (E) The comparison of ELOVL2 transcriptomic expression in LMS1 and LMS2. (F) The comparison of stromal score, immune score and ESTIMATE score calculated by ESTIMATE algorithm according to ELOVL2 expression. (G, H) Bubble diagram depicting the biological processes enriched by GO analysis and the signaling pathways enriched by KEGG analysis according to upregulated DEGs in ELOVL2high group, compared with ELOVL2low group. *P < 0.05, ***P < 0.001, and ****P < 0.0001.
Dedifferentiated liposarcoma (DDLPS) was often progressed from primary or recurrent well-differentiated liposarcomas (WDLPS), which constituted the most common pathological type of RPLS (36). Thus, we also investigated the heterogeneous expression of ELOVL2 in WDLPS and DDLPS. Interestingly, the expression of ELOVL2 was significantly higher in LMS2 than that in LMS1 (Figure 5E). Specifically, ELOVL2 exhibited a potential enrichment in DDLPS, as compared with WDLPS in both cohort-FD and cohort-GSE30929 (Supplementary Figures 5E, F), indicating the significant impacts of ELOVL2 in the dedifferentiation evolution of RPLS.
To further identify the association between immune status and ELOVL2 expression, a series of TIME profiles was also conducted. Accordingly, immune score, stromal score, ESTIMATE score (Figure 5F), infiltrating immune cells and TLS signature (Supplementary Figures 6A, B) were all calculated. As expected, we discovered that significant difference was found in the immune status between ELOVL2high and ELOVL2low subgroups, suggesting that RPLS patients with high expression of ELOVL2 might towards a status of immune desert. In addition, GO enrichment analysis indicated that high expression of ELOVL2 was involved in immune system process and cell surface reportor signaling. Similarly, KEGG enrichment analysis also validated cytokine-cytokine receptor interaction and cytokine signaling pathway (Figures 5G, H).
Accumulating evidence indicated that the increase of DNA methylation of ELOVL2 leaded to the decrease of its protein expression and polyunsaturated fatty acid synthesis, but the accumulation of short chain fatty acids, which is closely related to aging (37). Thus, we investigated the genome and epigenome landscape of ELOVL2 in cohort-TCGA (Supplementary Figure 6C), and 13 ELOVL2 related CpG sites were identified. However, only the cg20462795 exhibited significantly survival patterns (Supplementary Figure 6D), depicting that ELOVL2 associated epigenetic metabolic axis could be a novel therapeutic target in RPLS.
Transcription factor (TF) is one of the most common tool that involving in regulating gene expression (38). Thereafter, we further systematically screened and identified three ELOVL2 related TFs (TFDP1, TP73 and MYBL2) in TCGA-SARC cohort (Supplementary Figure 6E). Interestingly, high expression of them were significantly associated with decreased OS (Supplementary Figure 6F). Consistent with survival prediction occurring in our cohort, TFDP1 and MYBL2 exhibited similar survival patterns in both TCGA-SARC and GTEx databases (data not shown). Taken together, these results demonstrated that gene expression was determined by the synergistic regulation of both transcription and epigenetic factors.
ELOVL2 dominated lipid metabolism reprogramming and executive TIME affect prognosis in RPLS
Considering the dysregulation of immune status and metabolism remodeling involved in RPLS, we discovered that ELOVL2 played an essential role in tumor progression and even dedifferentiation, and served as the only one lipid metabolism related target biomarker. Moreover, PLCG1, ZC3HAV1L and NFAT5 were selected as immune hub genes through immune gene co-expression modules analysis and predicted great prognostic efficacy. Therefore, we re-analyzed the association between ELOVL2 and three hub genes.
Firstly, we validated these hub genes in cohort-FD and cohort-GSE30929 through survival analysis. Interestingly, the expression of PLCG1 was positive associated with ELOVL2 in all cohorts (Figure 6A). In addition, PLCG1 was significantly higher in LMS2 than that in LMS1, and exhibited similar survival patterns in all cohorts (Figure 6B). Next, the combined prognostic analysis of ELOVL2 and PLCG1 was also performed. Patient stratification based on these three groups presented that the Group I was associated with favorable prognosis, whereas the Group III was associated with dismal prognosis, and Group II was associated with medium prognosis (Supplementary Figure 7A). Subsequently, we developed another prognostic model for OS prediction based on ELOVL2 and PLCG1, which suggested a great prognostic efficacy (Figure 6C).
Figure 6 ELOVL2 dominated lipid metabolism reprogramming and executive TIME affect prognosis in RPLS. (A) Correlation among the transcriptomic expression of ELOVL2 and PLCG1 in three databases (cohort-TCGA, GSE30929 and FD). Spearman’s correlation coefficient was calculated. (B) Kaplan-Meier survival curves of OS for PLCG1 expression in cohort-FD. The comparison of PLCG1 transcriptomic expression in LMS1 and LMS2. (C) Kaplan-Meier survival curve of OS for combined riskscore of ELOVL2 and PLCG1 in cohort-TCGA. (D) Heatmap showing the transcriptomic expression of genes involving TIME according to ELOVL2 and PLCG1 expressions. The histogram on the right showing the results of differential expression analysis of each gene in indicated comparisons on the top of each column. (E) t-SNE plot showing of the overview of 11 cell clusters in the integrated single-cell transcriptomes of 54126 cells from four RPLS. Clusters are named as indicated cell subsets according to the specific gene expression patterns, color-coded according to ELOVL2 and PLCG1. (F) Representative mIHC images show the positivity of CD3, CD8, CD20, CD11b, CD68, and ELOVL2 in case 3 and case 1. Scale bar, 100 um. *P < 0.05, **P < 0.01, and ***P < 0.001.
In order to systematically comprehend the extent of relevancy between ELOVL2 and PLCG1 within TIME, we re-analyzed the expression profile in cohort-TCGA. Based on comparative analysis, we observed remarkable positive correlation between ELOVL2 and PLCG1. In addition, the infiltration of T cells, monocytics, dendritic cells, MDSC and M2-TAM were significantly lower with high expression of ELOVL2 and PLCG1. Moreover, the chemokine enrichment of CCL3, CCL4, CCL5 and CCL18, as well as PDCD1LG2 and HAVCR2 were significantly negative associated with high expression of ELOVL2 and PLCG1 (Figure 6D). These results demonstrated that highly expression of ELOVL2 and PLCG1 in TIME of RPLS was involved with a immune-exculded phenotype.
To fully characterize the specific cellular localization ELOVL2 and PLCG1 in RPLS, we first evaluated them in the human protein atlas and further validated them in situ single cell spatial phenotype analysis through single cell RNA-sequencing and immunohistochemistry (IHC) in cohort-FD. Interestingly, ELOVL2 was seemed enriching in fibroblasts (no images in database), while PLCG1 in T cells (Supplementary Figure 7B). At the single-cell level, ELOVL2 was expressed in tumor cells, cancer associated fibroblasts (CAFs) and smooth muscle cells (SMCs), while PLCG1 was presented in CD4+ and CD8+ T cells, as expected (Figure 6E; Supplementary Figure 7C). According to the analysis of bulk-RNA seq data, case 1/2 was defined as LMS2 and case 3/4 with LMS1, respectively. Next, mIHC was performed to validate previous findings in these four DDLPS. Consistent with RNA-seq data,
CD3+T cell, CD20+B cell, CD11b+DC cell and CD68+ macrophages presented a higher infiltration in case3 and case4, but less ELOVL2+ cells. On the contrary, less CD3+T cells, CD20+B cells and CD68+ macrophages, but more ELOVL2+ cells were found in case 1 and case 2 (Figure 6F; Supplementary Figure 7E). Moreover, the densities of CD3+ and cytotoxic CD8+ T cells in the area of tumor and invasive margin were quantified by IHC. We observed that the density of CD3+ and cytotoxic CD8+ T cells were significantly but negatively associated with ELOVL2+ cells (Supplementary Figure 7D). In concordance with previous findings, the mRNA expression level of CD3 and CD8 was negatively associated with ELOVL2 (Supplementary Figure 7F). Taken together, the dysregulation of lipid metabolism might remodel an immune desert of TIME further supported their crucial roles in the evolution of RPLS.
Discussion
RPLS is one of the most aggressive malignancies with a heterogeneous molecular profile, lipid metabolism dysregulation, limited medical efficacy and highly local recurrence rate (39). The combination of doxorubicin and ifosfamide is the first-line option in treating advanced RPLS, with limited clinical benefits and a median OS of just 8 to 14 months (40). Although immunotherapy has revolutionized oncology from the therapeutic point of view, its effectiveness in RPLS remains limited (41).
To the best of our knowledge, this is the first study aimed to systematically screen lipid metabolism related targets for RPLS. We conducted the profiles of somatic mutations and amplification in LMAGs, revealing potential targets that might be considered in RPLS. Considering the targets selected using the gene alteration profile might not be functionally significant, prognostic efficiency in OS and DFS and immune correlations were further evaluated to confirm the clinical relevance of the targets. ELOVL2 correlated with significantly prognosis and infiltration of APCs was identified. For instance, ELOVL2 was identified as an unique tissue-independent age-associated DNA methylation marker (42) and a specific superenhancer-associated gene implicated the LC-PUFA synthesis network as a critical metabolic dependency (43), and was associated with worsened patient survival in glioblastoma (44). In addition, ELOVL2 deficient indicated an increased infiltration of TH1/TH17 cells and a decrease of Treg (45). Thus, although basic and clinical investigation is further required, the potential of ELOVL2 to be successful targets was consolidated in RPLS.
Lipid metabolism and synthesis is based on the normal function of endoplasmic reticulum (ER). Previous results indicated that the accumulation of free fatty acids in the ER would eventually lead abnormal protein overloading and chronic ER stress (46). Meanwhile, the reduction of PUFA synthesis upon ELOVL2 decreasing can affect cellular metabolic homeostasis and mitochondrial energy metabolism (47). Unexpected, we found significantly activation of glucose metabolism pathway with ELOVL2 deficiency, such as gluconeogenesis, pentose phosphate pathway, pentose and glucuronate interconversions, fructose and mannose metabolism. This finds indicated that ELOVL2 deficiency may induce a switch in metabolism from the tri-carboxylic acid cycle to glycolysis, which eventually produces more reactive oxidative species (ROS) and causes oxidative stress. In addition, more mannose type O-glycan biosynthesis and lysine degradation was also found in the overexpression of ELOVL2. As the key regulator in PUFAs synthesis, more fatty acid elongation, but less arachidonic acid metabolism could be found in high group of ELOVL2, which was consistent with the recently published data (48). However, the relationship between abnormal expression phenotype of ELOVL2 and ER stress, mitochondrial dysfunction and cellular senescence still needs to be further studied.
Lipid metabolism reprogramming is a cretical marker of tumorigenesis and development. Previous study had indicated that a large number of lipid related genes underwent copy number amplification in the process of malignant development of RPLS (49). Therefore, it could be inferred that the amplification and high expression of lipid related gene-ELOVL2, was associated with malignant progression of RPLS with dismal prognosis. However, ELOVL2 displayed a heterogeneous role in the prognostic value in breast cancer (BC) (50) and glioma (44) (Supplementary Figure 8), respectively. In BC, ELOVL2 was hypermethylated and downregulated in the samples from tamoxifen resistance BC patients compared with those from tamoxifen-sensitive patients. Strikingly, in addition to having tumor suppressor activity, ELOVL2 was shown to recover tamoxifen sensitivity up to 70% in the MCF-7/tamoxifen resistance cells and in a xenograft mouse model (50). Furthermore, the depletion of ELOVL2 induced metastatic characteristics in BC cells via the SREBPs axis (51). In contrast, ELOVL2 depletion altered the phospholipid composition of the cell membrane by controlling fatty acid elongation, disrupting the structural characteristics of the cell membrane, and reducing EGFR signaling in glioblastoma cells (44). Taken together, we believed that abnormal lipid metabolism, including abnormal activation and inhibition, could both interfered with lipid metabolism reprogramming through different signal pathways in different types of cancers.
Since the tumor immune status is an another determinant of cancer associated efficacy, we further characterized the immune landscape in the different ELOVL2 subtypes. Interestingly, high ELOVL2 group is defined as immunological “cold” phenotype. In addition, the molecular signatures of this subtype is consistent with the immune status, indicating that patients with different immune signature may respond distinctly to immune therapy. To circumvent the poor immunogenicity of this subtype, targeted immune associated-biomarkers that reinvigorate the immune system by fascinating immune cell infiltrating may be a suitable option.
We demonstrated that RPLS patients could stratified into two LMSs with significant differences in molecular features and clinical prognosis. Patients with the LMS1 tumor had immune “hot” phenotype, whereas those with the LMS2 tumor had immune “cold” phenotype. Therefore, the lipid metabolism-associated molecular subtypes and risk model based on LMAGs were promising and complements the previous classification for RPLS.
With the prevalence of COVID-19 in the worldwide, the mRNA vaccine has highlighted its important strategic position, and greatly accelerating the development process of mRNA vaccine. Meanwhile, it also speeds up the research and development of mRNA cancer vaccine (52). The mRNA cancer vaccines represent promising novel method to treat malignancies with monotherapy or combination therapy (53). ELOVL2 is a potential target linking lipid metabolism to immune regulations for RPLS, specifically for patients with LMS2 tumors. However, given the function of ELOVL2 in lipid homeostasis and cellular homeostasis in normal cells, there is likely profound central tolerance against ELOVL2. On the contrast, ELOVL2 may be a target for small molecule inhibition for RPLS.
In summary, our study identified ELOVL2 as the potential effective target for RPLS, and patients with LMS2 tumor indicated an immune-excluded phenotype might benefit more from small molecule inhibition targerted ELOVL2.
Data availability statement
The transcriptome data presented in the study are deposited in the Sequence Read Archive (SRA) repository, accession number PRJNA987378.
Ethics statement
This study was reviewed and approved by the Institutional Review Board of Zhongshan Hospital (ID: B2022-586R), Shanghai, China. Informed consent was obtained from all study participants. All studies were performed in accordance with the Declaration of Helsinki.
Author contributions
Conception and design: ZW, LM, YZ, HT. Development of methodology: ZW, PT, PF, YH. Acquisition of data (acquired and managed patients, provided facilities, etc.): ZW, JW, YHZ, WL, YZ. Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): ZW, PT, TR, LM. Writing, review, and/or revision of the manuscript: ZW, LM. Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): ZW, LM, HT. Study supervision: LH, LM, YZ, HT. All authors read and approved the final manuscript.
Funding
This work was supported by grants from National Natural Science Foundation of China (81802302); Scientific Research Project of Shanghai Municipal Health Commission (20214Y0087, 20204Y0409); “Young Talents” Training Plan of Shanghai TCM-integrated Hospital (No. RCPY0063); Scientific Research Project of Shanghai TCM-integrated Hospital (No. 18-01-03); Scientific Research Project of Hongkou District Health Committee (No. 2003-02).
Acknowledgments
We wish to thank the Sinotech Genomics CO. LTD. Shanghai for their support of scRNA-seq library preparation. We thank Dr. Quanli Zou and his bioinformatics team for their reference codes. We thank Department of Pathology of Shanghai Skin Disease Hospital, Tongji University School of Medicine, for providing immunohistochemical staining and analysis platform.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1209396/full#supplementary-material
Glossary
References
1. Casadei L, de Faria F, Lopez-Aguiar A, Pollock RE, Grignol V. Targetable pathways in the treatment of retroperitoneal liposarcoma. Cancers (Basel) (2022) 14(6):1362. doi: 10.3390/cancers14061362
2. Hirata M, Asano N, Katayama K, Yoshida A, Tsuda Y, Sekimizu M, et al. Integrated exome and RNA sequencing of dedifferentiated liposarcoma. Nat Commun (2019) 10:5683. doi: 10.1038/s41467-019-13286-z
3. Zhu M, Shenasa E, Nielsen TO. Sarcomas: immune biomarker expression and checkpoint inhibitor trials. Cancer Treat Rev (2020) 91:102115. doi: 10.1016/j.ctrv.2020.102115
4. Delisle M, Gyorki D, Bonvalot S, Nessim C. Landmark series: a review of landmark studies in the treatment of primary localized retroperitoneal sarcoma. Ann Surg Oncol (2022) 29:7297–311. doi: 10.1245/s10434-022-12517-w
5. Gronchi A, Ferrari S, Quagliuolo V, Broto JM, Pousa AL, Grignani G, et al. Histotype-tailored neoadjuvant chemotherapy versus standard chemotherapy in patients with high-risk soft-tissue sarcomas (ISG-STS 1001): an international, open-label, randomised, controlled, phase 3, multicentre trial. Lancet Oncol (2017) 18:812–22. doi: 10.1016/S1470-2045(17)30334-0
6. Iwamoto H, Abe M, Yang Y, Cui D, Seki T, Nakamura M, et al. Cancer lipid metabolism confers antiangiogenic drug resistance. Cell Metab (2018) 28:104–17. doi: 10.1016/j.cmet.2018.05.005
7. Yu W, Lei Q, Yang L, Qin G, Liu S, Wang D, et al. Contradictory roles of lipid metabolism in immune response within the tumor microenvironment. J Hematol Oncol (2021) 14:187. doi: 10.1186/s13045-021-01200-4
8. Grunt TW. Interacting cancer machineries: cell signaling, lipid metabolism, and epigenetics. Trends Endocrinol Metab (2018) 29:86–98. doi: 10.1016/j.tem.2017.11.003
9. Martin-Perez M, Urdiroz-Urricelqui U, Bigas C, Benitah SA. The role of lipids in cancer progression and metastasis. Cell Metab (2022) 34:1675–99. doi: 10.1016/j.cmet.2022.09.023
10. Yin X, Xu R, Song J, Ruze R, Chen Y, Wang C, et al. Lipid metabolism in pancreatic cancer: emerging roles and potential targets. Cancer Commun (Lond) (2022) 42:1234–56. doi: 10.1002/cac2.12360
11. Gu X, Jiang C, Zhao J, Qiao Q, Wu M, Cai B. Identification of lipid metabolism-associated genes as prognostic biomarkers based on the immune microenvironment in hepatocellular carcinoma. Front Cell Dev Biol (2022) 10:883059. doi: 10.3389/fcell.2022.883059
12. Gong M, Liu X, Yang W, Song H, Zhao X, Ai X, et al. Identification of a lipid metabolism-associated gene signature predicting survival in breast cancer. Int J Gen Med (2021) 14:9503–13. doi: 10.2147/IJGM.S343426
13. Zhao J, Li G, Zhao G, Wang W, Shen Z, Yang Y, et al. Prognostic signature of lipid metabolism associated LncRNAs predict prognosis and treatment of lung adenocarcinoma. Front Oncol (2022) 12:986367. doi: 10.3389/fonc.2022.986367
14. Ye Z, Zou S, Niu Z, Xu Z, Hu Y. A novel risk model based on lipid metabolism-associated genes predicts prognosis and indicates immune microenvironment in breast cancer. Front Cell Dev Biol (2021) 9:691676. doi: 10.3389/fcell.2021.691676
15. Tawbi HA, Burgess M, Bolejack V, Van Tine BA, Schuetze SM, Hu J, et al. Pembrolizumab in advanced soft-tissue sarcoma and bone sarcoma (SARC028): a multicentre, two-cohort, single-arm, open-label, phase 2 trial. Lancet Oncol (2017) 18:1493–501. doi: 10.1016/S1470-2045(17)30624-1
16. Bleve A, Durante B, Sica A, Consonni FM. Lipid metabolism and cancer immunotherapy: immunosuppressive myeloid cells at the crossroad. Int J Mol Sci (2020) 21(16):5845. doi: 10.3390/ijms21165845
17. Tan PY, Chang CW, Duan K, Poidinger M, Ng KL, Chong YS, et al. E2F1 orchestrates transcriptomics and oxidative metabolism in wharton's jelly-derived mesenchymal stem cells from growth-restricted infants. PLoS One (2016) 11:e163035. doi: 10.1371/journal.pone.0163035
18. Liang Y, Huang X, Zhang Z, Deng K, An S, Gao X, et al. Spirulina supplementation improves lipid metabolism and autophagic activities in the liver and muscle of Hu lambs fed a high-energy diet. Arch Anim Nutr (2020) 74:476–95. doi: 10.1080/1745039X.2020.1820806
19. Blanco-Luquin I, Acha B, Urdanoz-Casado A, Sanchez-Ruiz DGJ, Vicuna-Urriza J, Roldan M, et al. Early epigenetic changes of alzheimer's disease in the human hippocampus. Epigenetics (2020) 15:1083–92. doi: 10.1080/15592294.2020.1748917
20. Osada H, Toda E, Homma K, Guzman NA, Nagai N, Ogawa M, et al. ADIPOR1 deficiency-induced suppression of retinal ELOVL2 and docosahexaenoic acid levels during photoreceptor degeneration and visual loss. Cell Death Dis (2021) 12:458. doi: 10.1038/s41419-021-03741-5
21. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. Omics (2012) 16:284–87. doi: 10.1089/omi.2011.0118
22. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47:D607–13. doi: 10.1093/nar/gky1131
23. 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:401–04. doi: 10.1158/2159-8290.CD-12-0095
24. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res (2020) 48:W509–14. doi: 10.1093/nar/gkaa407
25. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol (2019) 37:773–82. doi: 10.1038/s41587-019-0114-2
26. Petitprez F, de Reynies A, Keung EZ, Chen TW, Sun CM, Calderaro J, et al. B cells are associated with survival and immunotherapy response in sarcoma. Nature (2020) 577:556–60. doi: 10.1038/s41586-019-1906-8
27. Lin Z, Huang L, Li S, Gu J, Cui X, Zhou Y. Pan-cancer analysis of genomic properties and clinical outcome associated with tumor tertiary lymphoid structure. Sci Rep (2020) 10:21530. doi: 10.1038/s41598-020-78560-3
28. Haase F, Gloss BS, Tam P, Gold WA. WGCNA identifies translational and proteasome-ubiquitin dysfunction in Rett syndrome. Int J Mol Sci (2021) 22(18):9954. doi: 10.3390/ijms22189954
29. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol (2018) 36:411–20. doi: 10.1038/nbt.4096
30. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou YT, et al. The immune landscape of cancer. Immunity (2018) 48:812–30. doi: 10.1016/j.immuni.2018.03.023
31. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet (2019) 51:202–06. doi: 10.1038/s41588-018-0312-8
32. Kroemer G, Galassi C, Zitvogel L, Galluzzi L. Immunogenic cell stress and death. Nat Immunol (2022) 23:487–500. doi: 10.1038/s41590-022-01132-2
33. Wu C, Qin C, Long W, Wang X, Xiao K, Liu Q. Tumor antigens and immune subtypes of glioblastoma: the fundamentals of mRNA vaccine and individualized immunotherapy development. J Big Data (2022) 9:92. doi: 10.1186/s40537-022-00643-x
34. Bagaria SP, Gabriel E, Mann GN. Multiply recurrent retroperitoneal liposarcoma. J Surg Oncol (2018) 117:62–8. doi: 10.1002/jso.24929
35. Zhou Q, Yan X, Zhu H, Xin Z, Zhao J, Shen W, et al. Identification of three tumor antigens and immune subtypes for mRNA vaccine development in diffuse glioma. Theranostics (2021) 11:9775–90. doi: 10.7150/thno.61677
36. Lee A, Thway K, Huang PH, Jones RL. Clinical and molecular spectrum of liposarcoma. J Clin Oncol (2018) 36:151–59. doi: 10.1200/JCO.2017.74.9598
37. Chen D, Chao DL, Rocha L, Kolar M, Nguyen HV, Krawczyk M, et al. The lipid elongation enzyme ELOVL2 is a molecular regulator of aging in the retina. Aging Cell (2020) 19:e13100. doi: 10.1111/acel.13100
38. Donovan BT, Larson DR. Regulating gene expression through control of transcription factor multivalent interactions. Mol Cell (2022) 82:1974–75. doi: 10.1016/j.molcel.2022.05.015
39. Song Z, Wang S, Lu L, Xu J, Zhou Q, Lu W, et al. Lipidomics and transcriptomics differ liposarcoma differentiation characteristics that can be altered by pentose phosphate pathway intervention. Metabolites (2022) 12(12):1227. doi: 10.3390/metabo12121227
40. Judson I, Verweij J, Gelderblom H, Hartmann JT, Schoffski P, Blay JY, et al. Doxorubicin alone versus intensified doxorubicin plus ifosfamide for first-line treatment of advanced or metastatic soft-tissue sarcoma: a randomised controlled phase 3 trial. Lancet Oncol (2014) 15:415–23. doi: 10.1016/S1470-2045(14)70063-4
41. Song Z, Lu L, Gao Z, Zhou Q, Wang Z, Sun L, et al. Immunotherapy for liposarcoma: emerging opportunities and challenges. Future Oncol (2022) 18:3449–61. doi: 10.2217/fon-2021-1549
42. Li X, Wang J, Wang L, Gao Y, Feng G, Li G, et al. Lipid metabolism dysfunction induced by age-dependent DNA methylation accelerates aging. Signal Transduct Target Ther (2022) 7:162. doi: 10.1038/s41392-022-00964-6
43. Pan G, Cavalli M, Carlsson B, Skrtic S, Kumar C, Wadelius C. rs953413 regulates polyunsaturated fatty acid metabolism by modulating ELOVL2 expression. Iscience (2020) 23:100808. doi: 10.1016/j.isci.2019.100808
44. Gimple RC, Kidwell RL, Kim L, Sun T, Gromovsky AD, Wu Q, et al. Glioma stem cell-specific superenhancer promotes polyunsaturated fatty-acid synthesis to support EGFR signaling. Cancer Discovery (2019) 9:1248–67. doi: 10.1158/2159-8290.CD-19-0061
45. Chiurchiu V, Leuti A, Dalli J, Jacobsson A, Battistini L, Maccarrone M, et al. Proresolving lipid mediators resolvin D1, resolvin D2, and maresin 1 are critical in modulating T cell responses. Sci Transl Med (2016) 8:111r–353r. doi: 10.1126/scitranslmed.aaf7483
46. Hubel E, Fishman S, Holopainen M, Kakela R, Shaffer O, Houri I, et al. Repetitive amiodarone administration causes liver damage via adipose tissue ER stress-dependent lipolysis, leading to hepatotoxic free fatty acid accumulation. Am J Physiol Gastrointest Liver Physiol (2021) 321:G298–307. doi: 10.1152/ajpgi.00458.2020
47. Abdolmaleky HM, Sheng Y, Zhou JR. Bioactive nutraceuticals oligo-lactic acid and fermented soy extract alleviate cognitive decline in mice in part via anti-neuroinflammation and modulation of gut microbiota. Front Nutr (2023) 10:1116278. doi: 10.3389/fnut.2023.1116278
48. Hammouda S, Ghzaiel I, Khamlaoui W, Hammami S, Mhenni SY, Samet S, et al. Genetic variants in FADS1 and ELOVL2 increase level of arachidonic acid and the risk of alzheimer's disease in the Tunisian population. Prostaglandins Leukot Essent Fatty Acids (2020) 160:102159. doi: 10.1016/j.plefa.2020.102159
49. Hao Y, Li D, Xu Y, Ouyang J, Wang Y, Zhang Y, et al. Investigation of lipid metabolism dysregulation and the effects on immune microenvironments in pan-cancer using multiple omics data. BMC Bioinf (2019) 20:195. doi: 10.1186/s12859-019-2734-4
50. Jeong D, Ham J, Kim HW, Kim H, Ji HW, Yun SH, et al. ELOVL2: a novel tumor suppressor attenuating tamoxifen resistance in breast cancer. Am J Cancer Res (2021) 11:2568–89.
51. Kang YP, Yoon JH, Long NP, Koo GB, Noh HJ, Oh SJ, et al. Spheroid-induced epithelial-mesenchymal transition provokes global alterations of breast cancer lipidome: a multi-layered omics analysis. Front Oncol (2019) 9:145. doi: 10.3389/fonc.2019.00145
52. Lorentzen CL, Haanen JB, Met O, Svane IM. Clinical advances and ongoing trials on mRNA vaccines for cancer treatment. Lancet Oncol (2022) 23:e450–58. doi: 10.1016/S1470-2045(22)00372-2
Keywords: ELOVL2, retroperitoneal liposarcoma, TCGA, lipid metabolism, immune landscape
Citation: Wang Z, Tao P, Fan P, Wang J, Rong T, Hou Y, Zhou Y, Lu W, Hong L, Ma L, Zhang Y and Tong H (2023) Insight of a lipid metabolism prognostic model to identify immune landscape and potential target for retroperitoneal liposarcoma. Front. Immunol. 14:1209396. doi: 10.3389/fimmu.2023.1209396
Received: 20 April 2023; Accepted: 20 June 2023;
Published: 06 July 2023.
Edited by:
Thomas Dubensky, Tempest Therapeutics, Inc., United StatesReviewed by:
Chan Whiting, Actym Therapeutics, Inc., United StatesMarion Curtis, Mayo Clinic Arizona, United States
Copyright © 2023 Wang, Tao, Fan, Wang, Rong, Hou, Zhou, Lu, Hong, Ma, Zhang and Tong. 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: Liang Hong, ZG9jdG9yX2hsQDEyNi5jb20=; Lijie Ma, bWFsaWppZV8wOTE1QHNqdHUuZWR1LmNu; Yong Zhang, emhhbmcueW9uZzJAenMtaG9zcGl0YWwuc2guY24=; Hanxing Tong, dG9uZy5oYW54aW5nQHpzLWhvc3BpdGFsLnNoLmNu
†These authors have contributed equally to this work