Prognostic Value and Oncogenic Role of STING-Related Genes GAB3 and IL16 in Lung Adenocarcinoma: Implications for Immune Evasion and Treatment

doi:10.21203/rs.3.rs-5184237/v1

Download PDF

Article

Prognostic Value and Oncogenic Role of STING-Related Genes GAB3 and IL16 in Lung Adenocarcinoma: Implications for Immune Evasion and Treatment

https://doi.org/10.21203/rs.3.rs-5184237/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Lung adenocarcinoma (LUAD) is the most prevalent subtype of lung cancer (LC), and the stimulator of interferon genes (STING) is critical in inhibiting its progression. This study investigates the prognostic significance and molecular mechanisms of STING-related genes (STING-RGs) in LUAD. Differential expression analysis, weighted gene co-expression network analysis, as well as Cox regression (CR) identified GAB3 and IL16 as key prognostic genes. A LASSO-based risk model categorized LUAD patients into high-risk group (HRG) and low-risk group (LRG). HRGs exhibited lower GAB3 and IL16 expression and worse survival outcomes. A nomogram integrating risk scores (RS) and clinical factors effectively predicted patient survival. Functional enrichment, immune landscape, and mutation analyses revealed that HRGs were more likely to immune evasion, while LRGs responded better to targeted therapies. Mutation analysis showed lower survival in patients with high-risk scores (HRS) as well as high tumor mutational burden. Immunohistochemical staining confirmed that GAB3 was upregulated in LUAD tissues. In vitro experiments demonstrated that GAB3 overexpression promoted cancer cell proliferation and migration, while siRNA-mediated knockdown of GAB3 inhibited these processes, suggesting its role as an oncogene. In conclusion, GAB3 and IL16 are key prognostic markers, providing insights into STING-related immunotherapy strategies for LUAD.

Biological sciences/Cancer/Lung cancer/Non small cell lung cancer

Health sciences/Oncology/Cancer/Lung cancer/Non small cell lung cancer

Biological sciences/Immunology/Immune evasion

Lung adenocarcinoma

Stimulator of interferon genes

Prognostic

Immune Evasion

GAB3

Lung adenocarcinoma (LUAD) accounts for approximately 40% of all lung cancer (LC) cases, the most common cause of cancer-related death worldwide^1,2. The tumor microenvironme (TME) in LUAD is highly infiltrative and destructive, allowing cancer cells to invade blood vessels and lymph, leading to metastasis via the circulatory as well as lymphatic systems. The introduction of immunotherapies targeting immune checkpoints has significantly improved therapeutic outcomes and revolutionized the treatment of LUAD^3,4. Despite notable advancements in treatment, the five-year overall survival (OS) for LUAD is still below 20% ^5,6. Additionally, tumor heterogeneity, low response rates, and therapeutic resistance significantly hinder the efficacy of existing targeted therapies^7,8. Therefore, to enhance the therapy of LUAD, it is imperative to found certain molecular markers and establish precise therapeutic techniques.

The stimulator of interferon genes (STING) serves as a critical adaptor molecule triggering innate immune response following the recognition of intracellular pathogenic DNA, playing a critical role in controlling antiviral responses and detecting tumor formation^9–11. Research has identified two primary mechanisms through which STING influences tumor development. On the one hand, STING activation promotes type I interferon (IFN) synthesis and induces the maturation of dendritic cells (DCs) and macrophages, thereby enhancing the adaptive immune response and inhibiting tumor growth ^11–13. On the other hand, STING may promote tumor progression in the context of chronic inflammation by contributing to tumorigenesis or metastasis, a process linked to endoplasmic reticulum (ER) stress and T-cell survival, particularly CD8 + T cells. This mechanism is critical and independent of IFN signaling^14–15. Ze Lin et al. reported that low STING expression in LUAD tissues promotes cancer progression through the cyclic GMP-AMP synthase (cGAS)-STING pathway^16–17. Jianjun Wu et al. demonstrated that STING mutations result in lung inflammation and immune deficiency, characterized by reduced T-cell levels in mice, suggesting a link between this mechanism and ER stress¹⁴. However, few comprehensive studies exist, and no definitive conclusions have been drawn regarding the role of STING-related genes (STING-RGs) in LUAD progression.

Consequently, we aim to identify STING-related prognostic genes in LUAD. We explored the potential mechanisms by which STING regulates LUAD progression. Additionally, we constructed a new STING-related prognostic model and performed drug sensitivity analysis, offering a theoretical basis for clinical prognostication and LUAD treatment prediction.

2.1 Preparation of relevant datas

The Cancer Genome Atlas (TCGA) was queried for LUAD-related datasets, including gene expression and clinical information matrices for 526 LUAD samples and 59 control samples (training dataset). The data types included HTSeq-FPKM and HTSeq-Counts. For differential expression analysis, the HTSeq-Counts expression matrix was used, whereas the HTSeq-FPKM transformed TPM expression matrix (exp(log(fpkm) - log(sum(fpkm)) + log(1e6))) was employed for all other analyses. Gene symbols were employed to convert the ensemble IDs in the gene expression matrix, in conjunction with the corresponding annotation files for the respective species. Furthermore, the Gene Expression Omnibus provided the GSE13213 dataset, which is a validation dataset based on the GPL6480 platform and includes gene expression as well as clinical data matrices for 117 LUAD samples¹⁸. A query of the Molecular Signatures Database (MSigDB) with the term "STING" identified 47 STING-related genes (STING-RGs).

2.2 Analysis of the differential expression (DE) between the LUAD and control groups in the training dataset (TD)

Differential expression (DE) analysis between LUAD and control groups in training dataset was performed using DESeq2 (version 1.38.0) to identify differentially expressed genes (DEGs)¹⁹. The thresholds for differential expression analysis were set at |log2FC| > 0.5 and p < 0.05.

2.3 Acquisition of key modules and genes by gene co-expression network analysis

Dataset were found to differ between LUAD and control groups. The GSVA program (edition 1.42.0)'s single-sample gene set enrichment analysis (ssGSEA) was utilized to determine the 35 STING-RG scores for every LUAD sample in the TD²⁰. LUAD samples were then classified into HSG and LSG based on the median score. Kaplan-Meier (KM) analysis, using Wilcox test, was carried out on the two score groups to evaluate the impact of STING-RGs on LUAD patient survival.

