Figure S1 showed a flowchart of this study. The Lung adenocarcinoma (LUAD) mRNA expression profiles and related clinical data were collected from The Cancer Genome Atlas (TCGA) data portal. The mRNA expression data were transformed to values in transcripts per million (TPM). Gene Expression Omnibus (GEO) databases: GSE135222 and GSE126044 were downloaded from the GEO repository. Angiogenesis-associated genes (AGG) were obtained from GSEA Molecular Signatures Database (MSigDB-hallmark gene sets), a total of 36 genes in “HALLMARK_ANGIOGENESIS” pathway was extracted from the hallmark sub database.

Unsupervised clustering was used to separate patients into two groups, named angiogenesis-associated subtypes (AGSs)]. Next, a Kaplan-Meier (K-M) plot was developed to measure the OS differences. AGG that met the criterion of false discovery rate (FDR) <0.05 between AGS1 and AGS2 were identified as differentially expressed, and enrichment analysis was performed subsequently.

Infiltrating stromal cells and immune cells constitute the main part of normal cells in tumor tissue. They not only disrupt tumor signals in molecular research, but also play an important role in cancer biology. Based on “Estimating stromal cells and immune cells in tumors using expression data” (ESTIMATE) - a method that uses gene expression characteristics to infer the proportion of stromal cells and immune cells in tumor samples, we calculated the stromal and immune scores of AGSs.

Next, we obtained five immune expression signatures: wound healing (W.H), macrophage regulation (M.R), lymphocyte infiltration (Lym. Inf), IFN-γ response (IFN-γ), and TGF-β response (TGF-β); C1-C6 immune subtypes of each TCGA LUAD samples from Thorsson et al. (12). Moreover, we compared them between AGSs. Considering that immune infiltration was related to measures of DNA damage, we compared several measures of DNA damage, including fraction alteration, number of segments, aneuploidy score (AS), homologous recombination defects (HRD), and intratumor heterogeneity (ITH), which also achieved from Thorsson et al.

The whole-exome sequencing (WES) data obtained from the TCGA database was used to describe the mutation spectrum and interaction between AGSs, and then searched for specific hypermutation. Next, we calculated the tumor mutation burden (TMB) score with exons uniformly counted as 40M regions and compared it between two AGSs.

