Expression patterns of E2Fs identify tumor microenvironment features in human gastric cancer

Data acquisition of gastric cancer samples

Training data sets were retrieved using the R package “TCGAbiolinks” from The Cancer Genome Atlas Stomach Adenocarcinoma (TCGA-STAD) database, including transcriptome data,Fragments Per Kilobase of exon model per Million mapped fragments (FPKM), and single nucleotide variations (SNVs). The official corrected survival information (overall survival, OS) and clinical information (age, stage, gender, grade) were downloaded from the TCGA Pan-Cancer Clinical Data Resource (TCGA-CDR). A total of 351 tumor samples with both expression and survival information were retained for subsequent analysis.

The validation data sets, comprising six independent cohorts (GSE13861, GSE15459, GSE26899, GSE26901, GSE66229, GSE84433) were downloaded from the Gene Expression Omnibus (GEO).

The validation cohort for immunotherapy was acquired using the R package “IMvigor210CoreBiologies,” which included transcriptomic data and clinical information of 298 urothelial carcinoma patients who underwent anti-PD-1/PD-L1 therapy.

Unsupervised clustering of transcription factors and differentially expressed genes

The R package “ConsensusClusterPlus,” along with the unsupervised clustering method, was used to divide patients into distinct E2F clusters based on E2F expression levels, which was then repeated 1,000 times to verify the stability of the classification. Gene cluster count was determined using the same methods asE2F.

Identification of the immune characteristics of clusters

Single-sample gene-set enrichment analysis (ssGSEA) was employed to assess the relative abundance of immune cell infiltration in TME, which comprised 28 human immune cell types, including activated CD8 T cells, activated dendritic cells, and macrophages (Charoentong et al., 2017).

The ESTIMATE algorithm was adopted to calculate the immune score and stromal score of each tumor sample. The Wilcoxon test was then used to compare the differences in the scores between groups.

Selection and identification of differentially expressed genes between distinct subgroups

The R package “limma” and empirical Bayesian algorithm were adopted to filter out differentially expressed genes (DEG) between diverse E2Fclusters, which were defined as E2F-related genes, using the following screening criteria: fold change value ≥ 1.5 (—log2FC—>=0.585) and a significance threshold (FDR<0.05). These genes were further collected through a univariate Cox regression analysis to obtain a refined prognosis-related gene list. Then, the samples were divided into various gene clusters by applying an unsupervised clustering method according to the expression levels of genes ultimately screened out.

Construction of E2F scoring system

A principal component analysis (PCA) was employed to construct an E2F scoring system according to E2F prognostic-related genes through a dimension reduction method to transform many indicators into a couple of comprehensive indicators. Both principal components 1 and 2 were selected to act as E2F scores (Sotiriou et al., 2006). The advantage of this approach is that it retains the critical genetic characteristics while removing non-essential genetic information. The computed formula was as follows: ${\displaystyle \text{E2F score}=\sum \text{(PC1i + PC2i)}}$ where i represents the E2F prognostic-related gene expression level.

Analysis of functional mechanisms

Gene set variation analysis (GSVA) is a non-parametric unsupervised strategy to estimate changes in the activity of pathways and biological processes in samples, and was used in this study to determine underlying gene sets with E2F score changes. The hallmark gene set (msigdb.v7.4.symbols.gmt) was downloaded from the Molecular Signatures Database (MSigDB) for GSVA. Statistical significance was considered at P ≤ 0.05.

A Gene Ontology (GO) enrichment analysis involving molecular functions (MFs), biological pathways (BPs), and cellular components (CCs) and a Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis centering gene function were used to annotate the DEGs of distinct E2F phenotypes with p-value cutoff = 0.05 and pAdjustMethod = “BH” using the R package “clusterProfiler.”

Tumor drug sensitivity and efficacy prediction

Based on the Genomics of Drug Sensitivity in Cancer (GDSC) dataset, the sensitivity of chemotherapy drugs was measured by calculating the IC50 in the training set using the R package “oncoPredict” and the “calcPhenotype” algorithm.

The immunotherapy validation cohort was obtained using the R package “IMvigor210CoreBiologies,” which included transcriptomic data and clinical information of 298 urothelial carcinoma patients who underwent anti-PD-1/PD-L1 therapy.

Quantitative real-time PCR