To identify co-expressed gene modules as well as key modular genes associated with STING, co-expression networks were built via the weighted gene co-expression network analysis package, with STING-RG scores as the trait²¹. First, outlier samples were removed, followed by hierarchical cluster analysis. The optimal soft-threshold power (β) was then determined via a scale-free topology fit index of 0.9. Dissimilarity coefficients between genes were derived by calculating neighborhood and similarity measures, resulting in a systematic clustering tree with a minimum module size of 200. The association between STING-RG scores and module performance was then examined via Pearson correlation analysis. Modules exhibiting the strongest correlation with STING-RG scores, as depicted in a correlation heatmap, were designated as key modules. To illustrate the relationship between genes, modules, and traits, scatter plots of Gene Significance and Module Membership were generated. Genes with |MM| > 0.8 as well as |GS| > 0.2 were chosen as critical module genes.

2.4 Acquisition of differentially expressed STING-RGs (DE-STING-RGs) with prognostic significance

STING-RGs with prognostic significance were identified by performing the proportional hazards (PH) assumption test and univariate CR analysis on key module genes related with STING in disease samples from the training dataset (p < 0.05). Differentially expressed STING-RGs (DE-STING-RGs) with prognostic significance were then identified by intersecting DEGs and STING-RGs with prognostic significance. Additionally, the chromosomal distribution of DE-STING-RGs with prognostic significance was analyzed and visualized using Circos (http://circos.ca/)²².

2.5 Functional annotations for DE-STING-RGs with prognostic significance

To investigate the biological processes responsible for DE-STING-RGs with prognostic significance in the progression of LUAD, functional enrichment analysis of DE-STING-RGs with prognostic significance was performed using Metascape (https://www.metascape.org). A minimum overlap of 3, a minimum enrichment factor of 1.5, as well as a p-value threshold of less than 0.05 were the requirements for the screening process. Additionally, the Molecular Complex Detection (MCODE) algorithm was used to identify highly interconnected biological processes and annotate them in the network generated through functional enrichment analysis, using the PHYSICAL-CORE database. To build a PPI network, a range of 500 was selected for the maximum network size, with a minimum node size of 3.

2.6 The development, assessment, and verification of a risk model (RM)

The least absolute shrinkage and selection operator (LASSO) algorithm in the glmnet package (edition 4.1-4) was initially employed to screen prognostic genes and develop a RM based on DE-STING-RGs with prognostic significance²³. The formula was, in which Xi represented the relative expression of prognostic gene i, as well as coefi was the LASSO Cox coefficient for prognostic gene i. Subsequently, LUAD individuals in the TD were categorized into HRG and LRG based on the RS produced by the RM, via the optimal cut-off value. The validity of RM was assessed by conducting KM analysis using the survminer package (edition 0.4.9) for the two risk groups, with log-rank tests used to determine survival rate differences²⁴. Additionally, receiver operating characteristic (ROC) curves for LUAD patients at 1, 3, and 5 years were generated via survival ROC package (edition 1.0.3.1)²⁵. The Area Under the Curve (AUC) values were computed to evaluate RM’s effectiveness, with values exceeding 0.6 indicating commendable diagnostic efficacy. A heatmap was used to illustrate the expression patterns of prognostic genes in the two risk groups. Finally, the RM was validated using the validation dataset.

2.7 The associations with clinical features and RS

The clinical features of LUAD, including age (> 60 or ≤ 60 years), gender (female or male), tumor (T) stage (T1, T2, T3, T4), node (N) stage, as well as metastasis (M) (M0, M1, M2), were analyzed using the Wilcoxon test to evaluate their association with the RS¹⁸. The survival ROC function was used to generate ROC curves at 1, 3, as well as 5 years to evaluate the efficacy of RS and clinical features in predicting LUAD prognosis.

2.8 Evaluation of independent prognostic factors (IPFs) and nomogram construction

The clinical features and RS were incorporated into independent prognostic analyses¹⁸. Univariate CR analysis, the PH assumption test, as well as multivariate CR analysis were sequentially employed to find IPFs (p < 0.05). A nomogram for 1-, 3-, and 5-year survival in LUAD patients was then constructed on the basis of the independent prognostic factors. In addition, calibration curves were produced in order to evaluate the nomogram's validity for LUAD.

2.9 Enrichment analysis (EA) between two risk group

Differential expression (DE) analysis between the two risk groups in the TD was carried out using DESeq2, as well as log2FC values were produced. The log2FC values were ranked from largest to smallest, followed by gene set EA using clusterProfiler package (edition 4.7.1.3)^26–27. The reference gene set used was "cc2.cp.kegg.v7.4.symbols.gmt" from the MSigDB database. To further investigate variations in biological functions between the two risk groups, gene set variation analysis (GSVA) scores were generated for each pathway via the GSVA package, with the same reference gene set as in the gene set EA. Additionally, associations between the RS and SPs were explored.

2.10 Immune landscape analysis of LUAD

Immune cell infiltration levels in the two risk groups from TD were quantified using ssGSEA. The ssGSEA algorithm was then used to calculate enrichment scores (ESs) for 13 immunological functions in the two risk groups²⁸. Differences in immune function between the two groups were analyzed via Wilcoxon test (p < 0.05).

Immunotherapies targeting immune checkpoints eliminate tumor cells by modulating T-cell activity through co-suppressing or co-stimulating signaling pathways. An analysis of differences in immune checkpoints, using data from 48 published immune checkpoints, was applied to the two risk groups²⁸. To study the proportions of relevant immune cells, stromal cells, as well as tumor cells in the LUAD microenvironment, the StromalScore, ImmuneScore, and ESTIMATEScore of the biological samples were computed via ESTIMATE algorithm²⁸. Additionally, Spearman correlation analysis was used to evaluate associations between RS, StromalScore, ImmuneScore, and ESTIMATEScore.

2.11 Exploration of the the correlation with prognostic genes and the microenvironment of LUAD

To study the association between prognostic genes and single cells in a TME 16 single-cell RNA sequencing (scRNA-seq) datasets associated with LUAD were identified from the TME Single Cell Database (TISCH). The distribution and expression of prognostic genes in single cells were analyzed using the Seurat package (version 4.1.0)²⁹.

2.12 Analysis of immunotherapy responses between two risk group

The tumor immune response is a critical process in tumor development, and immune escape (IE) facilitates tumor proliferation, invasion, and metastasis. The TIDE method was utilized to determine the TIDE score, Dysfunction score,as well as Exclusion score for the two groups in order to evaluate the possibility of tumor IE in various risk categories of LUAD patients. An increased TIDE score indicated a higher probability of IE in LUAD patients during treatment. To further analyze the distinct impacts of immune checkpoints (IC), specifically PD-L1 as well as CTLA4, on cytotoxic T lymphocytes across different groups, Immunophenotype Score (IPS) data for IC inhibitors were obtained from the Cancer Immunity Database. Finally, ssGSEA scores for predicted immunotherapy pathways in the training dataset of LUAD samples were calculated using the GSVA package. Subsequently, Spearman analysis was employed to evaluate the relevance of RS with immunotherapy pathways as well as the immune cycle.

2.13 Analysis of mutated landscapes

To assess the relationship between tumor mutational burden (TMB) and survival, LUAD samples were grouped into high-TMB as well as low-TMB groups on the basis of their TMB levels. OS rates in the two risk groups as well as TMB groups was then compared using KM survival curves. The relationship between RS and TMB was assessed via Spearman correlation analysis, incorporating all samples in the training dataset. Somatic mutations in LUAD patients within the two risk groups were further visualized ans analyzed via maftools package to display tumor cell (TC) mutations³⁰.

2.14 Sensitivity analysis of LUAD to chemotherapeutic drugs

To assess the susceptibility of LUAD patients to chemotherapy drugs, drug sensitivity data were obtained from the GDSC database. The pRRophetic function was then used to compute the IC50 of each drug for the samples in the training dataset (version 0.5)³¹. Finally, the Wilcoxon test was employed to analyze variations in IC50 values between the two groups, identifying differential responses to chemotherapy drugs. The association between prognostic genes and drug efficacy was determined by Spearman correlation analysis, using the NCI-60 cell line as the reference.

2.15 Expression analysis of prognostic genes

To confirm prognostic gene expression, an analysis was conducted on the mRNA expression patterns of these genes in both tissues and cells. Specifically, the analysis focused on comparing gene expression in LUAD tissues with control tissues in the training dataset. Additionally, prognostic gene expression levels in 570 LUAD cell lines obtained from the CCLE online database were analyzed. Protein expression levels of prognostic genes in LUAD as well as control tissue samples were determined using immunohistochemistry images from the HPA online database.

2.16 Western blotting

Protein expression levels of protein expression levels of prognostic genes were confirmed using Western blot. Samples were treated with RIPA lysis buffer, as well as protein content was measured via a BCA test kit (Thermofisher, USA). 20 µg of protein was loaded into each lane. Proteins were transferred to a PVDF membrane (Millipore, USA) after separation by polyacrylamide gel electrophoresis. The membrane was blocked with 5% milk as well as incubated overnight at 4°C with primary antibodies diluted 1:1000: GAPDH (ab181602, Abcam, Cambridge, MA, USA), anti-IL16 (ab207181, Abcam, Cambridge, MA, USA), as well as anti-GAB3 (ab233411, Abcam, Cambridge, MA, USA). An improved chemiluminescence kit (Millipore, USA) was utilized to detect target proteins following a one-hour room temperature incubation with the matching horseradish peroxidase-linked secondary antibody.

2.17 Wound healing assay

GAB3 siRNA, overexpression A549 or NCL-H23 cells, as well as constructed empty vector (Ctrl) were all planted at a density of 1 × 10^5 cells/well on a 24-well culture plate. The cells were then grown for 24 hours at 37°C with 5% CO2. After removing the culture media, a 10 µl pipette tip was used to gently scrape the surface of the planted cells. The cells were gently washed twice with PBS before being given 1 milliliter of DMEM media containing 10% FBS. At 0 as well as 24 hours, pictures of the scratches were captured. There were three times of the experiment.The experiment was repeated three times.

2.18 CCK-8 assay

Constructed empty vector (Ctrl), GAB3 siRNA, as well as overexpression A549 or NCL-H23 cells were seeded in 96-well culture plates at a density of 1 × 10^4 cells/well as well as grown overnight in a 37°C, 5% CO2 incubator. Following a pair of PBS washes, 100 µl of DMEM media containing 10% FBS was introduced to the cells. In each well of the 96-well plate, 90 µl of medium containing 10 µl of CCK-8 solution was added after the medium was changed at 0, 24, and 48 hours. After that, the plates were incubated for two hours at 37°C. Using a microplate reader (BioTek Instruments, VT, USA), absorbance at 450 nm was determined for each well. Then, the ability of each group to proliferate was plotted.

2.19 Statistical analysis (SA)

For statistical analysis, R programming language (edition 4.2.3) as well as GraphPad 9.5 was used. The data was compared between different categories via the Wilcoxon test. Unless otherwise noted, p < 0.05 as well as p.adj < 0.05 was deemed statistically significant.

3.1 Totally 133 key module genes were identified

The training dataset (LUAD vs. control) identified 9,585 DEGs, encompassing 6,119 upregulated as well as 3,466 downregulated DEGs (Fig. 1A-B). The LUAD samples in the training dataset (TD) were categoried into high-score group (HSG) and low-score group (LSG) based on the STING-RG enrichment scores. Significant variations in survival between the two groups were found using KM survival analysis (p = 0.029). Consequently, STING-RG enrichment scores were used as a trait to construct co-expression networks (Fig. 1C). Subsequently, clustering analysis was carried out on the TD samples, as well as no outliers were excluded (Fig. 1D). To ensure biologically meaningful scale-free topology, an optimal β value of 5 was selected, considering a scale independence criterion of > 0.9 (Fig. 1E). 13 co-expression modules were found with a minimum module size of 200, among which the MEyellow module (Cor = 0.72, p < 0.05) exhibited the strongest positive correlation with STING-RG enrichment scores (Fig. 1F-G). A total of 133 key module genes were extracted from the MEyellow module on the basis of the standard |MM| > 0.8 and |GS| > 0.2 (Fig. 1H).

3.2 Totally 38 DE-STING-RGs with prognostic significance were screened from key module genes, some of which were involved in fatty acid metabolism

From the 133 key module genes, those that met the PH assumption test (p > 0.05) were encompassed in univariate CR analysis (p < 0.05), yielding 48 STING-RGs with prognostic significance (Fig. 2A). Subsequently, genes present in both the 9,585 DEGs and the 48 prognostically significant STING-RGs were identified, resulting in a total of 38 DE-STING-RGs with prognostic relevance (Fig. 2B). The majority of these genes were localized on human chromosome 7 and exhibited significant downregulation in LUAD samples (Fig. 2C). Functional enrichment analyses as well as protein-protein interaction (PPI) networks indicated ACVDAL, ACADM, ACADL, ACAA1, ACADSB, ACADS, and ACAT1 interact and cooperate in the fatty acid metabolism process (Fig. 2D-E).

3.3 GAB3 and IL16 were confirmed as prognostic genes, and RM based on prognostic genes had promising predictive capability for survival in LUAD

Prognostication is critically important for patients with LUAD. Two key genes, GAB3 and IL16, were identified from DE-STING-RGs with prognostic significance using LASSO regression analysis (λ = -3.507), and a RM was developed (Fig. 3A-B). To evaluate the reliability of the RM, the RS was computed using the equation: RS = (expression of GAB3 * -0.0599294226493178) + (expression of IL16 * -0.201038468078862). Based on the optimal RS cutoff (-1.335668), the LUAD samples were divided into LRG and HRG (Fig. 3C). As RS increased, a higher mortality rate was observed in LUAD samples, indicating a positive association between RS and mortality (Fig. 3D). Additionally, KM survival analysis showed that individuals in the low-risk LUAD group exhibited a remarkablly higher probability of OS (p < 0.05) (Fig. 3E). Moreover, AUC values for 1-, 3-, and 5-year survival exceeded 0.6, suggesting the RM effectively predicts survival outcomes in LUAD patients (Fig. 3F). Our analysis further revealed that high-risk LUAD samples were associated with low GAB3 and IL16 expression (Fig. 3G). Similar findings were replicated in the validation dataset, confirming the reliability of these results (Fig. 4A-E).

3.4 Significant association between clinical features and RS

The Wilcoxon test of risk score (RS) by sex, age, T stage (TS), N stage (NS), M stage (MS), and overall stage revealed significant differences for gender, stage (i-iii), TS (T1-T3), and NS (N2/N3-N0) (Fig. 5A-F). Moreover, AUC values of the ROC curves combining six clinical characteristics and RS were approximately 0.6 for 1, 2, and 3 years. Stage demonstrated superior predictive capability, with AUC values exceeding 0.6 for 1-, 2-, and 3-year survival (Fig. 5G-I).

3.5 A Nomogram integrating RS, TS, NS, and MS effectively predicted survival in patients with LUAD

RS, TS, NS, and MS were shown to be independent predictive factors (p < 0.05) (Fig. 6A-B). Therefore, a nomogram model was generated to forecast the likelihood of LUAD patient survival at 1, 3, as well as 5 years, incorporating RS, TS, NS, MS, and total points, with a higher total of independent prognostic factors corresponding to an increased likelihood of patient survival (Fig. 6C). The high precision of the nomogram in calculate LUAD patient survival was demonstrated through calibration curves (Fig. 6D).

3.6 Transcriptional and immune-related functional items were activated in the two risk groups, respectively

GSEA analysis revealed that the HRG showed ribosome activation, while the LRG showed activation of the intestinal immune network (IIN) for IgA production, graft-versus-host disease, autoimmune thyroid disease, as well as allograft rejection (Fig. 7A). Moreover, GSVA results indicated remarkable variations in 20 pathways between the two groups. In the HRG, aminoacyl-tRNA biosynthesis, RNA polymerase, base excision repair, and pyrimidine metabolism were activated. Conversely, the LRG demonstrated activation of the IIN for IgA production and hematopoietic cell lineage (Fig. 7B). In summary, transcriptional and immune-related pathways were activated in the respective risk groups. Notably, the two risk groups exhibited significant positive and negative correlations with IIN for IgA production (Fig. 7C).

3.7 Analysis of the immune landscape revealed a worse prognosis for high-risk LUAD

The ssGSEA algorithm identified 16 differentially infiltrating immune cells (ICs) between the two risk groups. All 16 cell types were heavily infiltrated in the LRG, including CD8 + T cells, activated DCs, as well as B cells (Fig. 8A). Additionally, 13 immune function enrichment scores were remarkably higher in LRG, revealing that high risk LUAD had diminished immune functions (Fig. 8B). In the training dataset, 39 immunological checkpoints were expressed, with 36 showing differential expression between the two groups. TMIGD2 expression was higher in the HRG, while the remaining 35 checkpoints were overexpressed in the LRG, suggesting that blocking TMIGD2 expression may benefit high-risk LUAD patients (Fig. 8C). Additionally, our findings revealed an inverse correlation between immunity score, stromal score, ESTIMATE score, and RS (Fig. 8D-G).

3.8 The correlation between prognostic genes and the microenvironment of the LUAD

Given the prognostic significance of the LUAD microenvironment in tumor development and progression, we analyzed the expression of GAB3 and IL16 in TME-associated cells using 16 datasets from the TISCH database. GAB3 expression was found to be high in CD4Tconv, CD8T, CD8Tex, as well as Tprolif cells in the NSCLC-GSE99254 and NSCLC-GSE151537 datasets. IL16 was highly expressed in CD4Tconv, CD8T, CD8Tex, Tprolif, as well as Treg cells in NSCLC-GSE99254 and NSCLC-GSE151537 datasets (Fig. 9A-B). We further analyzed the GSE99254 dataset, which was classified into six cell types. The number of each cell type, ranked in order, was CD4Tconv, CD8T, CD8Tex, Treg, Tprolif, and Mono/Macro (Fig. 9C). IL16 expression was higher in TME-associated cells compared to GAB3. These findings support the conclusion that the prognostic genes are closely related to LUAD (Fig. 9D).

3.9 Immunotherapy response that performed effectively in low-risk LUAD

The TIDE score, Exclusion score, as well as Dysfunction score was obtained via TIDE algorithm to assess immunotherapy performance in each LUAD patient across different groups. The results showed that the Dysfunction score was remarkably elevated in LRG, while the Exclusion score and TIDE score were substantially higher in HRG, suggesting high risk LUAD patients were more likely to immune escape. The response to immunotherapy in the two risk groups was 153:128 and 173:61, respectively, further demonstrating that immunotherapy was more effective for low-risk LUAD patients (Fig. 10A-B).

Additionally, analysis of immune checkpoint blockade (ICB) indicated remarkable variations between the two risk groups in their immunotherapeutic response to CTLA-4 as well as PD-1 in cytotoxic T lymphocytes. Scores for IPS, IPS-CTLA4 blocker, IPS-CTLA4 as well as PD1/PDL1/PDL2 blocker, along with IPS-PD1/PDL1/PDL2 blocker were all higher in LRG. This finding supports the efficacy of ICB targeting CTLA-4 and PD-1 in preventing cancer cells from evading the immune system in low-risk LUAD patients, thus contributing to their treatment (Fig. 10C).

Moreover, understanding TME heterogeneity could improve the precision of tumor immunotherapy. Hence, this study exmained the correlations between immunotherapy pathways,RS, and the immune cycle. Most signaling pathways and RS showed significant positive correlations, with the IFN-γ signature exhibiting the highest correlation with the RS (r = 0.330, p = 0.001). A strong positive association n was observed between the RS and three specific stages in the seven-step cancer immune cycle: cancer antigen presentation (step 2), induction and activation (step 3), as well as migration of ICs to the tumor (step 4). The correlation coefficients for these three stages were 0.450, 0.528, and 0.432, respectively (p = 0.001) (Fig. 10D-E).

Through the findings, we speculate low-risk LUAD patients have an activated immune response, making ICB more effective, while high-risk LUAD patients have a suppressed immune response, increasing the likelihood of cancer cell escape. However, as the RS increases, most immunotherapeutic pathways and immune cycles become activated, enabling sustained immunotherapy in high-risk patients (HRPs).

3.10 LUAD patients in the HRG as well as low-TMB group suffered from lower survival rates

Based on the combination of RS and TMB, LUAD patients were classified into four groups: High-TMB + High-risk, High-TMB + Low-risk, Low-TMB + High-risk, as well as Low-TMB + Low-risk. Patients with somatic hypermutations showed reduced survival rates (Fig. 11A). Correlation analysis further supported a positive association between RS and TMB (r = 0.19, p < 0.001) (Fig. 11B). A waterfall plot illustrated somatic mutations in patients across both risk groups, showing that somatic mutations were concentrated in TTN, TP53, CSMD3, MUC16, and PYR2. TTN showed a higher mutation frequency in HRPs, whereas TP53 showed a more consistent mutation pattern in low-risk LUAD patients. Overall, high-risk LUAD patients showed a higher likelihood of harboring mutations (Fig. 11C-D).

3.11 The high-risk LUAD patients displayed sensitivity to most chemotherapy treatments

Chemotherapy and targeted therapy are essential treatments for LUAD. The results suggested low-risk LUAD patients were more susceptible to several drugs, including Axitinib, Gefitinib, as well as Lenalidomide. Docetaxel, Paclitaxel, Vinorelbine, Erlotinib, and Lapatinib can be more effective in treating high-risk LUAD patients (Fig. 12). Correlation analyses showed that most anticancer drugs, such as Estramustine, Okadaic acid, and Dasatinib, were positively correlated with IL16, while Acetalax and GAB3 were negatively correlated (Supplementary Fig. 1).

3.12 Expression of prognostic genes

In LUAD samples from three surgically resected individuals, GAB3 protein expression was higher than in nearby normal tissues (Fig. 13A). Furthermore, GAB3 overexpression (GAB3-ov) appeared to promote the proliferation and migration of A549 as well as NCL-H23 cells. Conversely, GAB3 siRNA interference (siGAB3) inhibited the proliferation and migration of A549 and NCL-H23 cells (Fig. 13B-D). These results suggest that GAB3 may function as an oncogene.

STING affects cancer cell proliferation as well as invasion is linked to LUAD tumor progression through the cGAS-STING pathway, primarily by modulating the tumor immune microenvironment^11–14. However, there are limited comprehensive studies on investigating the role of STING-related genes (STING-RGs) in LUAD progression, and no definitive conclusions have been reached. In this study, we identified 38 DE-STING-RGs with prognostic significance in LUAD. These genes were associated with LUAD development and were all linked to poor prognostic outcomes. Functional enrichment analysis revealed that these 38 DE-STING-RGs were involved in the biosynthesis and metabolism of various amino acids, including tryptophan and branched-chain amino acids (like alanine, leucine, as well as isoleucine), as well as fatty acid metabolism, specifically oxidation and degradation. It is well established that there are distinct differences in carbohydrate, amino acid, as well as lipid metabolism between tumor cells (TCs) and normal cells (NCs)²⁶. Amino acid metabolism and its derivatives have been reported as key factors in the progression of LUAD². Zheng Wang et al. demonstrated that tryptophan metabolism is associated with tumor immunosuppression and immune escape, and that abnormal tryptophan metabolism can worsen the prognosis of LUAD patients³². Similarly, Kayo Ikeda et al. suggested that branched-chain amino acids promote tumor proliferation and stimulate tumor growth, which is associated with alterations in the tumor immune microenvironment^33–34. Furthermore, synthetic lipids are utilized by cancer cells as an energy source to drive angiogenesis, invasion, growth, and survival³⁵. According to Suyu Wang et al., fatty acid metabolism is essential for LUAD carcinogenesis, progression, and immune modulation^36–37. Our results agreed with those of these earlier investigations.

Based on the 38 DE-STING-RGs, we identified two prognostic genes in LUAD: GAB3 and IL16. Expression levels of both genes were significantly elevated in low- and high-risk samples compared to normal and LUAD samples, respectively. GAB3, or Grb2-associated binding protein 3, is an adaptor protein belonging to the Grb2-associated binders’ family of multi-substrate docking proteins. These proteins resemble insulin receptor substrate 1 and can bind phosphatidylinositol lipids in biological membranes^38–39. Grb2-associated binders facilitate immune cell activation by forming signalosomes and interacting with other signaling partners, acting as both positive and negative regulators⁴⁰. Specifically, GAB3 is involved in several growth factor and cytokine signaling pathways (SPs)⁴¹. It may also promote tumor growth in cancers such as ovarian cancer and gliomas^42–43. Lei Zhang et al. demonstrated that the LINC00324/miR-9-5p (miR-33b-5p)/GAB3 regulatory axis is critical in controlling tumor-associated macrophage risk as well as prognosis in LUAD patients^44–45.

In the present study, By analysing samples from three surgically resected LUAD patients, GAB3 protein expression was found to be higher in tumour tissues than in adjacent normal tissues. Further in vitro experiments showed that GAB3 overexpression significantly promoted the proliferation and migration of lung cancer cells, while GAB3 siRNA interference inhibited the proliferation and migration of these cells. These results suggest that GAB3 may function as a pro-oncogene in LUAD. The in vitro functional validation not only further supports our conclusions drawn from the raw letter analysis, but also highlights the significance of GAB3 as a potential therapeutic target.

Interleukin-16 (IL16) is a pro-inflammatory cytokine and chemoattractant for various immune cells^46–47. Lifeng Li et al. identified IL16 as a prognostic marker for LUAD and reported that IL16 expression levels in low-risk samples were remarkably higher than in high-risk samples, which agrees with our findings^48–49. IL16 also modulates T cell function and induces T cell immune responses, while IL16 deficiency can promote dendritic cell (DC) maturation⁵⁰. Interestingly, Jianjun Wu et al. showed that STING mutations regulate T cell function by inducing ER stress, thereby alleviating lung disease¹⁴. Notably, STING activation promotes the synthesis of IFNs, which induces dendritic cell maturation, enhances the adaptive immune response, and suppresses tumor growth¹². Therefore, we speculate that GAB3 and IL16 may influence LUAD progression by regulating T cell activation and immune response. Specifically, GAB3 was linked to macrophage activity, while IL16 was associated with dendritic cell maturation.

Subsequently, a RS model according to the two prognostic genes was constructed, as well as we found that LUAD risk was associated with multiple immune-related pathways, like B cell receptor, T cell receptor, chemokine, as well as JAK-STAT signaling pathways (SPs). The results agreed with Sijin Sun et al^51–52. The JAK-STAT pathway, which is controlled by cytokines, is critical in the innate immunity initiation, inflammation suppression, as well as regulation of immune responses, and it is also important in LUAD^53–54. Anqi Lin et al. identified the B cell receptor signaling pathway (SP) as a possible biomarker for predicting the effectiveness of immune checkpoint inhibitor therapy in LUAD patients⁵⁵. Notably, B cells activated through the B cell receptor SP can lead to the activation of tumor-specific CD8 + T cells, elevating anti-tumor immunity⁵⁶. Similarly, chemokines and the T cell receptor signaling pathway regulate immune cells within tumors, thereby influencing the therapeutic response and prognosis in LUAD patients^57–59. Additionally, we found that activated DCs, B cells, as well as CD8 + T cells were linked to LUAD risk. As previously mentioned, GAB3 and IL16 regulate T cell activation and immune response, with IL16 being linked to DC maturation. Single-cell localization results also showed that GAB3 and IL16 were primarily distributed in CD8 + and CD4 + T cells. These findings further support our hypothesis. Notably, IFN responses differed significantly between the risk groups, further confirming that IL16 promotes IFN synthesis and subsequently induces DC maturation. ESTIMATE analysis suggested that tumor purity was higher in the HRG, implying a greater inclination for tumor metastasis, consistent with the recurrent nature and poor prognosis of LUAD⁶.

Additionally, we found that Docetaxel, Paclitaxel, Vinorelbine, Erlotinib, and Lapatinib were potentially more effective in treating high-risk LUAD patients, while Axitinib, Gefitinib, and Lenalidomide showed greater sensitivity in low-risk patients. Furthermore, correlation analyses revealed that anticancer drugs such as Estramustine, Okadaic acid, and Dasatinib were positively correlated with IL16, while Acetalax was correlated with GAB3. These drugs have been widely reported in the treatment response of LUAD. The upregulated expression of GAB3 protein in clinical tumor samples, as well as its functional validation through in vitro experiments, further elucidates the oncogenic role of GAB3 in LUAD and underscores its potential as a therapeutic target.

In summary, this study, through bioinformatics analysis and experimental validation, identified two STING-related prognostic genes, GAB3 and IL16, which may influence LUAD progression by regulating T cell activation and immune responses. Moreover, we found that LUAD risk is associated with multiple immune-related pathways, with extensive roles played by T cells and dendritic cells. However, the precise influence of these prognostic genes on the mechanisms underlying LUAD progression remains to be confirmed, necessitating further in vivo and in vitro studies to provide more robust scientific evidence and support clinical interventions.

Ethics approval and consent to participate

This study was conducted in accordance with the ethical principles of the Declaration of Helsinki and approved by the Ethics Committee of the First Affiliated Hospital of Chengdu Medical College.Informed consent was obtained from patients for all tissue acquisitions.

Funding

This study was supported Scientific Research Funding Project of Ningxia Medical University(XT2024032), Natural Science Foundation of Ningxia Hui Autonomous Region in 2024(2024A1121), the Program of Sichuan Province Cadres Health Care (2024–2301) and Chengdu Medical College (CYZYB23-23).

Author Contribution

T J, XY Y, TT Z were responsible for performance of data analysis, and manuscript preparation; YJ C, YX X were responsible for performance of experiment; J C, K Y were responsible for manuscript writing and revision, and experimental design. Among them, J C and K Y serves as the corresponding author of the manuscript. All authors reviewed the manuscript.

Data Availability

The datasets analysed in this study are available in the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), including GSE13213. The Cancer Genome Atlas database, including TCGA-LUAD. Molecular Signatures Database (MSigDB), including STING-RGs.

Torre, L. A. et al. Global cancer statistics, 2012. CA Cancer J Clin 65, 87–108(2015).
Chang, W. et al. Identification of an amino acid metabolism-related gene signature for predicting prognosis in lung adenocarcinoma. Genes (Basel) 13,2295(2022).
Xu, J. Y. et al. Integrative proteomic characterization of human lung adenocarcinoma. Cell 182, 245–261.e217(2020).
Song, P. et al. Identification and validation of a novel signature based on NK cell marker genes to predict prognosis and immunotherapy response in lung adenocarcinoma by integrated analysis of single-cell and bulk RNA-sequencing. Front Immunol 13, 850745(2022).
Imielinski, M. et al. Mapping the hallmarks of lung adenocarcinoma with massively parallel sequencing. Cell 150, 1107–1120(2012).
Chang, J. T., Lee, Y. M. & Huang, R. S. The impact of the cancer genome atlas on lung cancer. Transl Res 166, 568–585(2015).
Consonni, D. et al. Lung cancer prognosis before and after recurrence in a population-based setting. J Natl Cancer Inst 107, djv059(2015).
Filosso, P. L. et al. Adenosquamous lung carcinomas: a histologic subtype with poor prognosis. Lung Cancer 74, 25–29(2011).
Ishikawa, H. & Barber, G. N. STING is an endoplasmic reticulum adaptor that facilitates innate immune signalling. Nature 455, 674–678(2008).
Ishikawa, H., Ma, Z. & Barber, G. N. STING regulates intracellular DNA-mediated, type I interferon-dependent innate immunity. Nature 461, 788–792(2009).
Li, T. & Chen, Z. J. The cGAS-cGAMP-STING pathway connects DNA damage to inflammation, senescence, and cancer. J Exp Med 215, 1287–1299(2018).
Ranoa, D. R. E. et al. STING promotes homeostasis via regulation of cell proliferation and chromosomal stability. Cancer Res 79, 1465–1479(2019).
Barber, G. N. STING: infection, inflammation and cancer. Nat Rev Immunol 15, 760–770(2015).
Wu, J. et al. STING-mediated disruption of calcium homeostasis chronically activates ER stress and primes T cell death. J Exp Med 216, 867–883(2019).
Merrild, D. M. et al. Erratum: Pit- and trench-forming osteoclasts: a distinction that matters. Bone Res 4, 16006(2016).
Lin, Z., Liu, Y., Lin, P., Li, J. & Gan, J. Clinical significance of STING expression and methylation in lung adenocarcinoma based on bioinformatics analysis. Sci Rep 12, 13951(2022).
Dong, Q. et al. Mutant RB1 enhances therapeutic efficacy of PARPis in lung adenocarcinoma by triggering the cGAS/STING pathway. JCI Insight 8, e165268(2023).
Yu, T. et al. CircRNAs in cancer metabolism: a review. J Hematol Oncol 12, 90(2019).
Escobar, L. F., Weaver, D. D., Bixler, D., Hodes, M. E. & Mitchell, M. Urorectal septum malformation sequence. Report of six cases and embryological analysis. Am J Dis Child 141, 1021–1024(1987).
Peng, H., Wang, Y. & Luo, W. Multifaceted role of branched-chain amino acid metabolism in cancer. Oncogene 39, 6747–6756(2020).
Ikeda, K. et al. Slc3a2 mediates branched-chain amino-acid-dependent maintenance of regulatory T cells. Cell Rep 21, 1824–1838(2017).
Amiri, M. et al. Diverse roles of fatty acid binding proteins (FABPs) in development and pathogenesis of cancers. Gene 676, 171–183(2018).
Wang, S. et al. Characterization of fatty acid metabolism in lung adenocarcinoma. Front Genet 13, 905508(2022).
Huang, D., Tang, E., Zhang, T. & Xu, G. Characteristics of fatty acid metabolism in lung adenocarcinoma to guide clinical treatment. Front Immunol 13, 916284(2022).
Gu, H. & Neel, B. G. The "Gab" in signal transduction. Trends Cell Biol 13, 122–130(2003).
Liu, Y. & Rohrschneider, L. R. The gift of Gab. FEBS Lett 515, 1–7(2002).
Wang, Z. et al. Gab2 and Gab3 redundantly suppress colitis by modulating macrophage and CD8(+) T-cell activation. Front Immunol 10, 486(2019).
Sármay, G., Angyal, A., Kertész, A., Maus, M. & Medgyesi, D. The multiple function of Grb2 associated binder (Gab) adaptor/scaffolding protein in immune cell signaling. Immunol Lett 104, 76–82(2006).
Berkel, C. & Cacan, E. GAB2 and GAB3 are expressed in a tumor stage-, grade- and histotype-dependent manner and are associated with shorter progression-free survival in ovarian cancer. J Cell Commun Signal 15, 57–70(2021).
Jia, P., Li, F., Gu, W., Zhang, W. & Cai, Y. Gab3 overexpression in human glioma mediates Akt activation and tumor cell proliferation. PLoS One 12, e0173473(2017).
Zhang, L. et al. Identification of a ceRNA network in lung adenocarcinoma based on integration analysis of tumor-associated macrophage signature genes. Front Cell Dev Biol 9, 629941(2021).
Wang, N., Zhou, X., Wang, X. & Zhu, X. Identification of Grb2-associated binding protein 3 expression to predict clinical outcomes and immunotherapeutic responses in lung adenocarcinoma. J Biochem Mol Toxicol 36, e23166(2022).
Glass, W. G., Sarisky, R. T. & Vecchio, A. M. Not-so-sweet sixteen: the role of IL-16 in infectious and immune-mediated inflammatory diseases. J Interferon Cytokine Res 26, 511–520(2006).
Cruikshank, W. W., Kornfeld, H. & Center, D. M. Interleukin-16. J Leukoc Biol 67, 757–766(2000).
Li, L. et al. Integrated analysis of dysregulated long non-coding RNAs/microRNAs/mRNAs in metastasis of lung adenocarcinoma. J Transl Med 16, 372(2018).
Feng, H. et al. Identification of signature genes and characterizations of tumor immune microenvironment and tumor purity in lung adenocarcinoma based on machine learning. Front Med (Lausanne) 9, 843749(2022).
Jia, R., Liu, S., Xu, J. & Liang, X. IL16 deficiency enhances Th1 and cytotoxic T lymphocyte response against influenza A virus infection. Biosci Trends 13, 516–522(2020).
Sun, S. et al. Development and validation of an immune-related prognostic signature in lung adenocarcinoma. Cancer Med 9, 5960–5975, doi:10.1002/cam4.3240 (2020).
Mendoza Diaz, M. I. et al. Dominant role of OH(-) and Ti(3+) defects on the electronic structure of TiO(2) thin films for water splitting. Dalton Trans 51, 15300–15311(2022).
Zhu, J., Li, Y. & Lv, X. IL4I1 enhances PD-L1 expression through JAK/STAT signaling pathway in lung adenocarcinoma. Immunogenetics 75, 17–25(2023).
Gao, Y. et al. PYCR1 knockdown inhibits the proliferation, migration, and invasion by affecting JAK/STAT signaling pathway in lung adenocarcinoma. Mol Carcinog 59, 503–511(2020).
Lin, A. et al. B cell receptor signaling pathway mutation as prognosis predictor of immune checkpoint inhibitors in lung adenocarcinoma by bioinformatic analysis. J Inflamm Res 15, 5541–5555(2022).
Ghosh, D., Jiang, W., Mukhopadhyay, D. & Mellins, E. D. New insights into B cells as antigen presenting cells. Curr Opin Immunol 70, 129–137(2021).
Tokunaga, R. et al. CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation - A target for novel cancer therapy. Cancer Treat Rev 63, 40–47(2018).
Ahrends, T. & Borst, J. The opposing roles of CD4(+) T cells in anti-tumour immunity. Immunology 154, 582–592(2018).
Jeong, S., Jang, N., Kim, M. & Choi, I. K. CD4(+) cytotoxic T cells: an emerging effector arm of anti-tumor immunity. BMB Rep 56, 140–144(2023).
Li, H. et al. Identification of lysosomal genes associated with prognosis in lung adenocarcinoma. Transl Lung Cancer Res 12, 1477–1495(2023).
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550(2014).
Hänzelmann, S., Castelo, R. & Guinney, J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14, 7(2013).
Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559(2008).
Li, Q. et al. Novel immunogenic cell death-related risk signature to predict prognosis and immune microenvironment in lung adenocarcinoma. J Cancer Res Clin Oncol 149, 307–323(2023).
Cao, R., Yuan, L., Ma, B., Wang, G. & Tian, Y. Immune-related long non-coding RNA signature identified prognosis and immunotherapeutic efficiency in bladder cancer (BLCA). Cancer Cell Int 20, 276(2020).
Liu, T. T. et al. Identification of CDK2-related immune forecast model and ceRNA in lung adenocarcinoma, a pan-cancer analysis. Front Cell Dev Biol 9, 682002(2021).
Zhang, S. et al. Development and validation of PET/CT-based nomogram for preoperative prediction of lymph node status in esophageal squamous cell carcinoma. Ann Surg Oncol 30, 7452–7460(2023).
Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2, 100141(2021).
Zhang, W. et al. Identification of cuproptosis and immune-related gene prognostic signature in lung adenocarcinoma. Front Immunol 14, 1179742(2023).
Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e3529(2021).
Wang, W. et al. Histone acetylation modification regulator-mediated tumor microenvironment infiltration characteristics and prognostic model of lung adenocarcinoma patients. J Thorac Dis 14, 3886–3902(2022).
Geeleher, P., Cox, N. & Huang, R. S. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9, e107468(2014).

No competing interests reported.

Download PDF

Reviewers invited by journal
11 Nov, 2024
Editor assigned by journal
10 Nov, 2024
Editor invited by journal
21 Oct, 2024
Submission checks completed at journal
18 Oct, 2024
First submitted to journal
30 Sep, 2024

You are reading this latest preprint version

Prognostic Value and Oncogenic Role of STING-Related Genes GAB3 and IL16 in Lung Adenocarcinoma: Implications for Immune Evasion and Treatment

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials and methods

2.1 Preparation of relevant datas

2.3 Acquisition of key modules and genes by gene co-expression network analysis

2.4 Acquisition of differentially expressed STING-RGs (DE-STING-RGs) with prognostic significance

2.5 Functional annotations for DE-STING-RGs with prognostic significance

2.6 The development, assessment, and verification of a risk model (RM)

2.7 The associations with clinical features and RS

2.8 Evaluation of independent prognostic factors (IPFs) and nomogram construction

2.9 Enrichment analysis (EA) between two risk group

2.10 Immune landscape analysis of LUAD

2.11 Exploration of the the correlation with prognostic genes and the microenvironment of LUAD

2.12 Analysis of immunotherapy responses between two risk group

2.13 Analysis of mutated landscapes

2.14 Sensitivity analysis of LUAD to chemotherapeutic drugs

2.15 Expression analysis of prognostic genes

2.16 Western blotting

2.17 Wound healing assay

2.18 CCK-8 assay

2.19 Statistical analysis (SA)

3. Results

3.1 Totally 133 key module genes were identified

3.4 Significant association between clinical features and RS

3.5 A Nomogram integrating RS, TS, NS, and MS effectively predicted survival in patients with LUAD

3.6 Transcriptional and immune-related functional items were activated in the two risk groups, respectively

3.7 Analysis of the immune landscape revealed a worse prognosis for high-risk LUAD

3.8 The correlation between prognostic genes and the microenvironment of the LUAD

3.9 Immunotherapy response that performed effectively in low-risk LUAD

3.10 LUAD patients in the HRG as well as low-TMB group suffered from lower survival rates

3.11 The high-risk LUAD patients displayed sensitivity to most chemotherapy treatments

3.12 Expression of prognostic genes

4. Discussion

Declarations

Ethics approval and consent to participate

Funding

Author Contribution

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Version 1