We performed ssGSEA to calculate scores of immune cells, whose markers were achieved from TISIDB (http://cis.hku.hk/TISIDB/), and compared their scores between AGSs. To assess the immune state for AGSs, we compared the expression of the following marker genes of myeloid-derived suppressor cells, which have been proven causes the immunosuppressive microenvironment: CD33, ITGAM, OLR1, S100A9 (13, 14).

We collected four major metabolic pathways from MSigDB and previous studies (15, 16), including glycolysis, fatty acid oxidation, pentose phosphate pathway, and glutamate metabolism ( Table S1 ).

Furthermore, we used the “pRRophetic” package (17) to predict the treatment response of each LUAD sample, which was determined by the half maximal inhibitory concentration (IC50) based on the Genomics of Drug Sensitivity in Cancer (GDSC) database.

The Univariate and multivariate cox regression analysis was performed to delineate the prognostic gene signature. Then, we analyzed the correlation of signature gene expression and measures of immune state, including lymphocyte infiltration and TGF-β pathway activation. Finally, JAG1 was selected as signature gene because its expression was significantly negatively correlated with lymphocyte infiltration and positively with TGF-β pathway.

To validate the function of JAG1, we analyzed the correlation between its expression and markers of basal cells, including endothelial cells: PECAM1, CD34, ICAM1, PTPRC; fibroblast cell: LUM, DCN, COL1A1; mesenchymal stem cells: ZEB1, MME, ANPEP, ITGB1, THY1, TCF4, and SOX2. Furthermore, we divided the patients of TCGA LUAD datasets into two groups based on JAG1 expression, samples with expression greater than the mean value were classified as JAG1-high group, otherwise as JAG1-low group. Then we compared their scores of glycolysis and fatty acid oxidation between JAG1-high and JAG1-low groups. Finally, we analyzed the JAG1 expression in two immunotherapy cohorts. We also performed a correlation analysis between JAG1 and immunotherapy signatures expression including CD274 (PD-L1) and CTLA4 (18).

The single-cell transcriptome data obtained from GSE131907, which included 208,506 cells derived from 58 lung adenocarcinomas from 44 patients, which covered primary tumor, lymph node and brain metastases, and pleural effusion in addition to normal lung tissues and lymph nodes. 12 of them were primary tumor samples. In order to reduce the heterogeneity among samples, we selected them for the following analysis. The single-cell transcriptome data were preprocessed using the “Seurat” package, and the cell annotation used results provided by the authors. We calculated percentages of cell types among samples, and analyzed their correlation with JAG1 positive cell rate. To explore the relationship between T cell status and JAG1 expression. We used ProjecTIL to parse human scRNA-seq T cell data in the context of murine TIL profiles (19), and described the atlas of T cells. Then, we researched the relationship between JAG1 positive ratio and subtype percentage. Furthermore, we counted the proportion of each subgroup of T cells in JAG1 positive and negative cells, and compared the T cell subgroup composition in them. Finally, we calculated the active scores of four metabolism pathways, then compared them between JAG1 positive and negative cells.

First, the positive rate of JAG1 at the single-cell level was not considerable overall, so we observed the distribution of JAG1-positive cells among various cells types by JAG1 staining. Considering that the clinical assay usually detects on tumor cells, we divided the samples into two groups by the proportion of positive cells in the epithelium. The cell-cell interaction strength between the two groups was analyzed by the “Cellchat” R package, which predicted major signaling inputs and outputs for cells and how those cells and signals coordinate for functions using network analysis and pattern recognition approaches (20). Further, based on the ligand receptor database, significant signal flow and ligand-receptor interaction in the positive group were found out.

All analyses were performed using R software (version 3.6.1). Differences between groups were evaluated using Wilcoxon rank-sum tests for continuous data and Fisher’s exact tests for categorical variables. Pearson’s test was used for the correlation analysis. The center point-based partitioning method with the pam() function in the “cluster” package was used to perform unsupervised clustering. ssGSEA scores were calculated using the GSVA package. K-M plots were constructed using the Kaplan–Meier “survival” package. All analyses were two-sided, and statistical significance was set at P < 0.05.

Firstly, we divided lung adenocarcinoma patients into two groups based on 36 angiogenesis-associated genes, which we named angiogenesis-associated subtypes (AGSs). 11 genes of AGGs showed significant differences between two AGSs, including JAG1, VCAN, POSTN, COL3A1, COL5A2, FSTL1, ITGAV, FGFR1, STC1, SPP1, and SERPINA5 (P _all < 0.05, Figure 1A ). Overall, AGS2 exhibited an abundant enrichment of AGGs. Importantly, we found that patients with AGS2 had a poorer prognosis than those with AGS1 ( Figure 1B ), and enrichment analysis showed AGS2 significantly active PI3K-Akt signaling pathway and epithelial mesenchymal transition ( Figure 1C , Figure S2A ). Additionally, we compared tumor purity and immune scores between the AGSs, and the results showed that patients in AGS2 had higher stromal scores (P < 0.05) but no difference in immune scores ( Figure 1D ). More importantly, we analyzed the five immune expression signatures between AGSs, and found that AGS2 had significantly lower lymphocyte infiltration, accompanied with higher TGF-β response, which indicated AGS1 is an immune “hot” phenotype, while AGS2 is an immune “cold” phenotype (P _all < 0.05, Figure 1E ). All above mentioned suggested that angiogenesis-associated genes enrichment may reprogram the “cold” phenotype of tumor immune microenvironment, thus leading to different survival outcomes.

To research the competence of tumor microenvironment, we calculated the ssGSEA score of 28 immune cells based on markers from the TISIDB database. The comparison results showed that AGS2 exhibited obviously lower activated CD8 T cells, and a variety of other immune cells also showed significant differences ( Figure 2A ). Considering the importance of myeloid-derived suppressor cells (MDSCs), which were known to suppress immune responses by inhibiting T cell proliferation and activation, we compared markers of MDSCs between AGSs. We found that CD33, ITGAM, OLR1, and S100A9 were significantly expressed higher in AGS2 (P = 2×10^-4, 2.86×10^-10, 1.47×10^-10, 8.2×10^-3, respectively, Figures 2B–E ).

The comparison of DNA damage measures, including AS, HRD, ITH, altered fraction, and number of segments, were all significantly stronger in AGS2 than in AGS1, representing a poorer prognosis of tumor patients with immunosuppressive state (P _all < 0.05, Figures 3A–E ). Next, we described the mutation spectra of AGSs, and revealed several high-frequency mutations. For example, FLG was specifically mutated in AGS2 ( Figure 3F ). Additionally, the mutation interactions were weaker in AGS2 ( Figure 3G ). There was no difference in TMB between AGSs ( Figure 3H ).

Considering that metabolic activation could affect the tumor microenvironment, we compared scores of four major metabolism pathways between AGSs. The result demonstrated that glycolysis is obviously activated in AGS2. Several studies have shown that tumor cells tend to select the metabolic mode of glycolysis, which leads to the accumulation of lactic acid in the microenvironment, thereby promoting the function of immunosuppressive cells such as TAMs, MDSCs and Treg (13). At the same time, fatty acid oxidation was suppressed in AGS2 (P _all < 0.05, Figure 4A ). Drug susceptibility prediction was performed using “pRRophetic” package based on the GDSC database. Screening of drugs sensitive to ARG2 group patients, and the result showed that patients in AGS2 were more sensitive to 15 drugs than ARG1, such as such as ABT 263 (navitoclax), cisplatin, dasatinib AP.24534 (ponatinib), and midostaurin (P _all < 0.05, Figure 4B ).

Multivariate regression analysis was performed to identify hub regulatory factors. JAG1 and STC1 showed significant effect on prognosis (hazard ratio (HR): 1.20, P = 0.028; HR: 1.15, P = 0.021, respectively, Figure 4C ). Kaplan–Meier plot of the comparison in the overall survival between JAG1-high and JAG1-low groups shown in Figure 4D . JAG1, a ligand of the Notch signaling pathway, regulates cell differentiation and proliferation in various cancers. To verify the relationship between JAG1 expression level and immune microenvironment, we performed the correlation analysis, which showed that the expression of JAG1 was obviously negative with lymphocyte infiltration, and positive with TGF-β response (R=-0.13, P = 5×10^-3; R=0.25, P = 3.246×10^-8, respectively, Figures 4E, F ). Furthermore, the EMT pathway was also positively correlated with JAG1 expression obviously (R=0.23, P=1.431×10^-7, Figure 4G ). However, the expression of STC1 did not show the same consequence ( Figures S3A, B ). Moreover, the JAG1 expression was a negative correlation with immune cell score, positive correlation with stromal score, although not significant (R = -0.04, P = 0.34; R = 0.06, P = 0.16, respectively, Figures S3C, D ). Assessing the level of JAG1 in multi-omics, we found that JAG1 can be detected from either copy number amplification or methylation level ( Figures S3E, F ).

Previously, we found more stromal cells in AGS2, and we explored the relationship between markers of various types of stromal cells and JAG1 expression. The result demonstrated that most markers positively interrelated with JAG1, especially mesenchymal stem cells, such as ITGB1 ( Figure 5A ). Next, we compared the activation of glycolysis and FAO pathway between JAG1-high and JAG1-low groups, which revealed that JAG1-high patients had stronger activated glycolysis than JAG1-low. (P < 0.05, Figures 5B, C ). In order to explore the efficacy of immunotherapy, we used two datasets of PD1/PD-L1 treatment cohorts to compare the expression of JAG1 between treatment-sensitive and resistant groups, accompanied with a correlation analysis of JAG1 expression with PD-L1 and CTLA4 expression. The results indicated that JAG1-high patients had higher expression of PD-L1 and CTLA4 than JAG1-low group (R=0.18, P = 5.563×10^-5; R=0.09, P = 4.149×10^-2, respectively, Figures 5F, G ). Moreover, JAG1-high patients were resistant to pure anti-PD-L1/anti-PD1 therapy ( Figures 5D, E ). The treatment resistance is likely a glycolysis-mediated immunosuppressive microenvironment resulting from angiogenesis.

To explore T cell status at the single-cell level accurately, we described the single-cell transcriptome atlas of lung adenocarcinoma from GSE131907. The cell annotation was downloaded from the original research. Firstly, we counted the compositional proportions of cell types among samples and sorted them by the proportion of T lymphocytes. Explore the relationship between the rate of JAG1+ cells and T lymphocytes, we found that the patients with a higher percentage of JAG1+ cells were mainly characterized by lower T cell infiltration ( Figures 6A, B ). Comparing the correlation between the rate of JAG1+ cells and the proportion of each cell type, we found that it was significantly positively correlated with malignant epithelial cells, and negatively with T lymphocytes, indicating that patients with high JAG1 expression had higher tumor purity, lower T cell infiltration ( Figure 6C ). Furthermore, we annotated subpopulations of T lymphocytes by ProjecTIL, and divided them into nine subgroups, including CD8+ terminally-exhausted (Tex) effector cluster, CD8+ precursor-exhausted (Tpex) cluster, CD8+ effector memory cluster, CD8+ early active cluster, CD8+ naïve cluster, CD4+ naïve cluster, CD4+ follicular-helper (Tfh) cluster, CD4+ Th1-like cells, and a cluster of regulatory T cells (Treg). The TSNE plot showed that the distribution densities of various T cell types in the JAG1+ and JAG- subsets were quite different ( Figure 6D ). At the same time, we calculated the fold change value of the proportion of T cell subsets in JAG1+ and JAG1- subsets, which manifested that Treg is enriched in the JAG+ subgroup, CD8+ early active cells are enriched in the JAG1- subgroup ( Figure 6E ). Finally, we evaluated the activation of metabolic pathways between JAG1+ and JAG1- groups, and found that glycolysis, pentose phosphate pathway, and glutamate metabolism activated special in JAG1+ cells, while fatty acid oxidation was activated in JAG1- cells ( Figure 6F ).

To dissect the interactions of the tumor microenvironment, we explored the intercellular communication between the JAG-positive and JAG1-negative groups. First, we stained the expression of JAG1 on umap plot, and the results showed that: JAG1+ cells were mainly distributed in epithelial cells, endothelial cells, fibroblasts and myeloid cells ( Figure 7A ). Next, we counted the proportion of JAG+ cells in epithelial cells in each sample and sorted them. Seven samples had a very low proportion of positive cells (<10%), and we defined them as JAG1-negative samples and the rest as JAG1-positive samples ( Figure 7B ). Then, we compared the cell-cell interaction strength between the two groups. Red represented the stronger interaction in the JAG1-positive samples in the heatmap. The results showed that the JAG1-positive samples had stronger cell-cell communication ( Figure 7C ). Furthermore, we classified the ligand-receptor interaction and explored the signal flow that was significantly activated between the two groups. The results showed that most signal flows in the JAG1 positive samples was stronger, including: SPP1 and SELL signal flow, which related to immunosuppressive microenvironment, and tumor angiogenesis-related pathway: NOTCH, VEGF, FGF signal flows ( Figure 7D ). Finally, we explored the significantly activated ligand-receptor interaction in JAG1-positive samples. The result showed that the HLA-A, HLA-B, HLA-C of epithelial cells, was significantly strongly interacted with CD8A of T cells in JAG1-negative samples, but not in JAG1-positive samples. Moreover, SPP1- (ITGAV+ITGB5), SPP1-(ITGAV+ITGB1), and SPP1-(ITGA8+ITGB1), which was associated with immunosuppressive, occurred in JAG1-positive group. The CD46-JAG1 ligand-receptor interaction occurred in epithelial cells and endothelial cells was obvious in JAG1-positive group, which suggested that targeting CD46 in tumor cells may decrease the expression of JAG1, thereby reducing the immunosuppressive phenotype ( Figure 7E ).