RNA was extracted from both normal and tumor tissues of 12 GC patients randomly selected from the tissue bank of the First Affiliated Hospital of Xi’an Jiaotong University. Written, informed consent was obtained from all patients. Ethical approval for the study was obtained from the Ethics Committee of The First Affiliated Hospital of Xi’an Jiaotong University (No. XJTU1AF2023LSK-451). The clinical information of the participating patients are listed in Table S1. TRIzol reagent was used to extract total RNA from normal tissues and tumor tissues, and spectrophotometry was used for quantification and quality control. Quantitative PCR analysis was performed with the SYBR-Green PCR master mix (TransGen, Beijing, China) on a CFX96 real-time PCR detection system (BioRad, United States). The primer sequences of E2F2, E2F8, and GAPDH are listed in Table S2. Fold differences of mRNA were calculated using the 2^−ΔΔCT method and normalized to GAPDH levels. The same methods were used to determine the mRNA level of genes E2F2 and E2F8 in cells.

Cell phenotype assay

Human stomach carcinoma cell line AGS was cultured in RPMI 1640 medium (Gibco, United States) with 10% fetal bovine serum (FBS, Sijiqing, Hangzhou, China) for in vitro experiments. AGS cells were then seeded into corresponding microplates for further detection 24 h after transfecting siRNA with RNATransMate (Sangon, Shanghai, China). Cell proliferation assays were performed using a Cell Counting Kit-8 (CCK8) colorimetric assay (Fude, Hangzhou, China) at a wavelength of 450 nm by a Spark microplate reader (TECAN, Switzerland). A colony formation experiment was also performed to determine cell proliferation. For transwell migration assays, transfected cells were seeded into 8-µm pore inserts (Falcon, United States). To induce cell migration, 1% FBS was added to the top chamber, and 10% FBS was added to the bottom chamber. After 48 h of migration, the cells were stained with 0.1% crystal violet (Biosharp, Anhui, China) and observed under a microscope. Cell counts were obtained using ImageJ software. The data was evaluated by normalizing to cells transfected with the negative control siRNA. The siRNA sequences of E2F2 and E2F8 are listed in Table S3. The overexpressing plasmids were purchased from Weizhen Biosciences (Jinan, China) and then used to transfect AGS cells to construct cell lines overexpressing E2F2 and E2F8. All plasmids were confirmed by DNA sequencing, and detailed plasmid sequence information is included in Table S4. After transfection by Lipo3000 (Invitrogen, Waltham, MA, USA), cells were used for proliferation and migration assays. The experiment and analysis methods are the same as above.

Statistical analysis

All statistical analyses were performed by the R software package, version 4.1.2 (R Core Team, 2021). Nonparametric tests, such as the Wilcoxon rank-sum test (comparison between two groups) or the Kruskal–Wallis test (comparison among multiple groups), were used for continuous variables (such as expression quantity and infiltration ratio). The survival analysis was performed using a Kaplan–Meier (K-M) procedure (log-rank test). The univariate Cox regression model was employed to calculate hazard ratios, while the multivariable Cox regression model was adopted to identify independent prognostic factors. All statistical tests were two-sided, with statistical significance set at P ≤ 0.05.

Expression profile and heritable variation of E2Fs in gastric cancer

This study identified the significant differential expression of total E2F transcription factors (E2F1–E2F8) between GC and normal samples from the TCGA dataset, suggesting the dysregulation of E2F expression plays a crucial role in the occurrence and development of GC (Fig. 1A).

A methylation analysis of these normal and tumor samples showed that the probes in the 1,500 bp range of the eight transcription factor promoter regions had minor changes. Even the highest E2F6 had a change of only 0.2 (Fig. 1B). The somatic mutation frequency of eight E2Fs in GC samples from the TCGA dataset were then calculated. Somatic mutation was observed at a low frequency. Only E2F7 reached 2.17%, and other E2Fs reached almost zero (Fig. 1C). Spearman correlation was calculated for the expression level of E2Fs, and there was a significantly positive correlation between the eight factors, indicating potential synergy (Fig. 1D). These results suggest a limited contribution of the E2F transcription factor family to the tumor immune microenvironment (TIME) at the genomic level, indicating the influence of E2Fs on TIME may mainly be reflected in the expression level.

Association of E2F with immunity and prognosis

Based on The Cancer Immunome Atlas, an ssGSEA analysis was used to profile the infiltration degree of 28 kinds of immune cells in GC samples, most of which were then verified to be correlated with the expression level of E2Fs (Fig. 2A).

