Lung adenocarcinoma datasets GSE31210 (N = 226) (Okayama et al., 2012), GSE37745 (N = 106) (Botling et al., 2013), GSE50081 (N = 180) (Der et al., 2014), and GSE72094 (N = 398) (Schabath et al., 2016) were downloaded from the NCBI gene expression omnibus (GEO) database. The samples with missing survival data were removed and all data were log2 (x+1) normalized and corrected. The ComBat method in the “SVA” package was used to eliminate batch effects among the four GEO datasets and to integrate the four GEO datasets for subsequent analyses. All datasets contained a total of 857 lung adenocarcinoma sample data containing clinical information such as survival data, gender, age, pathological stage, and smoking history. The lung adenocarcinoma transcriptome expression dataset was downloaded through the cancer genome atlas (TCGA) database containing a total of 515 lung adenocarcinoma samples (different samples from the same patient were homogenized by taking the mean value) and 59 normal tissue samples; all expression data were downloaded in the fragments per kilobase of exon per million mapped fragments format and converted to the transcripts per kilobase million (TPM) format with log2 (TPM+1) data normalization operation. The survival data of 500 tumor samples was known, of which the data of samples with a survival time of 0 was removed and all others were rounded to survival time in years. Mutation data and copy number variation (CNV) data for lung adenocarcinoma were downloaded from the TCGA Data Portal (https://tcga-data.nci.nih.gov/tcga/).

By reviewing the literature related to SUMOylation, we identified 33 genes recognized as critical regulators in SUMOylation (Chang and Yeh, 2020). Subsequently, we performed unsupervised consensus clustering using the “ConsensusClusterPlus” package based on the mRNA expression of the 33 SUMOylation regulators. We use the “km” clustering algorithm, the distance calculation algorithm is “euclidean,” and the random seed is set to “123,456". The optimal clusters were determined under the cumulative distribution curve (CDF), and the rationality of the clusters was further verified using principal component analysis (PCA). Consensus clustering was adopted with 1,000,000 replicates to optimize clustering results. The clustering method for unsupervised consensus clustering analysis based on the expression of differentially expressed key SUMOylation genes is consistent with the above.

We used the R package GSVA (Hanzelmann et al., 2013) to study the biological processes associated with different SUMOylation modification patterns. Hallmark gene sets (Liberzon et al., 2015) and well-defined gene sets of biological features were downloaded from MsigDB (http://www.gsea-msigdb.org/gsea/msigdb/). The gene set of tumor-associated biological pathways was were obtained from the compilation by Mariathasan et al. (Mariathasan et al., 2018) and included the following: angiogenesis, immune checkpoint, cell cycle regulators, Pan fibroblast transforming growth factor beta (Pan F TBRs), epithelial-mesenchymal transition processes involved in EMT1, EMT2, EMT3, and the cell cycle. The immune checkpoint genes (including immune co-stimulatory and immune co-inhibitory molecules) and MHC molecules were derived from the immune-related gene set compiled by Rooney et al. (Rooney et al., 2015).

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases are used to annotate genes and their RNA or protein products to determine their unique biological properties. GO annotation and KEGG annotation of SUMOylation-related genes was performed using the R package “clusterProfiler” with a cutoff value of false discovery rate (FDR<0.05).

The gene set for each immune cell-specific marker was compiled from a recent study, and Single Sample Gene Set Enrichment Analysis (ssGSEA) (Barbie et al., 2009) was used to quantify the relative abundance of 24 immune cell types in the TME and was expressed as an enrichment fraction (Bindea et al., 2013). Subsequently, the EPIC immune infiltration algorithm (Racle et al., 2017) and the TIMER immune algorithm (Li et al., 2017) were used to calculate the extent of infiltration of each immune cell in lung adenocarcinoma tissue. Tumor tissue with abundant immune cell infiltration is consistent with a higher immune score and lower tumor purity.

The R package “limma” was used to measure DEGs associated with different SUMOylation patterns in lung adenocarcinoma samples, with the significance filtering criteria for DEGs of |FC| > 1.5 and FDR <0.05. The interaction of DEGs was investigated using the STRING database and used to build PPI networks, and the screening network type was a physical subnetwork with a confidence score ≥0.7. We used the “CytoHubba” plugin (Chin et al., 2014) in “Cytoscape” software (Shannon et al., 2003) to calculate the number of neighboring nodes in the PPI network, and filtered genes with neighboring nodes ≥5 as SUMOylation-related key genes.

We performed a univariate cox regression analysis using DEGs associated with SUMOylation modification patterns. Next, we screened the significant prognostic correlates (p < 0.05) among them for the optimal variables for constructing SUMOylation scores using the least absolute shrinkage and selection operator (LASSO) regression method. The SUMO score was calculated as follows: RiskScore = Exp (Gene1)*β1+Exp (Gene2)*β2+. . . + Exp (Genen)*βn. Where Exp is the expression of the genes and β is the regression coefficient calculated by LASSO. Patients were then divided into the high-risk and low-risk groups according to the optimal cut-off value calculated using the surv-cutpoint function in the R package “survival”. Prognostic analyses were conducted using the R package “survminer”. Receiver operating characteristic (ROC) curves were used to evaluate the prognostic performance of the SUMOylation score and the nomogram model. The area under the curve (AUC) was calculated using the R package “timeROC,” while multiple regression analysis was used to verify the validity of the score as an independent prognostic factor for lung adenocarcinoma.

The immunotherapy cohort GSE126044 (Cho et al., 2020) contained a total of 16 lung adenocarcinoma recipients of nivolumab (anti-PD1 drug) treatment, five of which were responsive to (CR and PR) and 11 non-responsive to immunotherapy (PD and SD). (Tumor Immune Dysfunction and Exclusion) TIDE (http://tide.dfci.harvard.edu/) is applied to predict immunotherapy response for lung adenocarcinoma patients. Gene expression data from pretreatment tissue samples were collated and transformed into TPM format with log2 (TPM+1) normalization for further analysis. Chemotherapy drug sensitivity analysis was performed using the R package “pRRophetic” (Geeleher et al., 2014) to predict clinical response to chemotherapy.

Eighteen pairs of LUAD tissues and adjacent non-cancerous lung tissues were obtained from the First Affiliated Hospital of Soochow University after informed consent from patients. Pathological diagnostics for patients with LUAD were assessed according to the Revised International System for Staging Lung Cancer. All patients included in this study did not receive any radiation therapy, chemotherapy, or immunotherapy. The experiment was approved by the Academic Advisory Board of Soochow University.

Total RNA was isolated using Trizol reagent (Beyotime, China) and reverse transcribed into cDNA by Hiscript III Reverse Transcriptase (Vazyme, China). Real-Time PCR was performed using ChamQ SYBR qPCR Master Mix (Vazyme, China) on a LightCycler96 real-time PCR system (Roche, Switzerland). The specific primers we used are listed in Supplementary Table S1.

All statistical analyses in this study were performed using the R software 4.1.2 or GraphPad Prism 9. For quantitative data, normally distributed variables were analyzed using the Student’s t-test, and non-normally distributed variables were analyzed using the Wilcoxon rank-sum test. Statistical significant differences between three and more data sets were analyzed using the Kruskal–Wallis method for non-parametric statistical tests and one-way ANOVA for parametric statistical tests. We calculated correlation coefficients using Spearman and distance correlation analysis. Survival analyses were conducted using the Kaplan–Meier method and log-rank tests, while Cox proportional risk regression models were used to analyze the relationship between SUMOylation patterns and regulatory genes and prognosis. All statistical comparisons in this study were two-sided with α = 0.05 and the Benjamini-Hochberg method was employed to control the false discovery rate (FDR) for multiple hypothesis testing. *p < 0.05, **p < 0.01, ***p < 0.001.

The flow chart shows our research procedure (Figure 1). We investigated the role of 33 SUMOylation-modified regulatory genes in lung adenocarcinoma, namely, SUMO isoforms (SUMO1, SUMO2, SUMO3, SUMO4), SUMO-Activating Enzyme (SAE1, UBA2), SUMO-Conjugating Enzyme (UBE2I), SUMO E3 Ligases (BCL11A, CAPN3, CBX4, HDAC4, HDAC7, MUL1, NSMCE2, PIAS1, PIAS2, PIAS3, PIAS4, RANBP2, RANGAP1, RNF212, RWDD3, TOPORS, TRIM27, TRIM28, ZMIZ1), SUMO proteases (HINT1, SENP1, SENP2, SENP3, SENP6, SENP5, SENP7, USPL1). The reversible process of SUMOylation is regulated by the aforementioned regulatory genes (Figure 2A). We first performed GO and KEGG enrichment analysis of the 33 SUMOylation regulatory genes. GO enrichment revealed that the 33 SUMOylation-regulated genes were remarkably enriched in biological processes such as lysine modification, protein SUMOylation, regulation of small protein covalent protein modification (Supplementary Figure S1A), and molecular functions such as ubiquitinated protein convertase activity, SUMO convertase activity, among others (Supplementary Figure S1B). KEGG enrichment analysis revealed other SUMOylation-related pathways (Supplementary Figure S1C).

Subsequently, for the comprehensive understanding of the SUMOylation regulatory genes, we performed the Spearman correlation analysis of the 33 SUMOylation regulatory genes to assess their pattern of co-expression (Supplementary Figure S1D) and found a broad significant negative correlation between the expression of CAPN3 and PIAS1 and other regulatory genes, as well as a significant positive correlation between CAPN3 and PIAS1 expression. Notably, most of the SUMOylation regulatory genes were remarkably upregulated in tumor tissues, and overall, the expression of SUMOylation regulatory genes was remarkably different between the tumor and normal samples (Figure 2B). PCA based on paired tumor and normal samples identified the SUMOylation regulatory genes whose expression was significantly different between normal and tumor samples (Supplementary Figure S1E).

Moreover, CNVs were very common in the 33 regulatory genes; CNV amplifications were prevalent in PIAS3, NSMCE2, UBA2, CBX4, and SUMO2, while CNV deletions were prevalent in RWDD3, USPL1, SAE1, and CAPN3 (Figure 2C). The CNV occurred widely in genes on a variety of chromosomes (Figure 2D); however, they were concentrated mainly on chromosomes 1, 3, 4, 6, 12, 15, and 19. We further evaluated the prevalence of somatic mutations in SUMOylation regulatory genes. Of the 561 lung adenocarcinoma samples, 126 (22.46%) had alterations, mainly missense mutations and non-sense mutations, in SUMOylation regulatory genes. RANBP2 and BCL11A had the highest mutation frequency of 4%, followed by HDAC4, TOPORS, and SENP7 with 2% mutation frequency; most mutations occurred in the regulatory genes SUMO activator and E3 ligase, and less in SUMO isoforms (Figure 2E). Subsequently, we analyzed the mutational co-occurrence of the SUMOylation regulatory genes and found that SUMO2 co-occurred with mutations in both RANGAP1 and PIAS4 (Supplementary Figure S1F).

Furthermore, to determine the relationship between the 33 SUMOylation regulatory genes and the prognosis of lung adenocarcinoma patients, we performed a cox risk proportional regression model was used to determine the relationship between the 33 SUMOylation regulatory genes and the prognosis of lung adenocarcinoma patients, and the forest plot revealed that CAPN3, SENP6, SENP7, RANBP2, HINT1, BCL11A, TOPORS, and HDAC4 were protective factors for lung adenocarcinoma and were generally downregulated, while RANGAP1, SENP5, SAE1, SENP2, TRIM28, and TRIM27 were risk factors for lung adenocarcinoma and were upregulated in lung adenocarcinoma (Figure 2F). The multivariate Cox proportional Hazards regression analysis further established that SUMOylation regulatory genes were significantly associated with lung adenocarcinoma prognosis (Supplementary Figure S1G). In conclusion, we mapped the relationship between SUMOylation regulatory genes and lung adenocarcinoma prognosis (Figure 2G). In summary, we determined that SUMOylation regulatory genes differed significantly between normal and lung adenocarcinoma tissues and were also markedly correlated with the frequency of mutations and CNVs in lung adenocarcinoma. In addition, we demonstrated that altered expression and genetic variation of specific SUMOylation regulatory genes have a critical role in the development and prognosis of lung adenocarcinoma.

For exploring different SUMOylation patterns, we collated 857 lung adenocarcinoma patients in four datasets of the GEO database for unsupervised consistency clustering based on 33 SUMOylation regulatory genes. Based on the variation of the area under the CDF curve (Figures 3A, B; Supplementary Figure S2A), we determined the optimal number of clusters to be 2 (K = 2). We subsequently performed PCA of the two SUMOylation patterns based on mRNA expression data of the regulatory genes (Figure 3C) and found that the two different modification patterns could be completely distinguished by SUMOylation regulatory genes. Meanwhile, Kaplan–Meier survival analysis revealed that patients with the SUMO-C2 pattern showed a worse prognostic survival (log-Rank test p = 7.55e-07) (Figure 3D). The occurrence of the SUMO-C2 pattern was consistent with the downregulation of BCL11A, CAPN3, HDAC4, PIAS4, UBA2, UBE2I, SENP6, SENP7, and USPL1, and the upregulation of CBX4, NSMCE2, PIAS3, RANGAP1, TRIM28 and SENP2 (Figure 3E). The Heatmap further confirmed that the expression of the 33 SUMOylation-regulated genes was significantly different in the two SUMOylation patterns (Supplementary Figure S2B). The multivariate Cox regression results indicated, in agreement with Kaplan–Meier results, that the SUMOylation patterns can be an independent prognostic factor, and compared with SUMO-C1, SUMO-C2 was associated with poor prognostic survival, HR = 1.55, p < 0.001 (Figure 3F).

To further explore the differences in key biological pathways associated with the two SUMOylation patterns. The GSVA algorithm was used to calculate the Hallmark gene set scores and convert them into a scoring matrix of the gene set. We found that there were significant differences in GSVA scores of Hallmark gene sets consistent with SUMO-C1 and SUMO-C2 (Figure 3G), and SUMO-C2 were mainly enriched in biological pathways related to the cell cycle, glucose metabolism, and genetic material replication and repair, such as G2M checkpoint, E2F targets, DNA repair, Glycolysis, and PI3K-AKT-mTOR signaling pathways and was highly active in biological processes related to cancer progression such as hypoxia, angiogenesis, and epithelial-mesenchymal transition (Supplementary Figure S2C). In contrast, SUMO-C1 was mainly enriched in tumor immune-related processes, such as TNFA signaling via NF-κB, IL6-STAT3 signaling, complement, inflammatory response, and IL2-STAT5 signaling, allograft rejection signaling pathway. This suggests that the SUMOylation pattern based on SUMOylation regulatory genes may influence the immune microenvironment status. Biological mechanisms for the poor prognosis associated with SUMO-C2 maybe because of the active immune microenvironment of the SUMO-C1 and the active stromal microenvironment of the SUMO-C2. SUMO-C1 may present better antitumor effects and better prognostic survival in lung adenocarcinoma through high activation of antitumor immunity and inhibition of the pro-cancer progressive stromal microenvironment.

To explore the potential mechanisms by which SUMOylation patterns regulate the immune microenvironment, we used the ssGSEA algorithm to calculate the abundance of 24 different immune cells. A heatmap was used for visualizing the difference in immune cell infiltration (Figure 4A), and the Wilcoxon rank-sum test was used for comparing the discrepancy in immune cell infiltration consistent with SUMO-C1 and SUMO-C2 patterns (Figure 4B). Notably, SUMO-C1 showed a widespread active state of immune cells, while SUMO-C2 showed an inactive state of immune cells. The infiltration of B cells, T cells, T helper cells, Central Memory T cell (Tcm), Effector memory T cells (Tem), Follicular helper T cell (TFH), Th17 cells, CD8 T cells, cytotoxic cells, dendritic cells (DCs), and mast cells were significantly higher in the SUMO-C1 than SUMO-C2 group. In contrast, immune cells such as Th2 cells and Treg cells were significantly activated in SUMO-C2. Thereafter, we calculated immune cell infiltration levels in lung adenocarcinoma using the EPIC algorithm and TIMER algorithm to further support the differences in immune cell infiltration levels and compared the significance of the differences between the SUMO-C1 and SUMO-C2 groups using the Wilcoxon rank-sum test. We found that the SUMO-C1 group exhibited higher levels of immune cell infiltration, particularly CD8 T cells, which has been reported to have a killing effect on tumor cells thereby producing protective immunity against tumors (Supplementary Figures S3A, S3B). Furthermore, Spearman correlation between each SUMOylation regulatory gene and immune cell revealed that CAPN3 was widely associated with various immune cells such as CD8+T cells (r = 0.32), DCs (r = 0.343), iDCs (r = 0.432), Mast cells (r = 0.468), Tcm (r = 0.421), Tem (r = 0.275), TFH (r = 0.291), Th2 cells (r = 0.563), and Treg (r = 0.115), indicating that CAPN3 among SUMOylation regulatory genes is closely associated with immune cell infiltration in the tumor microenvironment of lung adenocarcinoma. Overall, the immune microenvironment was regulated by SUMOylation regulatory genes in lung adenocarcinoma mainly by modulating the levels of immune cells such as CD8+T cells, mast cells, T helper cells, Tcm, Tem, and Th2 cells (Figure 4C; Supplementary Figures S3A).

The above analysis supports that SUMO-C1 represents a type of TME with an immune characteristic of tumor suppression. Nevertheless, SUMO-C2 exhibited predominantly pro-cancer biological processes, as well as promoted the stromal microenvironment according to our previous GSVA enrichment analysis of the Hallmark gene set. Furthermore, we assessed the enrichment scores for angiogenesis, immune checkpoint, cell cycle regulation, Pan F TBRS, EMT, and DNA replication using the ssGSEA method. The ssGSEA results showed pathways that can promote immune escape such as EMT1, EMT2, and immune checkpoint pathways were remarkably activated in the SUMO-C2 group (Figure 4D). Moreover, the SUMO-C1 group showed a higher expression of co-stimulatory and MHC molecules and lower expression of co-inhibitory molecules (Figures 4E–G), such as CD274, PDCD1, and LAG3, which are now commonly used as immunotherapeutic checkpoints in lung adenocarcinoma. Finally, we assessed the immune score, stroma score, and tumor purity score using the Estimate package for lung adenocarcinoma expression data and found that SUMO-C1 had a better immune score as well as a higher stroma score (Figure 4H), while SUMO-C1 had a lower tumor purity score, which further explains its association with better prognosis. Altogether, the above analysis demonstrated that the SUMO-C2 SUMOylation pattern promoted the development of a silent anti-tumor immune microenvironment and an active tumorigenic stromal microenvironment.

To further explore the underlying biological processes modulated in different SUMOylation patterns, we identified 906 significant DEGs (|FoldChange|>1.5, FDR<0.05) as SUMOylation pattern-associated genes (SPAGs) using the R package limma. These DEGs were also significantly associated with lung adenocarcinoma prognosis, indicating the critical role of SUMOylation patterns in dictating prognosis. KEGG functional annotation showed significant enrichment of SUMOylation-related genes to biological pathways such as cell cycle-related processes (cell cycle), cell adhesion (ECM-receptor interaction), immune activation signaling pathways (complement and coagulation cascades, IL 17 signaling pathway) (Figure 5A), while GO functional annotations revealed enrichment of biological processes such as regulation of the cell cycle phase transition, cell growth, humoral immune response, cell-cell junction and response to hypoxia (Supplementary Figures S4A–C).

Next, we constructed a protein interaction network (PPI, network type = physical subnetwork, minimum required interaction score = 0.7) of SPAGs using the STRING database and calculated the relationship between nodes using the Cytohubba plugin. We then classified genes with neighboring nodes ≥5 of 62 genes as key SPAGs and mapped their interaction (Figure 5B; Supplementary Figure S4D). Of the key genes, CCNB1 had the most neighboring nodes, followed by PLK1, BUB1B, CDK1, and BRCA1. It has been recently reported that the key genes in the PPI network can represent the functional characteristics of all genes accurately. Differential expression analysis by the TCGA dataset showed that key genes were significantly differentially expressed in lung adenocarcinoma from normal tissue samples (Figure 5C). Meanwhile, key genes extensively indicated significant prognostic significance, and CBX7, ITGA8, ADRB2, and CAV1 were significant favorable factors for lung adenocarcinoma (Figure 5D).

Unsupervised consensus clustering was carried out based on the identified key SPAGs. Three genomic subtypes were obtained, namely, gene-T1, gene-T2, and gene-T3 by CDF curve (Figures 5E, F; Supplementary Figure S4E). PCA suggested significant differences among the three genomic subtypes identified based on key genes (Figure 5G). Kaplan–Meier survival analysis indicated that gene-T2 was associated with the worst, gene-T3 with the best, and gene-T1 with the intermediate prognostic outcome, indicating different clinical courses of the genomic subtypes (Figure 5H). Moreover, the heatmap revealed differential expression of the 62 key genes among the three genomic subtypes (Figure 5I). This is consistent with the results of the prognostic utility of key genes in lung adenocarcinoma, as the heatmap showed that CAV1, ADRB2, CBX7, and ITGAB8 which were favorable factors for lung adenocarcinoma were remarkably upregulated in gene-T3, and the remaining risk factors for lung adenocarcinoma were significantly upregulated in the gene-T2 subtype. We then combined prognostic features such as age, sex, pathological stage, and genomic subtype with survival information for multicox analysis, which revealed that all the above prognostic features including genomic subtype are independent prognostic factors in lung adenocarcinoma; gene-T2 vs. gene-T1 (HR = 1.53, p = 0.003) and gene-T3 vs. gene- T1 (HR = 0.66, p = 0.008) (Figure 5J). Finally, we determined that SUMOylation regulatory genes were significantly differentially expressed among the three genomic subtypes and their expression in the gene-T1 subtype was intermediate to that in the gene-T2 and gene-T3 subtypes (Supplementary Figure S4F), indicating that different genomic subtypes further reflect the differences in the SUMOylation patterns.

We aimed to determine the prognostic value and biological significance of SUMOylation modification patterns in lung adenocarcinoma and to develop a SUMOylation pattern-based prognostic model for the assessment of individual lung adenocarcinoma patients. For this, we first performed a univariate cox analysis of the 62 key SPAGs along with LASSO regression to screen the best variables. Eventually, eight SUMOylation key genes that were significantly differentially expressed between normal tissues and tissues in lung adenocarcinoma and significantly associated with prognosis in lung adenocarcinoma were screened, and a SUMOylation-based prognostic signature model was constructed based on the expression of these genes and their coefficients (Supplementary Figures S5A, B). The risk score was calculated formulas follows: risk score = ExpNDC80*(−0.1024) + ExpORC1*0.0027+ ExpCCNE1*0.0164+ ExpCENPN*0.0127+ ExpCENPL*(−0.0997) + ExpCBX7*(−0.0606) + ExpCCNB2*0.4661+ ExpADRB2*(−0.0568) (Figure 6A). All lung adenocarcinoma samples were divided into high-risk and low-risk groups using the Survival package to calculate the optimal cut-off values. Survival analysis revealed that patients with higher scores were significantly associated with worse prognosis (log-rank p = 2.22e-16) (Figure 6B); the 1-year ROC AUC was 0.67, 3-year AUC was 0.68, and 5-year AUC was 0.69, indicating that the risk score based on the SUMOylation model had relatively superior prognostic accuracy compared to the vast array of existing prognostic models (Figure 6C). Multivariate Cox regression analysis incorporating age, gender, pathological stage, smoking history, and risk score confirmed that the risk score was an independent prognostic factor for assessing patient outcomes. (low vs. high, HR = 0.42 (0.33, 0.54), p < 0.001) (Figure 6D). ADRB2 and CBX7 expression were significantly negatively correlated with the risk score (Figure 6E). More deaths were observed in the high-risk group.

Four datasets, GSE31210 (Figure 6F), GSE37745 (Figure 6G), GSE50081 (Figure 6H), and GSE72094 (Figure 6I), and the external dataset TCGA-LUAD (Figure 6J) were used to further validate the stability and prognostic value of SUMOylation scores. The results of survival analysis gave p-values <0.001 for log-rank tests in all data sets, suggesting that the SUMOylation model risk score can consistently forecast patient prognostic outcomes; the ROC analysis also showed a high diagnostic value of the risk score in predicting lung adenocarcinoma prognosis (Supplementary Figures S5C–G). Multivariate cox regression analysis of TCGA lung adenocarcinoma cohort incorporating various prognostic factors such as age further confirms that the risk score is an independent risk factor for lung adenocarcinoma (low vs. high, HR = 0.52 (0.38, 0.71), p < 0.001) (Supplementary Figure S5H). To better understand the direct relationship between SUMOylation patterns and genomic subtypes, we constructed a Sankey diagram and observed that the vast majority of SUMO-C2 and a small fraction of SUMO-C1 comprise genomic subtype-T1 and then almost all of them comprise the Low Riskscore group (Figure 6K).

To verify the expression levels of model genes, we collected tumor tissues from 18 pairs of lung adenocarcinoma patients and the corresponding normal tissues adjacent to the cancer and performed quantitative real-time PCR. We also analyzed the expression levels and prognostic value of the model genes using the TCGA lung adenocarcinoma cohort. NDC80, CENPL, ORC1, CENPN, CCNE1, and CCNB2 were significantly upregulated in the TCGA cohort as well as in the 18 pairs of lung adenocarcinoma patient samples, and high expression was significantly associated with poor prognosis in lung adenocarcinoma. However, CBX7 and ADRB2 were significantly downregulated in lung adenocarcinoma tumor tissues and may serve as favorable prognostic factors for lung adenocarcinoma (Supplementary Figures S6A–H, p < 0.05).

SUMOylation score indicates immune microenvironment status and immunotherapy response in lung adenocarcinoma.

Comparison by the Wilcoxon rank-sum test showed that SUMO-C2 was associated with a higher risk score (Figure 7A), gene-T1 with an intermediate, gene-T2 with the highest, and gene-T3 with the lowest risk score (Supplementary Figure S7A), implying that differences in high and low-risk scores could reflect prognostic differences attributable to different SUMOylation patterns as well as different genomic subtypes. To further determine the relationship between SUMOylation scores and immune infiltration and other biological processes, the differences in the tumorigenic microenvironment between the groups with high- and low-risk scores were investigated. Immune checkpoint, cell cycle regulators, EMT1, EMT2, the cell cycle, and DNA replication were significantly active in the high-risk group (Figure 7B). Assessment of the level of immune infiltration in the Estimate package revealed significantly lower immune scores and higher tumor purity in the high-risk than in the low-risk group (Figure 7C). Moreover, Spearman correlation analysis exhibited a significant negative correlation between SUMOylation score and immune score, and a significant positive correlation between SUMOylation score and tumor purity (Supplementary Figure S7B), reinforcing the accuracy of the assessment of active biological processes in different SUMOylation patterns through risk scores.

Subsequently, Spearman analysis was then used to investigate the correlation between risk score and different immune cell subpopulations, and risk score was overall significantly negatively correlated with the levels of various immune cells. The risk score was significantly positively correlated with Th2 cells (r = 0.68) and Treg cells (r = 0.2)., and was negatively correlated with immune cells such as CD8+T cells (r = −0.4), Tem (r = −0.33), Tcm (r = −0.39), TFH (r = −0.31), CD56 (bright) natural killer cells (r = −0.36), DCs (r = −0.31), (immature Dendritic Cells)iDC (r = −0.4) and mast cells (r = −0.6). In addition, there was a close correlation between the degree of infiltration of Th1 cells, CD8⁺ T cells, T cells, and cytotoxic cells, as well as that of macrophages and iDCs in lung adenocarcinoma. (Figure 7D). These findings suggest that SUMOylation regulatory genes are involved in immune cell infiltration and immune regulation of lung adenocarcinoma tumors. Of these, Th2, Treg, and CD8⁺ T cells, and mast cells are involved in immune dynamic regulation (Supplementary Figure S7C). We then assessed immune checkpoint expression between the high- and low-risk groups to further investigate the effect of SUMOylation on immunotherapy response and found that the expression of clinically common immune checkpoint co-inhibitory molecules such as CD274, PDCD1, CTLA4, PDCD1LG2, HAVCR2, LAG3, and TIGIT was increased in the high-risk group (Figure 7E), while co-stimulatory molecules, such as CD2, CD28, and CD40LG, were decreased in the high-risk group (Supplementary Figure S7D), indicating that SUMOylation scores broadly modulate immune checkpoint expression. In addition, SUMOylation also suppresses the activation of the immune microenvironment by reducing the expression of MHC molecules (Supplementary Figure S7E).

A subsequent analysis of the role of SUMOylation patterns on immunotherapy response using the immunotherapy cohort reported in the lung cancer dataset GSE126044 indicated that the responder group (CR/PR) had low-risk scores and the non-responder group had high-risk scores (Figure 7F), explaining the worse immunotherapy response in high-risk patients. Meanwhile, we evaluated the immunotherapy response of lung adenocarcinoma patients in both high and low SUMOylation score groups in the GEO cohort (High-risk vs. Low-risk, OR is 1.716, p-value < 0.001, Figure 7G) and TCGA cohort (High-risk vs. Low-risk, OR is 2.207, p-value < 0.001, Figure 7H) based on the TIDE algorithm and indicated that low-risk patients were more likely to benefit from immunotherapy. In addition to this, significantly high SUMOylation scores are observed in patients with recurrent lung adenocarcinoma (Figure 7I), and significantly high SUMOylation scores are observed in smokers (Figure 7J), suggesting that SUMOylation scores may be effective in comprehensively assessing clinical prognostic factors in patients with lung adenocarcinoma.

As the above results have confirmed that SUMOylation leads to remarkable changes in the tumor microenvironment in lung adenocarcinoma, therefore, we used the SUMOylation score to characterize the sensitivity of lung adenocarcinoma patients to various chemotherapeutic agents to guide the combination chemotherapy. Based on the mRNA levels in lung adenocarcinoma tissues, we calculated the IC50 of various chemotherapeutic agents for lung adenocarcinoma patients in different SUMOylation score groups using the pRRophetic package and found that the group with low SUMOylation score was more sensitive to AKT-inhibitor-VIII, bexarotene, erlotinib, GDC0941, MK-2206, among others than that with high SUMOylation score (Figure 8A). In contrast, the high SUMOylation score group was sensitive to more chemotherapeutic agents, including A-443654, camptothecin, CGP-60474, cisplatin, cytarabine, docetaxel, doxorubicin, elesclomol, paclitaxel, and RO-3306 (Figure 8B). These findings suggest that SUMOylation scores can accurately reflect significant biological differences between SUMOylation patterns and also correlate with clinical prognostic features.

Meanwhile, as we previously illustrated by multivariate Cox analysis that age, gender, pathological stage, smoking history, and SUMOylation score all can serve as independent prognostic factors for lung adenocarcinoma, we combined the above significant prognostic factors to construct column plots to more accurately and efficiently predict the overall survival (OS) of individual patients with lung adenocarcinoma (Figure 8C). Calibration curves showed that the OS at 1, 3, and 5 years predicted by the nomogram closely corresponded to the actual OS of lung adenocarcinoma patients (Figure 8D) and the time-dependent ROC showed high AUC of 0.78, 0.83, and 0.82 for 1-, 3-, and 5- years, respectively. Yang et al. reported that 0.9≥ AUC>0.8 indicates excellent discrimination (Shengping and Gilbert, 2017), indicating that the developed nomogram has a remarkably prognostic performance (Figure 8E). Furthermore, a comparison of the predictive performance of multiple clinical prognostic features revealed that the nomogram offered the best predictive utility relative to any other clinical feature as determined by the ROC curve (Figure 8F) as well as the DCA decision curve (Figure 8G). Finally, to further validate the accuracy of the multifactorial regression model Nomogram, we evaluated the AUC of the ROC curves of Nomogram using three independent lung adenocarcinoma cohorts and the 5-year AUC values were 0.79, 0.81, and 0.89, respectively (Supplementary Figures S7F–H7). Validation in three independent lung adenocarcinoma cohorts supports that the accuracy of our Nomogram is significantly better than a large number of reported prognostic models for lung adenocarcinoma (Mo et al., 2020; Song et al., 2020; Chen et al., 2021).