Development and validation of cuproptosis-related lncRNAs associated with pancreatic cancer immune microenvironment based on single-cell

Background: Cuproptosis, a novel mode of cell death associated with the tricarboxylic acid (TCA) cycle, is relevant to the development of cancer. However, the impact of single-cell-based Cuproptosis-associated lncRNAs on the Tumor immune microenvironment (TIME) of Pancreatic adenocarcinoma (PAAD) and its potential value for individualized immunotherapy has not been clarified.

Methods: 14 immune-related CRGs were screened by exploring the interaction between differentially expressed Immune-Related Genes (IRGs) and Cuproptosis-Related Genes (CRGs) in PAAD. Next, the expression amount and expression distribution of CRGs in single-cell samples were analyzed by focusing on 7-CRGs with significant expressions. On the one hand, MAP2K2, SOD1, and VEGFA, which were significantly differentially expressed between PAAD sites and normal tissues adjacent to them, were subjected to immunohistochemical validation and immune landscape analysis. On the other hand, from these 7-CRGs, prognostic signatures of lncRNAs were established by co-expression and LASSO-COX regression analysis, and their prognostic value and immune relevance were assessed. In addition, this study not only validated the hub CRGs and the lncRNAs constituting the signature in a PAAD animal model treated with immunotherapy-based combination therapy using immunohistochemistry and qRT-PCR but also explored the potential value of the combination of targeted, chemotherapy and immunotherapy.

Results: Based on the screening of 7-CRGs significantly expressed in a PAAD single-cell cohort and their co-expressed Cuproptosis-Related lncRNAs (CRIs), this study constructed a prognostic signature of 4-CRIs named CIR-score. A Nomogram integrating the CIR-score and clinical risk factors was constructed on this basis to predict the individualized survival of patients. Moreover, high and low-risk groups classified according to the median of signatures exhibited significant differences in clinical prognosis, immune landscape, bioenrichment, tumor burden, and drug sensitivity. And the immunohistochemical and qRT-PCR results of different mouse PAAD treatment strategies were consistent with the trend of inter-group variability in drug sensitivity of hub CRGs and CIR-score. The combination of immunotherapy, targeted therapy, and chemotherapy exhibited a better tumor suppression effect.

Conclusion: CIR-score, as a Cuproptosis-related TIME-specific prognostic signature based on PAAD single cells, not only predicts the prognosis and immune landscape of PAAD patients but also provides a new strategy for individualized immunotherapy-based combination therapy.

1 Introduction

Worldwide, PAAD ranks seventh in cancer-related deaths (1). As a highly fatal disease, PAAD mortality is expected to continue to rise in the coming decades, with more than 800,000 deaths expected by 2040 (2). A major feature of human PAAD is the lack of effective anti-tumor immunity. Multiple immunosuppressive and evasion mechanisms, established from the early steps of tumor transformation, effectively protect PAAD from immunity (3). Suppressive TIME, consisting of multiple immune cell populations and stromal components in tumor islets, is not only a fertile ground for PAAD development, progression, and metastasis but also a great obstacle to immunotherapy (4). Single-cell sequencing, as a technique for analysis at the individual cellular level, overcomes the problem of cellular heterogeneity (5). It not only better reflects how the different cell populations that makeup TIME regulate different tumor states from a more microscopic perspective, but also has important implications for unraveling the mechanisms of cancer progression and immunotherapeutic response (6, 7).

Cuproptosis refers to a novel form of cell death in the TCA cycle through the binding of copper to lipoacylase leading to subsequent protein aggregation, proteotoxic stress, and ultimately cell death (8, 9).CRGs and lipoacylated proteins abundance are highly correlated in human tumors, and cell lines with high levels of lipoacylated proteins are sensitive to copper-induced cell death (10). Considering the increased demand for copper by cancer cells, this suggests the potential value of copper ion carrier therapy in tumors with this metabolic profile. Previous studies have shown that serum copper ion concentrations are significantly elevated in patients with PAAD (11). Furthermore, in a mouse model of PAAD, not only did tumor growth rates increase significantly in those chronically exposed to elevated copper levels, but the copper chelator TTM delayed angiogenesis in precancerous lesions and the growth of advanced tumors to some extent (12).

lncRNAs are non-coding RNAs greater than 200 nucleotides in length and transcribed by RNA polymerase II. lncRNAs’ aberrant expression of oncogenic or tumor suppressive effects, high tissue and disease specificity, and ability to modulate TIME make them potential candidates for biomarker and cancer vaccine development (13, 14). lncRNAs not only participate in RNA regulatory mechanisms and control the expression of their downstream target genes, but can also mediate a range of multiple cellular processes through chromatin reprogramming, cis- or trans-regulation of neighboring genes, and post-transcriptional regulation of mRNA processing (15). In the current PAAD studies, proliferation and metastasis have been clearly shown to be affected by the negative regulation of the miR-34 transcriptional pathway by lncRNA HOTAIR (16), the negative regulation of HIF-1α pathway by the reduced expression of lncRNA ENST00000480739 (17), the miR-448 sponge pathway participated by lncRNA PVT1 (18), the inhibition of the oncogene E-cadherin pathway by lncRNA MALAT1 (19), and the regulation of the HOX gene pathway by lncRNA HOTTIP (20). Apart from affecting PAAD proliferation and metastasis, lncRNAs can also regulate TIME via macrophages (21), PD-1 on the surface of T cells (22), etc. In the field of PAAD, Chen H et al. have conducted some studies on the prognostic value of Cuproptosis-associated lncRNAs (23), but there is still a gap in the exploration of Cuproptosis-associated lncRNAs and their characteristic TIME based on PAAD single-cell data.

Given the limited evidence for effective early detection and therapeutic interventions for PAAD, a unified predictive approach needs to be incorporated into strategies for early detection and guiding treatment. In this study, CRGs with immune-related potential prognostic value were screened using differential expression levels and the degree of interaction with IRGs as indicators. The immune landscape of high- and low-expression cohorts of these genes were analyzed on the one hand, and the differences in expression and distribution of these genes in different cell clusters were explored individually at the single-cell level on the other. Subsequently, a prognostic signature of 4-CRIs named CIR-score was constructed based on the 9 most prognostic CRIs screened. This new Cuproptosis-related prognostic signature not merely reflects the tumor load, immune landscape, and immune checkpoint sensitivity of PAAD but also well predicts patient survival and drug sensitivity. Finally, a cohort of four PAAD mouse models with different treatment regimens was used to validate hub lncRNAs of signature and differentially expressed CRGs with prognostic value. The promising tumor suppressive effect demonstrated by immunotherapy combined with targeted therapy and chemotherapy also provides innovative therapeutic options for PAAD.

2 Materials and methods

2.1 Data acquisition and processing

