The genomic and transcriptomic landscapes of the Notch pathway genes

The Notch pathway genes may display similar genomic and transcriptomic features in the GIACs at different anatomic locations. In the TCGA database, mutational commonalities among GIACs were observed, such as the mutations in CREBBP, FBXW7, CNTN1/6, NCOR2, and EP300 (Fig. 2A). Moreover, the frequencies of CNVs among GIACs were highly consistent (Fig. 2B). Most mutational events were not associated with the changes of mRNA expression (Fig. 2C [COAD]; Supplemental Fig. S1A [ESAD, STAD, and READ]), while CNVs were largely linked with transcriptional changes (Fig. 2D [COAD]; Supplemental Fig. S1B [ESAD, STAD, and READ]).

Fig. 2

Genomic alterations of 60 Notch pathway genes in gastrointestinal adenocarcinomas. A-B. Mutational events (A) and copy number variations (B) of the 60 Notch pathway genes in GIACs. C-D. Associations of the mRNA expression with mutational events (C) and copy number variations (D) of the 60 Notch pathway genes in COAD. Abbreviations: COAD = colon adenocarcinoma, GIAC = gastrointestinal adenocarcinoma

Among the 60 Notch pathway genes, only 5 showed a significant difference in mRNA expression among GIACs (Supplemental Table S2). Of these, lower-GIACs expressed higher levels of RBPJL, CNTN6, and DLL3, and lower levels of PTCRA and DTX4 (Supplemental Fig. S2). For the key genes involved in the Notch pathway, higher expression levels of NOTCH1, NOTCH2, and NOTCH3 are observed in STAD compared to normal tissues (Supplemental Fig. S3). In addition, the pairwise correlations between the mRNA level of each Notch pathway gene were also similar among GIACs (Supplemental Fig. S4). Most expression levels were positively correlated with each other, while the correlations were relatively weaker in the “HES” genes (HES1, HES2, HES4, HES5, and HES7). Taken together, the similar genomic and transcriptomic features of the Notch pathway genes among GIACs might suggest the analogous role of the Notch signaling in promoting carcinogenesis along the gastrointestinal tract.

The Notch pathway defines an aggressive subtype of GIACs with immune suppression

Given the shared features of the Notch pathway genes in GIACs, we sought to perform a joint analysis of all GIAC samples ranging from ESAD to READ to investigate the impact of the Notch pathway in GIACs. Unsupervised clustering based on the mRNA levels of the 60 Notch pathway genes classified all GIAC samples into three clusters (determined by delta area plot, Supplemental Fig. S5) spanning anatomic boundaries (Fig. 3A). The cluster 3 with lower expression of the Notch pathway genes, compared to the other two clusters, had a higher fraction of lower-GI adenocarcinomas (P = 0.035), more senior patients (P = 0.004), and a slightly earlier TNM stage (P = 0.044, Supplemental Table S3). As for the pan-GIAC subtype suggested by Liu et al. [30], the EBV, MSI, and GS subtypes were, respectively, enriched in the cluster 1, 2, and 3, while the proportions of CIN and HM-SNV subtypes were comparable among three clusters (Supplemental Table S3). In addition, a poor prognosis was observed in the cluster 1 with higher expression of the Notch pathway genes (log-rank P for trend = 0.044, Fig. 3B).

Fig. 3

Prognostic and molecular correlates of the Notch pathway genes in GIACs. A. Heatmap illustrating the clinicopathological features and mRNA expression of 60 Notch pathway genes in the GIACs of the TCGA dataset. B-E. Associations of the clusters with OS (B), immune cell infiltration (C), and expression levels of immune checkpoints (D), and ssGSEA scores of multiple signatures (E). Abbreviations: GIAC = gastrointestinal adenocarcinoma, OS = overall survival, TCGA = The Cancer Genome Atlas

