Association of aging-related genes with prognosis and immune infiltration in pancreatic adenocarcinoma
Shengbai Xue1† 
Weiyu Ge1† 
Kexuan Wang2† Tiebo Mao1 
Xiaofei Zhang1 Haiyan Xu1 Yongchao Wang1 Jiayu Yao1 Shumin Li1 
Ming Yue1 
Jingyu Ma1 
Yanling Wang1 
Daiyuan Shentu1 
Jiujie Cui1* 
Liwei Wang1*- 1Department of Oncology, Renji Hospital, School of Medicine, Shanghai Jiaotong University, Shanghai, China
- 2Department of Nursing, School of Nursing, Xuzhou Medical University, Xuzhou, China
Pancreatic adenocarcinoma (PAAD) is one of the deadliest malignancies. Aging is described as the degeneration of physiological function, which is complexly correlated with cancer. It is significant to explore the influences of aging-related genes (ARGs) on PAAD. Based on The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) datasets, we used univariate Cox regression analysis and acquired eight differentially expressed ARGs with prognostic values. Two molecular subtypes were identified based on these ARGs to depict PAAD patients’ overall survival (OS) and immune microenvironments preliminarily. Cluster 1 had a poor OS as well as a worse immune microenvironment. Through least absolute shrinkage and selection operator (LASSO) regression analysis, we constructed a seven-ARG risk signature based on the TCGA dataset and verified it in Gene Expression Omnibus (GEO) and International Cancer Genome Consortium (ICGC) to predict the prognoses, immune microenvironments, signal pathways, tumor mutations, and drug sensitivity of PAAD patients. The high-risk group possessed an unfavorable OS compared with that of the low-risk group. We also verified the independence and clinical availability of the risk signature by Cox regression analyses and the establishment of a nomogram, respectively. The higher risk score was associated with several clinical factors such as higher grade and advanced tumor stage as well as lower immunoscore and cluster 1. The negative associations of risk scores with immune, stroma, and estimate scores proved the terrible immune microenvironment in the high-risk group. Relationships between risk score and immune checkpoint gene expression as well as signal pathways provided several therapeutic targets. PAAD patients in the low-risk group possessed lower tumor mutations as well as a higher susceptibility to axitinib and vorinostat. The high-risk group bore a higher TMB and cisplatin and dasatinib may be better options. We used immunohistochemistry and qPCR to confirm the expression of key ARGs with their influences on OS. In conclusion, we identified two ARG-mediated molecular subtypes and a novel seven-ARG risk signature to predict prognoses, immune microenvironments, signal pathways, tumor mutations, and drug sensitivity of PAAD patients.
Introduction
As a lethal malignant tumor, pancreatic adenocarcinoma (PAAD) remains the fourth-leading cause of cancer-related death globally with overall survival (OS) of less than 10% (Mizrahi et al., 2020). Surgery is considered the only potential way to cure PAAD (Okasha et al., 2017). However, late detection results in most patients being failed to tolerate surgery. Recently, immunotherapy has become a novel treatment option. Including immune checkpoint inhibitors (ICIs), tumor vaccines, and monoclonal antibodies, immunotherapy has shown unique efficacy in PAAD (Schizas et al., 2020). However, PAAD was described as a “cold tumor” with a terrible immune microenvironment (Clark et al., 2007). Therefore, it is necessary to explore novel biomarkers for predicting survival outcomes as well as immune microenvironment infiltrations of PAAD patients.
Aging, considered a gradual break of physiological integrity, has been a crucial risk factor for cancer development (López-Otín et al., 2013). Older people take up a large proportion of cancer patients, and cancer has been the first cause of death in people aged 60–79 years (Siegel et al., 2022). The interaction between aging and cancer has been studied widely. Accumulation of DNA damage, oxidative stress, and senescence-associated secretory phenotype secretion associated with aging contribute to initial tumorigenesis (Berben et al., 2021). The decline of immune functions caused by aging also leads to the increasing incidence of cancer (Camous et al., 2012). In PAAD, a previous study has found the aging of normal fibroblasts promoted the progression of cancer (Sarsour et al., 2020). However, the relationships between aging-related genes (ARGs) and prognoses as well as immune infiltrations of PAAD patients remain unclear. In previous studies, ferroptosis-related signatures (Jiang et al., 2021), immune-related signatures (Wang et al., 2020), and glycolysis-related signatures (Song et al., 2021) have been constructed to predict the OS of PAAD patients. Therefore, in this study, our purpose was to identify molecular subtypes as well as establish a risk signature based on ARGs to predict prognoses and immune microenvironments of PAAD patients.
Materials and methods
Public datasets
RNA sequencing data and relevant clinical features of PAAD patients (Table 1) were obtained from The Cancer Genome Atlas (TCGA), Genotype-Tissue Expression project (GTEx), Gene Expression Omnibus (GEO), and International Cancer Genome Consortium (ICGC). The expression data from TCGA, GEO, and ICGC was processed by Perl software and normalized by the R package “sva.” Patients with incomplete clinical characteristics or unavailable expression profiles were excluded. The list of ARGs (Supplementary Table S1) was achieved from the Human Ageing Genomic Resources (HAGR). Immunohistochemical results were obtained from Human Protein Atlas (HPA).
TABLE 1. Clinical characteristics of pancreatic adenocarcinoma patients from TCGA, GEO, and ICGC.
Identification of differentially expressed aging-related genes with prognostic ability
By R package “limma,” we merged and normalized the data from TCGA and GTEx by averaging the multiple expressions of the same genes. We also used it to perform differential analysis between tumor samples from the TCGA dataset (n = 178) and normal samples from TCGA (n = 4) and GTEx (n = 167) datasets with the criteria of |log2FC| ≥ 2 as well as false discovery rate (FDR) <0.05. Taking the intersection of ARGs with differentially expressed genes (DEGs), selected differentially expressed ARGs were involved in univariate Cox regression analysis for the identification of prognostic ARGs (p < 0.05).
Analyses of molecular subtypes
To identify the potential functions of prognostic ARGs, 177 patients from the TCGA dataset were classified into two clusters using the “ConsensusClusterPlus” R package (Seiler et al., 2010). The survival conditions of two subtypes were compared by the R package “survival.” Furthermore, we estimated the immune infiltration of each molecular subtype. The main categories of immune cells in PC were obtained using “CIBERSORT R script v1.03.” A comparison of expression of programmed cell death-ligand 1 (PD-L1) was conducted between two clusters. We also investigated the correlation of prognostic ARGs with PD-L1. Finally, we explored the expressions of the eight ARGs and clinical features in two molecular subtypes.
Construction and validation of an aging-related gene risk signature
PAAD patients from TCGA cohort (n = 177) were involved in training set, while GSE57495 (n = 63) and ICGC-PACA-AU (n = 80) were selected as validation sets. RNA expression data from different datasets were normalized by the R package “sva.” To develop and validate an ARG risk signature, the least absolute shrinkage and selection operator (LASSO) regression analysis was adopted to explore key ARGs most tightly associated with OS. The risk score was computed by the formula:
According to the median risk score, PAAD patients from the training set were separated into low-risk and high-risk subgroups. Furthermore, we performed the survival analysis by R package “survival” to find whether there existed a significant difference between the two groups. To confirm the stability of ARG risk signature at 1-, 2- and 3-year, receiver operating characteristic (ROC) curves were implemented. Afterward, analyses in validation sets validated our aforementioned results.
Association of risk score with molecular subtypes, immunoscore, and clinicopathological characteristics
Univariate and multivariate Cox regression analyses were taken to verify the independent prognostic value of our risk signature. Prognostic factors with p < 0.05 in both analyses were considered independent. For clinical use to predict outcomes of PAAD patients, a nomogram with a risk score and other clinicopathological signatures was established. We drew 1-, 2-, and 3- year calibration curves to confirm its prognostic precision. Moreover, we explored the relationships between risk score and molecular subtypes, ImmuneScore as well as clinical characteristics based on the TCGA cohort. ImmuneScore of each patient was acquired by the R package “estimate.”
Identification of immune microenvironment affected by aging-related genes
Several acknowledged methods were used to explore the relevance between risk score and immune microenvironment, including XCELL, TIMER, QUANTISEQ, MCPOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT. In addition, we evaluated the correlations of risk score and prognostic ARGs with immune checkpoint genes by the R package “limma.”
Correlation of risk signature with signal pathways, tumor mutation, and chemosensitivity
For investigation of the relevance between signal pathways and risk score, we performed Gene Set Variation Analysis (GSVA), and the result was visualized by R packages “reshape2” and “ggplot”. Moreover, we used GeneMANIA (http://genemania.org/) to identify functionally similar genes of key ARGs and their relationships. Mutations and genetic alterations of key ARGs were visualized by cBioPortal. With the help of the R package “maftools,” we compared the mutations of the top 20 genes between two risk groups. Tumor mutation burden (TMB) scores of PAAD patients were computed by somatic mutation analysis. We also evaluated the correlation of risk score with TMB. Finally, we analyzed the survival outcomes of patients from two risk groups with different TMBs. To assess the ability of risk signature to predict chemosensitivity of PAAD, we investigated therapeutic half inhibitory centration (IC50) in two risk groups using the R package “pRRophetic”.
Expression level validation and survival analysis of aging-related genes expression
Immunohistochemical results of ARGs involved in the risk signature were obtained from HPA for validation of ARGs expression in normal and tumor tissue. We furtherly investigated their influences on survival outcomes of PAAD patients based on the TCGA dataset. By RT-qPCR (Real Time Quantitative Polymerase Chain Reaction), we verified the expressions of two key ARGs (BSCL2: HR = 0.653, Coef = −0.13, and TOP2A: HR = 1.104, Coef = 0.07) in normal pancreatic cell line HPNE and two pancreatic adenocarcinoma cell lines PANC1 and CFPAC1. RT-qPCR was performed by a LightCycler®480 real-time PCR system. The ARGs’ primer sequences were shown below: BSCL2-forward 5′-3′ GTCTGTCTTCCTCTATGGCTCC, BSCL2-reverse 5′-3′ CCTTAGTCAGCGAGACATTGGC; TOP2A-forward 5′-3′ ACCATTGCAGCCTGTAAATGA and TOP2A-reverse 5′-3′ GGGCGGAGCAAAATATGTTCC. The expression level of each gene was calculated using the 2−ΔCT method.
Statistical analysis
R software (4.1.1) and relevant packages were used for data normalization and statistical analyses. The significant standard was listed as p < 0.05. Independent prognostic factor was identified using Cox regression analyses. ROC curves determined the precision of our risk signature. The difference in OS was analyzed by log-rank tests and visualized by Kaplan-Meier curves. Our research was completed in accordance with the Declaration of Helsinki as revised in 2013.
Results
Identification of prognostic aging-related genes expressing differentially in pancreatic adenocarcinoma
The workflow was generalized in Figure 1. A total of 1,589 DEGs between normal and tumor samples were selected based on both TCGA and GTEx datasets and illustrated in the volcano map (Figure 2A). Furthermore, we acquired 24 ARGs (Supplementary Table S2) from these DEGs, and a heatmap exhibited their differential expression of them (Figure 2B). We implemented a univariate Cox regression analysis and eight ARGs with prognostic values were obtained (Figure 2C, p < 0.05). PTGS2, TOP2A, IGFBP3, HOXB7, and PLAU were dangerous, while BSCL2, UCP2, and FGFR1 were protective.
FIGURE 1. Workflow of the study.
FIGURE 2. Identification of ARGs with predictive ability. (A) Volcano map of differentially expressed genes in TCGA and GTEx. (B) Heatmap of differentially expressed ARGs. (C) Forest plot presenting the hazard ratios of the prognostic ARGs. ARGs, aging-related genes; TCGA, the cancer genome atlas; GTEx, genotype-tissue expression.
Survival and immune microenvironment of molecular subtype
Patients with PAAD from the TCGA cohort were classified into two molecular clusters by consensus clustering with the standard of highest clustering stability (Figure 3A, k = 2). OS was lower in cluster 1 (n = 58) compared with cluster 2 (n = 119) (Figure 3B, p < 0.05). Moreover, we evaluated the immune infiltration of each cluster (Figure 3C). Naive B cells (p < 0.01), plasma cells (p < 0.05), CD8 T-cells (p = 0.001), and monocytes (p < 0.05) were higher infiltrated in cluster 2, while cluster 1 possessed higher infiltrations of M0 macrophages (p < 0.001) and activated dendritic cells (p < 0.05). Expression of PD-L1 was higher in cluster 1 (Figure 3D, p < 0.05). TOP2A (p < 0.05) and PLAU (p < 0.05) were correlated with PD-L1 positively, while negative correlation was found between HOXB7 (p < 0.05) with PD-L1 (Figure 3E). In addition, a heatmap depicted the expressions of prognostic aging-related genes (Supplementary Table S3) and clinical features in two clusters (Figure 3F). We furtherly compared the clinical features between two molecular subtypes and found that cluster 1 possessed a higher grade as well as advanced stage and TNM stage of the tumor (Figure 3G).
FIGURE 3. Molecular subtypes of PAAD patients mediated by ARGs. (A) Concordance matrix of subtypes. (B) Survival analysis demonstrates that cluster 1 has poor survival. (C) Vioplot visualizes immune microenvironments of two clusters. (D) PD-L1 expressed higher in cluster 1. (E) Correlation of PD-L1 with ARGs. (F) Heatmap depicts the expressions of prognostic aging-related genes and clinical features of two clusters. (G) Comparison of clinical features in two clusters. PAAD, pancreatic adenocarcinoma; ARGs, aging-related genes; PD-L1, programmed cell death-ligand 1; *p < 0.05.
Development and verification of an aging-related gene risk signature
Eight prognostic ARGs were incorporated into the LASSO regression analysis based on the TCGA database (Figures 4A,B). Finally, seven ARGs were identified for the construction of the risk signature (Supplementary Table S4). The risk score was calculated by the sum of each ARG expression multiply each ARG coefficient, and PAAD patients from the TCGA cohort were classified into two risk subgroups by median risk score (Figure 4C). The scatter plot showed the survival rates of patients according to risk scores (Figure 4D), and the heatmap further illustrated the different expressions of seven ARGs between two groups (Figure 4E). High-risk group possessed a poor OS (Figure 4F, p < 0.001). Time-dependent ROC curves verified the precision of the risk signature at 1, 2, and 3 years (Figure 4G). All of the aforementioned results were validated in GSE57495 (Figures 5A–E) and ICGC-PACA-AU (Figures 5F–J). Taken together, our ARG risk signature could precisely predict PAAD patients’ survival outcomes.
FIGURE 4. Construction of the seven-ARG risk signature. (A,B) LASSO regression analyses for the identification of seven ARGs with the highest prognostic value. (C) PAAD patients from the TCGA cohort are sorted by the median risk score. (D) Survival conditions of PAAD patients. (E) Heatmap of seven ARGs involved in risk signature. (F) Kaplan-Meier survival curve displays that the low-risk group has a better survival outcome. (G) Time-dependent ROC curves verify the predictive precision of the risk signature. ARGs, aging-related genes; LASSO, Least absolute shrinkage, and selection operator; PAAD, pancreatic adenocarcinoma; TCGA, the cancer genome atlas; ROC, receiver operating characteristic.
FIGURE 5. Validation of the seven-ARG risk signature. (A) Grouping of PAAD patients from GSE57495 based on the same median risk score. (B) Survival status of PAAD patients. (C) Heatmap of prognostic ARGs in GSE57495. (D) PAAD patients from the high-risk group have terrible survival conditions. (E) Time-dependent ROC curves based on GSE57495. (F–J) Same validation in ICGC-PCAC-AU. ARG, aging-related gene; PAAD, pancreatic adenocarcinoma; ROC, receiver operating characteristic; ICGC, international cancer genome consortium.
Independence analysis of aging-related gene risk signature and establishment of a nomogram
We explored whether our risk signature had independently predictive value for PAAD patients from the TCGA cohort. Univariate Cox regression analysis (Figure 6A, Hazard ratio (HR): 2.902, 95% Confidence interval (CI): 1.542–5.406, p < 0.001) and multivariate Cox regression analysis (Figure 6B, HR: 3.298, 95% CI: 1.615–6.735, p = 0.001) proved it together. Multiple ROC curves verified its predictive stability compared with other clinical factors (Figure 6C). For clinical use, we constructed a nomogram combining risk scores with age, gender, grade, and stage (Figure 6D). Patients at stage III-IV (n = 7 of 177) were excluded due to the heterogeneity they brought to the results. Calibration curves were drawn and the precision of our nomogram was well confirmed at 1-, 2-, and 3 years (Figures 6E–G).
FIGURE 6. Independence and clinical use of the risk signature. (A) Univariate Cox regression analysis based on TCGA database. (B)Multivariate Cox regression analysis based on TCGA database. (C) Multiple ROC curves depict the predictive precision of risk score and clinicopathological features. (D) Nomogram combing risk score and clinical factors. (E–G) Calibration curves at 1, 2, and 3 years. TCGA, the cancer genome atlas; ROC, receiver operating characteristic.
Correlation of risk score with molecular subtype, ImmuneScore, and clinicopathological characteristics
The heatmap visualized the associations between ARG expression, risk score, cluster subtype, ImmuneScore, and clinical factor based on the TCGA cohort (Figure 7A). Significant correlations were found between risk score and ImmuneScore as well as the cluster. PAAD patients with high-risk scores took a huge proportion in cluster 1 and possessed lower ImmuneScore. For validation, we explored PAAD patients’ risk scores in different groups. Consistent with aforementioned results, patients from cluster 1 (Figure 7B, p < 0.001) and low ImmuneScore subgroup (Figure 7C, p < 0.001) retained high risk scores. Interestingly, patients in grades 3-4 (Figure 7D, p < 0.01) or Tumor (T) stage 3-4 (Figure 7E, p < 0.05) also possessed high-risk scores, which hadn’t been illustrated in the heatmap.
FIGURE 7. Connections between risk score, molecular subtypes, ImmuneScore, and clinicopathological characteristics. (A) Heatmap of ARGs expressed in different subgroups. (B–E) Risk score in different clusters, ImmuneScore, grade, and T stage. ARGs, aging-related genes; ***p < 0.001.
Immune microenvironment affected by aging-related genes
Since PAAD patients from the high-risk group achieved lower ImmuneScore, we analyzed their immune microenvironment in detail based on the TCGA cohort. As confirmed by different tests, the risk score was positively associated with M0 macrophage, M1 macrophage, neutrophil, and cancer-associated fibroblast (CAF), while the correlations of NK cell, CD8 T-cell, CD4 T-cell, and B cell with risk score were negative (Figure 8A). Immune score, stroma score, and microenvironment score were also negatively relevant to risk score, demonstrating the terrible immune microenvironment in the high-risk group. Furthermore, we conducted analyses of immune checkpoint gene expression. Expressions of TNFSF4, TNFSF9, CD80, CD70, CD44, CD274 (PD-L1), and CD276 were positively correlated with a risk score, while other significant correlations were negative (Figure 8B). These significant correlations could provide options for ICIs in clinical treatment for PAAD.
FIGURE 8. Association of risk score with immune microenvironment. (A) Association of immune infiltration with risk score. (B) Correlation of immune checkpoint genes expression with risk score. *p < 0.05, **p < 0.01, ***p < 0.001.
Analyses of pathways, tumor mutation, and half inhibitory centration
For the exploration of the potential pathways related to the risk score, we performed GSVA. Several pathways related to the development of PAAD such as WNT, VEGF, MTOR, and MAPK pathways were positively associated with a risk score, and the p53 pathway was most closely correlated with a risk score (Figure 9A). A total of 20 functionally similar genes co-expressing with key ARGs were found by GeneMANIA (Figure 9B). Mutations of seven key ARGs were investigated in 849 PAAD patients and found in 62 (7%) of them (Figure 9C). PTGS2 mutation accounted for the largest and amplification was the main form of genetic alteration. We then analyzed the mutations of the top 20 genes in two risk groups. The oncoplots illustrated that the high-risk group retained higher genetic mutation (Figures 9D,E). Furthermore, we identified that the high-risk group bore a higher TMB (Figure 9F). Survival analysis demonstrated that patients with lower TMBs got better OS compared with those with higher TMBs in either risk group (Figure 9G). Finally, we explored the ability of our risk signature to predict the chemosensitivities of several drugs for PAAD. On the basis of the prediction made by the pRRophetic algorithm, patients from the high-risk group may be more sensitive to cisplatin (Figure 9H, p < 0.001) and dasatinib (Figure 9I, p < 0.001), while axitinib (Figure 9J, p < 0.001) and vorinostat (Figure 9K, p < 0.001) may be better options for the low-risk group.
FIGURE 9. Analyses of pathways, tumor mutation, and IC50 in two risk groups. (A) GSVA performed in risk signature based on TCGA. (B) Twenty functionally similar genes co-expressing with key ARGs. (C) Mutation of key ARGs in PAAD patients. (D) Mutation of top 20 genes in the low-risk group. (E) Mutation of top 20 genes in the high-risk group. (F) TMBs of PAAD patients in two risk groups. (G) Survival analysis illustrates that PAAD patients with higher TMB have worse OS in each risk group. Chemosensitivities of (H) Cisplatin, (I) Dasatinib, (J) Axitinib, and (K) Vorinostat in two risk groups. IC50, half inhibitory centration; GSVA, gene set variation analysis; TCGA, the cancer genome atlas; ARGs, aging-related genes; TMBs, tumor mutation burdens; PAAD, pancreatic adenocarcinoma; OS, overall survival; *p < 0.05, **p < 0.01, ***p < 0.001.
Protein and mRNA expression level of aging-related genes and survival analysis
For validation of ARG expression in normal and tumor tissues, we acquired immunohistochemical results from HPA. Except for IGFBP3 which hadn’t been found in HPA, as we predicted, protein expressions of PTGS2 (Figures 10A,B, p < 0.01) and TOP2A (Figures 10C,D, p < 0.001) were higher in pancreatic tumor tissue, which was consistent with their correlation with poor survival of PAAD patients. BSCL2 (Figures 10E,F, p < 0.01), FGFR1 (Figures 10G,H, p < 0.01), and UCP2 (Figures 10I,J, p < 0.001) were expressed higher in normal pancreatic tissue, and their higher expressions were relevant to better survival outcomes. It was worth noting that UCP2 was expressed higher in tumor samples based on TCGA and GTEx, which was contrary to its immunohistochemical result. The exact result should be verified in a large cohort further. PLAU (Figure 10K,L, p < 0.01) was expressed higher both in normal and tumor tissues, but its association with worse OS was consistent with what we found before. Finally, we verified the expression of the key ARGs in normal and tumor pancreatic cell lines. Compared with HPNE, BSCL2 was expressed lower in PANC1 and CFPAC1 (Figure 10M, p < 0.001), while TOP2A was expressed higher (Figure 10N, p < 0.001).
FIGURE 10. Immunohistochemical and RT-qPCR results of key ARGs expressions with their influences on OS. (A–B) PTGS2. (C–D) TOP2A. (E–F) BSCL2. (G–H) FGFR1. (I–J) UCP2. (K–L) PLAU. (M) BSCL2 is expressed higher in the normal cell line. (N) TOP2A is expressed higher in tumor cell lines. RT-qPCR, real time quantitative polymerase chain reaction; ARGs, aging-related genes; OS, overall survival; *p < 0.05, **p < 0.01, ***p < 0.001.
Discussion
Aging, as an inevitable decline of physiological function, has shown tight correlations with the incidence and death rate of cancer. As the main risk factor of cancer, aging resulted in the generation of inflammatory mediators, which were considered a crucial part of the tumor microenvironment and a trigger of tumor (Bottazzi et al., 2018). The decline of the immune system in elder individuals also weakens their resistance to tumorigenesis (Franceschi and Campisi, 2014). In terms of immune microenvironment, a previous study has proved that M2 tumor-associated macrophages were infiltrated in the microenvironment of aged mice, which led to immunosuppression and promotion of tumors (Jackaman et al., 2013). Moreover, increased expression of PD-L1, as well as immunosenescence of effector immune cells, has also been found in aged mice (Golomb et al., 2015). However, cellular senescence also plays an antitumor role in malignancy, which might be induced as a therapeutic (Collado et al., 2005; te Poele et al., 2002). The complex correlation of aging with tumors prompts us to explore the relationship between ARGs and cancer. Yue et al. (2021) have identified an ARG signature to predict immune infiltration as well as the prognosis of colorectal cancer. Nevertheless, the survival outcome and immune microenvironment characterization influenced by ARGs have not been completely elucidated in PAAD.
In this study, combing the TCGA dataset with the GTEx dataset, we first obtained eight differentially expressed ARGs with prognostic values and then identified two molecular subtypes mediated by them. There existed distinct differences in OS and immune infiltrations between the two clusters. Compared to patients in cluster 2, patients in cluster 1 had less effector immune cell infiltrations and worse OS. Moreover, we found that cluster 1 had a higher grade and advanced stage of the tumor, which may explain its poor survival.
We also constructed and validated a risk signature to investigate ARGs’ predictive ability of prognosis and immune microenvironments of PAAD patients from the TCGA cohort. LASSO regression analysis was used to screen robust prognostic biomarkers to establish an effectively prognostic ARG risk signature. The most important seven key ARGs (PTGS2, BSCL2, TOP2A, IGFBP3, UCP2, FGFR1, PLAU) were finally identified. These key aging-related genes play different roles in aging. Senescence surveillance, immune-mediated destruction of senescence cells, is partly dependent on PTGS2 and its product PGE2 (Gonçalves et al., 2021). Seipin is a protein encoded by BSCL2. Seipin plays a crucial role in the synthesis of lipid droplets, which interact with peroxisomes, key organelles in aging, and central degenerative disease (Sánchez-Iglesias et al., 2021). TOP2A encodes DNA topoisomerases to maintain the stability of DNA replication, ensuring normal cell proliferation (Wang, 2002). IGFBP3 and UCP2 are both related to muscle phenotypes associated with sarcopenia, an age-related loss of muscle function (Pratt et al., 2019). Aging suppressor αKlotho binds to FGF23 through FGFR1, and the ternary complexes are involved in the regulation of phosphate and vitamin D homeostasis (Chen et al., 2018). PLAU, closely connected with lots of physiologic and pathologic processes, such as inflammatory, p53 signal pathway, invasion, cell proliferation, and apoptosis, is directly associated with aging and age-related diseases such as Alzheimer’s disease and complications of diabetes (Cardoso et al., 2018). Senescent cells expressed and secreted PLAU to mediate cell proliferation and apoptosis (Hildenbrand et al., 2008).
Some biomarkers involved in our gene signature have been studied in pancreatic cancer. A previous study proved that TOP2A accelerated the development of pancreatic cancer through the activation of the β-catenin pathway (Pei et al., 2018). Similarly, DGCR5 was proved that can promote pancreatic cancer by targeting miR-3163/TOP2A (Liu et al., 2021). Moreover, Mody et al. (2019) suggested increasing expression of TOP2A was found in PD-L1 positive tumors. In our study, correlation evaluation also demonstrated that TOP2A expression was positively relevant to PD-L1, implying that TOP2A could be a dangerous biomarker and a potential target of treatment for PAAD.
PAAD patients with higher risk scores had worse OS compared with patients in the low-risk group. ROC curves verified the prognostic precision of ARG risk signature for 1-, 2-, and 3-year OS. After controlling confounding factors, our risk score was proved to independently predict PAAD patients’ OS. For clinical application, we established a nomogram combining risk scores with other clinical features, which facilitated the availability of the ARG risk signature. Moreover, the risk score was remarkably relevant to the ImmuneScore, grade, and T stage of PAAD patients as well as their molecular subtypes. Taken together, the ARG risk signature had a reliably predictive value for the prognosis of PAAD patients.
The immune microenvironment is a crucial component of cancer as well as a therapeutic target for PAAD. The immune microenvironment is primarily infiltrated by stromal and immune cells, and evaluations of them are associated with clinical signatures and prognosis in PAAD (Quail and Joyce, 2013; Pitt et al., 2016). In our study, negative correlations were found between risk score and effector immune cells (such as CD4+ T-cells, CD8+ T-cells, and NK cells) infiltrations. Moreover, a higher risk score was relevant to lower immune, stroma, and estimate scores. These results certified that immune microenvironments of PAAD patients with higher risk scores were worse, which possibly explained their unfavorable OS. CAFs were important components of the immune microenvironment and several promoted while others suppressed the development of PAAD (Geng et al., 2021). We found that CAFs were positively correlated with a risk score. Therefore, it is necessary to classify the subtypes of CAFs infiltrated in the immune microenvironments of PAAD patients, and perform subgroup analyses further to verify their relationships with a risk score. Overall, our ARG risk signature could be used for immune microenvironment stratification of patients with PAAD.
In addition, the risk score was also relevant to the expressions of several immune checkpoint genes positively. Inhibitor of PD-L1 was recognized as effective antitumor immunotherapy for different kinds of cancer, including pancreatic cancer. Moreover, CD44 was expressed higher in several solid tumors, such as pancreatic cancer, which was considered an important target for therapeutic (Mattheolabakis et al., 2015). CD276 also retained high expression in pancreatic cancer and was associated with worse OS (Inamura et al., 2018). Therefore, our research provided potential immunotherapeutic targets for PAAD patients based on their risk score, and targeting aging could serve as a valuable regulation for immunotherapy in PAAD.
Aging is correlated with several biological processes of various tumors, including PAAD (Campisi, 2013; Campisi and d’Adda di Fagagna, 2007). We further explored potential pathways related to risk score positively by GSVA. Several molecules promoted tumorigenesis, progression, and chemoresistance of pancreatic cancer through activation of the wnt/β-catenin pathway (Zhou et al., 2020; Chen et al., 2021). Oxymatrine played a role in antiangiogenesis in PAAD by inhibiting the VEGF pathway (Chen et al., 2013). Moreover, the mTOR pathway was also activated in pancreatic cancer to the promotion of progression and drug resistance (Mann et al., 2016). The tightest link was found between the p53 pathway and risk score, and inactivation of p53 was found in almost every tumor (Joerger and Fersht, 2016). Zhang et al. (2018) have proved that blockade of cell proliferation and acceleration of cell senescence could be triggered by CCNB1 silencing through activation of the p53 pathway. These correlations between risk scores and pathways could offer several potential therapeutic targets to PAAD patients, while more details should be validated by fundamental experiments. In addition, we found 20 functionally similar genes of key ARGs. Their potential functions and pathways needed to be explored further.
The accumulation of gene mutation leads to tumorigenesis. Thus, we first investigated the mutations of key ARGs in PAAD patients and found the mutations were rare. Moreover, we found significant differences in genetic alterations among the two risk groups. Higher genetic mutations were found in patients from the high-risk group, such as KRAS, TP53, SMAD4, and CDKN2A, which were major drivers of PAAD (Jones et al., 2008; Waddell et al., 2015; Makohon-Moore and Iacobuzio-Donahue, 2016). Patients with high-risk scores bore higher TMBs. The correlation of higher TMB with better efficacy of immunotherapy has been proved (Samstein et al., 2019). Thus, patients from the high-risk group may benefit from immunotherapy, which should be confirmed in clinical further. In addition to the poor survival in the high-risk group, higher TMB was also related to worse OS. However, TMB could not be an effectively prognostic biomarker on its own (Addeo et al., 2021). More clinical trials should be implemented to explore the relationship between TMB and OS.
PAAD shows resistance to chemotherapy due to complex factors (Springfeld et al., 2019). Therefore, we explored the potential sensitive drugs for PAAD patients from different risk subgroups. Since patients with different risk scores were predicted to be sensitive to different drugs, our risk signature could offer several therapeutic options for them. The final effect still needed to be verified in clinical trials.
Finally, we used immunohistochemical results from HPA to confirm the expressions of key ARGs. Their expressions in normal and tumor samples were in keeping with their influences on OS. However, UCP2 showed a contrary tendency in expression compared with that in the TCGA dataset. Moreover, the result of IGFBP3 was not found in HPA, and PLAU didn’t express differentially in normal and tumor samples. In the future, we hope to validate their expressions ourselves. The results of RT-qPCR verified the different mRNA levels of key aging-related genes in normal and tumor cells, which was consistent with our previous analyses.
There were some limitations existing in our study. On the one hand, more experiments in vivo and vitro were needed for verification of the results we identified by bioinformatic analyses. We hope to conduct functional validation of the prognostic ARGs we found in the future. On the other hand, PAAD samples from TCGA, ICGC, and GEO were deficient, thus we needed more samples from a large cohort to verify our risk signature further.
Conclusion
We identified two molecular subtypes and established a novel ARG risk signature based on ARGs. According to the risk score generated from our signature, we predicted the prognoses, immune microenvironments, tumor mutations, and drug sensitivity of PAAD patients. Our study may offer a novel understanding of the molecular mechanism in PAAD as well as new targets and options for the clinical treatment of PAAD.
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
Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author contributions
LW and JC conceived and designed the study. LW provided administrative and funding support. SX, WG, and KW collected and analyzed the data. TM, XZ, HX, YoW, JY, SL, MY, JM, YaW, and DS provided study materials. SX, WG, and KW drafted the manuscript. All authors reviewed and made the final approval of the manuscript.
Funding
This work was supported by the Innovation Group Project of Shanghai Municipal Health Commission (2019CXJQ03), National Natural Science Foundation of China (81874048, 82171824), Shanghai Municipal Commission of Health and Family Planning (2018ZHYL0223), Fostering Fund of Renji Hospital Affiliated to Shanghai Jiao Tong University School of Medicine (PYIV-17-001), Clinical Research Plan of SHDC (No. SHDC2020CR1035B), Scientific and Technological Innovation Project of Science and Technology Commission of Shanghai Municipality (21JC1404300), Project from CSCO Clinical Oncology Research Foundation (Y-2019AZZD-0513), Shanghai Key Clinical Speciality (Oncology), Shanghai Leading Talents Project, and Innovative Research Team of High-level Local Universities in Shanghai.
Acknowledgments
We wish to appreciate medical staff contributing to organizing and updating the TCGA, GTEx, GEO, ICGC, HAGR, HPA, and cBioPortal datasets. We thank for everyone who participates in this study.
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/fcell.2022.942225/full#supplementary-material
Abbreviations
PAAD, pancreatic adenocarcinoma; OS, overall survival; ICIs, immune checkpoint inhibitors; ARGs, aging-related genes; TCGA, the cancer genome atlas; GTEx, genotype-tissue expression; GEO, gene expression omnibus; ICGC, international cancer genome consortium; HAGR, human aging genomic resources; HPA, human protein atlas; DEGs, differentially expressed genes; PD-L1, programmed cell death-ligand 1; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; GSVA, gene set variation analysis; TMB, tumor mutation burden; IC50, half inhibitory centration; RT-qPCR, real time quantitative polymerase chain reaction; HR, hazard ratio; CI, confidence interval; CAF, cancer associated fibroblast
References
Addeo, A., Friedlaender, A., Banna, G. L., and Weiss, G. J. (2021). TMB or not TMB as a biomarker: that is the question. Crit. Rev. Oncol. Hematol. 163, 103374. doi:10.1016/j.critrevonc.2021.103374
Berben, L., Floris, G., Wildiers, H., and Hatse, S. (2021). Cancer and aging: Two tightly interconnected biological processes. Cancers (Basel) 13 (6), 1400. doi:10.3390/cancers13061400
Bottazzi, B., Riboli, E., and Mantovani, A. (2018). Aging, inflammation and cancer. Semin. Immunol. 40, 74–82. doi:10.1016/j.smim.2018.10.011
Camous, X., Pera, A., Solana, R., and Larbi, A. (2012). NK cells in healthy aging and age-associated diseases. J. Biomed. Biotechnol. 2012, 195956. doi:10.1155/2012/195956
Campisi, J. (2013). Aging, cellular senescence, and cancer. Annu. Rev. Physiol. 75, 685–705. doi:10.1146/annurev-physiol-030212-183653
Campisi, J., and D'Adda Di Fagagna, F. (2007). Cellular senescence: when bad things happen to good cells. Nat. Rev. Mol. Cell Biol. 8, 729–740. doi:10.1038/nrm2233
Cardoso, A. L., Fernandes, A., Aguilar-Pimentel, J. A., De Angelis, M. H., Guedes, J. R., Brito, M. A., et al. (2018). Towards frailty biomarkers: candidates from genes and pathways regulated in aging and age-related diseases. Ageing Res. Rev. 47, 214–277. doi:10.1016/j.arr.2018.07.004
Chen, G., Liu, Y., Goetz, R., Fu, L., Jayaraman, S., Hu, M. C., et al. (2018). α-Klotho is a non-enzymatic molecular scaffold for FGF23 hormone signalling. Nature 553, 461–466. doi:10.1038/nature25451
Chen, H., Zhang, J., Luo, J., Lai, F., Wang, Z., Tong, H., et al. (2013). Antiangiogenic effects of oxymatrine on pancreatic cancer by inhibition of the NF-κB-mediated VEGF signaling pathway. Oncol. Rep. 30, 589–595. doi:10.3892/or.2013.2529
Chen, H., Zhao, J., Jiang, N., Wang, Z., and Liu, C. (2021). Hyperglycemia promotes pancreatic cancer initiation and progression by activating the wnt/β-catenin signaling pathway. Anticancer. Agents Med. Chem. 21, 2592–2602. doi:10.2174/1871520621666210201095613
Clark, C. E., Hingorani, S. R., Mick, R., Combs, C., Tuveson, D. A., Vonderheide, R. H., et al. (2007). Dynamics of the immune reaction to pancreatic cancer from inception to invasion. Cancer Res. 67, 9518–9527. doi:10.1158/0008-5472.CAN-07-0175
Collado, M., Gil, J., Efeyan, A., Guerra, C., Schuhmacher, A. J., Barradas, M., et al. (2005). Tumour biology: senescence in premalignant tumours. Nature 436, 642. doi:10.1038/436642a
Franceschi, C., and Campisi, J. (2014). Chronic inflammation (inflammaging) and its potential contribution to age-associated diseases. J. Gerontol. A Biol. Sci. Med. Sci. 69 (Suppl. 1), S4–S9. doi:10.1093/gerona/glu057
Geng, X., Chen, H., Zhao, L., Hu, J., Yang, W., Li, G., et al. (2021). Cancer-associated fibroblast (CAF) heterogeneity and targeting therapy of CAFs in pancreatic cancer. Front. Cell Dev. Biol. 9, 655152. doi:10.3389/fcell.2021.655152
Golomb, L., Sagiv, A., Pateras, I. S., Maly, A., Krizhanovsky, V., Gorgoulis, V. G., et al. (2015). Age-associated inflammation connects RAS-induced senescence to stem cell dysfunction and epidermal malignancy. Cell Death Differ. 22, 1764–1774. doi:10.1038/cdd.2015.21
Gonçalves, S., Yin, K., Ito, Y., Chan, A., Olan, I., Gough, S., et al. (2021). COX2 regulates senescence secretome composition and senescence surveillance through PGE(2). Cell Rep. 34, 108860. doi:10.1016/j.celrep.2021.108860
Hildenbrand, R., Gandhari, M., Stroebel, P., Marx, A., Allgayer, H., Arens, N., et al. (2008). The urokinase-system--role of cell proliferation and apoptosis. Histol. Histopathol. 23, 227–236. doi:10.14670/HH-23.227
Inamura, K., Takazawa, Y., Inoue, Y., Yokouchi, Y., Kobayashi, M., Saiura, A., et al. (2018). Tumor B7-H3 (CD276) expression and survival in pancreatic cancer. J. Clin. Med. 7, E172. doi:10.3390/jcm7070172
Jackaman, C., Radley-Crabb, H. G., Soffe, Z., Shavlakadze, T., Grounds, M. D., Nelson, D. J., et al. (2013). Targeting macrophages rescues age-related immune deficiencies in C57BL/6J geriatric mice. Aging Cell 12, 345–357. doi:10.1111/acel.12062
Jiang, P., Yang, F., Zou, C., Bao, T., Wu, M., Yang, D., et al. (2021). The construction and analysis of a ferroptosis-related gene prognostic signature for pancreatic cancer. Aging (Albany NY) 13, 10396–10414. doi:10.18632/aging.202801
Joerger, A. C., and Fersht, A. R. (2016). The p53 pathway: Origins, inactivation in cancer, and emerging therapeutic approaches. Annu. Rev. Biochem. 85, 375–404. doi:10.1146/annurev-biochem-060815-014710
Jones, S., Zhang, X., Parsons, D. W., Lin, J. C., Leary, R. J., Angenendt, P., et al. (2008). Core signaling pathways in human pancreatic cancers revealed by global genomic analyses. Science 321, 1801–1806. doi:10.1126/science.1164368
Liu, S. L., Cai, C., Yang, Z. Y., Wu, Z. Y., Wu, X. S., Wang, X. F., et al. (2021). DGCR5 is activated by PAX5 and promotes pancreatic cancer via targeting miR-3163/TOP2A and activating Wnt/β-catenin pathway. Int. J. Biol. Sci. 17, 498–513. doi:10.7150/ijbs.55636
López-Otín, C., Blasco, M. A., Partridge, L., Serrano, M., and Kroemer, G. (2013). The hallmarks of aging. Cell 153, 1194–1217. doi:10.1016/j.cell.2013.05.039
Makohon-Moore, A., and Iacobuzio-Donahue, C. A. (2016). Pancreatic cancer biology and genetics from an evolutionary perspective. Nat. Rev. Cancer 16, 553–565. doi:10.1038/nrc.2016.66
Mann, K. M., Ying, H., Juan, J., Jenkins, N. A., and Copeland, N. G. (2016). KRAS-related proteins in pancreatic cancer. Pharmacol. Ther. 168, 29–42. doi:10.1016/j.pharmthera.2016.09.003
Mattheolabakis, G., Milane, L., Singh, A., and Amiji, M. M. (2015). Hyaluronic acid targeting of CD44 for cancer therapy: From receptor biology to nanomedicine. J. Drug Target. 23, 605–618. doi:10.3109/1061186X.2015.1052072
Mizrahi, J. D., Surana, R., Valle, J. W., and Shroff, R. T. (2020). Pancreatic cancer. Lancet 395, 2008–2020. doi:10.1016/S0140-6736(20)30974-0
Mody, K., Starr, J., Saul, M., Poorman, K., Weinberg, B. A., Salem, M. E., et al. (2019). Patterns and genomic correlates of PD-L1 expression in patients with biliary tract cancers. J. Gastrointest. Oncol. 10, 1099–1109. doi:10.21037/jgo.2019.08.08
Okasha, H., Elkholy, S., El-Sayed, R., Wifi, M. N., El-Nady, M., El-Nabawi, W., et al. (2017). Real time endoscopic ultrasound elastography and strain ratio in the diagnosis of solid pancreatic lesions. World J. Gastroenterol. 23, 5962–5968. doi:10.3748/wjg.v23.i32.5962
Pei, Y. F., Yin, X. M., and Liu, X. Q. (2018). TOP2A induces malignant character of pancreatic cancer through activating β-catenin signaling pathway. Biochim. Biophys. Acta. Mol. Basis Dis. 1864, 197–207. doi:10.1016/j.bbadis.2017.10.019
Pitt, J. M., Marabelle, A., Eggermont, A., Soria, J. C., Kroemer, G., and Zitvogel, L. (2016). Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann. Oncol. 27, 1482–1492. doi:10.1093/annonc/mdw168
Pratt, J., Boreham, C., Ennis, S., Ryan, A. W., and De Vito, G. (2019). Genetic associations with aging muscle: a systematic review. Cells 9, E12. doi:10.3390/cells9010012
Quail, D. F., and Joyce, J. A. (2013). Microenvironmental regulation of tumor progression and metastasis. Nat. Med. 19, 1423–1437. doi:10.1038/nm.3394
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, 202–206. doi:10.1038/s41588-018-0312-8
Sánchez-Iglesias, S., FERNáNDEZ-Pombo, A., Cobelo-GóMEZ, S., Hermida-Ameijeiras, Á., ALARCóN-MARTíNEZ, H., Domingo-JIMéNEZ, R., et al. (2021). Celia's encephalopathy (BSCL2-Gene-Related): current understanding. J. Clin. Med. 10, 1435. doi:10.3390/jcm10071435
Sarsour, E. H., Son, J. M., Kalen, A. L., Xiao, W., Du, J., Alexander, M. S., et al. (2020). Arachidonate 12-lipoxygenase and 12-hydroxyeicosatetraenoic acid contribute to stromal aging-induced progression of pancreatic cancer. J. Biol. Chem. 295, 6946–6957. doi:10.1074/jbc.RA120.012798
Schizas, D., Charalampakis, N., Kole, C., Economopoulou, P., Koustas, E., Gkotsis, E., et al. (2020). Immunotherapy for pancreatic cancer: a 2020 update. Cancer Treat. Rev. 86, 102016. doi:10.1016/j.ctrv.2020.102016
Seiler, M., Huang, C. C., Szalma, S., and Bhanot, G. (2010). ConsensusCluster: a software tool for unsupervised cluster discovery in numerical data. Omics 14, 109–113. doi:10.1089/omi.2009.0083
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2022). Cancer statistics, 2022. Ca. Cancer J. Clin. 72 (1), 7–33. doi:10.3322/caac.21708
Song, W., He, X., Gong, P., Yang, Y., Huang, S., Zeng, Y., et al. (2021). Glycolysis-related gene expression profiling screen for prognostic risk signature of pancreatic ductal adenocarcinoma. Front. Genet. 12, 639246. doi:10.3389/fgene.2021.639246
Springfeld, C., JäGER, D., BüCHLER, M. W., Strobel, O., Hackert, T., Palmer, D. H., et al. (2019). Chemotherapy for pancreatic cancer. Presse Med. 48, e159–e174. doi:10.1016/j.lpm.2019.02.025
Te Poele, R. H., Okorokov, A. L., Jardine, L., Cummings, J., and Joel, S. P. (2002). DNA damage is able to induce senescence in tumor cells in vitro and in vivo. Cancer Res. 62, 1876–1883.
Waddell, N., Pajic, M., Patch, A. M., Chang, D. K., Kassahn, K. S., Bailey, P., et al. (2015). Whole genomes redefine the mutational landscape of pancreatic cancer. Nature 518, 495–501. doi:10.1038/nature14169
Wang, J. C. (2002). Cellular roles of DNA topoisomerases: a molecular perspective. Nat. Rev. Mol. Cell Biol. 3, 430–440. doi:10.1038/nrm831
Wang, W., Yan, L., Guan, X., Dong, B., Zhao, M., Wu, J., et al. (2020). Identification of an immune-related signature for predicting prognosis in patients with pancreatic ductal adenocarcinoma. Front. Oncol. 10, 618215. doi:10.3389/fonc.2020.618215
Yue, T., Chen, S., Zhu, J., Guo, S., Huang, Z., Wang, P., et al. (2021). The aging-related risk signature in colorectal cancer. Aging (Albany NY) 13, 7330–7349. doi:10.18632/aging.202589
Zhang, H., Zhang, X., Li, X., Meng, W. B., Bai, Z. T., Rui, S. Z., et al. (2018). Effect of CCNB1 silencing on cell cycle, senescence, and apoptosis through the p53 signaling pathway in pancreatic cancer. J. Cell. Physiol. 234, 619–631. doi:10.1002/jcp.26816
Zhou, C., Yi, C., Yi, Y., Qin, W., Yan, Y., Dong, X., et al. (2020). LncRNA PVT1 promotes gemcitabine resistance of pancreatic cancer via activating Wnt/β-catenin and autophagy pathway through modulating the miR-619-5p/Pygo2 and miR-619-5p/ATG14 axes. Mol. Cancer 19, 118. doi:10.1186/s12943-020-01237-y
Keywords: pancreatic adenocarcinoma, aging, molecular subtype, prognostic signature, immune microenviroment
Citation: Xue S, Ge W, Wang K, Mao T, Zhang X, Xu H, Wang Y, Yao J, Li S, Yue M, Ma J, Wang Y, Shentu D, Cui J and Wang L (2022) Association of aging-related genes with prognosis and immune infiltration in pancreatic adenocarcinoma. Front. Cell Dev. Biol. 10:942225. doi: 10.3389/fcell.2022.942225
Received: 12 May 2022; Accepted: 30 June 2022;
Published: 08 August 2022.
Edited by:
Zuoren Yu, Tongji University, ChinaCopyright © 2022 Xue, Ge, Wang, Mao, Zhang, Xu, Wang, Yao, Li, Yue, Ma, Wang, Shentu, Cui and Wang. 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: Jiujie Cui, cuijiujie@126.com; Liwei Wang, liweiwang@shsmu.edu.cn
†These authors have contributed equally to this work and share first authorship