The overall flow of this study is shown in the figure (Figure 1). A total of 53 CRGs were summarized from previous studies (10, 24). This method of summarizing CRGs has been used by Cai Z et al. in the field of renal clear cell carcinoma (25). A total of 2350 IRGs were screened from the Gene Set Enrichment Analysis database (GSEA, https://www.gsea-msigdb.org/gsea/index.jsp) based on 10 tumor immune-related datasets. A total of 2483 IRGs were downloaded from the Immport database (Updated. July 2020, https://www.immport.org). 178 PAAD samples and 4 normal tissue samples were obtained from The Cancer Genome Atlas Program (TCGA, https://portal.gdc.cancer.gov/). The samples were collated for quantitative gene expression data and clinical data. Transcriptomic data of 167 normal pancreatic tissues were obtained from the Genotype-Tissue Expression Project database (GTEx, https://gtexportal.org/). The cohorts of TCGA-PAAD and GTEx were matched. Single-cell data were obtained from GSE212966 and GSE205049 in the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/).

FIGURE 1

Figure 1 The flowchart of our study process.

2.2 Screening and interaction analysis of CRGs and IRGs

The 877 IRGs obtained by taking the intersection of two databases, GSEA and Immport, were subjected to analysis of variation and Univariate COX in the TCGA-PAAD cohort (p=0.05). 71 differentially expressed IRGs with prognostic values were obtained. The String database (https://cn.string-db.org) and Cytoscape software (version 3.9.0) were used to explore the interaction between the two groups of 53 CRGs and 71 IRGs, and only genes with intergroup interactions were retained. Immunohistochemical results of differentially expressed CRGs were obtained through the Human Protein Atlas database (HPA, https://www.proteinatlas.org/).

2.3 Immune landscape of differentially expressed CRGs

The TIME evaluation (StromalScore, ImmuneScore, and ESTIMATEScore) of the three immune-related CRGs with significant differential expression was evaluated using the “ESTIMATE” package. “CIBERSORT” was used to obtain the results of immune cell infiltration in different groups with high and low expression of CRGs. The correlation between different immune cells and the expression of different CRGs was obtained by Spearman correlation analysis, and the correlation coefficients and p-values were calculated. The “ggplot2” package was used to draw a Lollipop plot of the correlation between genes and immune cells. The “limma” package and the “corrplot” package were used to plot the correlation heat map between different CRGs and immune checkpoints.

2.4 Expression and distribution of CRGs in single cells

Single-cell RNA sequencing (scRNA-seq) count matrices for all 6 PAAD samples were downloaded from GSE212966. These 6 matrices were combined into 1 Seurat object using the CreateSeuratObject function (“Seurat” R package, version 4.3.0). Single cells with more than 2500 or less than 200 genes detected per cell or with a percentage of mitochondria-derived UMI counts greater than 5% were considered low-quality cells and filtered out for further analysis. Single-cell data from 6 PAAD samples were normalized on the basis of the first 2000 high-variable genes (HVG). Subsequently, similar cells were identified using the FindNeighbors function in the “Seurat” package. The FindClusters function (resolution=0.4) was used to identify the major cell clusters. Based on this, UMAP non-linear dimensionality reduction is performed on these cell clusters. The results of the dimensionality reduction were visualized using DimPlot. Cell clusters were identified using the “SingleR” package, and annotation of cell clusters was done based on marker genes obtained from previous studies. Finally, the accuracy of the annotation is verified using the plotScoreHeatmap function. Using the FeaturePlot function, the expression and distribution of the target CRGs are marked on the UMAP map. In addition, GSM6567169 and GSM6567170 with incomplete data and GSM6567167 with too low quality were deleted, and the same analysis as above was done for 3 single cell samples (GSM6567165, GSM6567166, and GSM6567171) with normal tissues adjacent to PAAD from GSE212966 as controls. Furthermore, single-cell sequencing data from the 9-case PAAD samples and their neighboring 9-case normal tissue samples of GSE205049 were included in the study. They were combined into one Seurat object for cancer tissues and one Seurat object for normal tissues using the same data analysis methods described above. The GSE205049 dataset was used as a complement and support to the results of the analysis of the top CRGs based on the samples from GSE212966.

2.5 Identification of Cuproptosis-associated lncRNAs signature with prognostic significance in PAAD

7 CRGs that were both co-expressed with IRGs and significantly expressed in PAAD single-cell data were screened. Wilcoxon test was used to obtain lncRNAs with co-expression relationships with these 7 CRGs (correlation coefficient >0.4, p<0.001). “ggalluvial” and “ggplot2” packages were used to draw a correlation Sankey diagram. Next, to further filtrate prognosis-related lncRNAs, we collected 9 CRIs by univariate Cox regression analysis (p< 0.05). Least absolute shrinkage and selection operator (LASSO) regression analysis and lambda spectra were used to screen prognosis-associated lncRNAs to prevent overfitting when constructing prognostic risk models. The area under the ROC curve (AUC) was calculated for these models. When the AUC value reached the maximum value, the model was indicated as the best candidate. Based on the 4 CRIs identified, a new PAAD risk scoring model was identified. This model was named The Cuproptosis-related Immune Risk Score (CIR-score) in this study, and the risk score for each PAAD patient can be derived using the following equation (β: coefficients, Exp: gene expression level):.

$$CIRscore=\Sigma ({\beta }_i\times Ex{p}_i)$$

In addition, to demonstrate the superiority of the signature in distinguishing between high and low-risk groups of patients, the “scatterplot3d” package was used to perform a principal component analysis (PCA) to visualize the distribution pattern of the CIR-score in the PAAD sample.

2.6 Independent prognostic and clinical correlation analysis

The CIR-score was analyzed by univariate and multifactorial Cox regression. The 1-, 3-, and 5-year ROC curves were plotted using the “timeROC” package. The ROC curves of CIR-score and clinical indicators were obtained in the same way as above. The Concordance Index (C-index) was calculated for the CIR-score and clinical indicators, respectively, to evaluate the predictive power of the signature. The prognosis of individual PAAD patients was evaluated by plotting the CIR-score-related Nomogram with the “regplot” package. The calibration curves of the Nomogram were plotted using the Calibrate function. Whereafter, the TCGA-PAAD cohort were randomly grouped into a train set and a test set. The randomness of the grouping was demonstrated by performing a round-robin for clinical traits and obtaining p values (p>0.05) for the differences in clinical traits between the train and test sets. The “pheatmap” package was used to map the differences in CIR-score expression patterns between the TCGA-PAAD cohort, the training set, and the test set for the high and low-risk groups. “survival” and “survminer” packages were used to compare the survival differences between high and low-risk groups. Overall survival (OS) and progression-free survival (PFS) of PAAD patients were analyzed separately using the Kaplan-Meier (KM) method. Besides, clinical traits were grouped and cycled for each group to compare survival differences between high and low-risk groups with different clinical characteristics.

2.7 Immune microenvironment landscape analysis

The “GSEABase” package was used to perform a single sample Gene Set Enrichment Analysis (ssGSEA) of immune-related functions in the high- and low-risk groups. The results of the “CIBERSORT” package were used to compare the differences in the infiltration of all 22 immune cell types in the high and low-risk groups. Spearman correlation analysis was used to calculate correlation coefficients and p-values by the immune cell cycle and to plot the correlation between the 4 immune cells with the most significant differences and the CIR-score. Considering that the TISIDB database (http://cis.hku.hk/TISIDB/) is an integrated repository portal for tumor-immune system interactions. To explore the correlation between CIR-score and immune checkpoints, this study summarized the specific immune checkpoint genes of these four immune cells from the TISIDB database. The “corrplot” package was used to visualize the correlations separately.

2.8 Functional analysis and tumor mutation burden analysis

The “limma” package was used to screen the differentially expressed genes (DEGs) in high and low-risk groups (∣log2FC∣>1, FDR<0.05). GO and KEGG analyses of DEGs were performed by enrichGO function and enrichKEGG function. Among them, GO analysis included 3 domains of biological process (BP), cellular component (CC), and molecular function (MF). Next, tumor mutation burden (TMB) values were obtained for each sample based on the mutation data downloaded and collated from the TCGA-PAAD cohort. The TMBs of the high- and low-risk groups were compared, and the differences were visualized as violin plots using the ggviolin function. The samples were divided into high and low TMB groups according to the median, and the survdiff function and ggsurvplot function were used to analyze the difference in OS between the high and low mutation load groups. On this basis, survival curves for TMB combined with patient risk were plotted in the same way to further explore the correlation between CIR-score and TMB. Moreover, considering Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) is a database that allows estimating multiple published transcriptomic biomarkers based on tumor pre-treatment expression profiles to predict patient response. Its scoring of PAAD samples was used to analyze differences in potential immune escape potential between high and low-risk groups.

2.9 Drug sensitivity analysis

Immunotherapy sensitivity scores for the TCGA-PAAD cohort were downloaded from The Cancer Immunome Atlas (TCIA, https://www.tcia.at/home) database. The TCIA database provides comprehensive immune next-generation sequencing data (NGS) from TCGA and other data sources for 20 solid cancer genomic analysis results. The “ggpubr” package was used to compare the immunotherapy sensitivity of different risk groups by combining the expression of CRGs, hub lncRNAs, and CIR-score size for each sample. The “pRophetic” package was used to calculate the difference in Half Maximal Inhibitory Concentration (IC50) of the different drugs between the different groups, the lower the IC50, the more sensitive the patient was to the drug.

2.10 Validation of hub CRGs and lncRNAs in the immune, targeted combination chemotherapy mouse PAAD model

2.10.1 Mice and drugs

6-8 week female mice (C57BL/6). Mouse pancreatic cancer PANC-02 cells. Nanoalbumin paclitaxel (Jiangsu Hengrui Pharmaceutical Co., Ltd., Specification: 100 mg, Lot No.: 220310AF), Gemcitabine, Recombinant Human Vascular Endothelial Inhibitor Injection (Shandong Xiangsheng Biopharmaceutical Co., Ltd., Specification: 15 mg/3 ml, Lot No.: 202109053), PD-1 monoclonal antibody (Jiangsu Hengrui Pharmaceutical Co., Ltd., Specification: 200 mg, Lot No. 202009043A). The above drugs were dissolved and diluted with saline to the desired concentration.

2.10.2 PAAD model preparation

The cultured mouse pancreatic cancer cells PANC-02 cell suspension was collected at a concentration of 5×10⁷ cells/ml and inoculated subcutaneously in the right axilla of mice at 0.1 ml each. The mice transplanted tumors were measured with vernier calipers to measure the diameter of the transplanted tumors, and the animals were randomly grouped when the tumors grew to 100 mm³. At the same time, each group of mice started to administer the drug, with 6 mice in each group.

2.10.3 Grouping and drug administration regimen

The experimental groupings in this study were a, b, c and d. Group-a was the model control group; group-b was the chemotherapy group; group-c was the Recombinant human endostatin + immunotherapy group (PD-1 Monoclonal Antibody group); and group-d was the combination of Recombinant human endostatin + immunotherapy (PD-1 Monoclonal Antibody) + chemotherapy group. The model control mice in group-a were administered with saline intravenous injection, 3 times on days 1, 4, and 7 after grouping. The chemotherapy mice in group-b were treated with Nanoalbumin Paclitaxel injection (20 mg/kg) and Gemcitabine (50 mg/kg) for 3 times on days 1, 4, and 7 after grouping. The Recombinant human endostatin + immunotherapy mice in group-c were treated with Recombinant human endostatin (20 mg/kg) for 15 consecutive days and PD-1 Monoclonal Antibody (200 μg/20g) for 3 times on days 1, 4, and 7 after grouping. The Recombinant human endostatin + immunotherapy + chemotherapy mice in group-d were treated with a combination of Recombinant human endostatin (20 mg/kg) for 15 consecutive days, PD-1 Monoclonal Antibody (200 μg/20g) for 3 times on days 1, 4, and 7 after grouping and Nanoalbumin Paclitaxel (20 mg/kg) + Gemcitabine (50 mg/kg) for 3 times on days 1, 4, and 7 after grouping. More detailed grouping and dosing regimen are shown in a table (Supplementary Table S1).

2.10.4 In vivo imaging

The tumor marker probe HGC675 (BioActs, LOT: AQN00037) was injected into the tail vein at 0W, 1W, and 2W of administration, respectively, 200 μl/each, and fluorescence in vivo imaging was performed at the small animal in vivo 3D imaging system IVIS^® Spectrum (PerkinElmer) 4 h after injection. Excitation wavelength: 635 nm, emission wavelength: 670 nm.

2.10.5 Immunohistochemical detection of hub CRGs

After the tumor tissues of each group of mice were rinsed, fixed, and dehydrated, the tissues were paraffin-embedded and the wax blocks were fixed on a microtome for continuous sectioning (4 μm thick). Routine immunohistochemistry was performed to detect the expression of 3 hub CRGs in PAAD tumor tissues under 4 different dosing regimens. Antibody information: rabbit anti SOD1 (proteintech, 10269-1-AP) at a dilution ratio of 1:500. rabbit anti VEGFA (proteintech, 19003-1-AP) at a dilution of 1:100. Dilution ratio 1:1000 for mouse anti MEK2 (proteintech, 67410-1-Ig). The kit used MaxVision-mouse/rabbit (Fuzhou Maishin Biotechnology Co., KIT-5010). DAB was used for color development, hematoxylin contrast staining, and observation under light microscopy. Image-Pro Plus (version 6.0) was used to analyze the immunohistochemical results.

2.10.6 Tissue detection of hub lncRNAs

Total RNA was extracted from tissues using TRIzol reagent (Invitgen, MA, USA). extracted RNA was subjected to quantitative polymerase chain reaction using the One Step TB Green™ PrimeScript™ RT-PCR Kit II (SYBR Green). qRT-PCR was performed with SYBR Green PCR Master Mix (Vazyme) in a fluorescent quantitative PCR cycler (ABI Step one plus Real time-PCR system, USA). GAPDH was used as an internal control, and the results of each sample were normalized to GAPDH expression. The sequences of the primers are shown in Supplementary Table S2.

2.11 Statistical analysis

All data analyses were performed using R software (version 4.2.1, http://www.R-project.org) and GraphPad Prism software (version 8.0.2). Where not explicitly mentioned, thresholds for statistical significance were considered as P< 0.05, with P values being bilateral.

3 Results

3.1 Development of immune-related CRGs

First, a total of 4833 IRGs were summarized from the Immport database and the 10 GSEA tumor immune-related subsets. 310 of them with significant differences constituted the IRGs cohort by comparing the mRNA expression levels between tumor tissues in the TCGA-PAAD cohort and adjacent normal tissues in the GTEx cohort (Figure 2A, Supplementary Table S3). Then, 53 CRGs were summarized from the study of Tsvetkov P and Ge EJ, constituting the CRGs cohort (Supplementary Table S4). The within-group relationships between the IRGs cohort and the CRGs cohort were eliminated, and only the between-group interaction relationships were retained. The results showed that 14 CRGs interacted significantly with 35 IRGs (Figure 2B, Supplementary Table S5). Among these 14 immune-related CRGs, MAP2K2, SOD1, and VEGFA were most significantly differentially expressed between the TCGA-PAAD cohort and the GTEx cohort (Figure 2C-E). Next, to analyze the reliability of the differences, immunohistochemical results from the HPA database were selected for validation (Figure 2F). Among them, the immunohistochemical results of MAP2K2 and VEGFA were consistent with the mRNA differential expression results. MAP2K2 staining in normal pancreatic tissues (Staining: High, Intensity: Strong, Quantity:>75%) was significantly higher than that in PAAD tissues (Staining: Low, Intensity: Weak, Quantity: 75%-25%). VEGFA staining in normal pancreatic tissue (Staining: High, Intensity: Strong, Quantity:>75%) was also significantly higher than in PAAD tissue (Staining: Low, Intensity: Weak, Quantity: >75%). Interestingly, however, the immunohistochemical results of SOD1 were not completely consistent with the differential expression analysis. SOD1 was lowly expressed in exocrine cells of normal pancreatic tissues (Staining: Low, Intensity: Weak, Quantity:>75%); moderately expressed in endocrine cells (Staining: Medium. Intensity: Moderate, Quantity:>75%); but also generally low in PAAD tissues (Staining: Low, Intensity: Weak, Quantity: >75%). We speculate that the inconsistency is related to the overall low expression of SOD1.

FIGURE 2

Figure 2 Development of immune-related CRGs. (A) Intersection of IRGs from the Immport database, 10 immune-related subsets of IRGs from GSEA and DEGs from the TCGA-PAAD cohort (B) 14 CRGs have significant interactions with the 35 IRGs screened (Green nodes: CRGs; Red/Blue nodes: IRGs with high/low expression in PAAD; arrows: direction of action; nodes size: HR value) (C–E) Differential expression of MAP2K2, SOD1, and VEGFA in TCGA-PAAD cohort and GTEx cohort ("***": P<0.001) (F) Immunohistochemical results.

3.2 Immune landscape of differentially expressed CRGs

Based on the screening of 3 genes that were differentially expressed at the PAAD transcriptome level, their immune microenvironment was further analyzed in the TCGA-PAAD cohort. The high and low expression groups of MAP2K2, SOD1, and VEGFA were first scored for TME (Figure 3A-C). The results showed that the Stromal Score, Immune Score, and ESTIMATE Score of the low expression group of these 3 genes were higher than the high expression group. This indicates that for MAP2K2, SOD1, and VEGFA, the stromal cell component and immune cell component were higher in the low expression group than in the high expression group, suggesting that immunosuppressive microenvironment was prevalent in the high expression group. Interestingly, the expression of MAP2K2 and VEGFA did not match the TME scoring. To further understand the specific composition of TME in the high- and low-expression groups, we analyzed the immune cell infiltration of the three CRGs (Figure 3D-F). The results showed that T cells CD4 memory resting and T cells CD4 memory activated were mainly distributed in the MAP2K2 low expression group; T cells regulatory (Tregs) and NK cells activated were mainly distributed in the MAP2K2 high expression group. Plasma cells, T cells CD8, T cells CD4 memory activated, and Mast cells resting were mainly distributed in the low VEGFA expression group; Macrophages M0 and T cells CD4 memory resting were mainly distributed in the high VEGFA expression group. And the difference of SOD1 was not significant. To judge the accuracy of immune cell infiltration analysis, the correlation between MAP2K2, SOD1, and VEGFA with immune cells was analyzed separately (Figure 3G-I), and the results were consistent with previous. Among them, Macrophages, as the immune cell component associated with all three CRGs, Macrophages M0 was positively correlated with the expression of MAP2K2 (Cor=0.1557, P=0.0420) and VEGFA (Cor=0.2755, P=0.0003), and Macrophages M1 was negatively correlated with the expression of SOD1(Cor=-0.1538, P=0.0446). The 47 immune check loci were weakly to moderately correlated with MAP2K2, SOD1, and VEGFA (Figure 3J-L).

FIGURE 3

Figure 3 Immune landscape of differentially expressed CRGs. TME evaluation (A–C) and immune cell infiltration (D–F) of MAP2K2, SOD1, and VEGFA high and low expression groups ("***": P<0.001, "**": P<0.01, "*": P<0.05). Correlation of MAP2K2, SOD1, and VEGFA expression profiles with immune cells (G–I) and 47 immune checkpoints (J–L).

3.3 Expression patterns of immune-related CRGs at single-cell level

All 6 cancer tissue samples from 6 PAAD patients of GSE212966 participated in this study. After quality control (200< nFeatureRNA< 2500, per cent.mt< 5), a total of 9431 single cells and 28435 genes were included in the study for further analysis (Supplementary Figure S1A, 1B). Normalization was performed on the basis of the first 2000 highly variable genes, of which the first 10 genes are presented in the figure (Supplementary Figure S1C). The key genes and their importance in the first 2 principal components exhibited satisfactory between-group variability (Supplementary Figure S1D). As shown, these cells could be divided into 10 major cell clusters (Figure 4A). Among them, the first 5 marker genes of each cell cluster had significant intergroup differences (Supplementary Figure S1E), and these marker genes were the basis of the annotation. By annotation, these 10 cell clusters were clearly grouped into 7 categories (Figure 4B). In addition to Epithelial cells, which contain cancer cells, Tissue stem cells and 5 types of immune cells (T cells, Neutrophils, B cells, Macrophages, and NK cells) are also present. By comparing the expression patterns of CRGs in single cells, 19 of these genes were found to be significantly expressed in PAAD samples. Among them, COX17, FDX1, GLS, MT1E, MT1F, MT1G, MT1X, MT2A, and PDE3B were mainly clustered in NK cell clusters (Supplementary Figure S2A). VEGFA was mainly clustered in Neutrophils clusters (Supplementary Figure S2B). ATOX1 was mainly clustered in Macrophages (Supplementary Figure S2C). MAP2K2 was mainly clustered in Epithelial cells (Supplementary Figure S2D). SOD1, MT-CO2, UBE2D2, and UBE2D3 were commonly and significantly expressed in PAAD monocytes (Supplementary Figure S2E). Next, to further explore the differentially expressed 3 CRGs at the transcriptome level, we performed an in-depth analysis of their specific expression levels and distribution. Among them, MAP2K2 was mainly expressed in Epithelial cells, B cells, and T cells clusters (Figure 4C), SOD1 was mainly expressed in Endothelial cells and T cells clusters (Figure 4D), and VEGFA was mainly expressed in Endothelial cells, Neutrophils and T cell clusters (Figure 4E). The annotation accuracy was checked and the distribution was verified by quantifying the expression (Figure 4F, G), which was consistent with the counts results. In addition to this, the paraneoplastic normal 3-case samples demonstrated satisfactory quality under the same QC treatments described above (Supplementary Figure S3). The patterns of 3 differentially expressed immune-associated CRGs, MAP2K2, SOD1, and VEGFA, in 3 single-cell samples of normal tissues adjacent to cancer were further analyzed as controls (Supplementary Figure S4A-S4B). As seen in the figure, MAP2K2 was mainly expressed in NK cells and CMP clusters (Supplementary Figure S4C), SOD1 was commonly and significantly expressed in all clusters (Supplementary Figure S4D), and VEGFA was mainly expressed in monocytes and CMP clusters (Supplementary Figure S4E). These results were validated using the same method as above (Supplementary Figure S4F-S4G). The single-cell results of PAAD showed more pronounced infiltration of B cells, macrophages, NK cells, and T cells compared to normal tissues, which suggests that the inflammatory features of PAAD seem to change with disease progression. Furthermore, single-cell sequencing data of the 9-case PAAD samples and their adjacent 9-case normal tissue samples from GSE205049 were used as external validation and support for the above single-cell analysis results. The analysis results showed that the 9 PAAD samples could be annotated into 5 cell clusters (Supplementary Figure S5A, 5B), whereas the 9 normal para-cancerous samples could be annotated into 7 cell clusters (Supplementary Figure S6A, 6B). MAP2K2 was expressed predominantly in the NK cell, T cell, and monocyte clusters (Supplementary Figure S5C). SOD1 was commonly and significantly expressed across all the clusters (Supplementary Figure S5D). VEGFA was expressed mainly in the monocyte and NK cell clusters (Supplementary Figure S5E). Compared to normal tissues (Supplementary Figure S6C-6E), the single-cell results of PAAD showed more dense and more pronounced immune cells infiltration (T cells, B cells, NK cells and monocytes) under the same filtration conditions. This is consistent with the results of the GSE212966 analysis. Annotation accuracy heatmaps and expression distribution violin plots demonstrated satisfactory reliability of the conclusions (Supplementary Figures S5F, 5G, 6F, 6G). And the differential distribution of CRGs at the single cell level in PAAD and its adjacent normal tissues suggests that the occurrence of Cuproptosis is not only associated with the cancer cells themselves but also with the type of immune cells infiltrating the inflammatory TIME in PAAD. However, the exact mechanism needs to be further explored.

FIGURE 4

Figure 4 Expression patterns of immune-related CRGs at single-cell level. (A) 10 major cell clusters in the scRNA-PAAD cohort (B) 7 cell types identified based on cellular markers: Epithelial cells, Tissue stem cells, T cells, Neutrophils, B cells Macrophages, and NK cells (C-E) Expression patterns of the differentially expressed 3 CRGs in cancer tissues: MAP2K2, SOD1, and VEGFA (F) Accuracy of cell cluster annotation identification (G) Specific expression levels of the 3 CRGs in each cell cluster.

3.4 Establishment of Cuproptosis-related signature to assess individual prognosis

A total of 37 CRIs were identified in PAAD based on the 7 CRGs (Supplementary Table S6) that had interactions with IRGs (|PearsonR| > 0.4 and p< 0.001). And most of the CRGs were positively correlated with CRIs (Figure 5A, Supplementary Table S7). 9 differentially expressed CRIs with prognostic value were obtained by Lasso-Cox regression (Figure 5B, C) and one-way regression analysis (Figure 5D), namely, LIF-AS1, AC007292.2, AC015660.1, MEG9, LINC01133, FAM27E3, AC092171.5, LINC01091, and AC007292.1. Based on these 9 CRIs and multi-factor regression analysis, 4 key CRIs were extracted to construct the signature, including AC007292.2, AC015660.1, LINC01091, and MEG9 (Figure 5E). We named this signature of 4-CRIs as CIR-score (Table 1). In addition, we respectively performed PCA downscaling analysis of the TCGA-PAAD cohort according to CRGs, CRIs, and risk lncRNAs of CIR-score (Figure 5F–H). The results showed that the descending results conditional on the CIR-score were significantly better than the other two groups, validating the reliability of the CIR-score.

FIGURE 5

Figure 5 Establishment of Cuproptosis-related signature to assess individual prognosis. (A) Co-expression relationship of 7 CRGs with CRIs (B, C) Lasso-Cox regression analysis (D) Univariate COX regression analysis demonstrating 9 CRIs and their P value and Hazard Ratio (E) Risk lncRNAs of CIR-score and 7 CRGs correlation (F-H) PCA analysis according to CRGs, CRIs, and CIR-score.

TABLE 1

Table 1 Hub lncRNAs and their correlation coefficients of CIR-score.

3.5 Validation of the independent predictive power of the CIR-score to assess individual prognosis

Whether the CIR-score is an indicator with independent prognostic value was explored by performing univariate and multifactorial Cox regression analyses in the training set of TCGA-PAAD (Figure 6A, B). The results showed that in the univariate Cox regression analysis, age (p=0.015, Hazard Ratio=1.026), Grade (p=0.028, Hazard Ratio=1.383), and CIR-score (p<0.001, Hazard Ratio=1.276) of PAAD patients were significantly associated with prognosis correlation. The ROC curve showed good results in evaluating the survival of PAAD patients at 1, 3, and 5 years (Figure 6C). Additionally, the AUC of CIR-score was the highest compared to the AUC of age, gender, tumor grade, and clinical stage at 0.697 (Figure 6D). C-index curves showed the same results as the former (Figure 6E). This evidence suggests that the CIR-score has a better ability to predict prognosis compared with other clinical characteristics. On this basis, we constructed a Nomogram based on the CIR-score and other clinically relevant indicators to assess the 1-, 2-, and 3-year survival of individual patients (Figure 6F). Its calibration curve also confirmed its good predictive power (Figure 6G).

FIGURE 6

Figure 6 Validation of the independent predictive power of the CIR-score to assess individual prognosis. (A, B) Univariate and multifactorial Cox analysis considering the CIR-score and clinical characteristics (C) ROC curves for 1-, 3- and 5-year overall survival (D) ROC curves for the CIR-score and other clinical indicators (E) CIR-score and other clinical indicators c- index curves (F) A Nomogram that predicts 1-, 2-, and 3-year survival in patients with PAAD ("***": P<0.001, "**": P<0.01, "*": P<0.05) (G) Calibration curves for this Nomogram.

3.6 Prognostic power and clinical relevance of CIR-score

Based on the CIR-score, samples from the TCGA-PAAD cohort, train set, and test set were scored separately and divided into high- and low-risk groups based on the median (Supplementary Table S8). No between-group differences in clinical risk factors between the train and test sets demonstrated the reliability of the randomized groupings (Supplementary Table S9). The standardized distribution of CIR-score, scatter plots of survival time and survival status, and heat maps of expression of risk lncRNAs were plotted for the high- and low-risk groups, respectively (Figure 7A-C). The results showed that the survival time and survival rate were significantly higher in the low-risk group than in the high-risk group. The differences in the expression levels of the 4 risk lncRNAs of the CIR-score were also consistent with each other in the high- and low-risk groups. Furthermore, survival analysis and Kaplan-Meier curves were performed on samples from the TCGA-PAAD cohort, train set, and test set respectively (Figure 7D). As expected, PAAD cases with higher CIR-score obtained significantly worse OS in the TCGA-PAAD cohort (p< 0.001), train set (p< 0.001), and test set (p=0.006). This is consistent with Progression Free Survival (PFS) in the TCGA-PAAD cohort (p= 0.004), and test set (p= 0.027) (Figure 7E). However, the variability of PFS in the train set was not significant (p= 0.071). We speculate that this is related to the size of the sample size. Then, the TCGA-PAAD cohort was subjected to CIR-score scoring, and the sample was divided into 2 groups of high and low risk according to the median. The TCGA-PAAD cohort was subgrouped according to different prognosis-related clinical traits, and survival analysis was performed for each group of samples (Figure 7F). Results showed that Survival probability was significantly higher in the low-risk group than in the high-risk group in the Grade 1-2 cohort (p= 0.005) and the Grade 3-4 cohort (p< 0.001). In the Stage I-II cohort (p< 0.001), higher risk scores were consistent with survival. In addition, in cohorts differentiated according to TNM staging, risk scores were consistent with survival in the T1-2 cohort (p= 0.029), T3-4 cohort (p< 0.001), and N1 cohort (p< 0.001). Nevertheless, in the Stage III-IV cohort (p= 0.055) and N0 cohort (p= 0.058), the differences were not significant due to insufficient sample size or too many missing values.

FIGURE 7

Figure 7 Prognostic power and clinical relevance of CIR-score. (A-C) Standardized distribution of scores, distribution of survival time and survival status, and differences in expression of risk lncRNAs in the TCGA-PAAD cohort, train set, and test set with high- and low-risk (D) Survival analysis of these 3 groups (E) Differences in progression-free survival of these 3 groups (F) Survival analysis of the high- and low-risk groups for different clinical indicators.

3.7 Differences in the TIME landscape between high and low-risk groups

Correlation analysis of immune activity was performed for the high- and low-risk groups according to the CIR-score (Figure 8A). The results showed that the immune function scores of the low-risk group were generally higher than those of the high-risk group, suggesting a significant state of immunosuppression in the high-risk group. Next, the fraction of all 22 immune cells in the high- and low-risk groups was explored (Figure 8B). The most significant intergroup variability was found for Macrophages M0, NK cells resting, T cells CD8, and T cells regulatory (Tregs). Among them, Macrophages M0 and NK cells resting were more significantly infiltrated in the high-risk group. T cells CD8 and T cells regulatory were more significantly infiltrated in the low-risk group. That is, Macrophages M0 (R = 0.43, p = 4.4e-05) and NK cells resting (R = 0.27, p = 0.012) were positively correlated with CIR-score, and T cells CD8 (R = -0.26, p = 0.015) and T cells regulatory (R = - 0.23, p = 0.035) showed a negative correlation with CIR-score (Figure 8C-F). Besides, by summarizing the immune checkpoint genes of different immune cells (26), we further explored the correlation of CIR-score with these four immune microenvironment cells (Figure 8G-J). The results showed that the CIR-score was moderately correlated with the immune checkpoints of Macrophage, Natural killer cell, Activated CD8 T cell, and Regulatory T cell.

FIGURE 8

Figure 8 Differences in the TIME landscape between high and low-risk groups. (A) Differences between immune-related activities in high and low-risk groups ("*": P<0.05) (B) Immune cell infiltration in tumor microenvironment in high and low-risk groups ("**": P<0.01) (C–F) Correlation analysis between immune cells (Macrophages M0, NK cells resting, T cells CD8, T cells regulatory) and CIR-score (G–J) Correlation analysis between immune checkpoints (Macrophage, NK cells, Activated CD8 T cells, Regulatory T cells) and CIR-score.

3.8 Enrichment analysis, tumor burden, and immune checkpoint sensitivity prediction

To further analyze the intergroup differences between high and low-risk groups, DEGs of high and low-risk groups were screened. Then, GO and KEGG enrichment analyses were performed on these DEGs. In GO analysis (Supplementary Table S10), DEGs were mainly enriched in signaling as well as hormone-related biological processes (BP), cytological components (CC) such as presynapse, and molecular biological functions (MF) related to multiple channel activities (Figure 9A). In KEGG analysis (Supplementary Table S11), Neuroactive ligand-receptor interaction, Pathways of neurodegeneration-multiple diseases, and Dopaminergic synapse related pathways were represented in a larger proportion (Figure 9B). Therefore, we hypothesized that DEGs mainly function in cell signaling and hormonal regulation. Next, we assessed Tumor Mutation Burden (TMB) in the high- and low-risk groups (Figure 9C, D). The high-risk group had a higher TMB compared to the low-risk group (p=0.0047). And the high TMB group had a significantly worse prognosis than the low-risk group (p=0.008). By further analysis of the high-low risk group versus the high-low TMB group (Figure 9E), it was obtained that the CIR-score had a stronger potential to predict survival than TMB (p<0.001). Additionally, to explore the potential clinical efficacy of immunotherapy in different subgroups, the Tumor Immune Dysfunction and Exclusion score (TIDE) was applied (Figure 9F). The TIDE prediction score was significantly higher in the low-risk group than in the high-risk group, indicating that patients in the low-risk group were more likely to experience immune escape and patients were less likely to benefit from Immune Checkpoint Inhibitor (ICI) therapy (p<0.01). However, given the uncertainty of TIDE for tumor scores other than the melanoma dataset and the non-small cell lung cancer dataset, we still need to make further predictions about the likelihood of drug treatment in the high- and low-risk groups.

FIGURE 9

Figure 9 Enrichment analysis, tumor burden, and immune checkpoint sensitivity prediction. (A) Results of GO enrichment analysis (Biological Process, Cellular Components, and Molecular Function) (B) Results of KEGG enrichment analysis (Pathways) (C) TMB differences between high and low-risk groups (D) High and low TMB groups for survival analysis (E) Survival analysis of high TMB-high risk group, high TMB-low risk group, low TMB-high risk group, and low TMB-low risk group (F) TIDE assessment of high and low-risk groups ("**": P<0.01).

3.9 Drug sensitivity intergroup differences

To deeply explore the different focus of drug therapy in patients in high and low-risk groups. The Immune cells Proportion Score (IPS) differences between the high and low expression groups of CRGs and CRIs and the high and low-risk groups of CIR-score were evaluated based on the expression levels of hub CRGs and lncRNAs, and the differences in immunotherapy sensitivity between the different groups were compared. Higher IPS was associated with higher immunogenicity. The results indicated that although the differences in PD-1 were not significant in any of the 3 CRGs high and low expression groups. Nevertheless, for hub lncRNAs, the AC007292.2 low expression group (Figure 10A), AC015660.1 high expression group (Figure 10B), LINC01091 low expression group (Figure 10C), MEG9 low expression group (Figure 10D), and CIR-score high-risk group (Figure 10E) possessed significantly higher immunogenicity. Moreover, the difference in PD-1 monoclonal antibody sensitivity between the CIR-score high- and low-risk groups was far more significant than that of single hub lncRNA, which demonstrated that CIR-score is more significant for guiding immunotherapy dosing strategies than single hub lncRNA and single CRG. In addition, considering that Nanoalbumin Paclitaxel (Abraxane^®) in combination with Gemcitabine is the preferred chemotherapy for non-surgical PAAD patients (27), we further probed into the intergroup differences in the IC50 of these two agents. A lower IC50 was associated with higher drug sensitivity. Paclitaxel was found to be more sensitive in the VEGFA high expression group, AC015660.1 high expression group, LINC01091 low expression group, and CIR-score high-risk group (Figure 10F). Gemcitabine was more sensitive in the LINC01091 low expression group, and MEG9 low expression group (Figure 10G). The above results suggest that patients with lower CIR-score are more likely to be resistant to immunotherapy and chemotherapy.

FIGURE 10

Figure 10 Drug sensitivity intergroup differences. (A-D) Immunogenicity differences of the 4 hub lncRNAs: AC007292.2, AC015660.1, LINC01091, and MEG9 with different expression (E) Immunogenicity differences between the high and low-risk groups (F) IC50 differences of Paclitaxel in VEGFA, AC015660.1, LINC01091, and CIR-score with different levels (G) IC50 differences of Gemcitabine in LINC01091 and MEG9 with different levels.

3.10 Validation of hub CRGs and lncRNAs in immune combined targeted and chemotherapy mouse PAAD model

The clinical benefit of monotherapy is limited to a small population. Currently, the efficacy of immunotherapy (PD-1 monoclonal antibody) alone plus chemotherapy (28) and immunotherapy combined with other molecules is suboptimal (29). In contrast, chemotherapy combined with recombinant human vascular endothelial inhibitor (Endostar ^®) has achieved some achievements in clinical treatment (30). Therefore, to explore new strategies for combination therapy, 4 PAAD mouse treatment cohorts of model control group, chemotherapy group (Nanoalbumin Paclitaxel + Gemcitabine), recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody group, and recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody + chemotherapy group were included in this study. Although we assessed the efficacy of immunotherapy in patients by screening immune-related CRGs as raw material for the construction of CIR-score and characterization, it is necessary to directly compare the response rates of hub CRGs and lncRNAs under different treatment regimens in a mouse PAAD model. Considering that MAP2K2 is positively correlated with AC007292.2 (Coef. = -2.898), LINC01091 (Coef. = -0.779), and MEG9 (Coef. = -0.852) (Figure 5E), indicating that the lower MAP2K2 expression had a higher CIR-score. Similarly, SOD1 was mildly negatively correlated with AC007292.2, LINC01091, and MEG9, and mildly positively correlated with AC015660.1 (Coef. = 0.497). VEGFA is positively correlated with AC007292.2, AC015660.1, LINC01091, and MEG9, but the correlation coefficient with AC015660.1 is the highest. That is, the higher the expression of SOD1 and VEGFA the higher the CIR-score. Therefore, the immunohistochemical results showed that the recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody + chemotherapy group had the lowest CIR-score, the model control group had the highest CIR-score, and the recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody group was slightly higher than the chemotherapy alone group (Figure 11A). In contrast, with almost all current drugs, cancer tissue becomes resistant over time as patients use them (31). Compared to those already on the drug, cancer patients who did not experience treatment were more sensitive to the same drug. And the high-scoring group was more sensitive to immunotherapy (Figure 10E). Thus, the differences in MAP2K2, SOD1, and VEGFA expression in different treatment cohorts were consistent with intergroup differences in immunotherapy sensitivity. Furthermore, MAP2K2 expression was significantly higher in the chemotherapy group than in the non-chemotherapy group, and SOD1 and VEGFA expression was significantly lower in the chemotherapy group than in the non-chemotherapy group. That is, the scores in the chemotherapy group were significantly lower than those in the non-chemotherapy group. This is also consistent with the result that the higher-scoring group was more sensitive to chemotherapy (Figure 10F, G). On the other hand, qRT-PCR results showed that the expression of MEG9 was significantly higher in the immunotherapy cohort in the presence of PD-1 monoclonal antibody treatment than in the non-immunotherapy cohort (Figure 11B), which is also consistent with the analysis of sensitivity differences between immunotherapy groups. Unfortunately, we did not retrieve the sequences of AC007292.2, AC015660.1, and LINC01091 in mice, and therefore did not validate them by qRT-PCR in mouse models. These 3 lncRNAs will be planned to be validated on human PAAD tissues in the future. The above results further demonstrate that CIR-score is qualified to predict the efficacy of immunotherapy in PAAD patients and provide guidance for individualized and precise application of immunotherapy. Interestingly, in vivo imaging showed that the Total Radiant Efficiency of recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody + chemotherapy group was significantly lower than the other 3 groups after 1 and 2 weeks of administration (Figure 11C, D). It indicates that immunotherapy combined with targeted and chemotherapy significantly improved the proliferation and progression of PAAD, which provides a new possibility for a multidrug combination therapy strategy for PAAD with immunotherapy as the main treatment.

FIGURE 11

Figure 11 Validation of hub CRGs and lncRNAs in immune combined targeted and chemotherapy mouse PAAD model. (A) Immunohistochemical results of 3 differentially expressed immune-related CRGs in mouse PAAD model (a: Model Control Group; b: Chemotherapy Group; c: Recombinant human endostatin + PD-1 Monoclonal Antibody Group; d: Recombinant human endostatin + PD-1 Monoclonal Antibody + Chemotherapy Group) (B) Hub lncRNA MEG9 relative expression differences in the 4 treatment cohorts ("***": P<0.001, "**": P<0.01, "*": P<0.05) (C) In vivo imaging of PANC-02 transplanted tumors in different treatment cohorts of PAAD mice at 0, 1, and 2 weeks (D) Differences in total radiation efficiency in different treatment cohorts at 0, 1, and 2 weeks.

4 Discussion

Cuproptosis is a novel cell death mechanism different from the known ones. Copper ions can mediate cell death through a mitochondria-dependent increase in energy metabolism and cytotoxicity induced by reactive oxygen species (ROS) accumulation (10). Copper deficiency can lead to early embryonic death or congenital malformations. This is due to copper’s role as a cofactor for mitochondrial cytochrome c oxidase, which is required to meet the energy requirements of rapidly dividing cells (24). Similarly, the copper requirement of cancer cells is higher compared to non-dividing cells. First, in studies of PAAD (11), gastric cancer (32), hepatocellular carcinoma (33), gallbladder cancer (34), lung cancer (35), thyroid cancer (36), and prostate cancer (37), serum copper ion levels were found to be significantly higher in tumor patients compared to normal patients. Higher serum copper ion concentrations in lung cancer (35) and breast cancer (38) were also associated with poorer clinical prognosis. Second, copper also reflects TIME status at the immune cell level. Copper ion concentration is inversely correlated with the degree of infiltration of T cells and macrophages in TIME. In mesothelioma, reducing bioavailable copper slows tumor growth, normalizes blood vessels, and promotes T-cell infiltration (39). In hepatocellular carcinoma, copper with disulfiram inhibits T-cell infiltration by inhibiting PARP1 activity and enhancing GSK3β Ser9 site phosphorylation, which in turn upregulates PD-L1 expression (40). In breast cancer, the copper-depleting compound tetrathiomolybdate (TM) resulted in reduced collagen deposition, decreased levels of myeloid-derived suppressor cells, and increased infiltration of CD4+ T cells (41). In addition, downregulation of the copper-containing metabolism MURR1 structural domain 1 (COMMD1) in cancer cells enhances the inflammatory response and creates favorable conditions for macrophage recruitment (42). Therefore, the present study screened for immune-related CRGs by exploring the interaction between IRGs and CRGs. The results revealed that among the 14 CRGs that interacted significantly with 35 IRGs, VEGFA, SOD1, and MAP2K2 had the most significant differential expression between PAAD and adjacent normal tissues. Among them, VEGFA encodes vascular endothelial growth factor A, which induces angiogenesis, vascular endothelial cell growth, and increased vascular permeability. Although VEGFA is upregulated in many tumors and is a target for many cancer therapies. However, approximately 90% of PAAD are ductal adenocarcinomas with KRAS mutations as the primary driver, which are typically characterized by fibrous tissue hyperplasia and reduced vascularity, which in turn leads to immunosuppression and resistance to chemotherapy and immunotherapy (43). Thus, the low expression of VEGFA in PAAD tissues is consistent with the low vascular density characteristic of cancer tissues. Copper-zinc superoxide dismutase encoded by SOD1 is one of the predominant intracellular antioxidant enzymes and an important target for cancer drug design. Overexpression of SOD1 in PAAD protects cancer cells from oxidative stress (44). Hence, SOD1 is an important means of regulating superoxide anion radical (O2–) and hydrogen peroxide (H2O2) levels in cancer cells. However, most of the reported SOD1 inhibitors such as ATN-224 (45), DDC (46), and LD100 (47) have shortcomings like inefficient inhibition or lack of selectivity for other copper-containing proteins. MAP2K2, one of the core drivers of cancer, encodes mitogen-activated protein kinase kinase 2 (MEK2) which plays a key role in mitogenic growth factor signaling, and knockdown of MEK2 inhibits the invasive ability of PAAD cancer cells (48). Interestingly, however, the expression level of MAP2K2 in PAAD tissues was not exactly consistent with previous studies. We speculate that this correlates with the cell subpopulation category in which MAP2K2 is highly expressed. To further explain the high expression of MAP2K2 in PAAD paracancerous tissues and to investigate the differences in the spatial distribution of CRGs expression levels, the present study further explored the expression patterns of all 14 CRGs with co-expression relationships with IRGs at the single-cell level.

The transformation of normal pancreatic follicular cells to PAAD is often accompanied by the appearance of chronic inflammation, which induces and chemotacticizes a large number of immune cells distributed in the dense pancreatic interstitial tissue. These immune cells undergo altered functions in TIME, which eventually produce amplified immunosuppressive signals leading to immune escape of tumor cells and promoting tumorigenesis (49). Considering the complexity and highly suppressive nature of PAAD TIME, single-cell sequencing, as an emerging technology, provides a new perspective for the study of TIME. On the one hand, single-cell sequencing technology can better reveal the heterogeneity of tumor cells at the cellular and molecular levels, capture different tumor states, and provide new perspectives for understanding tumor drug resistance (50). On the other hand, cell clusters and cell states with similar gene expression profiles can be more thoroughly analyzed from a multi-omics perspective (genomics, transcriptomics, epigenomics, and proteomics) (51). To further visualize the differences in expression of CRGs at the cellular level between PAAD and its adjacent normal tissues, and to further explore the association of CRGs with different immune cell subpopulations in TIME, all screened CRGs were analyzed in depth in this study in single-cell sequencing samples. Among them, SOD1, MT-CO2, UBE2D2, and UBE2D3 were commonly and significantly expressed in each single cell subpopulation of PAAD, and the cluster with the highest number of aggregated genes was the NK cells cluster (Supplementary Figure S2). Overall, for the 3 immune-related CRGs that were significantly differentially expressed in the transcriptomic data, the single-cell data were consistent with the results of transcriptomic data analysis. Among them, SOD1 was significantly more expressed in each cell cluster in PAAD tissues than normal. Previous studies have shown that SOD1 overexpression not only protects cancer cells from oxidative stress but also suppresses pro-inflammatory immune responses in colitis by preventing oxidative stress (52). We speculate that similarly, SOD1 overexpression in TIME could inhibit the immune response by suppressing immune cells and their secreted cytokines and enzymes. Second, VEGFA expression in neutrophils, NK cells, and epithelial cells was significantly higher in PAAD tissues than normal. However, in CMP cell clusters the cancer VEGFA expression was significantly lower than normal. Moreover, MAP2K2, as an oncogene, plays a key role in mitogenic growth factor signaling. But single-cell sequencing analysis showed that MAP2K2 was significantly more expressed in NK and CMP cell subsets in normal pancreatic tissues adjacent to cancer than in PAAD tissues. We speculate that the main reason for the difference in MAP2K2 transcriptome data does not lie in the cancer cells in the tissues, but in the NK cells in the TIME.

lncRNAs are promising circulating biomarkers in cancer diagnosis and prognosis, and exceptional candidates for further therapeutic exploration. However, considering that individual genes or lncRNAs are not sufficient to achieve the desired prediction, a new form of indicator is needed to provide new strategies for individualized therapy. lncRNA signature is a risk indicator that can systematically predict cancer survival prognosis, immune landscape, and drug sensitivity by combining multiple lncRNA molecules as variables assigned with their respective coefficients. It can circumvent the limitations of individual lncRNAs in predicting prognosis, but also avoid the defects of interactions among lncRNAs. Currently, studies on the prognostic signature of Cuproptosis-associated lncRNAs have been conducted in tumors such as bladder cancer (53), head and neck squamous cell carcinoma (54), and gastric cancer (55). However, most studies have used co-expression analysis for screening of lncRNAs, but few studies have systematically explored the genes from which lncRNAs originate, as well as their distribution and expression differences at the single-cell level. Furthermore, the role of Cuproptosis-associated lncRNAs, in the regulation of PAAD and its immunotherapy is unclear. Therefore, this study screened for 7-CRGs based on single-cell sequencing, which were both co-expressed with IRGs and significantly expressed in single cells. The immune-related 4-CRGs with the most prognostic value were identified by co-expression and regression analysis: MEG9, LINC01091, AC015660.1, and AC007292.2. Among them, MEG9 could protect endothelial cells from DNA damage-induced cell death and predict cancer prognosis through the PI3K-AKT signaling pathway (56). In addition, MEG9 was significantly associated with hepatitis B virus (HBV) infection status in hepatocellular carcinoma (57), and epidermal growth factor receptor (EGFR) status in adenocarcinoma patients with non-small-cell lung cancer (NSCLC) (58). In the PAAD field, MEG9 is involved in the construction of m6A-related prognostic signature as hub lncRNA (59). LINC01091 can coordinate the microRNA-128-3p/ELF4/CDX2 axis and thus promote gastric cancer growth and metastasis by way of exosomes on one hand (60) and can act as hub lncRNA in DNA methylation-related prognostic signatures of prostate cancer (61). AC015660.1 is involved in the angiogenesis of PAAD (62) and inflammatory correlation of gastric cancer (63) as hub lncRNA of prognostic signatures. In contrast, AC007292.2 has no previous relevant studies. To verify the prognostic value, this study further discussed the survival, clinical relevance, bioenrichment, and gene mutations in the high- and low-risk groups according to the CIR-score median. The results showed that the CIR-score could more clearly classify PAAD patients into high- and low-risk clusters compared to CRGs or CRIs alone. And CIR-score also showed a significant positive correlation with risk factors such as pathological grade, clinical grade, and TMB.

Furthermore, differences in the immune landscape between the high and low-expression groups of CRGs and the high and low-risk groups of the CIR-score raised additional concerns. Immune cells, tumor-associated fibroblasts (CAFs), vascular system, and extracellular matrix (ECM) are the 4 major components that constitute the highly immunosuppressive microenvironment of PAAD (64). 2 of the major types of immune cells are most associated with immunosuppressive TIME in PAAD: tumor-associated macrophages (TAMs) and tumor-infiltrating T cells. Previous studies have confirmed that TAMs not only promote tumor initiation but also serve as central drivers of immunosuppressive TIME by expressing cell surface receptors, secreting cytokines, chemokines, and enzymes that regulate the recruitment and function of multiple immune cell subtypes (65). In this study, we found that macrophage M0 was positively correlated with the expression of MAP2K2 and VEGFA; macrophage M1 and monocytes were negatively correlated with the expression of SOD1. The expression of CRGs in macrophage clusters of PAAD single-cell samples was significantly higher than in normal tissues. M0 is the state in which macrophages are not activated. In contrast, activated macrophages have the typical activation phenotype M1 and the selective activation phenotype M2. M1 subtypes not only promote antitumor immune responses but also increase radiotherapy sensitivity (66, 67). Higher numbers of M2 subtypes are associated with larger PAAD tumor volumes, early liver and local recurrence, accelerated lymphatic metastasis, and shorter survival (68). However, activation of TAMs inhibits TLR-mediated M1 polarization (69). To address this feature Guiducci C et al. investigated a TLR9 ligand that induced a shift from the M2 phenotype to M1 in macrophages and resulted in tumor shrinkage (70). Thus, targeting macrophages provides a direction for improving TIME in PAAD. On the other hand, this study found significantly lower infiltration of regulatory T cells (Tregs) in the high-risk group than in the low-risk group. Zhang Y et al. suggested that depletion of Tregs accelerates pancreatic carcinogenesis, leads to differentiation of inflammatory fibroblast subpopulations, and consequently leads to compensatory immunosuppression in advanced disease (71). The more pronounced Tregs depletion and worse OS in the high-risk group are consistent with previous studies. This cross-talk between Tregs and fibroblasts in PAAD reveals a potential new therapeutic approach to alleviate immunosuppression in pancreatic cancer.

Immunotherapy is considered to be the fourth major oncology treatment after surgery, chemotherapy and radiotherapy. Although good efficacy has been achieved in some advanced tumors, the overall effect of immunotherapy in PAAD is poor. Compared with the poorer response to PD-1 inhibition alone, combination therapy with targeted chemotherapy and other agents has shown a more durable response, providing new ideas for the clinical treatment of PAAD (72). Recombinant human vascular endothelial inhibitor alleviates immunosuppression by normalizing tumor vasculature (73) and specifically replicates in PAAD cancer cells and kills them (74). PD-1 monoclonal antibody in combination with recombinant human vascular endothelial inhibitor and chemotherapy has shown good efficacy and safety in patients with advanced non-small cell lung cancer (NSCLC) (75). Combining recombinant human vascular endothelial inhibitor with first-line standard chemotherapy also improved the quality of life in patients with small cell lung cancer (30). However, no such combination therapy has been reported in the field of PAAD. Considering that Nanoalbumin Paclitaxel plus Gemcitabine is the chemotherapy of choice for PAAD patients who are non-surgical candidates but in good physical condition (27). And the CIR-score high and low-scoring groups and their associated high and low-expression groups of CRGs and CRIs exhibited significant drug sensitivity differences. To explore a more optimal combination treatment strategy for PAAD and to validate the reliability of the CIR-score, this study divided the mouse model into 4 treatment cohorts: model control group, chemotherapy group (Nanoalbumin Paclitaxel + Gemcitabine), recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody group, and recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody + chemotherapy group. MAP2K2, SOD1, and VEGFA, as the most prognostically valuable CRGs screened and co-expressed with CIR-score, were differentially expressed in the chemotherapy and non-chemotherapy groups in agreement with the pharmacoresensitivity prediction results (Figure 10). The qRT-PCR results of MEG9 in the hub lncRNAs of CIR-score were also consistent with the intergroup differences in immunotherapy sensitivity. These results further demonstrate that the CIR-score is qualified to predict immunotherapy efficacy in PAAD patients and provide guidance for individualized and precise application of immunotherapeutic agents. In addition, the recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody + chemotherapy group in the in vivo imaging results demonstrated a significant advantage in delaying tumor progression. Therefore, not only the CIR-score can predict the effect of immunotherapy in patients and reflect the difference of drug sensitivity in patients; but also immunotherapy combined with chemotherapy and targeted therapy can be a direction of clinical drug treatment. The combination of these two aspects suggests a focused multi-modality combination therapy for specific patients’ signature scores for the purpose of achieving optimal therapeutic outcomes, providing a new strategy for individualized treatment of PAAD.

There are also some limitations in this study. Firstly, unfortunately, we were unable to quantify all 4 hub lncRNAs by qRT-PCR in a mouse model because we did not retrieve the sequences of AC007292.2, AC015660.1, and LINC01091 in mice. These 3 IncRNAs will be planned to be further validated on human PAAD tissues in the future. Secondly, because of relying mainly on public databases and animal experiments, this study still needs to further investigate and validate of the value of prognostic CIR-score in clinical PAAD cases, the value of reflecting the immune landscape, and the value of guiding individual clinical medication use. For the future study and validation of CIR-score in human PAAD tissues, it is possible to assess the expression of CIR-score hub lncRNAs in biopsied tumor tissues of PAAD patients, and then bring the expression results into the CIR-score formula constructed in this study to assess the risk level of patients. The risk level of each patient can also be used as a reference index for the assessment of the patient’s immunotherapy sensitivity. In addition, the Nomogram constructed in this study can be used to individualize the survival rate of each patient based on the routine clinical risk factors collected from the patients and the CIR-score calculated in the clinic.

5 Conclusion

In conclusion, the construction and validation of a Cuproptosis-associated PAAD lncRNAs signature based on single-cell sequencing carried out in this study not only reveals the organizational principles that shape the TIME complexity of PAAD but also provides certain ideas for immunotherapy-based combination treatment strategies.

Data availability statement

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

Ethics statement

The animal study was approved by Animal Ethics Committee of Jiangsu University. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

YF provided direction and guidance throughout the preparation of this manuscript. YS wrote and edited the manuscript. ZG, CM, and LY collected and prepared the related papers. RH and ZG reviewed and made significant revisions to the manuscript. All authors contributed to the article and approved the submitted version.

Funding

China International Medical Exchange Fund(Z-2014-06-2104), Zhenjiang 169 Talent Project Scientific Research Project (2021-169-XX-19), Zhenjiang Key Research and Development Fund (SH2022032), and Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX23_3766).

Acknowledgments

The authors sincerely appreciate the constructive suggestions and comments of the reviewers on this article and the clinical samples collected by TCGA and GEO databases.

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

References

1. Stoffel EM, Brand RE, Goggins M. Pancreatic cancer: changing epidemiology and new approaches to risk assessment, early detection, and prevention. Gastroenterology (2023) 164:752–65. doi: 10.1053/j.gastro.2023.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Cabasag CJ, Ferlay J, Laversanne M, Vignat J, Weber A, Soerjomataram I, et al. Pancreatic cancer: an increasing global public health concern. Gut (2022) 71(8):1686–7. doi: 10.1136/gutjnl-2021-326311

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Clark CE, Hingorani SR, Mick R, Combs C, Tuveson DA, Vonderheide RH. Dynamics of the immune reaction to pancreatic cancer from inception to invasion. Cancer Res (2007) 67(19):9518–27. doi: 10.1158/0008-5472.CAN-07-0175

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Hingorani SR. Epithelial and stromal co-evolution and complicity in pancreatic cancer. Nat Rev Cancer (2023) 23(2):57–77. doi: 10.1038/s41568-022-00530-w

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Baslan T, Hicks J. Unravelling biology and shifting paradigms in cancer with single-cell sequencing. Nat Rev Cancer (2017) 17(9):557–69. doi: 10.1038/nrc.2017.58

PubMed Abstract | CrossRef Full Text | Google Scholar

6. He J, Shen J, Luo W, Han Z, Xie F, Pang T, et al. Research progress on application of single-cell tcr/bcr sequencing technology to the tumor immune microenvironment, autoimmune diseases, and infectious diseases. Front Immunol (2022) 13:969808. doi: 10.3389/fimmu.2022.969808

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Vegliante R, Pastushenko I, Blanpain C. Deciphering functional tumor states at single-cell resolution. EMBO J (2022) 41(2):e109221. doi: 10.15252/embj.2021109221

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Chen L, Min J, Wang F. Copper homeostasis and cuproptosis in health and disease. Signal Transduction Targeted Ther (2022) 7(1):378. doi: 10.1038/s41392-022-01229-y

CrossRef Full Text | Google Scholar

9. Tong X, Tang R, Xiao M, Xu J, Wang W, Zhang B, et al. Targeting cell death pathways for cancer therapy: recent developments in necroptosis, pyroptosis, ferroptosis, and cuproptosis research. J Hematol Oncol (2022) 15(1):174. doi: 10.1186/s13045-022-01392-3

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M, et al. Copper induces cell death by targeting lipoylated tca cycle proteins. Science (2022) 375(6586):1254–61. doi: 10.1126/science.abf0529

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Lener MR, Scott RJ, Wiechowska-Kozłowska A, Serrano-Fernández P, Baszuk P, Jaworska-Bieniek K, et al. Serum concentrations of selenium and copper in patients diagnosed with pancreatic cancer. Cancer Res Treat (2016) 48(3):1056–64. doi: 10.4143/crt.2015.282

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ishida S, Andreux P, Poitry-Yamate C, Auwerx J, Hanahan D. Bioavailable copper modulates oxidative phosphorylation and growth of tumors. Proc Natl Acad Sci U.S.A. (2013) 110(48):19507–12. doi: 10.1073/pnas.1318431110

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Sharma GG, Okada Y, Von Hoff D, Goel A. Non-coding rna biomarkers in pancreatic ductal adenocarcinoma. Semin In Cancer Biol (2021) 75:153–68. doi: 10.1016/j.semcancer.2020.10.001

CrossRef Full Text | Google Scholar

14. Barczak W, Carr SM, Liu G, Munro S, Nicastri A, Lee LN, et al. Long non-coding rna-derived peptides are immunogenic and drive a potent anti-tumour response. Nat Commun (2023) 14(1):1078. doi: 10.1038/s41467-023-36826-0

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Guttman M, Rinn JL. Modular regulatory principles of large non-coding rnas. Nature (2012) 482(7385):339–46. doi: 10.1038/nature10887

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Li C-H, Xiao Z, Tong JH-M, To K-F, Fang X, Cheng AS, et al. Ezh2 coupled with hotair to silence microrna-34a by the induction of heterochromatin formation in human pancreatic ductal adenocarcinoma. Int J Cancer (2017) 140(1):120–9. doi: 10.1002/ijc.30414

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Sun YW, Chen YF, Li J, Huo YM, Liu DJ, Hua R, et al. A novel long non-coding rna enst00000480739 suppresses tumour cell invasion by regulating os-9 and hif-1α in pancreatic ductal adenocarcinoma. Br J Cancer (2014) 111(11):2131–41. doi: 10.1038/bjc.2014.520

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhao L, Kong H, Sun H, Chen Z, Chen B, Zhou M. Lncrna-pvt1 promotes pancreatic cancer cells proliferation and migration through acting as a molecular sponge to regulate mir-448. J Cell Physiol (2018) 233(5):4044–55. doi: 10.1002/jcp.26072

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Cheng Y, Imanirad P, Jutooru I, Hedrick E, Jin U-H, Rodrigues Hoffman A, et al. Role of metastasis-associated lung adenocarcinoma transcript-1 (Malat-1) in pancreatic cancer. PloS One (2018) 13(2):e0192264. doi: 10.1371/journal.pone.0192264

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Li Z, Zhao X, Zhou Y, Liu Y, Zhou Q, Ye H, et al. The long non-coding rna hottip promotes progression and gemcitabine resistance by regulating hoxa13 in pancreatic cancer. J Transl Med (2015) 13:84. doi: 10.1186/s12967-015-0442-z

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Liu Y, Shi M, He X, Cao Y, Liu P, Li F, et al. Lncrna-pacerr induces pro-tumour macrophages via interacting with mir-671-3p and M6a-reader igf2bp2 in pancreatic ductal adenocarcinoma. J Hematol Oncol (2022) 15(1):52. doi: 10.1186/s13045-022-01272-w

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Yao J, Gao R, Luo M, Li D, Guo L, Yu Z, et al. Exosomal linc00460/mir-503-5p/anln positive feedback loop aggravates pancreatic cancer progression through regulating T cell-mediated cytotoxicity and pd-1 checkpoint. Cancer Cell Int (2022) 22(1):390. doi: 10.1186/s12935-022-02741-5

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Chen H, Yu Y, Zhou L, Chen J, Li Z, Tan X. Cuproptosis-related lncrnas signature as biomarker of prognosis and immune infiltration in pancreatic cancer. Front Genet (2023) 14:1049454. doi: 10.3389/fgene.2023.1049454

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Ge EJ, Bush AI, Casini A, Cobine PA, Cross JR, DeNicola GM, et al. Connecting copper and cancer: from transition metal signalling to metalloplasia. Nat Rev Cancer (2022) 22(2):102–13. doi: 10.1038/s41568-021-00417-2

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cai Z, Ye He, Yu Z, Hu J, Xiao Z, Zu X, et al. Cuproptosis-related modification patterns depict the tumor microenvironment, precision immunotherapy, and prognosis of kidney renal clear cell carcinoma. Front Immunol (2022) 13:933241. doi: 10.3389/fimmu.2022.933241

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Ru B, Wong CN, Tong Y, Zhong JY, Zhong SSW, Wu WC, et al. Tisidb: an integrated repository portal for tumor-immune system interactions. Bioinformatics (2019) 35(20):4200–2. doi: 10.1093/bioinformatics/btz210

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Kamisawa T, Wood LD, Itoi T, Takaori K. Pancreatic cancer. Lancet (2016) 388(10039):73–85. doi: 10.1016/S0140-6736(16)00141-0

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Padrón LJ, Maurer DM, O'Hara MH, O'Reilly EM, Wolff RA, Wainberg ZA, et al. Sotigalimab and/or nivolumab with chemotherapy in first-line metastatic pancreatic cancer: clinical and immunologic analyses from the randomized phase 2 prince trial. Nat Med (2022) 28(6):1167–77. doi: 10.1038/s41591-022-01829-9

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang H, Duan X-R, Xing L-Y, Jia Y-X, Zhou J-Y, Ma J-J. Exploring novel systemic therapies for pancreatic cancer: A review of emerging anti-pd-1/pd-L1 combination therapy. Neoplasma (2022) 69(5). doi: 10.4149/neo_2022_220306N246

CrossRef Full Text | Google Scholar

30. Lu S, Li L, Luo Y, Zhang L, Wu G, Chen Z, et al. A multicenter, open-label, randomized phase ii controlled study of rh-endostatin (Endostar) in combination with chemotherapy in previously untreated extensive-stage small-cell lung cancer. J Thorac Oncol (2015) 10(1):206–11. doi: 10.1097/JTO.0000000000000343

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang L, Ye B, Chen Z, Chen Z-S. Progress in the studies on the molecular mechanisms associated with multidrug resistance in cancers. Acta Pharm Sin B (2023) 13(3):982–97. doi: 10.1016/j.apsb.2022.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yaman M, Kaya G, Yekeler H. Distribution of trace metal concentrations in paired cancerous and non-cancerous human stomach tissues. World J Gastroenterol (2007) 13(4):612–8. doi: 10.3748/wjg.v13.i4.612

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Fang A-P, Chen P-Y, Wang X-Y, Liu Z-Y, Zhang D-M, Luo Y, et al. Serum copper and zinc levels at diagnosis and hepatocellular carcinoma survival in the guangdong liver cancer cohort. Int J Cancer (2019) 144(11):2823–32. doi: 10.1002/ijc.31991

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Basu S, Singh MK, Singh TB, Bhartiya SK, Singh SP, Shukla VK. Heavy and trace metals in carcinoma of the gallbladder. World J Surg (2013) 37(11):2641–6. doi: 10.1007/s00268-013-2164-9

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Wang W, Wang X, Luo J, Chen X, Ma K, He H, et al. Serum copper level and the copper-to-zinc ratio could be useful in the prediction of lung cancer and its prognosis: A case-control study in northeast China. Nutr Cancer (2021) 73(10):1908–15. doi: 10.1080/01635581.2020.1817957

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Kosova F, Cetin B, Akinci M, Aslan S, Seki A, Pirhan Y, et al. Serum copper levels in benign and Malignant thyroid diseases. Bratisl Lek Listy (2012) 113(12):718–20. doi: 10.4149/BLL_2012_162

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Safi R, Nelson ER, Chitneni SK, Franz KJ, George DJ, Zalutsky MR, et al. Copper signaling axis as a target for prostate cancer therapeutics. Cancer Res (2014) 74(20):5819–31. doi: 10.1158/0008-5472.CAN-13-3527

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Choi R, Kim M-J, Sohn I, Kim S, Kim I, Ryu JM, et al. Serum trace elements and their associations with breast cancer subgroups in Korean breast cancer patients. Nutrients (2018) 11(1). doi: 10.3390/nu11010037

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Crowe A, Jackaman C, Beddoes KM, Ricciardo B, Nelson DJ. Rapid copper acquisition by developing murine mesothelioma: decreasing bioavailable copper slows tumor growth, normalizes vessels and promotes T cell infiltration. PloS One (2013) 8(8):e73684. doi: 10.1371/journal.pone.0073684

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhou B, Guo L, Zhang B, Liu S, Zhang K, Yan J, et al. Disulfiram combined with copper induces immunosuppression via pd-L1 stabilization in hepatocellular carcinoma. Am J Cancer Res (2019) 9(11):2442–55.

PubMed Abstract | Google Scholar

41. Liu YL, Bager CL, Willumsen N, Ramchandani D, Kornhauser N, Ling L, et al. Tetrathiomolybdate (Tm)-associated copper depletion influences collagen remodeling and immune response in the pre-metastatic niche of breast cancer. NPJ Breast Cancer (2021) 7(1):108. doi: 10.1038/s41523-021-00313-w

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Yeh DW, Chen YS, Lai CY, Liu YL, Lu CH, Lo JF, et al. Downregulation of commd1 by mir-205 promotes a positive feedback loop for amplifying inflammatory- and stemness-associated properties of cancer cells. Cell Death Differ (2016) 23(5):841–52. doi: 10.1038/cdd.2015.147

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Buscail L, Bournet B, Cordelier P. Role of oncogenic kras in the diagnosis, prognosis and treatment of pancreatic cancer. Nat Rev Gastroenterol Hepatol (2020) 17(3):153–68. doi: 10.1038/s41575-019-0245-4

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Kodydkova J, Vavrova L, Stankova B, Macasek J, Krechler T, Zak A. Antioxidant status and oxidative stress markers in pancreatic cancer and chronic pancreatitis. Pancreas (2013) 42(4):614–21. doi: 10.1097/MPA.0b013e318288360a

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Doñate F, Juarez JC, Burnett ME, Manuia MM, Guan X, Shaw DE, et al. Identification of biomarkers for the antiangiogenic and antitumour activity of the superoxide dismutase 1 (Sod1) inhibitor tetrathiomolybdate (Atn-224). Br J Cancer (2008) 98(4):776–83. doi: 10.1038/sj.bjc.6604226

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Faustova M, Nikolskaya E, Sokol M, Zabolotsky A, Mollaev M, Zhunina O, et al. High-effective reactive oxygen species inducer based on mn-tetraphenylporphyrin loaded plga nanoparticles in binary catalyst therapy. Free Radic Biol Med (2019) 143:522–33. doi: 10.1016/j.freeradbiomed.2019.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Li X, Chen Y, Zhao J, Shi J, Wang M, Qiu S, et al. The specific inhibition of sod1 selectively promotes apoptosis of cancer cells via regulation of the ros signaling network. Oxid Med Cell Longev (2019) 2019:9706792. doi: 10.1155/2019/9706792

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Zhou L, Tan X, Kamohara H, Wang W, Wang B, Liu J, et al. Mek1 and mek2 isoforms regulate distinct functions in pancreatic cancer cells. Oncol Rep (2010) 24(1):251–5. doi: 10.3892/or_00000853

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Protti MP, De Monte L. Cross-talk within the tumor microenvironment mediates th2-type inflammation in pancreatic cancer. Oncoimmunology (2012) 1(1):89–91. doi: 10.4161/onci.1.1.17939

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Bai X, Quek C. Unravelling tumour microenvironment in melanoma at single-cell level and challenges to checkpoint immunotherapy. Genes (2022) 13(10). doi: 10.3390/genes13101757

CrossRef Full Text | Google Scholar

51. Zhang J, Song C, Tian Y, Yang X. Single-cell rna sequencing in lung cancer: revealing phenotype shaping of stromal cells in the microenvironment. Front Immunol (2021) 12:802080. doi: 10.3389/fimmu.2021.802080

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Hwang J, Jin J, Jeon S, Moon SH, Park MY, Yum D-Y, et al. Sod1 suppresses pro-inflammatory immune responses by protecting against oxidative stress in colitis. Redox Biol (2020) 37:101760. doi: 10.1016/j.redox.2020.101760

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Bai Y, Zhang Q, Liu F, Quan J. A novel cuproptosis-related lncrna signature predicts the prognosis and immune landscape in bladder cancer. . Front Immunol (2022) 13:1027449. doi: 10.3389/fimmu.2022.1027449

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Sun Q, Qin X, Zhao J, Gao T, Xu Y, Chen G, et al. Cuproptosis-related lncrna signatures as a prognostic model for head and neck squamous cell carcinoma. Apoptosis (2023) 28(1-2):247–62. doi: 10.1007/s10495-022-01790-5

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Song X, Hou L, Zhao Y, Guan Q, Li Z. Metal-dependent programmed cell death-related lncrna prognostic signatures and natural drug sensitivity prediction for gastric cancer. Front Pharmacol (2022) 13:1039499. doi: 10.3389/fphar.2022.1039499

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Zhou H, Gao Y, Li X, Shang S, Wang P, Zhi H, et al. Identifying and characterizing lincrna genomic clusters reveals its cooperative functions in human cancer. J Transl Med (2021) 19(1):509. doi: 10.1186/s12967-021-03179-5

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Yılmaz Susluer S, Kayabasi C, Ozmen Yelken B, Asik A, Celik D, Balci Okcanoglu T, et al. Analysis of long non-coding rna (Lncrna) expression in hepatitis B patients. Bosn J Basic Med Sci (2018) 18(2):150–61. doi: 10.17305/bjbms.2018.2800

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Hong Y, Choi H-M, Cheong HS, Shin HD, Choi CM, Kim WJ. Epigenome-wide association analysis of differentially methylated signals in blood samples of patients with non-small-cell lung cancer. J Clin Med (2019) 8(9). doi: 10.3390/jcm8091307

CrossRef Full Text | Google Scholar

59. Cao P-W, Liu L, Li Z-H, Cao F, Liu F-B. Prognostic value of drug targets predicted using deep bioinformatic analysis of M6a-associated lncrna-based pancreatic cancer model characteristics and its tumour microenvironment. Front Genet (2022) 13:853471. doi: 10.3389/fgene.2022.853471

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Wang Q, Zhang C, Cao S, Zhao H, Jiang R, Li Y. Tumor-derived exosomes orchestrate the microrna-128-3p/elf4/cdx2 axis to facilitate the growth and metastasis of gastric cancer via delivery of linc01091. Cell Biol Toxicol (2022) 39:519–36. doi: 10.1007/s10565-022-09728-y

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Nikas JB, Nikas EG. Genome-wide DNA methylation model for the diagnosis of prostate cancer. ACS Omega (2019) 4(12):14895–901. doi: 10.1021/acsomega.9b01613

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Cao G, Chang Y, Yang G, Jiang Y, Han K. A novel risk score model based on four angiogenesis long non-coding rnas for prognosis evaluation of pancreatic adenocarcinoma. Aging (Albany NY) (2022) 14(22):9090–102. doi: 10.18632/aging.204387

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Zhang S, Li X, Tang C, Kuang W. Inflammation-related long non-coding rna signature predicts the prognosis of gastric carcinoma. Front Genet (2021) 12:736766. doi: 10.3389/fgene.2021.736766

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Hosein AN, Brekken RA, Maitra A. Pancreatic cancer stroma: an update on therapeutic targeting strategies. Nat Rev Gastroenterol Hepatol (2020) 17(8):487–505. doi: 10.1038/s41575-020-0300-1

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Xiang X, Wang J, Lu D, Xu X. Targeting tumor-associated macrophages to synergize tumor immunotherapy. Signal Transduction Targeted Ther (2021) 6(1):75. doi: 10.1038/s41392-021-00484-9

CrossRef Full Text | Google Scholar

66. Tong F, Mao X, Zhang S, Xie H, Yan B, Wang B, et al. Hpv + Hnscc-derived exosomal mir-9 induces macrophage M1 polarization and increases tumor radiosensitivity. Cancer Lett (2020) 478:34–44. doi: 10.1016/j.canlet.2020.02.037

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Li Z, Fu W-J, Chen X-Q, Wang S, Deng R-S, Tang X-P, et al. Autophagy-based unconventional secretion of hmgb1 in glioblastoma promotes chemosensitivity to temozolomide through macrophage M1-like polarization. J Exp Clin Cancer Res (2022) 41(1):74. doi: 10.1186/s13046-022-02291-8

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Kurahara H, Shinchi H, Mataki Y, Maemura K, Noma H, Kubo F, et al. Significance of M2-polarized tumor-associated macrophage in pancreatic cancer. J Surg Res (2011) 167(2):e211–e9. doi: 10.1016/j.jss.2009.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Ubil E, Caskey L, Holtzhausen A, Hunter D, Story C, Earp HS. Tumor-secreted pros1 inhibits macrophage M1 polarization to reduce antitumor immune response. J Clin Invest (2018) 128(6):2356–69. doi: 10.1172/JCI97354

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Guiducci C, Vicari AP, Sangaletti S, Trinchieri G, Colombo MP. Redirecting in vivo elicited tumor infiltrating macrophages and dendritic cells towards tumor rejection. Cancer Res (2005) 65(8):3437–46. doi: 10.1158/0008-5472.CAN-04-4262

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Zhang Y, Lazarus J, Steele NG, Yan W, Lee H-J, Nwosu ZC, et al. Regulatory T-cell depletion alters the tumor microenvironment and accelerates pancreatic carcinogenesis. Cancer Discovery (2020) 10(3):422–39. doi: 10.1158/2159-8290.CD-19-0958

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Nagaraju GP, Malla RR, Basha R, Motofei IG. Contemporary clinical trials in pancreatic cancer immunotherapy targeting pd-1 and pd-L1. Semin In Cancer Biol (2022) 86(Pt 3):616–21. doi: 10.1016/j.semcancer.2021.11.003

CrossRef Full Text | Google Scholar

73. Chen C, Tang Y, Huang H, Jia L, Feng L, Zhao J, et al. Relieving immunosuppression by endo@Plt targeting anti-angiogenesis to improve the efficacy of immunotherapies. Chem Commun (Camb) (2022) 58(19):3202–5. doi: 10.1039/d2cc00205a

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Shan Y-f, Fang Y-f, Wang X-q, Jin R, Zhang Q-y, Andersson R. Experimental studies on treatment of pancreatic cancer with double-regulated duplicative adenovirus adtphre-hendo carrying human endostatin gene. Pancreatology (2013) 13(4):393–400. doi: 10.1016/j.pan.2013.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Pu X, Wang Q, Liu L, Chen B, Li K, Zhou Y, et al. Rh-endostatin plus camrelizumab and chemotherapy in first-line treatment of advanced non-small cell lung cancer: A multicenter retrospective study. Cancer Med (2022) 12:7724–33. doi: 10.1002/cam4.5526

CrossRef Full Text | Google Scholar