As for immune characteristics, the cluster 1 tended to have more “resting”/ “anti-inflammatory” rather than “activated”/ “pro-inflammatory” tumor-infiltrating immune cells (TIICs, Fig. 3C). For instance, more resting memory CD4+ T cells, resting natural killer (NK) cells, M0 macrophages, and resting mast cells, and fewer memory B cells and plasma cells were observed in the cluster 1. In addition, the difference between the cluster 1 and the clusters 2/3 of the infiltrating level of anti-inflammatory M2 macrophage is numerically higher compared to pro-inflammatory M1 macrophage. Moreover, nearly the mRNA levels of all immune checkpoints were higher in the cluster 1 (Fig. 3D).

To delineate the molecular features of these clusters, we compared the ssGSEA scores of several cancer-related signatures among these three clusters. The cluster 1 had higher scores of the signatures concerning Notch, cell cycle, tyrosine receptor kinases, VEGF, RAS/RAF, mTOR, and TGF-beta (Fig. 3E) which attenuates tumor response to PD-1 blockade [18]. As for metabolism signatures, a higher score of glycolysis rather than tricarboxylic acid cycle (TCA) was observed in the cluster 1 (Fig. 3E), suggesting the association between Notch activation and anaerobic metabolism in GIACs. Low-glucose and high-lactate environment is immunosuppressive [31], and a high pre-treatment level of serum lactate dehydrogenase is linked with poor immunotherapy efficacy [32]. Taken together, these results of TIICs, immune checkpoints, and molecular correlates suggest the associations of Notch activation with immune suppression and potentially poor response to immunotherapy in GIACs.

Identification of a 6-mRNA risk stratification signature for GIACs

Of the Notch pathway genes, 16 had a consistent trend of prognostic power in at least three of the four types of GIAC in the TCGA dataset (Fig. 4A). Using the lambda.min cutoff threshold, LASSO penalized Cox regression analysis on the 16 genes identified a 6-mRNA signature including APH1B, CNTN1, DLL4, HES7, HEYL, and SNW1 (hereafter referred to as “mRNA-score,” Fig. 4B), which exhibited comparable prognostic power in all four types of GIAC (Fig. 4A). We calculated the prognostic effect of the mRNA-score at every cutoff value ranging from 10th percentile to 90th percentile and observed a robust association between a high mRNA-score and poor OS in the training set (Fig. 4C). These findings were validated in multiple additional cohorts, including GSE19417 (ESAD, Fig. 4D), GSE84437 (STAD, Fig. 4E), and GSE40967 (COAD, Fig. 4F). Patients were divided into three groups by tertiles, and the patients with a high mRNA-score had poor prognosis in the training (log-rank P for trend < 0.001, Fig. 4G) and the three independent validation cohorts (GSE19417: log-rank P for trend = 0.025, Fig. 4H; GSE84437: log-rank P for trend = 0.001, Fig. 4I; GSE40967: log-rank P for trend = 0.007; Fig. 4J).

Fig. 4

Development and validation of the mRNA-score. A. Prognostic effects of the 16 Notch pathway genes with a consistent trend of prognostic power in at least three of the four types of GIACs in the TCGA dataset. B. Coefficients of the six genes involved in the prognostic model. C-F. Association between the mRNA-score and OS in the training set (TCGA [C]) and 3 independent validation sets (GSE19417 [D]; GSE84437 [E]; and GSE40967 [F]). G-J. The OS of the three groups determined by tertiles in the training set (TCGA [G]) and 3 independent validation sets (GSE19417 [H]; GSE84437 [I]; and GSE40967 [J]). Abbreviations: GIAC = gastrointestinal adenocarcinoma, OS = overall survival, TCGA = The Cancer Genome Atlas

To investigate whether the prognostic effect of the mRNA-score is independent of other variables, we firstly analyzed the association between clusters and key clinicopathological variables and then performed subgroup analysis and multivariable analysis. The mRNA-score was higher in male patients and was positively correlated with T and N stages rather than M stage (Fig. 5A), indicating the association of the Notch signaling with local advancement instead of distant metastasis. The mRNA-score showed comparable associations with OS in the different subgroups classified by age, sex, tumor type, TNM stage, and pan-GIAC subtype (Fig. 5B). In addition, the prognostic effect of the mRNA-score was independent after adjusting for these potential covariates (multivariable HR = 1.62, 95% CI 1.15–2.28, P = 0.005, Fig. 5C).