The samples were then divided into higher- and lower-expression groups according to the median expression level of each transcription factor and the clinical and transcriptome data of the TCGA-STAD cohorts were transformed into K-M survival curves to estimate the prognostic value of the eight E2Fs (Fig. S1). Univariate Cox regression demonstrated that E2F2 and E2F8 showed a protective influence on survival, but the other E2Fs showed no independent contribution to survival (Fig. 2B). These results indicate that the eight E2Fs played a role in the regulation of TIME through a collaborative model and further affected clinical outcomes.

Identification of consensus clusters of E2F and corresponding clinical features

An unsupervised consensus clustering analysis identified two subgroups (Cluster1 and Cluster2) of GC patients associated with the expression level of the eight E2Fs supported by the cumulative distribution functions (CDF; Fig. 3A and Fig. S2). PCA was also conducted to verify the independence and diversity of the transcriptional profile in the two clusters (Fig. 3B). The differential expression levels of the eight E2Fs between the two E2FClusters are shown in the heat map in Fig. 3C. Poor tumor differentiation (Grade 3/4) dominated Cluster2, and patients with GC molecular subtype were characterized by Cluster2’s expression pattern, while patients with HM-indel and HM-SNV subtypes aligned with Cluster1’s expression pattern. Previous studies regarding molecular subtypes suggest that GC subtypes are significantly associated with adverse clinical outcomes, whereas MSI (HM-indel and HM-SNV) contribute to improved survival (Sohn et al., 2017; Cristescu et al., 2015).

K-M curves confirmed that between the two types of E2FClusters identified, Cluster1, with 229 samples, had better survival than Cluster2, with 122 samples (Fig. 3D). A multivariate Cox analysis combining multiple clinical features revealed that the expression patterns of the transcription factors identified were independent prognostic factors (Fig. 3E). These results indicate that the Cluster2 E2F expression pattern is an independent risk factor for GC prognosis.

Immune infiltration and genetic differences in E2F clusters

Cluster2 had elevated immune scores along with stromal scores using the “ESTIMATE” algorithm to assess immune cell infiltration (Fig. 4A). Previous studies have found that immune cells can infiltrate tumors with immune rejection phenotype, but they are confined to the fibroblast-and collagen-rich peritumoral stroma and cannot penetrate into the core of tumor cells (Mariathasan et al., 2018). Activated stromal cells in the TME are known to inhibit T cells (Lakins et al., 2018), so stromal activation in Cluster2 may inhibit the anti-tumor effect of immune cells. A hallmark pathway analysis demonstrated that stromal activity in Cluster2 was significantly enhanced, including epithelial-mesenchymal transition (EMT), transforming growth factor-β (TGF-β), angiogenesis, and many more, confirming that stromal activation may help explain the negative outcomes of patients in Cluster2 (Fig. 4B).

To further investigate the potential biological roles of genes characterized by E2Fs mode, a total of 1,340 differentially expressed genes (DEGs) were screened between the two clusters (Fig. 4C). A KEGG analysis indicated that these genes are involved in the cell cycle, DNA replication, mismatch repair, and other cell progressions, confirming that E2F transcription factors are involved in the immune regulation of TME (Fig. 4D). According to the unsupervised clustering of these E2F DEGs, the 351 GC patients were classified into three E2F-related gene clusters (E2F gene clusters 1, 2, and 3), which were strongly correlated with the expression of the eight E2Fs. The K-M curves indicated that these gene clusters were also correlated with survival status: E2F gene cluster 1 had the best OS when compared to clusters 2 and 3 (Figs. 4F–4G and Fig. S3). A univariate Cox analysis was performed on the DEGs and 58 were found to be correlated with prognosis (Fig. 4E).

Establishment and validation of the E2F scoring system

PCA was used to calculate an E2F score based on the 58 genes associated with prognosis to quantify the E2F transcription factor expression patterns in individual gastric tumors. Cluster1 scored significantly higher than Cluster2 (Fig. 5A). Also, E2F gene cluster 1 scored the highest among the three gene subtypes. Based on the median PCA score, these 351 patients were stratified into two groups (high and low E2F score). The OS of the low-E2F score group was worse than the high-E2F score group (Fig. 5B). A multivariate Cox analysis showed that E2F score was an independent prognostic factor, indicating that the E2F score model could characterize the prognosis of the sample with decent accuracy (Fig. 5C). These associations were then validated in six independent GEO datasets (GSE13861, GSE15459, GSE26899, GSE26901, GSE66229, and GSE84433). All the validation sets revealed consistent findings: patients with a high E2F score had a distinctly better prognosis than patients with a low E2F score, which reflected the stability and credibility of the E2F scoring system (Fig. 5D).

Molecular mechanism of the E2F scoring system

