Chao Ma
Feng Li
Zhuoyu Gu- Department of Thoracic Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
Background: Cuproptosis is a newly discovered non-apoptotic form of cell death that may be related to the development of tumors. Nonetheless, the potential role of cuproptosis-related lncRNAs in tumor immunity formation and patient-tailored treatment optimization of lung adenocarcinoma (LUAD) is still unclear.
Methods: RNA sequencing and survival data of LUAD patients were downloaded from The Cancer Genome Atlas (TCGA) database for model training. The patients with LUAD in GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081 were used for validation. The proofed cuproptosis-related genes were extracted from the previous studies. The Pearson correlation was applied to select cuproptosis-related lncRNAs. We chose differentially expressed cuproptosis-related lncRNAs in the tumor and normal tissues and allowed them to go to a Cox regression and a LASSO regression for a lncRNA signature that predicts the LUAD prognosis. Kaplan–Meier estimator, Cox model, ROC, tAUC, PCA, nomogram predictor, decision curve analysis, and real-time PCR were further deployed to confirm the model’s accuracy. We examined this model’s link to other regulated cell death forms. Applying TMB, immune-related signatures, and TIDE demonstrated the immunotherapeutic capabilities of signatures. We evaluated the relationship of our signature to anticancer drug sensitivity. GSEA, immune infiltration analysis, and function experiments further investigated the functional mechanisms of the signature and the role of immune cells in the prognostic power of the signature.
Results: An eight-lncRNA signature (TSPOAP1-AS1, AC107464.3, AC006449.7, LINC00324, COLCA1, HAGLR, MIR4435-2HG, and NKILA) was built and demonstrated owning prognostic power by applied to the validation cohort. Each signature gene was confirmed differentially expressed in the real world by real-time PCR. The eight-lncRNA signature correlated with 2321/3681 (63.05%) apoptosis-related genes, 11/20 (55.00%) necroptosis-related genes, 34/50 (68.00%) pyroptosis-related genes, and 222/380 (58.42%) ferroptosis-related genes. Immunotherapy analysis suggested that our signature may have utility in predicting immunotherapy efficacy in patients with LUAD. Mast cells were identified as key players that support the predicting capacity of the eight-lncRNA signature through the immune infiltrating analysis.
Conclusion: In this study, an eight-lncRNA signature linked to cuproptosis was identified, which may improve LUAD management strategies. This signature may possess the ability to predict the effect of LUAD immunotherapy. In addition, infiltrating mast cells may affect the signature’s prognostic power.
Introduction
Lung cancer is one of the most common malignant tumors affecting human health and quality of life, and it has the highest mortality rate among all cancers (Siegel et al., 2021). Nearly 80%–85% of lung cancer cases are non-small cell lung cancers (NSCLC) (Siegel et al., 2021). The 5-year survival rate for lung cancer in the United States is approximately 20%. Lung adenocarcinoma (LUAD), accounting for nearly 40% of NSCLC, has surpassed lung squamous cell carcinoma (LUSC) in incidence as the most common histological subtype of lung cancer (Siegel et al., 2021). It has been shown that early diagnosis and the advent of new treatments have prolonged the survival of people with LUAD (Spella and Stathopoulos, 2021). Due to LUAD’s insidious onset, patients with this disease have lymph node metastasis at screening, and the postoperative recurrence rate is high, leading to poor outcomes (Spella and Stathopoulos, 2021). Despite significant breakthroughs in diagnostic and remedial technologies, there is still a low long-term survival rate in LUAD patients (Spella and Stathopoulos, 2021). To improve the clinical outcome, new therapeutic targets need to be identified. It is, therefore, urgent and imperative that a new prognostic model for a feasible targeted biomarker be developed.
Programmed cell death (PCD) is crucial for tissue homeostasis and animal development. PCDs, such as apoptosis, necrosis, pyroptosis, and ferroptosis, are demonstrated to be critical for tumorigenesis, progression, and metastasis (Liu et al., 2022). One recently study published in “Science” have revealed cuproptosis as a new type of PCD, distinct from ferroptosis and apoptosis (Tsvetkov et al., 2022). Tsvetkov et al. proposed that cuproptosis requires mitochondrial respiration but that it is less affected by ATP produced by glycolysis (Tsvetkov et al., 2022). Their results show that copper seems to be closely linked to the TCA cycle, implying a close connection between cuproptosis and mitochondrial metabolism (Tsvetkov et al., 2022). It remains unclear how cuproptosis contributes to tumor initiation, progression, and immunity (Chen et al., 2022). It has been demonstrated that cuproptosis-related genes are expressed aberrantly and impact outcomes in clear cell renal cell carcinoma and hepatocellular carcinoma (Bian et al., 2022; Zhang et al., 2022). A recent study examined the significance of cuproptosis-related genes in immune infiltration and melanoma prognosis (Lv et al., 2022). Understanding the heterogeneity of LUAD may be improved by determining the molecular profiles of genes associated with cuproptosis. Long non-coding RNAs (lncRNAs), which are transcripts >200 nucleotides in length but not translated into protein, regulate gene expression at multiple levels and participate in various biological processes, especially cell death mechanisms (Alishahi et al., 2019). Therefore, the relationship between lncRNA and cuproptosis has potential significance for the clinical research of LUAD. Because RCD is more immuno-targeted than others (Liu et al., 2022), we need to understand how cuproptosis is initiated, propagated, and executed in LUAD cells, which have important implications for possible combined diagnostic and therapeutic interventions.
This work aims to fill the gap in studying cuproptosis-related lncRNA and LUAD prognosis. Using experimentally confirmed cuproptosis-related genes, we constructed a lncRNA signature to predict LUAD outcomes and further demonstrated its prognostic value. Moreover, we examined the correlation between the signature and other RCDs, and how each gene was expressed in the real world. More importantly, we investigated this model’s immunotherapeutic potential and anticancer drug selection value. For functional exploration, we performed gene set enrichment analysis (GSEA) to reveal the pathways in which the signatures reside. The final analysis of tumor-infiltrating immune cells provided the basis for further investigation of lncRNA signatures.
Materials and methods
Data selection and preprocessing
The project TCGA-LUAD contains LUAD samples and their accordingly clinical data. We downloaded the RNA-seq data (type: HTSeq - Counts) and clinical phenotype of patients in the project TCGA-LUAD on the GDC Xena Hub (https://gdc.xenahubs.net). All the data available in GDC Xena Hub is sourced exclusively from the official GDC pipeline (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/). It is worth noting that GDC Xena Hub carries out log2 (count+1) processing exclusively. To get more LUAD samples, we searched the Gene Expression Omnibus (GEO) (Clough and Barrett, 2016) database (https://www.ncbi.nlm.nih.gov/geo/) using the keyword “lung adenocarcinoma” and applying the following screening criteria: 1) total RNA expression data is available; 2) total LUAD cases that contain survival data are greater than 80. Finally, datasets named GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081 were selected, and their gene expression matrix and clinical parameters were retrieved from GEO online portal. These GEO datasets are all generated using Platform GPL570, and proper protocols were followed for data collection and generation. Additional details can be found at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570. For the datasets selected above, only the tumor sample that contained gene expression and survival data were collected. Eligible patients in TCGA-LUAD were collected and served as the training cohort, and the downloaded data were directly used in our analysis without preprocessing. For the preprocessing of GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081, we used the R package inSilicoMerging to merge them (Taminau et al., 2012), and then we adopted the method that developed by Johnson et al. (Johnson et al., 2007) to remove the batch effect and finally obtained the data matrix treated as a validation cohort.
Screening of experimental proofed cuproptosis-related genes
We screened cuproptosis-related genes according to the following criteria: 1) Search the PubMed (https://pubmed.ncbi.nlm.nih.gov/) using the keyword “copper induced cell death pathway cancer”; 2) Studies unrelated to the search topic or with only Silico analysis were manually excluded; 3) The genes involved in wet experiment confirmed copper-induced cell death pathways in each study were collected, and unique ones were taken.
Construction and validation of the prognosis model
LncRNA expression data were categorized depending on the annotations provided by the GENCODE project (Derrien et al., 2012). Cuproptosis-related lncRNAs were identified by comparing cuproptosis-related genes and lncRNAs using the limma R package with the following criteria: | Pearson correlation coefficient | > 0.3 and p < 0.05. For the selection of differentially expressed ones form the cuproptosis-related lncRNAs in normal and tumor tissue, we adopted the limma R package and applied false discovery rate <0.05.
Then, these screened lncRNAs were subjected to univariate Cox regression model for choosing the ones with potential prognostic power. The “survival” R package helped to build Cox models, and the cuproptosis-related lncRNA with a p-value less than 0.05 were selected. As a precaution against overfitting, we used the R package “glmnet” to combine Cox proportional hazards regression and the least absolute shrinkage and selection operator (LASSO) for the selected lncRNAs. We used 10-fold cross-validation to pick the optimal lambda (Tibshirani, 1997; Sauerbrei et al., 2007; Friedman et al., 2010; Goeman, 2010). The following formula was used to calculate each patient’s risk score:
with βi indicating the coefficient, Expi representing the relative expression levels of each lncRNA standardized by z-score, and n indexing each hub lncRNA in the signature.
Patients in each cohort and subtypes were divided into low-risk groups and high-risk groups based on the best cutoff defined by the “surv_cutpoint ()” function of the “survminer” R package. The prognostic value of the signature was assessed in LUAD patients in the training and validation cohorts in Kaplan–Meier curve, univariate Cox analysis, and multivariable Cox analysis, the receiver operating characteristic (ROC) curve (Cao and Lopez-de-Ullibarri, 2019), time-dependent AUC (tAUC) analysis, and principal component analysis (PCA). The ROC curves and the tAUC analysis were made possible with the help of the “timeROC” and “survival” R packages. The “scatterplot3D” R package was used to evaluate the distribution of patients with different risk scores by PCA. Moreover, using the R packages “survival”, “survminer,” “rms” and “regplot”, we built nomograms predicting overall survival at 1, 3, and 5 years and applied calibration curves to evaluate model predictions versus actual outcomes consistency between. The R package “ggDCA” carried out the decision curve analysis to assess the nomogram accuracy for LUAD outcomes.
Validation of the signature genes’ expression profile using the real-time PCR
The LUAD cell line A549 and the human lung fibroblasts cell line WI-38 were provided by Jiangsu KeyGEN BioTECH Co., Ltd. (Nanjing, China). A549 cells were grown in 90% RPMI-1640 (Thermo, United States) containing 10% fetal bovine serum (FBS), while WI-38 cells were grown in a 90% DMEM medium (Thermo, United States) containing 10% FBS. In addition, all cell media were supplemented with penicillin (100 U/mL) and streptomycin (100 U/mL) at 37°C in a 5% CO2 incubator. We isolated total RNA with TRIzol reagent (Invitrogen, United States). First-strand cDNA was synthesized from 1 μg total RNA using a PrimeScript™ 1st strand cDNA Synthesis Kit (Takara Bio, Japan). The One Step TB Green™ PrimeScript™ RT-PCR Kit II was used for real-time PCR (Takara Bio, Japan). The cycling conditions were 30 s of polymerase activation at 95°C followed by 40 cycles at 95°C for 5 s and 60°C for 30 s. GAPDH was used as an internal loading control. The relative level was calculated by the relative quantification 2−ΔΔCT method.
GSEA
For GSEA, we obtained the GSEA software (version 4.2.3) from the GSEA website (http://www.gsea-msigdb.org/gsea/downloads.jsp) (Subramanian et al., 2005). Patients were divided into low-risk groups and high-risk groups based on the best cutoff defined by the “surv_cutpoint ()” function of the “survminer” R package. We downloaded “c2.cp.kegg.v7.4.symbols.gmt” from the Molecular Signatures Database (Liberzon et al., 2011) (http://www.gsea-msigdb.org/gsea/downloads.jsp) to equip the GSEA software to evaluate related molecular mechanisms and pathways based on gene expression profiles and groupings. The parameter was as follows: the minimum gene set to 5; the maximum gene set to 5000; one thousand resamplings; p-value < 0.05 and FDR <0.25 were considered statistically significant.
Correlations between the cuproptosis-related lncRNA signature and apoptosis, necroptosis, pyroptosis, and ferroptosis
For better knowing the interactions between our cuproptosis-related lncRNA signature and other types of “cell death”, we adopted the Pearson analysis. Apoptosis, necroptosis, and pyroptosis-related genes were extracted from the GeneCard and GSEA online databases, respectively, by applying the following steps: 1) search the GeneCard using the corresponding keyword; 2) search the GSEA using the corresponding keyword: 3) merge the above results and take the unique genes. FerrDb is the first database dedicated to ferroptosis regulators and ferroptosis-related diseases (Zhou and Bao, 2020). Ferroptosis-related genes were obtained from FerrDb (http://www.datjar.com:40013/bt2104/).
Correlations between the lncRNA signature and immunotherapy
As part of the “gsva” R package, the “single-sample gene set enrichment analysis (ssGSEA)” function was used to assess 13 immune-related pathways. The “maftools” was used to identify the top 20 mutated genes and visualize mutations and their frequencies in all samples of the training cohort. Patients were divided into low-risk groups and high-risk groups based on the best cutoff defined by the “surv_cutpoint ()” function of the “survminer” R package. We applied the chi-square test to assess differences in gene mutation frequencies in low and high risk LUAD populations.
Tumor mutational burden (TMB) is an emerging therapeutic measure of immunotherapy sensitivity. TMB is defined as the frequency of certain mutations within a tumor’s genes (Chalmers et al., 2017). The TMB rank score of each case with LUAD was obtained as previously described (Chalmers et al., 2017). We deployed the Pearson coefficient together with the Wilcoxon rank-sum to calculate the connections between the risk score and TMB. Then, we selected CD274 (Fabrizio et al., 2018), CTLA4 (Rowshanravan et al., 2018), HAVCR2 (Holderried et al., 2019), IDO1 (Platten et al., 2014), LAG3 (Andrews et al., 2017), PDCD1 (Ribas and Wolchok, 2018), CD8A (Raskov et al., 2021), CXCL10 (Shi et al., 2021), CXCL9 (Liang et al., 2021), GZMA (Inoue et al., 2016), GZMB (Hurkmans et al., 2020), IFNG (Jorgovanovic et al., 2020), PRF1 (Fan et al., 2021), TBX2 (Li et al., 2021), and TNF (Freeman et al., 2021) from previous studies as immune relevant signatures, and used the Pearson and Wilcoxon rank-sum analysis to measure the correlations between our signature and the above listed immune signatures. To further test whether our lncRNA signature could guide immunotherapy, we deployed the Kaplan–Meier estimator to test the prognostic sensitivity of individual immune relevant signatures in high and low-risk patients, respectively. Tumor Immune Dysfunction and Exclusion (TIDE) is a computational framework that integrates T cell dysfunction expression signatures and T cell exclusion to model tumor immune evasion. Using TIDE, tumor immune evasion can be modeled in two ways and can be used to predict immunotherapy outcomes (Jiang et al., 2018; Fu et al., 2020; Chen Y. et al., 2021). More importantly, we tested if our signature correlated with the TIDE.
PDCD1 and PRF1 expression detection and colony-formation assays
A549 cells were stably infected with a virus carrying either MIR4435-2HG, NKILA and their respective controls and seeded in 6-well plates (2 × 105 cells/well) containing culture medium with 10% FBS for indicated time points.
48 h after transfection with lentivirus, the expressions of PDCD1 and PRF1 in the MIR4435-2HG and NKILA overexpression groups and their control groups were detected by real-time PCR. PDCD1 primer sequence is as follows: Forward: GTGTCACACAACTGCCCAAC, Reverse: CTGCCCTTCTCTCTGTCACC. PRF1 primer sequence is as follows: Forward: ACTCACAGGCAGCCAACTTT, Reverse: GGGTGCCGTAGTTGGAGATA.
For the colony-formation assay, MIR4435-2HG-overexpressing and NKILA-overexpressing A549 cells (1000 cells/well) were seeded in 6-well plates and cultured under standard conditions. After 7∼10 days, cells were fixed with methanol and stained with 1% crystal violet. The number of colonies was manually counted using free ImageJ software (version 1.53e, National Institutes of Health, Bethesda, MD, USA).
Drug sensitivity and the signature
Half-maximal inhibitory concentration (IC50) data for anticancer drugs were obtained from the Genomics of Cancer Drug Sensitivity (GDSC) database (https://www.cancerrxgene.org/) (Geeleher et al., 2014a). The drug distribution in the high and low risk score groups was analyzed and visualized using the “pRRophetic” R package (Geeleher et al., 2014b). The difference and correlation in the IC50 between the different risk groups was tested by the Wilcoxon test and the Spearman’s rank correlation coefficient, respectively.
Relationships between the TICs and our signature
TICs play an essential role in the development and progression of cancer (Ye et al., 2019). We invoked the “CIBERSORT” R package developed by Newman and colleagues (Newman et al., 2015) to assess the abundance of 22 TICs using gene expression data from the training cohort. To better illustrate a clear portrait of the lncRNA signature, we first used the Pearson test and Wilcoxon rank-sum analysis to assess the signatures’ correlations with 22 TICs. Next, we applied the univariate Cox model and Kaplan-Meier curve to determine the TICs that had prognostic power. Finally, combining the above results, we inferred candidate TICs may contribute to the prognostic power of the lncRNA signature.
Statistical analysis
We employed a Wilcox test to compare gene expression levels between tumor and non-tumor tissues. To evaluate the predictive power of our survival model for patient survival, we utilized the tAUC. The overall survival rates of each group were compared using a Kaplan-Meier analysis in conjunction with a log-rank test. To identify molecules associated with prognosis, Cox regression of overall survival was performed using univariate data. We analyzed all statistics using R software, and we performed two-sided t-tests for all statistical tests unless otherwise stated. Statistical significance was defined as p < 0.05.
Results
Patient characteristics
A total of 1054 LUAD cases with available survival data were included in the present study. The flowchart of this study is shown in Figure 1. We included 500 patients from the TCGA project for the risk model construction. We downloaded five validation datasets (GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081) from GEO. According to our selection criteria, data from the five GEO datasets contained 554 pieces of usable LUAD information and was used to validate the gene signature we built. The patient information used in this study is listed in Table 1. We merged five validation datasets into one cohort to better process our study. Batch effects are nonbiological differences between two or more datasets. We removed the bias caused by batch effects and made the transcriptional profiles more similar based on our methods stated before (Figures 2A, B).
FIGURE 1. Research design and analysis process. TCGA, The Cancer Genome Atlas; LUAD, lung adenocarcinoma; HR, hazard ratio; ROC, receiver operating characteristic; TP, true positive rate; FP, false positive rate; AUC, area under the ROC curve; PC, principal component; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; IC50, half maximal inhibitory concentration; TICs, tumor-infiltrating immune cells.
TABLE 1. Clinical characteristics of patients enrolled in the study.
FIGURE 2. Data preprocessing and risk model development. (A) UMAP plot showing validation cohort merged before batch effect removing, displaying that the samples of each dataset are separated from each other. (B) UMAP plot showing validation cohort merged after batch effect removing, displaying that samples of each dataset are clustered and intertwined with each other. (C) The heatmap for identifying differentially expressed lncRNAs in normal and tumor tissues. (D) The volcano plot for identifying differentially expressed lncRNAs in normal and tumor tissues. (E) LASSO coefficient profiles of prognostic lncRNAs enrolled. (F) LASSO regression with ten-fold cross-validation obtained eight prognostic genes using the minimum Lambda (0.006898922). UMAP, Uniform Manifold Approximation and Projection; LASSO, least absolute shrinkage and selection operator; FDR, false discovery rate.
A cuproptosis-related eight-lncRNA signature generated
We got 60 cuproptosis-related genes following our criteria, shown in Table 2. We used the Pearson algorithm to detect lncRNAs that significantly associated with cuproptosis-related genes in the training cohort samples, and finally 8277 unique lncRNAs were determined (Supplementary Table S1; Supplementary Figure S1). From the 8277 lncRNAs, the “limma” R language package identified 3261 of them differentially expressed in normal and tumor tissues (Figures 2C, D). Then, we used the expression profiles of the 3261 lncRNAs in the training cohort to assess for candidates with prognostic ability. And 9 lncRNAs were identified having potential prognostic ability (Supplementary Table S2). By utilizing the “glmnet” R package (version 4.1.2, New Zealand), LASSO Cox regression analysis was carried out to narrow the scope of applicant genes and develop a predictive model. Finally, the optimal lambda was selected as 0.006898922, at which point 8 lncRNAs and their coefficients were preserved (Figures 2E, F; Table 3). In addition, we displayed the correlations between each of the 8 lncRNAs and the 60 cuproptosis-related genes, which are shown in Supplementary Figure S2A.
TABLE 2. The cuproptosis-related genes that extracted from the previous research.
TABLE 3. Prognostic lncRNAs obtained from the LASSO Cox regression model and their primer sequences.
The eight-lncRNA signature proven to have stable prognostic capacity
In Supplementary Figures S2B, C, the upper parts showed the patients sorted by increasing risk score, the scatter plot in the middle showed the survival status of the LUADs, and the heatmap in the lower part showed the relative expression levels of eight hub lncRNAs. Based on a risk score calculator we developed, a risk score was calculated for each LUAD case enrolled. Kaplan-Meier estimator suggested that high risk patients had a worse survival in the TCGA cohort (log-rank test, p < 0.001, Figure 3A) and in the validaiton cohort (log-rank test, p < 0.001, Figure 3B) than low risk patients. In addition, the tumor stage subtype analysis suggested that our risk score can precisely predict the patients’ outcomes, where the higher risk score forebode worse disease ends (Figures 3A, B). In addition, in Supplementary Figure S3A, we displayed each eight lncRNA’s prognosis ability in the form of Kaplan-Meier curves using the two cohorts’ data, showing that the MIR4435-2HG and NKILA performed stable unfavorable impact on LUAD patients, while TSPOAP1-AS1, AC107464.3, AC006449.7, LINC00324, COLCA1, and HAGLR helped the prognosis improvement of LUADs.
FIGURE 3. Validation of the eight-lncRNA signature using the involved studied cohorts and real-time PCR experiment. (A,B) Kaplan-Meier analysis was performed in the whole and subtypes of training and validation cohorts. Patients in each cohort and subtypes were divided into low-risk groups and high-risk groups based on the best cutoff defined by the “surv_cutpoint ()” function of the “survminer” R package. The log-rank test with a p-value < 0.05 suggests the survival difference is significant. The bottom part displays the number of patients at risk. (C) Univariate and multivariable Cox proportional hazards models. *: the variables involved in the studied cohorts, explains as follows: Gender: male vs. female; Race: white vs. non-white; Ethnicity: Hispanic or Latino vs. non-Hispanic or Latino; Prior malignancy: yes vs. no; Tumor origin: upper lobe lung vs. non-upper lobe lung; Smoking history: ever vs. never. (D) ROC curves. The ROC curves valued the accuracy for LUAD outcome prediction of our signature at 1-, 3-, and 5-year, respectively. (E) The tAUC analyses compared our signature’s prognostic ability with other available clinical characteristics. The larger the AUC, the stronger the model’s predictive ability. (F) Principal component analysis scatter plot. (G) The nomogram, a quantitative model for predicting clinical prognosis, to predict 1-year, 3-year, and 5-year OS in the LUAD patients of the TCGA-LUAD cohort using seven factors, including age, grade, tumor stage, risk score, smoking history, tissue origin, and prior malignancy. (H) The calibration curves indicated that the nomogram accurately predicted the 1-, 3-, and 5-year OS of LUAD patients in the TCGA cohort. (I) Decision curve analysis. The curves assessed the nomogram prediction accuracy for LUAD outcomes at 1-, 3-, and 5-year, respectively. (J) The expression levels of the 8 signature genes in A549 (n = 9) and WI-38 (n = 9) cell lines detected by real-time PCR. Data were means ± STD. ****: p-value < 0.0001; HR, hazard ratio; L95%, 95% confidence interval lower; H95%, 95% confidence interval higher; HR, hazard ratio; ROC, receiver operating characteristic; AUC, area under the ROC curve; tAUC, time-dependent AUC; TCGA, The Cancer Genome Atlas; LUAD, lung adenocarcinoma; PCA, Principal components analysis; OS, overall survival; Exp, relative expression.
Subsequently, we assessed the independent prognostic value of our gene signature. Univariate and multivariate Cox regression analyses were performed to evaluate whether signature risk score was independent of other clinical parameters (age, gender, race, ethnicity, tumor stage, tumor origin, etc.) as prognostic factors for LUAD patients (Figure 3C). Univariate Cox regression analysis showed that in training and validation cohorts, risk scores were significantly correlated with overall survival. Multivariate Cox regression analysis showed that risk scores were independently associated with survival in the training cohort (HR = 2.93, 95% CI: 1.92–4.47, p = 6.45e-07) and validation cohort (HR = 2.64, 95% CI: 1.51–4.64, p = 7.14e-04). In addition, the age of validation cohort also showed independent prognostic value, however, its significance did not show a consistence in the two cohorts. Moreover, each gene’s univariate Cox regression was shown in the chart exhibited in Supplementary Figure S3B. These results further confirmed the high predictive accuracy of our gene signature, suggesting that the model could be independently used to predict the prognosis of LUAD patients.
Then we applied ROC analysis (Figure 3D) and time-dependent AUC (Figure 3E) to verify the prognostic ability of our lncRNA signature. The ROC curve for assessing the predictive strength of the gene prognostic signature exhibited an AUC of 0.720, 0.627, and 0.635 at 1, 3, and 5 years, respectively, suggesting that the risk signature would have a tolerable predictive capacity in the training database. Concurrently, this signature displayed similar results for the validation cohort in ROC curve analysis at 1, 3, and 5 years, with AUCs of 0.583, 0.621, and 0.612, respectively. According to the time-dependent AUC performed in the training cohort (Figure 3E), our risk score was close compared to tumor stage, which is regarded as the prognosis gold standard. The AUC was generally greater than 0.7 when we combined the risk score and tumor stage. Consistently, the AUC performance of the risk score combined with the tumor stage outperformed all remaining factors at all time points, suggesting that our risk score is an excellent complement to the tumor stage (Figure 3E). PCA results showed significant heterogeneity between high-risk and low-risk patients in the training and validation cohorts (Figure 3F), which certified the superior discrimination of the risk score model. Moreover, we constructed a nomogram using seven factors, including age, grade, tumor, etc., a quantitative model for predicting clinical outcomes, predicting 1-, 3-, and 5-year overall survival in the TCGA-LUAD cohort population (Figure 3G). The calibration curves confirmed that the nomogram accurately predicted the 1-, 3-, and 5-year overall survival of LUAD patients in the TCGA cohort (Figure 3H). Decision curve analysis assessed the nomogram accuracy for LUAD outcomes at 1-, 3-, and 5-year, respectively, and the results confirmed that the nomogram is superior to other factors (Figure 3I).
The expression level of each gene of the signature was validated by real-time PCR
To better understand the real-world expression pattern of each gene in the gene signature, we applied the real-time PCR to detect the above genes in human LUAD cell lines (n = 9) and human normal cell lines (n = 9) difference in expression. Table 3 shows the primer sequences for the 8 genes, TSPOAP1-AS1, AC107464.3, AC006449.7, LINC00324, COLCA1, HAGLR, MIR4435-2HG, and NKILA. Notably, all genes were differentially expressed in tumor and normal cell lines (Figure 3J). Specifically, genes of MIR4435-2HG and NKILA were highly expressed in LUAD cell lines, while the remaining were under-expressed in LUAD cell lines. The upregulated genes in LUAD like MIR4435-2HG and NKILA were also showed having unfavorable prognosis powerful in Supplementary Figure S3, and the downregulated genes in LUAD like TSPOAP1-AS1, AC107464.3, AC006449.7, LINC00324, COLCA1, and HAGLR displayed owning protectable function for the LUAD outcomes, which further proved the validity of the gene signatures we found and provided clues for deeper research.
GSEA determined the mechanisms of the prognosis signature
The GSEA analysis performed based on risk scores discovered significantly enriched gene sets in the lncRNA signature, which the top ten items were primarily related to allograft rejection, fatty acid metabolism, inositol phosphate metabolism, arachidonic acid metabolism, PPAR signaling pathway, intestinal immune network for IgA production, vascular smooth muscle contraction, alpha-Linolenic acid metabolism, and non-homologous end-joining (Figure 4; Supplementary Table S3).
FIGURE 4. GSEA analysis with the KEGG gene set as the background identified relevant pathways of our signature. The significance threshold of this analysis was set as: p-value < 0.05, and FDR <0.25. GSEA, Gene Set Enrichment Analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate.
The eight-lncRNA signature’s relationships with apoptosis, necroptosis, pyroptosis, and ferroptosis
We got apoptosis, necroptosis, pyroptosis, and ferroptosis genes following our criteria, shown in Supplementary Table S4. The Pearson coefficient examined the relationships between our prognosis model and apoptosis, necroptosis, pyroptosis, and ferroptosis-related genes, respectively (Supplementary Table S4; Supplementary Figure S4). The analysis showed that SELENBP1, METTL7A, ANLN, SARM1, GGA2, GAPDH, MOAP1, CTSH, RBP5, and CEP55 were the top ten correlated apoptosis-related genes, and overall, 2321/3681 (63.05%) genes significantly linked with the lncRNA signature. RIPK3, FADD, TRPM7, CYLD, TLR3, IPMK, ZBP1, TNF, FASLG, and PGLYRP1 were discovered as the top necroptosis-related that correlated with the eight-lncRNA signature. As a whole, there were 11/20 (55.00%) necroptosis-related genes correlated with the signature pronouncedly. Moreover, the Pearson test found the top pyroptosis-related genes that correlated with our signature are NLRP1, CARD8, GSDMB, CYCS, IRF2, GSDMD, IL1A, NLRP6, CASP1, and CHMP2B. In total, 34/50 (68.00%) pyroptosis-related genes correlated with our signature. The examination found SLC2A1, RRM2, AURKA, MDM4, PLA2G6, DUOX1, GLS2, NRAS, HIF1A, and SIRT3 were top ferroptosis-related genes that correlated with our signature. To sum up, 222/380 (58.42%) ferroptosis-related genes correlated with our signature.
The role of risk score participating in immunotherapy
As for the related immune functions, the scores for the CCR, Check−point, Cytolytic_activity, HLA, Inflammation−promoting, T_cell_co−inhibition, T_cell_co−stimulation, and Type_II_IFN_Reponse were lower in the high-risk group than in the low-risk group (Figure 5A). These findings revealed that the eight-lncRNA signature might be associated with the immunological status of LUADs. We thoroughly explored the mutation characteristics of all tumor samples in the TCGA-LUAD cohort. As top 20 mutated shown in Figure 5B, gene TP53 mutated most frequently approximately accounting for 53.3% in the cohort, followed by TTN (50.6%) and MUC16 (43.8%). Among the alterations, missense mutation was the most common variant classification. Interestingly, we noticed that the mutation distribution of these 20 genes in the high and low risk groups was statistically significant (p < 0.05).
FIGURE 5. Determination of the relationship between the eight-lncRNA signature and immunotherapy. (A) The immune function distribution differences in the high and low risk score populations of the LUAD. (B) The waterfall plot shows the top 20 genes mutated in the LUAD, and their difference in the high and low risk groups. (C) The correlation between TMB and the signature tested by the Pearson coefficient. (D) The TMB difference in the high and low-risk patients tested by the Wilcoxon rank-sum. (E) The Pearson coefficient evaluated the correlations between the signature and the immune relevant signatures. (F) The Wilcoxon rank-sum revealed the distribution differences of the immune relevant signatures in high and low-risk patients. (G) The Kaplan–Meier estimator measured the immune relevant signatures’ prognosis sensitive zone by testing in high and low-risk groups, respectively. (H) The correlation between TIDE and the signature tested by the Pearson coefficient. (I) The TIDE difference in the high and low-risk patients tested by the Wilcoxon rank-sum. p-value < 0.05 was considered statistically significant. (J) Colony-forming abilities of OE-MIR4435-2HG manipulated A549 cells were measured by performing colony-formation assays. The left two panels: representative photomicrographs. Data presented as the mean ± SD of three independent experiments. (K) Colony-forming abilities of OE-NKILA manipulated A549 cells were measured by performing colony-formation assays. The left two panels: representative photomicrographs. Data presented as the mean ± SD of three independent experiments. (L) Effect of MIR4435-2HG overexpression on PDCD1 and PRF1 expression. Data presented as the mean ± SD of three independent experiments. (M) Effect of NKILA overexpression on PDCD1 and PRF1 expression. Data presented as the mean ± SD of three independent experiments. TMB, Tumor mutational burden; TIDE, Tumor Immune Dysfunction and Exclusion. ns: p-value > 0.05; *: p-value < 0.05; **: p-value < 0.01; ***: p-value < 0.001; ****: p-value < 0.0001.
We explored the TMB difference among the risk score, and the Pearson coefficient found that the TMB was positively correlated with the risk score (Figure 5C). Moreover, the Wilcoxon test displayed that the higher risk score group had a higher TMB (Figure 5D). According to clinical trials and preclinical studies, the immune checkpoint blockade offers more durable clinical benefits, including treatment responses and long-term survival, to patients with higher TMB (Samstein et al., 2019; Stenzinger et al., 2019). Our results here demonstrated that our high risk LUADs might expect more responses from immunotherapy to a certain extent.
We noticed that 10/15 of the checkpoint relevant signatures, GZMB, CD274, GZMA, CD8A, LAG3, PDCD1, CTLA4, PRF1, TNF, and TBX2 correlated with the eight-lncRNA signature by the Pearson coefficient (Figure 5E). The Wilcoxon examination found TNF, TBX2, PRF1, PDCD1, LAG3, HAVCR2, GZMA, CXCL9, CTLA4, and CD8A, were expressed differently between high and low-risk groups significantly (Figure 5F). Incorporating the Wilcoxon and Pearson analyses, eight genes, including CD8A, CTLA4, GZMA, LAG3, PDCD1, PRF1, TBX2, and TNF were closely connected to the eight-lncRNA signature. Subsequently, we focused on the eight identified checkpoint relevant signatures. We tested these signatures’ prognosis roles in the high and low-risk groups to see the “comfort risk score zone” for potential immunotherapy. As shown in Figure 5G, CD8A, CTLA4, GZMA, LAG3, TBX2, and TNF did not impact the prognosis in neither high nor low-risk LUADs, hinting that these checkpoint genes were incapable of utilizing the risk score for influencing the LUAD prognosis. Notably, PDCD1 could protect the prognosis more in the low risk-groups than that in the high-risk group, and PRF1 impacted the LUAD survival probability more in the high-risk group than that in its low-risk group, which implied PDCD1 and PRF1 target therapies might maximize their influences in their specific risk score zone, respectively. These results suggested that our risk score system could potentially guide immunotherapy choices based on each checkpoint’s “comfort risk score zone”, however, more clinical data are needed to support our conclusions.
Following that, we examined the potential clinical efficacy of immunotherapy based on the risk score subgroups using the TIDE. As a surrogate biomarker, TIDE scores can provide insight into whether a NSCLC patient will respond to immune checkpoint blockades, including anti-PD1 and anti-CTLA4, if therapy is initiated. Higher TIDE prediction scores represent a higher likelihood of immune evasion, suggesting that patients are less likely to benefit from immunotherapy (Jiang et al., 2018; Fu et al., 2020; Chen Y. et al., 2021). According to our results, high-risk patients had lower TIDE scores than low-risk patients, indicating that immunotherapy was more beneficial for them (Figures 5H, I), which corresponded with our “TMB difference” findings above. The TIDE results were consistent with the “comfort risk score zone” of PRF1 we mentioned and could potentially guide the planning of the anti-PD1 (anti-PDCD1) strategy for LUAD treatment based on the risk score.
The bioinformatics analysis results of PDCD1 and PRF1 are exciting to us. The role of two of the eight genes in our model, MIR4435-2HG and NKILA, negatively affecting prognosis in LUAD was also of interest. Using a colony formation assay, we further examined the effect of MIR4435-2HG and NKILA on long-term growth (7∼10 days) of A549 cells. As shown in Figures 5J, K, the numbers of colonies of the A549 cell line increased following the overexpression of MIR4435-2HG and NKILA, respectively. These results indicated that MIR4435-2HG and NKILA might be deeply involved in LUAD progression by inducing the growth and motility of cells. The real-time PCR results (Figures 5L, M) showed that the immune checkpoints PDCD1 and PRF1 were downregulated after overexpression of MIR4435-2HG or NKILA, suggesting that PDCD1 and PRF1 had a mutual regulatory relationship with MIR4435-2HG or NKILA, respectively.
Risk score and anticancer drug sensitivity
We screened out 8 chemotherapy drugs (cisplatin, paclitaxel, doxorubicin, etoposide, gemcitabine, docetaxel, methotrexate, and bleomycin), 2 epidermal growth factor receptor (EGFR) inhibitors (gefitinib, erlotinib), and 1 immunosuppressive drug (rapamycin) from the GDSC database, which had been used in LUAD clinical practice. By adopting “pRRophetic” package, we found 7 chemotherapy drugs (cisplatin, paclitaxel, doxorubicin, etoposide, gemcitabine, docetaxel, and methotrexate) and 1 EGFR inhibitor (erlotinib) showed significantly sensitivity differences in the high and low-risk patients, while, Bleomycin, gefitinib, and rapamycin exhibited no correlation with our risk score model. As specifically shown in Figure 6, cisplatin, paclitaxel, doxorubicin, etoposide, gemcitabine, docetaxel, and methotrexate received more feedback from the high-risk patients, otherwise, erlotinib got a more intensive response in the low-risk LUADs, which hinted the eight-lncRNA signature’s potential usage in the clinical practice.
FIGURE 6. The disparities in our selected anticancer drugs’ responses between groups with low and high risk scores. The difference and correlation in the IC50 between the different risk groups was analyzed by the Wilcoxon test (lower plots) and the Spearman’s rank correlation coefficient, respectively (upper plots). IC50: half maximal inhibitory concentration.
Mast cells’ vital roles in the eight-lncRNA signature’s prediction ability
Supplementary Figure S5A details the distribution of the 22 TICs in each patient and high and low-risk groups. Supplementary Figure S5B shows the mutual internal relationships of the 22 TICs in the LUADs. As shown in Figure 7A, the Wilcoxon rank sum test identified 9 TICs related to the risk score. The Pearson coefficient (Figure 7B; Supplementary Table S5) found 12 TICs closely linked to our signature. In summary (Figure 7C), a total of 9 TICs are significantly related to the lncRNA signature, which included Mast cells resting, Macrophages M0, B cells memory, Mast cells activated, T cells CD4 memory activated, T cells CD4 memory resting, Neutrophils, Dendritic cells activated, and Monocytes. Specifically, our signature was positively correlated with Macrophages M0, Mast cells activated, T cells CD4 memory activated, Neutrophils, Dendritic cells activated, and negatively correlated with the rest.
FIGURE 7. Comprehensive analysis determined the relationship between 22 TICs and the eight-lncRNA signature. (A) Patients were divided into low-risk groups and high-risk groups based on the best cutoff defined by the “surv_cutpoint ()” function of the “survminer” R package. The Wilcoxon rank-sum was applied to detect the difference in TIC distribution between the high and low-risk groups. (B) The Pearson coefficient examined the correlations between TICs and our signature. Here, we only plotted the TIC correlation with a p-value < 0.05. (C) The results of Wilcoxon’s rank-sum test and the Pearson coefficient were intersected to determine stable and critical TICs. (D) We performed Cox analysis to determine the survival predicting power for the 22 TICs. (E) Kaplan-Meier estimators were plotted to check the TICs that can differ the survival possibility between low and high-risk groups. Only showed Kaplan-Meier curves with p-values < 0.05. (F) The distribution of T-cell inflamed scores of the high- and low -risk groups is displayed in the form of a box plot. Significant differences were detected by the Wilcoxon test. HR, hazard ratio; TIC, tumor-infiltrating immune cell; p-value < 0.05 is considered significant; ns: p-value > 0.05; *: p-value < 0.05; **: p-value < 0.01; ***: p-value < 0.001; ****: p-value < 0.0001.
We then constructed the Cox analysis (Figure 7D) and Kaplan–Meier curve (Figure 7E) to assess the survival predictive power of each of the 22 TICs. Cox analysis can tell that Mast cells resting (HR = 0.764, 95% CI = 0.645–0.905, p = 1.84e-03) and Mast cells activated (HR = 1.141, 95% CI = 1.019–1.279, p = 2.25e-02) significantly affected the LUAD prognosis. Furthermore, Mast cells resting, Mast cells activated, and Dendritic cells resting significantly impacted the survival possibility of LUAD in the constructed Kaplan-Meier analysis (Figure 7E; Supplementary Table S6). On the whole, the survival predictive analysis suggested the Mast cells resting and Mast cells activated pronounced closely connected to LUAD outcomes. T-cell related score were obtained from the 22 TICs and combined as T-cell inflamed score. Figure 7F revealed a disparity in the T-cell inflamed score between the high-risk and low-risk score groups, which was unexpected. However, no statistical difference was observed after conducting a correlation test between the T-cell inflamed and risk scores (Supplementary Figure S5C). The KM analysis of the T-cell inflamed score indicated that it had no impact on the prognosis of LUAD (Supplementary Figure S5D).
To sum, our study in this section revealed that mast cells, which include the mast cells resting and mast cells activated, significantly correlated with the eight-lncRNA signature and bounded up with the LUAD outcomes, insinuating that mast cells infiltration might play an important part in our discovered signature.
Discussion
Programmed cell death has dramatically expanded the anticancer arsenal. Since ferroptosis, necroptosis, and pyroptosis were discovered, we have better-understood carcinogenesis mechanisms and clinical assessment (Lelievre et al., 2020; Ge et al., 2022). Copper plays an essential cofactor in multiple biological processes during tumor metastasis (Lelievre et al., 2020; Ge et al., 2022). A recent scientific study proposes a new form of programmed cell death: cuproptosis. Cuproptosis leads to protein aggregation, iron-sulfur cluster protein loss, and proteotoxic stress, ultimately leading to cell death (Tsvetkov et al., 2022). Studies have shown that cuproptosis affects the immune microenvironment of various tumors, such as sarcoma (Chu et al., 2022) and melanoma (Lv et al., 2022), thereby affecting the treatment effect and prognosis of patients. Multiple studies have shown that modulating lncRNAs may become a new therapeutic approach for human cancers (Bhan et al., 2017; Peng et al., 2017). However, no studies have focused on systematically assessing cuproptosis-related lncRNA signatures and their association with overall survival in LUAD patients. In the present research, we innovatively adopt the novel cuproptosis concept to establish a cuproptosis-related eight-lncRNA signature predicting LUAD outcomes by digging TCGA and GEO databases (with a total sample size exceeding 1000 cases). A key novelty was the use of comprehensive bioinformatics analysis, such as LASSO regression, Kaplan-Meier curves, Cox models, ROC curves, tAUC, nomogram, and decision curve analysis, as well as the validation in a large independent cohort. Most importantly, we detailed the expression pattern of the signature genes in real world, the role of risk score participating in immunotherapy and distinguishing ability for some anticancer-drug sensitive. The immune infiltration analysis speculated that mast cells may contribute to the predictive power of the signature.
Our signature contains eight lncRNAs (Table 3), which were TSPOAP1-AS1, AC107464.3, AC006449.7, LINC00324, COLCA1, HAGLR, MIR4435-2HG, and NKILA. In our research, MIR4435-2HG and NKILA adversely affected LUAD outcomes, while the remaining genes displayed favorable impacts (Supplementary Figure S3). We performed real-time PCR validation and found that all genes were significantly differentially expressed in LUAD and normal cell lines. In addition, the impact of MIR4435-2HG and NKILA on the extended growth (7∼10 days) of A549 cells was verified through a colony formation assay. Overexpression of MIR4435-2HG and NKILA resulted in an increase in the number of A549 cell line colonies, suggesting their potential involvement in promoting cell growth and motility and their potential significance in the progression of LUAD. There is a strong expression of MIR4435-2HG in lung cancer tissues and tumor cell lines (Zhong et al., 2022). In addition, MIR4435-2HG is involved in the mechanism of action of anticancer drugs, including cisplatin for the treatment of non-small cell lung cancer and colon cancer, and carboplatin for the treatment of triple-negative breast cancer (Zhong et al., 2022). Clinically, NKILA overexpression in tumor-specific CTL and TH1 cells correlates with shorter patient survival and increased apoptosis (Huang et al., 2018; Dai et al., 2021). The NKILA lncRNA stimulates activation-induced cell death in T cells, promoting tumor immunoevasion (Huang et al., 2018). In our wet-experiments, we also confirmed that the MIR4435-2HG and NKILA were up-expressed in LUAD cell lines, which are consistently with their prognosis function found (Supplementary Figure S3) in this study and the previous findings shown above. We uncovered that TSPOAP1-AS1, AC107464.3, AC006449.7, LINC00324, COLCA1, and HAGLR were down-expressed in LUAD cell lines and displayed them having favorable impacts on the LUAD outcomes, which also confirmed the gene signature components’ power. In exceptional cases, the relationship between the differential expression pattern of specific genes and its impact on the prognosis of LUAD is not static, which is related to their regulatory status in the disease, and more research is needed to explore. We found that 8 genes were differentially expressed in the real world, further corroborating the gene signature’s validity and coming forward for further work.
The way of regulated cell death of cuproptosis, apoptosis, necroptosis, pyroptosis, and ferroptosis varies, but they are also somewhat related (Chen X. et al., 2021; Tsvetkov et al., 2022). Apoptosis is mediated by proteolytic enzymes called caspases, which trigger cell death by cleaving specific proteins in the cytoplasm and nucleus. Beyond classical apoptosis, several forms of regulated cell death have been identified. In terms of the exact mechanism by which copper ions cause cell death, several hypotheses have been proposed, including the induction of apoptosis, caspase-independent cell death, the induction of reactive oxygen species, and inhibition of the ubiquitin-proteasome system (Elmore, 2007; Li et al., 2022). Necroptosis is an alternative mode of regulated cell death mimicking features of apoptosis and necrosis (Dhuriya and Sharma, 2018). Ferroptosis is a type of programmed cell death dependent on iron and characterized by the accumulation of lipid peroxides, and is genetically and biochemically distinct from other forms of regulated cell death (Ma et al., 2021). Pharmacological or genetic inhibition of apoptosis (using the caspase inhibitor Z-VAD-FMK or double knockout of BAK and BAX), ferroptosis (using ferrostatin-1 and liproxstatin-1), and necroptosis (using necrostatin-1) failed to suppress cell death induced by the ES–Cu complex in multiple cancer cell lines (Tang et al., 2022). Pyroptosis represents a form of cell death that is triggered by proinflammatory signals and associated with inflammation. Inhibition of cuproptosis has been shown to be involved in mediating inflammation, resulting in the excessive survival and proliferation of a variety of immune cells, such as fibroblast-like synoviocytes, effector T cells, and macrophages (Zhao et al., 2022). From our research, these regulated cell deaths seem to be somewhat related, such as our cuproptosis-related signature correlated with 2321/3681 (63.05%) apoptosis-related genes, 11/20 (55.00%) necroptosis-related genes, 34/50 (68.00%) pyroptosis-related genes, and 222/380 (58.42%) ferroptosis-related genes. Our discovery may provide potential explanations and inspirations for further research of cell death-related tumor mechanisms.
Continued innovation and advances in cancer immunotherapy have extended survival for some patients with deadly cancers (Vanneman and Dranoff, 2012; Waldman et al., 2020). The approach is revolutionizing the field of oncology as more and more cancer patients become eligible for immuno-based therapies (Vanneman and Dranoff, 2012; Waldman et al., 2020). Targeted strategies inhibit tumor progression by interfering with crucial molecular pathways, while immunotherapy produces durable and effective tumor destruction by stimulating the host’s own response (Vanneman and Dranoff, 2012; Waldman et al., 2020). To maximize the benefits of immunotherapy, it is crucial to determine whether a specific biomarker is suitable for the host (Suresh et al., 2018). This study gives hints about which immunotherapy targets to use or under what circumstances to apply. We first found that our risk score was associated with TMB, suggesting that our signature appeared to guide immunotherapy. Next, we followed the trail and found eight checkpoints, including CD8A, CTLA4, GZMA, LAG3, PDCD1, PRF1, TBX2, and TNF related to our eight-lncRNA signature score. We analyzed the sensitivity of each target in different risk score intervals and found that the targeting of PDCD1 was more effective in low-risk patients and triggering for PRF1 can benefit the high-risk group more. Our findings hinted that PDCD1 and PRF1 target therapies might maximize their influences in their specific lncRNA signature risk score zone, respectively. Interestingly, Voli and colleagues stated a strong correlation between PD-L1 expression and intracellular copper levels, noting that copper dysregulators inhibit PD-L1 in vitro and in vivo, resulting in increased tumor-infiltrating lymphocytes (Voli et al., 2018). In a recent study, Zeng demonstrated that FOXO1, a gene significantly associated with PDCD1, downregulated copper homeostasis is a novel indicator of breast cancer prognosis and immune response (Zeng et al., 2022). The TIDE database related analysis confirmed the above points of our view. Additionally, the real-time PCR results (Figures 5L, M) showed that the immune checkpoints PDCD1 and PRF1 were downregulated after overexpression of MIR4435-2HG or NKILA, suggesting that PDCD1 and PRF1 had a mutual regulatory relationship with MIR4435-2HG or NKILA, respectively. The above evidence affirmed that our risk score could bring hope to precisely targeted therapy.
In the section on TICs analysis of our work, we concluded that the mast cells might help our signature bring off secure prognostic power. Mast cell is known for its involvement in allergic disorders, but in recent years, accumulating evidence has shed light on its involvement in cancer, including LUAD (Stoyanov et al., 2012; Bo et al., 2018; Salamon et al., 2020). Mast cells are involved in disease processes characterized by inflammation and tissue remodeling, and their activated forms promote lung health through innate and adaptive immune responses to respiratory pathogens (Virk et al., 2016). Mast cells promote angiogenesis and secrete growth factors, including VEGF, evidence that mast cells are associated with poor prognosis in NSCLC (Virk et al., 2016). Since the immune microenvironment plays a vital part in the malignancy progression, mast cells, as a critical stromal component of the immune system, are undoubtedly a key regulator for maintaining tissue homeostasis. Also, our study examined the correlation between the T-cell inflamed score and the risk score and observed significant variations in the T-cell inflamed score between the high and low-risk score groups. These findings suggest that T-cell research may emerge as a promising avenue for future LUAD studies. There is no dispute that the current findings are far from sufficient, and further research should pay close attention to the parts of mast cells in tumor microenvironment remodeling and tumor immunity.
This study has certain limitations. We generated this eight-lncRNA signature from publicly accessible data. Although it has been confirmed to have stable prognosis ability through applied to another large independent cohort, its clinical applicability needs further confirmation with more parameters. Furthermore, there are still no wet-lab experimental facts to hold up the eight lncRNAs’ parts in cuproptosis-related mechanisms. Therefore, more research, which focuses in vivo and in vitro, is urgently needed to reveal more clues that support the signature’s promising future.
Conclusion
The present research constructed a novel and capable cuproptosis-related eight-lncRNA prognostic signature for LUAD. By applying the signature to independent cohort, its stability and broad applicability was validated. The signature owns the potential ability to undertake the role of precise immunotherapy and anticancer drug selection. The immune infiltration analysis hinted that the mast cells might help the signature to maintain its predictive power. Our work may promote the evolution of diagnosis and treatment of LUAD.
Novelty & impact statements
In the present research, we adopt the novel concept to establish a cuproptosis-related eight-lncRNA signature predicting LUAD outcomes by digging TCGA and GEO databases (with a total sample size exceeding 1000 cases). Most importantly, we detailed the expression pattern of the signature genes in the real world, the role of risk score participating in immunotherapy, and the model’s distinguishing power for anticancer drugs. The mast cells may contribute to the predictive power of the signature.
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: The following publicly available datasets were used in this study: TCGA-LUAD: https://gdc.xenahubs.net; GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081: https://www.ncbi.nlm.nih.gov/geo/; The data of real-time PCR can be reached by contacting the corresponding authors.
Author contributions
CM participated in the study design and drafting of the manuscript. FL and ZG were responsible for data analysis. YY, and YQ supervised the entire research project. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by funds from the Leading Talent Cultivation Project of Henan Health Science and Technology Innovation Talents (Grant No. YXKC2021016, YXKC2022014), the Henan Provincial College and University Research Project (Grant No. 22A320072), the Henan Provincial Science and Technology Development Project (Grant Nos. SBGJ202102122, 222102310239), the Medical Science and Technology Research Project of Henan Province: Joint Construction Project (Grant No. LHGJ20220282), and the Natural Science Foundation of Henan Province: Youth Project (Grant No. 232300420236).
Acknowledgments
CM thanks the post-doctoral research station of the First Affiliated Hospital of Zhengzhou University.
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/fphar.2023.1146840/full#supplementary-material
Abbreviations
LUAD, lung adenocarcinoma; LASSO, least absolute shrinkage and selection operator Cox regression model; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic; AUC, Area under the ROC curve; tAUC, Time-dependent AUC; GSEA, Gene Set Enrichment Analysis; TICs, tumor-infiltrating immune cells; HR, hazard ratio; CI, confidence interval; FDR, false discovery rate; TIDE, Tumor Immune Dysfunction and Exclusion; PC, principal component; KEGG, Kyoto Encyclopedia of Genes and Genomes; IC50, half maximal inhibitory concentration; UMAP, Uniform Manifold Approximation and Projection; PCA, Principal components analysis; OS, overall survival; TMB, Tumor mutational burden.
References
Alishahi, M., Ghaedrahmati, F., Kolagar, T. A., Winlow, W., Nikkar, N., Farzaneh, M., et al. (2019). Long non-coding RNAs and cell death following ischemic stroke. Metab. Brain Dis. 34 (5), 1243–1251. doi:10.1007/s11011-019-00423-2
Andrews, L. P., Marciscano, A. E., Drake, C. G., and Vignali, D. A. (2017). LAG3 (CD223) as a cancer immunotherapy target. Immunol. Rev. 276 (1), 80–96. doi:10.1111/imr.12519
Bhan, A., Soleimani, M., and Mandal, S. S. (2017). Long noncoding RNA and cancer: a new paradigm. Cancer Res. 77 (15), 3965–3981. doi:10.1158/0008-5472.CAN-16-2634
Bian, Z., Fan, R., and Xie, L. (2022). A novel cuproptosis-related prognostic gene signature and validation of differential expression in clear cell renal cell carcinoma. Genes (Basel). 13 (5), 851. doi:10.3390/genes13050851
Bo, X., Wang, J., Suo, T., Ni, X., Liu, H., Shen, S., et al. (2018). Tumor-infiltrating mast cells predict prognosis and gemcitabine-based adjuvant chemotherapeutic benefit in biliary tract cancer patients. BMC Cancer 18 (1), 313. doi:10.1186/s12885-018-4220-1
Cao, R., and Lopez-de-Ullibarri, I. (2019). ROC curves for the statistical analysis of microarray data. Methods Mol. Biol. 1986, 245–253. doi:10.1007/978-1-4939-9442-7_11
Chalmers, Z. R., Connelly, C. F., Fabrizio, D., Gay, L., Ali, S. M., Ennis, R., et al. (2017). Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 9 (1), 34. doi:10.1186/s13073-017-0424-2
Chen, L., Min, J., and Wang, F. (2022). Copper homeostasis and cuproptosis in health and disease. Signal Transduct. Target Ther. 7 (1), 378. doi:10.1038/s41392-022-01229-y
Chen, X., Kang, R., Kroemer, G., and Tang, D. (2021b). Broadening horizons: the role of ferroptosis in cancer. Nat. Rev. Clin. Oncol. 18 (5), 280–296. doi:10.1038/s41571-020-00462-0
Chen, Y., Li, Z. Y., Zhou, G. Q., and Sun, Y. (2021a). An immune-related gene prognostic index for head and neck squamous cell carcinoma. Clin. Cancer Res. 27 (1), 330–341. doi:10.1158/1078-0432.ccr-20-2166
Chu, B., Zheng, H., Zheng, X., Feng, X., and Hong, Z. (2022). Cuproptosis-associated lncRNAs discern prognosis and immune microenvironment in sarcoma victims. Front. Cell Dev. Biol. 10, 989882. doi:10.3389/fcell.2022.989882
Clough, E., and Barrett, T. (2016). The gene expression Omnibus database. Methods Mol. Biol. 1418, 93–110. doi:10.1007/978-1-4939-3578-9_5
Dai, S., Liu, T., Liu, Y. Y., He, Y., Liu, T., Xu, Z., et al. (2021). Long non-coding RNAs in lung cancer: the role in tumor microenvironment. Front. Cell Dev. Biol. 9, 795874. doi:10.3389/fcell.2021.795874
Derrien, T., Johnson, R., Bussotti, G., Tanzer, A., Djebali, S., Tilgner, H., et al. (2012). The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 22 (9), 1775–1789. doi:10.1101/gr.132159.111
Dhuriya, Y. K., and Sharma, D. (2018). Necroptosis: a regulated inflammatory mode of cell death. J. Neuroinflammation 15 (1), 199. doi:10.1186/s12974-018-1235-0
Elmore, S. (2007). Apoptosis: a review of programmed cell death. Toxicol. Pathol. 35 (4), 495–516. doi:10.1080/01926230701320337
Fabrizio, F. P., Trombetta, D., Rossi, A., Sparaneo, A., Castellana, S., and Muscarella, L. A. (2018). Gene code cd274/PD-L1: from molecular basis toward cancer immunotherapy. Ther. Adv. Med. Oncol. 10, 1758835918815598. doi:10.1177/1758835918815598
Fan, C., Hu, H., Shen, Y., Wang, Q., Mao, Y., Ye, B., et al. (2021). PRF1 is a prognostic marker and correlated with immune infiltration in head and neck squamous cell carcinoma. Transl. Oncol. 14 (4), 101042. doi:10.1016/j.tranon.2021.101042
Freeman, A. J., Kearney, C. J., Silke, J., and Oliaro, J. (2021). Unleashing TNF cytotoxicity to enhance cancer immunotherapy. Trends Immunol. 42 (12), 1128–1142. doi:10.1016/j.it.2021.10.003
Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33 (1), 1–22. doi:10.18637/jss.v033.i01
Fu, J., Li, K., Zhang, W., Wan, C., Zhang, J., Jiang, P., et al. (2020). Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 12 (1), 21. doi:10.1186/s13073-020-0721-z
Ge, E. J., Bush, A. I., Casini, A., Cobine, P. A., Cross, J. R., DeNicola, G. M., et al. (2022). Connecting copper and cancer: from transition metal signalling to metalloplasia. Nat. Rev. Cancer 22 (2), 102–113. doi:10.1038/s41568-021-00417-2
Geeleher, P., Cox, N., and Huang, R. S. (2014b). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468
Geeleher, P., Cox, N. J., and Huang, R. S. (2014a). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol. 15 (3), R47. doi:10.1186/gb-2014-15-3-r47
Goeman, J. J. (2010). L1 penalized estimation in the Cox proportional hazards model. Biom J. 52 (1), 70–84. doi:10.1002/bimj.200900028
Holderried, T. A. W., de Vos, L., Bawden, E. G., Vogt, T. J., Dietrich, J., Zarbl, R., et al. (2019). Molecular and immune correlates of TIM-3 (HAVCR2) and galectin 9 (LGALS9) mRNA expression and DNA methylation in melanoma. Clin. Epigenetics 11 (1), 161. doi:10.1186/s13148-019-0752-8
Huang, D., Chen, J., Yang, L., Ouyang, Q., Li, J., Lao, L., et al. (2018). NKILA lncRNA promotes tumor immune evasion by sensitizing T cells to activation-induced cell death. Nat. Immunol. 19 (10), 1112–1125. doi:10.1038/s41590-018-0207-y
Hurkmans, D. P., Basak, E. A., Schepers, N., Oomen-De Hoop, E., Van der Leest, C. H., El Bouazzaoui, S., et al. (2020). Granzyme B is correlated with clinical outcome after PD-1 blockade in patients with stage IV non-small-cell lung cancer. J. Immunother. Cancer 8 (1), e000586. doi:10.1136/jitc-2020-000586
Inoue, H., Park, J. H., Kiyotani, K., Zewde, M., Miyashita, A., Jinnin, M., et al. (2016). Intratumoral expression levels of PD-L1, GZMA, and HLA-A along with oligoclonal T cell expansion associate with response to nivolumab in metastatic melanoma. Oncoimmunology 5 (9), e1204507. doi:10.1080/2162402X.2016.1204507
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1
Johnson, W. E., Li, C., and Rabinovic, A. (2007). Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8 (1), 118–127. doi:10.1093/biostatistics/kxj037
Jorgovanovic, D., Song, M., Wang, L., and Zhang, Y. (2020). Roles of IFN-gamma in tumor progression and regression: a review. Biomark. Res. 8, 49. doi:10.1186/s40364-020-00228-x
Lelievre, P., Sancey, L., Coll, J. L., Deniaud, A., and Busser, B. (2020). The multifaceted roles of copper in cancer: a trace metal element with dysregulated metabolism, but also a target or a bullet for therapy. Cancers (Basel). 12 (12), 3594. doi:10.3390/cancers12123594
Li, C., Pan, J., Jiang, Y., Yu, Y., Jin, Z., and Chen, X. (2021). Characteristics of the immune cell infiltration landscape in gastric cancer to assistant immunotherapy. Front. Genet. 12, 793628. doi:10.3389/fgene.2021.793628
Li, S. R., Bu, L. L., and Cai, L. (2022). Cuproptosis: lipoylated TCA cycle proteins-mediated novel cell death pathway. Signal Transduct. Target Ther. 7 (1), 158. doi:10.1038/s41392-022-01014-x
Liang, Y. K., Deng, Z. K., Chen, M. T., Qiu, S. Q., Xiao, Y. S., Qi, Y. Z., et al. (2021). CXCL9 is a potential biomarker of immune infiltration associated with favorable prognosis in ER-negative breast cancer. Front. Oncol. 11, 710286. doi:10.3389/fonc.2021.710286
Liberzon, A., Subramanian, A., Pinchback, R., Thorvaldsdottir, H., Tamayo, P., and Mesirov, J. P. (2011). Molecular signatures database (MSigDB) 3.0. Bioinformatics 27 (12), 1739–1740. doi:10.1093/bioinformatics/btr260
Liu, J., Hong, M., Li, Y., Chen, D., Wu, Y., and Hu, Y. (2022). Programmed cell death tunes tumor immunity. Front. Immunol. 13, 847345. doi:10.3389/fimmu.2022.847345
Lv, H., Liu, X., Zeng, X., Liu, Y., Zhang, C., Zhang, Q., et al. (2022). Comprehensive analysis of cuproptosis-related genes in immune infiltration and prognosis in melanoma. Front. Pharmacol. 13, 930041. doi:10.3389/fphar.2022.930041
Ma, C., Li, F., and Luo, H. (2021). Prognostic and immune implications of a novel ferroptosis-related ten-gene signature in lung adenocarcinoma. Ann. Transl. Med. 9 (13), 1058. doi:10.21037/atm-20-7936
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12 (5), 453–457. doi:10.1038/nmeth.3337
Peng, W. X., Koirala, P., and Mo, Y. Y. (2017). LncRNA-mediated regulation of cell signaling in cancer. Oncogene 36 (41), 5661–5667. doi:10.1038/onc.2017.184
Platten, M., von Knebel Doeberitz, N., Oezen, I., Wick, W., and Ochs, K. (2014). Cancer immunotherapy by targeting Ido1/TDO and their downstream effectors. Front. Immunol. 5, 673. doi:10.3389/fimmu.2014.00673
Raskov, H., Orhan, A., Christensen, J. P., and Gogenur, I. (2021). Cytotoxic CD8(+) T cells in cancer and cancer immunotherapy. Br. J. Cancer 124 (2), 359–367. doi:10.1038/s41416-020-01048-4
Ribas, A., and Wolchok, J. D. (2018). Cancer immunotherapy using checkpoint blockade. Science 359 (6382), 1350–1355. doi:10.1126/science.aar4060
Rowshanravan, B., Halliday, N., and Sansom, D. M. (2018). CTLA-4: a moving target in immunotherapy. Blood 131 (1), 58–67. doi:10.1182/blood-2017-06-741033
Salamon, P., Mekori, Y. A., and Shefler, I. (2020). Lung cancer-derived extracellular vesicles: a possible mediator of mast cell activation in the tumor microenvironment. Cancer Immunol. Immunother. 69 (3), 373–381. doi:10.1007/s00262-019-02459-w
Samstein, R. M., Lee, C. H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y., et al. (2019). Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat. Genet. 51 (2), 202–206. doi:10.1038/s41588-018-0312-8
Sauerbrei, W., Royston, P., and Binder, H. (2007). Selection of important variables and determination of functional form for continuous predictors in multivariable model building. Stat. Med. 26 (30), 5512–5528. doi:10.1002/sim.3148
Shi, Z., Shen, J., Qiu, J., Zhao, Q., Hua, K., and Wang, H. (2021). CXCL10 potentiates immune checkpoint blockade therapy in homologous recombination-deficient tumors. Theranostics 11 (15), 7175–7187. doi:10.7150/thno.59056
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2021). Cancer statistics, 2021. CA Cancer J. Clin. 71 (1), 7–33. doi:10.3322/caac.21654
Spella, M., and Stathopoulos, G. T. (2021). Immune resistance in lung adenocarcinoma. Cancers (Basel) 13 (3), 384. doi:10.3390/cancers13030384
Stenzinger, A., Allen, J. D., Maas, J., Stewart, M. D., Merino, D. M., Wempe, M. M., et al. (2019). Tumor mutational burden standardization initiatives: recommendations for consistent tumor mutational burden assessment in clinical samples to guide immunotherapy treatment decisions. Genes Chromosom. Cancer 58 (8), 578–588. doi:10.1002/gcc.22733
Stoyanov, E., Uddin, M., Mankuta, D., Dubinett, S. M., and Levi-Schaffer, F. (2012). Mast cells and histamine enhance the proliferation of non-small cell lung cancer cells. Lung Cancer 75 (1), 38–44. doi:10.1016/j.lungcan.2011.05.029
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis:a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102
Suresh, K., Naidoo, J., Lin, C. T., and Danoff, S. (2018). Immune checkpoint immunotherapy for non-small cell lung cancer: benefits and pulmonary toxicities. Chest 154 (6), 1416–1423. doi:10.1016/j.chest.2018.08.1048
Taminau, J., Meganck, S., Lazar, C., Steenhoff, D., Coletta, A., Molter, C., et al. (2012). Unlocking the potential of publicly available microarray data using inSilicoDb and inSilicoMerging R/Bioconductor packages. BMC Bioinforma. 13, 335. doi:10.1186/1471-2105-13-335
Tang, D., Chen, X., and Kroemer, G. (2022). Cuproptosis: a copper-triggered modality of mitochondrial cell death. Cell Res. 32 (5), 417–418. doi:10.1038/s41422-022-00653-7
Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Stat. Med. 16 (4), 385–395. doi:10.1002/(sici)1097-0258(19970228)16:4<385:aid-sim380>3.0.co;2-3
Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 375 (6586), 1254–1261. doi:10.1126/science.abf0529
Vanneman, M., and Dranoff, G. (2012). Combining immunotherapy and targeted therapies in cancer treatment. Nat. Rev. Cancer 12 (4), 237–251. doi:10.1038/nrc3237
Virk, H., Arthur, G., and Bradding, P. (2016). Mast cells and their activation in lung disease. Transl. Res. 174, 60–76. doi:10.1016/j.trsl.2016.01.005
Voli, F., Lerra, L., Kimpton, K., Saletta, F., Cirillo, G., Kavallaris, M., et al. (2018). Harnessing copper in cancer to enhance anti-tumor immune response. Ann. Oncol. 29, x35. doi:10.1093/annonc/mdy487.030
Waldman, A. D., Fritz, J. M., and Lenardo, M. J. (2020). A guide to cancer immunotherapy: from T cell basic science to clinical practice. Nat. Rev. Immunol. 20 (11), 651–668. doi:10.1038/s41577-020-0306-5
Ye, L., Zhang, T., Kang, Z., Guo, G., Sun, Y., Lin, K., et al. (2019). Tumor-infiltrating immune cells act as a marker for prognosis in colorectal cancer. Front. Immunol. 10, 2368. doi:10.3389/fimmu.2019.02368
Zeng, R., Peng, B., and Peng, E. (2022). Downregulated copper homeostasis-related gene FOXO1 as a novel indicator for the prognosis and immune response of breast cancer. J. Immunol. Res. 2022, 9140461. doi:10.1155/2022/9140461
Zhang, Z., Zeng, X., Wu, Y., Liu, Y., Zhang, X., and Song, Z. (2022). Cuproptosis-related risk score predicts prognosis and characterizes the tumor microenvironment in hepatocellular carcinoma. Front. Immunol. 13, 925618. doi:10.3389/fimmu.2022.925618
Zhao, J., Guo, S., Schrodi, S. J., and He, D. (2022). Cuproptosis and cuproptosis-related genes in rheumatoid arthritis: implication, prospects, and perspectives. Front. Immunol. 13, 930278. doi:10.3389/fimmu.2022.930278
Zhong, C., Xie, Z., Zeng, L. H., Yuan, C., and Duan, S. (2022). MIR4435-2HG is a potential pan-cancer biomarker for diagnosis and prognosis. Front. Immunol. 13, 855078. doi:10.3389/fimmu.2022.855078
Keywords: lung adenocarcinoma, lncRNA, signature, cuproptosis, copper, tumor microenvironment
Citation: Ma C, Li F, Gu Z, Yang Y and Qi Y (2023) A novel defined risk signature of cuproptosis-related long non-coding RNA for predicting prognosis, immune infiltration, and immunotherapy response in lung adenocarcinoma. Front. Pharmacol. 14:1146840. doi: 10.3389/fphar.2023.1146840
Received: 17 January 2023; Accepted: 10 August 2023;
Published: 21 August 2023.
Edited by:
Zhaoshi Bai, The Affiliated Cancer Hospital of Nanjing Medical University, ChinaCopyright © 2023 Ma, Li, Gu, Yang and Qi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yang Yang, fccyangy1@zzu.edu.cn; Yu Qi, qiyu@zzu.edu.cn
†ORCID: Chao Ma, https://orcid.org/0000-0003-1444-4668