Fig. 5

Prognostic performance of a nomogram involving the mRNA-score in GIACs. A. Associations of the mRNA-score with sex, age, TNM stage, and pan-GIAC subtype. B-C. Prognostic performance of the mRNA-score assessed by subgroup analysis (B), and univariable and multivariable analyses (C). D. A nomogram involving the mRNA-score with other key parameters. E. The calibration curves illustrating the performance of the model-score in predicting 1-, 3-, and 5-year OS. F. Area under the curve of the model score in predicting OS. G. Association between the mRNA-score and OS in the TCGA dataset. H. The OS of the three groups determined by tertiles in the TCGA dataset. Abbreviations: GIAC = gastrointestinal adenocarcinoma, OS = overall survival, TCGA = The Cancer Genome Atlas

A nomogram involving the mRNA-score together with other key clinicopathological variables was constructed using the TCGA dataset based on the above multivariable model (Fig. 5D). The predicted survival and actual survival were highly consistent (Fig. 5E). The time-dependent AUROC was above 0.70 across time points (1 year: 0.743, 3 year: 0.778, 5 year: 0.706, Fig. 5F). The association between a high model score and short OS was remarkable at every cutoff (HR < 0.31, P < 0.001, Fig. 5G). All patients were divided by the tertiles of the model score, and the median OS was not reached, 1686 days and 675 days for the patients with low, medium, and high model score, respectively (Fig. 5H). Taken together, these results demonstrate the potential use of the Notch signaling for prognostication in GIACs.

mRNA-miRNA-lncRNA network and an 11-lncRNA risk stratification signature for GIACs

Based on the mRNA clusters, we further investigated the Notch pathway gene-related mRNA-miRNA-lncRNA network. In total, 6,607 mRNAs, 1,141 miRNAs, and 2,320 lncRNAs were identified with different expression (P < 0.05) between the mRNA clusters 1 and 3. According to the databases including miRTarBase, miRDB, and TargetScan, 110 miRNA-lncRNA pairs and 83 miRNA-mRNA pairs were observed (mRNA-miRNA-lncRNA network: Fig. 6A).

Fig. 6

The Notch pathway-related mRNA-miRNA-lncR NA network and the prognostic correlation of the Notch pathway-related lncRNAs. A. The Notch pathway-related mRNA-miRNA-lncRNA network and its schematic diagram. B. Heatmap illustrating the clinicopathological features and the expression of 110 lncRNAs in the GIACs of the TCGA dataset. C. Association of the clusters with OS. D. Coefficients of the 11 lncRNAs involved in the lncRNA-score. E–F. The OS of the three groups determined by tertiles in the training set (TCGA) and the independent validation set (GSE40967). Abbreviations: GIAC = gastrointestinal adenocarcinoma, OS = overall survival, TCGA = The Cancer Genome Atlas

There were 110 lncRNAs and 9 miRNAs in the mRNA-miRNA-lncRNA network. Given the small number of miRNAs, we focused on lncRNAs in the following analysis. Based on these 110 lncRNAs, all samples were classified into four clusters (determined by delta area plot, Supplemental Fig. S6A) by unsupervised clustering. Of these, one cluster had an extremely small sample size (n = 6, Supplemental Fig. S6B), which was merged with its closest cluster. The lncRNA cluster 2 with high lncRNA expression (Fig. 6B) had a higher fraction of the mRNA cluster 3 with low mRNA expression (P = 0.026), a lower mRNA-score (P < 0.001), and longer OS (P < 0.001) among the three clusters (Fig. 6C).

Of the 110 lncRNAs, 17 showed similar prognostic effects in the GIACs at different anatomic locations in the training set (the TCGA dataset). Based on these 17 lncRNAs, an 11-lncRNA signature was identified by LASSO (hereafter referred to as “lncRNA-score,” Fig. 6D). A high lncRNA-score was associated with poor OS in the training set (log-rank P for trend < 0.001, Fig. 6E). Among the validation sets (GSE19417, GSE84437, and GSE40967), only GSE40967 has available lncRNA data, in which the prognostic effect of the lncRNA-score was validated (log-rank P for trend = 0.004, Fig. 6F).