A gene-set enrichment analysis was used to assess molecular pathway enrichment to help explain the mechanisms behind the E2F score. The pathway analysis results showed that the group of genes with low E2F scores were enriched in matrix activation pathways involving the calcium signaling pathway, adhesion molecules, extracellular matrix (ECM) receptor interaction, and the MAPK signaling pathway (Fig. 6A). A correlation analysis between E2F score and hallmark pathway analysis indicated that E2F score was negatively correlated with the matrix activation pathways, such as the NOTCH pathway and IL2 STAT5 pathway, and positively correlated with immune-related pathways, such as DNA repair (Fig. 6B). In terms of immune cell infiltration, levels of activated CD4 T cell infiltration in the high-E2F score group were significantly higher than in the low-E2F score group, while the levels of T helper cell infiltration and other immune cell infiltration were significantly lower in the high-E2F score group than in the low-E2F score group (Figs. 6C–6D). Additionally, a higher tumor mutational burden (TMB) was observed in the high-E2F score group (Fig. 6E). TMB leads to tumor-specific and potentially highly immunogenic neoantigens, which can be recognized by the major histocompatibility complex (MHC) to enhance the immune response and improve prognosis. The distribution of these different classification systems and subtypes are shown in the alluvial diagram in Fig. 6F.

Treatment response prediction and validation of the E2F scoring system

The E2F scoring system reflected the regulation and modulation of the TIME, so this study further explored the promising role of E2F score in predicting therapeutic response and monitoring the effectiveness of treatments. Combined with drug information and expression profiles in the training set samples of GDSC, the Spearman correlation between the IC50 of several common GC drugs and E2F score was calculated. The results showed that the IC50 values of 5-Fluorouracil, Cisplatin, Docetaxel, and Paclitaxel were significantly negatively correlated with E2F score (Fig. 7A). Additionally, in the immunotherapy cohort database, a high E2F score was correlated with a better OS when compared with a low E2F score (Fig. 7B). When assessing patient response to immunotherapy, responders had significantly higher average E2F scores than non-responders (Fig. 7C). There were also significant differences in the distribution of high and low E2F scores among patients with distinct responses, and most non-responders had low E2F scores (Fig. 7D). These results showed that the higher the E2F score, the higher the sensitivity to the drugs, and the better the patient prognosis, confirming that E2F score is an underlying indicator for predicting chemotherapy and immunotherapy responses.

In vitro validation

Given the independent prognostic role of E2F2 and E2F8 in the database, these two transcription factors were chosen for validation in the collected samples. The qPCR assay results showed the mRNA levels of E2F2 and E2F8 were significantly upregulated in the gastric cancer tissues compared with the corresponding normal tissues, which was consistent with both the transcriptomic database and the findings of previous studies (Fig. 8A).

To validate the prognostic implications of E2F2 and E2F8 in gastric cancer, E2F2 and E2F8 were transiently knocked down with siRNA in AGS cell lines, and the cells were phenotypically and functionally characterized. Transfection effectively suppressed the RNA levels of E2F2 and E2F8 compared to the siRNA negative control (Fig. 8B). Cell proliferation was assessed at 24 h, 48 h, and 72 h after transfection. Compared with the negative control, the knockdown of E2F2 resulted in the significant loss of tumor cell proliferation, while the knockdown of E2F8 increased tumor cell proliferation (Fig. 8C). Colony formation experiment was also performed to determine cell viability. E2F2-deficient cells formed significantly fewer colonies, but the opposite was observed in E2F8-deficient cells (Fig. 8D). The in vitro transwell migration assay results showed that the low expression of E2F2 inhibited migration, and the low expression of E2F8 induced migration (Fig. 8E). We further overexpressed E2F2 and E2F8 by transfecting AGS cells with corresponding plasmids, using the empty vector as control. The efficacy of overexpressed plasmids was validated by qPCR and the RNA levels of E2F2 and E2F8 significantly increased (Fig. 8F). Compared with the negative control, the overexpression of E2F2 accelerated tumor cell proliferation, while the overexpression of E2F8 inhibited tumor cell proliferation (Fig. 8G). As for Colony formation experiment, AGS cells overexpressing E2F2 displayed stronger cell viability than negative control, while overexpression of E2F8 reduced cell viability (Fig. 8H). In vitro cell migration experiments demonstrated that overexpressed E2F2 increased cell migration, and overexpressed E2F8 suppressed migration (Fig. 8I). These results suggested that E2F2 and E2F8 affect tumor cell migration and proliferation, supporting their prognostic impact.