The AVEN mRNA expression in pan-cancer was analyzed by the GSCA web tool (http://bioinfo.life.hust.edu.cn/GSCA). GEPIA (Li et al., 2021) (http://gepia.cancer-pku.cn/) web tool was applied for the survival analysis in COAD, KIRC, KIRP, LUAD, LUSC and THCA. In order to conduct a more detailed investigation into the role of AVEN in LUAD, RNA-seq data and clinical data in LUAD patients derived from TCGA database were obtained from the UCSC Xena website (http://xena.ucsc.edu/). Log2(FPKM+1) value obtained from RNA-seq data was converted to TPM (Transcripts Per Million) value. Subsequently, patients with the highest 25% of AVEN expression were categorized into AVEN^high group, and patients with the lowest 25% of AVEN expression were classified into the AVEN^low groups. Overall survival analysis was performed in AVEN^high and AVEN^low groups.

To compare the characteristics between AVEN^high and AVEN^low patients, clinical data was downloaded from the UCSC Xena website, including age, TNM classification, gender, radiation therapy status, race, AVEN expression, and smoking status. The R package moonBook was exploited to visualize characteristics of patients between AVEN^high and AVEN^low groups.

The whole exome data of LUAD patients was downloaded from the cBioportal (Cerami et al., 2012; Gao et al., 2013) to examine genetic alterations (https://www.cbioportal.org/). Driver genes such as TP53, EGFR, KRAS, ERBB2, BRAF, ALK, RET, FGFR3, NTRK3 and ROS1 were selected, and their alteration pattern was examined in AVEN^high and AVEN^low LUAD patients. The R package ComplexHeatmap was used to generate an oncoprint plot.

Spearman’s rank correlation test, which works with rank-order variables instead of raw data value of the variables, was used to obtain differentially expressed genes (DEGs) between AVEN^high and AVEN^low group. Genes with an absolute R-value>0.4 were used for pathway analysis by ConsensusPathDB (Kamburov et al., 2009) (http://cpdb.molgen.mpg.de/MCPDB). Significantly altered pathways were selected with the criteria of p < 0.05 and were visualized using SRplot (https://www.bioinformatics.com.cn/srplot). Genes in Cell cycle and VEGFA-VEGFR2 pathways were further visualized by using the R package ComplexHeatmap. GSEA analysis was performed to compare the pathway enrichment between AVEN^high and AVEN^low groups by using these gene sets as references: “SHEDDEN_LUNG_CANCER_POOR_SURVIVAL_A6″, “HallMARK_MTORC1_SIGNALINF”, “HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION”, and “VEGF_A_UP.V1_DN”.

Immune score, stromal score, and estimate score representing immune infiltration, stromal cell level, and purity of tumor, respectively were obtained from Estimate website (https://bioinformatics.mdanderson.org/estimate/index.html) and compared between AVEN^high and AVEN^low patients. Furthermore, multiple deconvolution tools including Timer, CIBERSORT, EPIC, xCEll, Quanti-seq and MCP-counter (Li et al., 2017; Racle et al., 2017; Sturm et al., 2019) were utilized to examine various types of immune cells in tumor tissue. Immune infiltration data was obtained from Timer (https://cistrome.shinyapps.io/timer/) and was visualized by boxplot, using the R package via ggplot2.

To evaluate the prediction value of AVEN in immunotherapy response, Tumor Immune Dysfunction and Exclusion (TIDE) score (Jiang et al., 2018) was calculated (http://tide.dfci.harvard.edu/). Consequently, the expression levels of immune checkpoint genes and functional genes associated with cytotoxic T cells were analyzed in AVEN^high and AVEN^low LUAD. Furthermore, we extended our investigation to encompass additional immunotherapy response markers indicative of B cells, testing these markers in AVEN^high and AVEN^low LUAD patients.

To determine the correlation between each gene from the DEGs and the overall survival of LUAD patients, Univariate Cox Regression model was employed. AVEN-derived genes with p-value < 0.01 were regarded as AVEN-derived prognostic factors. Followed by Random Survival Forest analysis, the relative importance of each gene was calculated. Genes with relative importance >0.5 were used in the Multivariate Cox Regression model. Step forward Cox regression was utilized to optimize the model. The AVEN-derived genomic model was formulated as follows (Abd ElHafeez et al., 2021): R i s k   s c o r e = h 0 t × exp ( K R T 6 A × 0.0002919 + S L C 16 A 3 × 0.0045 + C T S L × 0.0008009 + L D H A × 0.007940 + C D C 42 E P 2 × 0.0147 )

The h0(t) represents the baseline hazard at time t, which denotes the hazard of an individual when all predictor variables are set to 0. Subsequently, based on the calculated risk scores, patients were categorized into high-risk and low-risk groups using the mean risk score. The survival and survminer packages were utilized to determine the survival state of LUAD patients in both the TCGA and GEO cohorts (GSE50081, GSE31210). ROC curves were generated to test the specificity and sensitivity of the AVEN-derived prognostic gene model by using the survivalROC package.

PrognoScan (Mizuno et al., 2009) was employed to assess the prognostic significance of AVEN across multiple LUAD data cohorts. Additionally, AVEN protein abundance was examined using the cProSite (Wang et al., 2023).

Statistical analysis was performed with R software (v4.2.1) and its suitable packages. In this study, group comparisons were performed using Student’s t-test, and the interaction between variables were examined using the Spearman correlation test.

The role of AVEN in cancer is systematically studied according to the workflow shown in Figure 1. First, the GSCA web tool was exploited (Liu et al., 2023) to identify the AVEN mRNA level in various cancer types. AVEN was highly expressed in six types of cancer compared to normal tissue: colon adenoma (COAD), Kidney renal clear cell carcinoma (KIRC), Kidney renal papillary cell carcinoma (KIRP), Lung adenocarcinoma (LUAD), Lung squamous cell carcinoma (LUSC), and Thyroid carcinoma (THCA) (Figure 2A). On the other hand, a slight downregulation of AVEN was observed in two types of cancer, Cholangiocarcinoma (CHOL), and Cervical squamous cell carcinoma (CESC). Besides mRNA, AVEN protein level was analyzed using the cProSite database, which showed a higher AVEN protein abundance in Breast cancer, Colon cancer, Kidney cancer, Liver cancer, Lung adenocarcinoma, Lung squamous cell carcinoma, and ovarian cancer compared to its adjacent normal tissue (Supplementary Figure S1). Given our interest in the oncogenic role of AVEN, we performed GEPIA analysis to determine the overall survival in the six tumor types with elevated AVEN mRNA expression: COAD, KIRC, KIRP, LUAD, LUSC and THCA (Supplementary Figure S2). Notably, AVEN showed a significant prognostic effect only in LUAD. To further validate this finding and investigate the oncogenic features of AVEN, we acquired the mRNA sequencing data and clinical information of LUAD patients from the TCGA dataset. Based on the expression level of AVEN, LUAD patients were classified into two groups: AVEN^high (top 25%) and AVEN^low (bottom 25%), and overall survival was investigated using Kaplan-Meier analysis (Figure 2B). Consistent with GEPIA results, high AVEN expression was markedly associated with poor survival in LUAD. Before studying details of underlying molecular and cellular mechanism of AVEN, patient characteristics of AVEN^high and AVEN^low gruop were analyzed (Table 1). As shown in Table 1, AVEN expression did not show significant differences in age, race, and smoking but only in gender. Interestingly, the value of the T stage and N stage showed significant differences between the two groups. TNM classification is a system that defines tumor size (T), regional lymph amount (N), and spread of cancer (M) in a patient’s body (Rosen and Sapra, 2022). Thus, we wondered whether AVEN contributes to tumor progression according to TNM classification. Notably, AVEN showed a gradual increase with tumor progression when AVEN expression was seen in normal tissue and different T stages (Figure 2C). AVEN expression also increased in the N1 stage compared to N0 (Figure 2D). We further analyzed AVEN expression in different M stages, however, there were no significant differences in the number of patients in different M stages between AVEN^high and AVEN^low patients (Table 1). Consistently, AVEN expression was similar between M0 and M1 stages (Figure 2E). Since genetic alterations can lead to the activation of various signaling pathways that promote the growth and survival of cancer cells (Skoulidis and Heymach, 2019), we next examined genetic alterations in driver genes that were well defined in LUAD. However, no significant associations between AVEN expression and genetic mutations in EGFR, KRAS, ALK, ROS1, BRAF, etc., Were observed (Figure 2F).

We next investigated molecular and cellular pathways associated with AVEN. First, DEGs between AVEN^high and AVEN^low groups were analyzed by using Spearman’s rank correlation test with the criteria of the absolute R-value over 0.4. A total of 838 genes were obtained (Supplementary Table S1), followed by pathway enrichment analysis via the ConsensusPathDB website tool. We found that those DEGs were mainly involved in the biological process pathways (such as metabolism of proteins, and membrane trafficking), oncogenic pathways (such as the cell cycle and VEGFA-VEGFR2 Signaling Pathway), and immune regulation process (such as neutrophil degranulation, and innate immune system) (Figure 3A).

Given the fact that VEGF plays an important role in tumor progression and angiogenesis (Jiang et al., 2020), we suggested that AVEN might enhance tumor aggressiveness by promoting the cell cycle and angiogenesis. To better demonstrate a correlation between AVEN and the gene expression profile in the cell cycle or VEGFA-VEGFR2 Signaling Pathway, heatmaps were introduced for visualization. Genes in these pathways were distinctly upregulated by AVEN overexpression (Figure 3B). Gene set enriched analysis (GSEA) was applied to further investigate the oncogenic role of AVEN. Genes in MTORC1 signaling, epithelial-mesenchymal-transition (EMT), and VEGF pathways were remarkably enriched in AVEN^high patients (Figure 3C). MTOR pathway is a key pathway to regulate cell growth. EMT and VEGFA also play a critical role in tumor progression (Brabletz et al., 2018; Zou et al., 2020). Furthermore, we found that AVEN overexpression was associated with a poor survival genomic signature (Figure 3D) in NSCLC (NES = 3.43, p-value = 0.0), consistent with results shown in Figure 2B.

The tumor microenvironment (TME) consists of various cells such as normal epithelial, vascular cells, stromal cells, and immune cells. Immune cells, such as T cells, B cells, natural killer cells, dendritic cells, macrophages, and others are attracted to the site of a tumor by various signals, including chemokines and cytokines produced by tumor cells and stromal cells and regulate tumor cells proliferation/apoptosis (Galli et al., 2020). Considering our finding of AVEN-associated signaling pathways such as innate immune system, neutrophil degranulation, and EMT (Figures 3A, C), it suggested unique TME features associated with AVEN expression.

Thus, we first analyzed immune cells and stromal cells using the ESTIMATE algorithm in AVEN^high and AVEN^low groups. Despite the fact that there was no significant difference in stromal score, AVEN^high patients showed lower immune scores compared to AVEN^low patients, which suggests the functional possibility of AVEN to inhibit the abundance of immune cells in the TME. Based on the stromal score and immune score, the Estimate score was calculated to infer the purity of tumor tissue (Yoshihara et al., 2013), and the result showed AVEN^high and AVEN^low groups shared a similar tumor purity (Figure 4A). To explore the diverse types of immune cells, multiple deconvolution tools were employed, including Timer, CIBERSORT, EPIC, xCEll, Quanti-seq, and MCP-counter (Sturm et al., 2019). These tools utilize different methodologies to analyze the RNA-seq data and provide insights into the composition and abundance of immune cell populations (Im and Kim, 2023). The analysis using various deconvolution tools revealed that B cells were significantly more abundant in AVEN^low groups. However, different patterns were observed for other immune cell types depending on the specific tool used. For instance, Timer and EPIC indicated that CD4 T cells were significantly higher in AVEN^low groups, while Quanti-Seq suggested higher levels in AVEN^high groups. Additionally, CIBERSORT analysis showed that activating memory CD4 T cells were more abundant in AVEN^high groups, whereas resting memory CD4 T cells were higher in AVEN^low groups (Figure 4B). EPIC, MCP-counter and xCEll subtracted cancer associated fibroblast (CAF) from tumor. Specifically, EPIC and MCP-counter indicated that CAF were significantly more abundant in AVEN^high groups. However, in contrast to these findings, xCEll tool showed lower levels of CAF in AVEN^high groups (Figure 4B).

Immune checkpoint inhibitors (ICIs) in cancer treatments have demonstrated significant potential in enhancing antitumor immune responses and improving patient outcomes (Bondhopadhyay et al., 2020; Robert, 2020). However, their responsiveness remains limited, necessitating the identification of predictive biomarkers (Bai et al., 2020; Błach et al., 2021). Therefore, we investigated whether AVEN could serve as a predictive biomarker for ICIs response using the TIDE algorithm, in addition to assessing the expression of immune checkpoint genes. To explore the AVEN-associated immunotherapy response, RNA-seq data from 260 patients in AVEN^high and AVEN^low groups were processed in TIDE tool. We found that among the 158 non-responders, 62% (n = 98) belonged to the AVEN^high group, while 38% (n = 60) were in the AVEN^low group. Specifically, among the AVEN^high patients (n = 130), only 24% (n = 32) were predicted to be responders based on the TIDE tool (Figure 5A, Supplementary Table S2). Subsequently, expression of immune checkpoint genes was examined. Immune checkpoint genes such as CD274, CTLA4, and LAG3 (Qin et al., 2019; Wang et al., 2019; Hu et al., 2021) had no significant differences between AVEN^high and AVEN^low groups (Figure 5B). Genes associated with anti-tumor CD8 T cells such as CD8A and GZMA (van der Leun et al., 2020) were also examined, but no differences were observed (Figure 5C). Immune checkpoints and CD8 T cell function are regarded as the essential indicators for immunotherapy respond, however the immunotherapy resistance is far more complicated than that. For example, a lower mutation burden (TMB), dysfunctional MHCs complex, or immunosuppressive tumor microenvironment (Lei et al., 2020) can result in a lower responsiveness in immunotherapy treatment without altering the expression level of immune checkpoint and CD8 T cell markers. Moreover, in our current study, the observation of diminished B cell infiltration (Figure 4B) in AVEN^high patients have spurred us to propose additional possibilities that B cell mediates immunotherapy resistance. With the knowledge that B cells also express the receptors of PD1/PDL1/CTLA4 (Kim et al., 2021) and have been established as a favorable prognostic marker in NSCLC (Germain et al., 2014), there exists a theoretical basis for B cells to respond to immune checkpoint inhibitors. We extended our analysis to assess the expression of B cell markers in LUAD. As depicted in Figure 5D, there was a significant reduction in mRNA levels of CD19, CD20, and CD22 in patients with elevated AVEN expression, further suggesting a distinct immune phenotype in AVEN^high patients, characterized by reduced B cell infiltration.

As AVEN displayed significant predictive value in terms of LUAD survival within the TCGA database, we were motivated to investigate its prognostic utility across other cohorts. Employing PrognoScan, we delved into the influence of AVEN expression on different datasets. While a substantial correlation with LUAD survival was observed in some cohorts, it is important to note that certain cohorts exhibited limitations in significant prognostic effect of AVEN (Supplementary Table S3). Since AVEN-derived DEGs showed remarkable effect in pivotal oncogenic pathways, we further investigated the prognosis effect of 838 DEGs by using Univariate Cox Regression (Ma et al., 2022; Yu et al., 2022). 305 of survival-related genes with p-value < 0.01 (Supplementary Table S4) were shown to contribute to LUAD patient survival and determined as prognostic factors, which were then used as input for the Random Survival Forest analysis (Figure 6A). Following the Random Survival Forest analysis, the relative importance of the above prognostic genes was calculated and ranked. Eight genes, KRT6A, SLC16A3, AHNAK2, CTSL, FAM83A, LDHA, CDC42EP2, and SPHK1, were selected with a relative importance value of over 0.5 (Figure 6B), and Multivariate-Cox regression model was developed with those genes. In order to optimize survival prediction model, step-forward Multivariate Cox analysis was used to further screen the valuable prognostic genes containing KRT6A, SLC16A3, CTSL, LDHA, and CDC42EP2. According to the model, the risk score of each patient was identified as follows: h 0 t × exp K R T 6 A × 0.0002919 + S L C 16 A 3 × 0.0045 + C T S L × 0.0008009 + L D H A × 0.007940 + C D C 42 E P 2 × 0.0147 . In accordance with the mean risk score, LUAD patients were classified into two subpopulations, high and low risk score. Kaplan-Meier survival analysis showed LUAD individuals with low-risk scores had a greater benefit in survival compared with a high-risk score population, which revealed that the AVEN-derived prognostic model was able to estimate patients’ prognosis according to risk score (Figure 6C). Additionally, the time-dependent ROC curves were performed and confirmed that AVEN derived prognostic model showed potency to predict patient prognosis (Figure 6D).

Two independent external datasets, GSE50081 and GSE31210, were used to evaluate AVEN-derived prognostic model. As shown in Figures 7A, B, LUAD patients with low-risk score had better survival, which indicated that the AVEN-derived prognostic model is effective and sufficient in predicting LUAD survival.