Associations of the mRNA-score with drug sensitivity and immunotherapy efficacy

Given the linkages of the Notch pathway with cell stemness and immunological regulation [33, 34], we hypothesized that the mRNA-score representing the Notch pathway may be associated with drug sensitivity and immunotherapy efficacy.

Based on the IC50 data of tumor cell lines in the GDSC database, we analyzed the associations of the mRNA-score with the sensitivities to chemotherapeutic drugs and targeted agents (Fig. 7A). As for the commonly used drugs, a high mRNA-score was linked with high sensitivities to cisplatin, oxaliplatin, paclitaxel, docetaxel, and irinotecan in the ESAD cell lines and paclitaxel, docetaxel, and lapatinib in the STAD cell lines.

Fig. 7

Associations of the mRNA-score with drug sensitivity and immunotherapy efficacy. A-B. Associations of the mRNA-score with drug sensitivity in the Genomics of Drug Sensitivity database (A), immune cell infiltration in the GIAC samples of the TCGA dataset (B), response to neoadjuvant chemoradiation plus atezolizumab in the patients with resectable ESAD (C), and histological grade and response to pembrolizumab in the patients with advanced STAD (D). E–G. Associations of the mRNA-score immune cell infiltration in the tumor samples other than GIACs of the TCGA dataset (E), OS on atezolizumab in patients with advanced UC (F), and OS on nivolumab in patients with advanced ccRCC (G). Abbreviations: ESAD = esophageal adenocarcinoma, GIAC = gastrointestinal adenocarcinoma, OS = overall survival, STAD = stomach adenocarcinoma, TCGA = The Cancer Genome Atlas, TRG = tumor regression grade. Other abbreviations of TCGA cancer codes are shown in Supplemental Table S4

In the TCGA dataset, a high mRNA-score was correlated with more “resting”/ “anti-inflammatory” rather than “activated”/ “pro-inflammatory” TIICs (Fig. 7B), similar to the result in Fig. 3C, which suggests that a high mRNA-score may predict a poor response to immunotherapy. In the PERFECT trial involving 32 patients with resectable ESAD who received neoadjuvant chemoradiation plus atezolizumab, the patients with a good pathological response (Mandard, tumor regression grade 1/2 [TRG1/2]) had a high mRNA-score (P = 0.048, Fig. 7C), and the AUC for predicting TRG1/2 was 0.71 (95% CI 0.53–0.90, P = 0.047, Fig. 7C). In another trial including 45 patients with advanced STAD who underwent pembrolizumab monotherapy, a high mRNA-score was associated with poor response (P = 0.049, Fig. 7D) rather than poor histological grade (P = 0.86, Fig. 7D). The AUC for predicting objective response was 0.69 (95% CI 0.54–0.83, P = 0.047), and the one for predicting disease control was 0.68 (95% CI 0.52–0.84, P = 0.044, Fig. 7D).

Moreover, we analyzed the association between the mRNA-score and the fraction of TIICs in other cancer types (Fig. 7E). The “immune-suppressive” relationship was observed in BLCA, BRCA, CESC, HNSC, KIRC, LUAD, LUSC, OV, PCPG, SARC, SKCM, TGCT, THCA, UCEC, and UCS (full name, see Supplemental Table S4). Given these, we further explored the association between the mRNA-score and immunotherapy efficacy in the other two cohorts, IMvigor210 (advanced urothelial cancer, n = 299, atezolizumab monotherapy) and CheckMate-009/010/025 (advanced clear cell renal cell carcinoma, n = 181, nivolumab monotherapy). A low mRNA-score was associated with a favorable OS on immunotherapy in both cohorts (IMvigor210: log-rank P for trend = 0.006, Fig. 7E; CheckMate-009/010/025: log-rank P for trend = 0.038, Fig. 7F). The above results further validated, to some extent, the associations of the Notch pathway with immune microenvironment and immunotherapy efficacy.

In summary, the established mRNA-score based on the Notch pathway genes might improve the patient selection for chemotherapy, targeted therapy, and anti-PD-(L)1 immunotherapy in pan-GIACs and potentially other solid tumors.