In total, 52 patients who were diagnosed with stage T1-2N0 primary laryngeal cancer were retrospectively and continuously collected from January, 2017 to January, 2020. All patients received laryngectomy without neo- or adjuvant chemotherapy or other systemic therapies, and all patients were followed up from the time of diagnosis until January 2022. A diagnosis of recurrence required clinical confirmation by computed tomography and magnetic resonance imaging exams. The collected formalin-fixed, paraffin-embedded (FFPE) samples were sent to Nanjing Simcere Diagnostics Co., Ltd. (Nanjing, China) for targeted DNA-based next-generation sequencing (NGS) analysis by TruSight Oncology 500 (TSO500) and NanoString RNA gene expression assays.

Considering the median time between primary diagnosis and local recurrence was 1.7 years (range 0.6–8.9 years) in 38 relapsed patients of total 303 patients with T1 glottic cancer (22), and the postoperative local recurrence rate of stage T1 supraglottic cancers is slightly higher than glottic cancers (1). Together with our clinical observations, we chosed 24 months as the relapse threshold in resected T1-2N0 laryngeal cancer to discuss the causes and mechanisms of rapid postoperative recurrence in our study. Thus, during the postoperative follow-up period, patients with local or distant recurrence within 24 months were classified into the relapsed group (n=23), and those without local recurrence more than 2 years after surgery were classified into the non-relapsed group (n=19). After excluding unqualified and unavaliable follow-up specimens, 42 postoperative FFPE samples were included in the molecular characterization analysis. 21 relapsed samples with qualified quality control were analyzed for immunologic features according to specifically defined criteria. The patient selection flowchart is shown in Figure 1 . This study was approved by the (blinded for peer review).

For clinical FFPE samples used for TSO500 evaluation, a minimum of 20% tumor content is needed. Genomic DNA (gDNA) was extracted from tissues using a tissue sample DNA extraction kit (Kai Shuo). All extraction procedures were carried out according to the manufacturer’s protocols. The DNA samples were captured and enriched by a two-step specific capture probe, and the DNA library and corresponding cDNA library were standardized by the library homogenization method, purified by magnetic beads, and sequenced using the Illumina NextSeq 550Dx platform. The sequencing results were analyzed by the TSO500 Docker pipeline. TSO500 NGS libraries were prepared using the TSO500 DNA Kit. Library preparation was performed manually according to the manufacturer’s protocol, with 24 samples per batch.

Prior to library normalization, the NGS libraries enriched by hybridization capture were quantified using the Qubit dsDNA HS Assay Kit. NGS sequencing was performed on the NextSeq 550 instrument (Illumina) with 8 libraries per sequencing. Base calls obtained through Illumina NovaSeq were converted to FASTQ files. The fastp tool (v.2.20.0) was used for adapter pruning and to filter low-quality substrates (23). The internal database was used to screen for germline variations. Copy number variation (CNV) involves amplification and deletion identified by CRAFT copy number callers from the TSO500 pipeline. Manta (Version 1.6.0) is used to detect large-scale structural variation (SV) in RNA libraries, and only fusions with at least three unique reads, one of which being a split read crossing a fusion break point, are considered candidate fusions (24). After screening for germline variations (internal database) and using the high COSMIC database to count mutations, the number of eligible somatic mutations (coding or high confidence regions, more than 50 times coverage, more than 5% VAF single nucleotide variations (SNVs) and INDELs were measured per megabase (Mb) by NGS (25). TMB measurements consider only SNV and INDELs in the coding region. Raw sequencing data and the TMB were evaluated using Illumina TSO500 support software. Based on the top 25% of the TMB scores, the 42 patients were divided into TMB-H and TMB-Low (TMB-L) groups. The TSO500 Docker pipeline was used to determine the sample microsatellite status in the form of a microsatellite instability (MSI) score. Samples with MSI scores ≥ 20 were considered microsatellite instable, and the rest were considered microsatellite stable.

Total RNA was isolated from 5-μm FFPE slices using the Qiagen RNeasy FFPE Kit, and 100 ng of RNA was hybridized to a version of the NanoString PanCancer code set for reading on the nCounter platform (NanoString). The expression of 289 immune-related genes, including housekeeping genes, was assessed ( sTable 1 ). Using nSolver 2.6 software, housekeeping genes were used to normalize the expression values as recommended by the manufacturer. According to the manufacturer’s specification, the genes are divided into 14 immune cell types (T cells, B cells, mast cells, dendritic cells, macrophages, neutrophils, cytotoxic cells, exhausted CD8 cells, NK-CD56 cells, CD8 T cells, CD45 cells Th1 cells, NK cells and Treg cells) ( sTable 2 ). Meagene scores were calculated based on the geometric mean expression levels of member genes (26). We also studied 10 previously published prognostic and immunotherapy response meta-genes. We calculated these metagene scores using the methods described in the respective publications ( sTable 3 ).

Differentially expressed genes (DEGs) were selected based on the NOTCH1 mutation status of relapsed laryngeal cancer patients through the “DEseq2” software package, with log2 |fold change| >1 and False Discovery Rate (FDR) <0.05. Heatmaps of DEGs were made using the “HeatMap” package. We then used Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to analyze the possible biological processes in which the overlapping DEGs are involved. Also, genomic enrichment analysis (GSEA) analysis was used to further explore gene expression differences in laryngeal cancer patients.

Continuous data with a normal distribution were presented as means and standard deviations and that with a screwed distribution were presented as medians and ranges. Categorical data were presented as frequencies or percentages. Differences between independent groups were statistically measured using the Wilcoxon Mann-Whitney test. Survival curves were drawn using the Kaplan-Meier method. Relapse-free survival (RFS) time was defined as the time between the diagnosis and confirmed disease relapse of the patient. Overall survival (OS) time was defined as the time between the diagnosis and death of the patient. A p-values < 0.05 was considered statistically significant. All statistical analyses were performed using the Graphpad (V 10.0), R (V. 4.1.0) and R Bioconductor packages (https://www.r-project.org).

Data from 52 patients with stage T1-2N0 laryngeal cancer were retrospectively and continuously collected from January, 2017 to January, 2020. Two patients with a history of other malignant tumor and four patients with unavailable follow-up were excluded. Besides, in the process of DNA-sequencing, two samples without sufficient proportion of tumors and another two samples with severe DNA degradation were also excluded. Finally, there were 42 patients (37 males and 5 females) analyzed in the study, with a median age of 67 years (40 to 87 years). The proportion of patients with stage T1 was 83% (35/42). Other than the TMB value, there were no significant differences between the two groups terms of in sex, age, smoking history, drinking history, clinical stage or anterior involvement. The clinicopathological features of the patients are shown in Table 1 .

We performed targeted NGS of 523 cancer-relevant genes ( sTable 4 ) using treatment-naïve tumor biopsies from 42 laryngeal cancer patients. Figure 2A shows that the number of mutations per patient ranged from 2 to 31, with a median of 7.0. The mean number of mutations in the relapsed group was higher than that in the non-relapsed group (9.0 vs. 6.0, p=0.045), indicating that more somatic mutations were present in the former group. Figure 2B depicts the genetic changes in the entire cohort, a total of 469 genomic alterations were detected in 211 distinct cancer-relevant genes, and 51 mutated genes (>5% patients) within the 523-gene panel are shown. The genes found to be mutated in more than five patients (>10%) included TP53 (78.5%), FAT1 (26%), LRP1B (19%), CDKN2A (17%), TET2 (17%), NOTCH1 (12%) and NRG1 (12%). All samples presented microsatellite stable status.

Among these genes, TP53, a gene related to DNA repair, displayed the highest mutation frequency in the cohort, and the most common mutation type was mis-sense mutation ( Supplementary Figure 1 ). Patients in the relapsed group tended to have a higher TP53 mutation rate (87% vs. 68%). The second most commonly mutated gene was FAT1, a tumor suppressor preventing cancer development, FAT1 promotes the migration of HCC cells by regulating the expression levels of EMT-associated genes (27). The major mutation type in our cohort was frame-shift mutation ( Supplementary Figure 1 ). Interestingly, pairwise mutual exclusivity and cooccurrence analysis revealed that FAT1/CDKN2A, and FAT1/SMAD family member 4 (SMAD4) were significantly co-occurring (p<0.05), and FAT1 typically displayed exclusivity with LRB1P (p<0.1) ( Figure 3A ). LRP1B, which encoding endocytic LDL-family receptor, is among the top 10 significantly mutated genes in human cancer (28). LRP1B mutation is associated with a higher TMB and better survival outcomes in patients with melanoma and NSCLC (29). In our study, the only type of mutation found in LRP1B was missense mutation ( Supplementary Figure 1 ), and LRB1P mutation tended to co-occur with GNAS complex locus (GNAS), kelch like ECH associated protein 1 (KEAP1), lysine methyltransferase 2D (KMT2D), RAN binding protein 2 (RANBP2) and SRY-box transcription factor 17 (SOX17) mutation (p<0.1) ( Figure 3A ). In addition, cell cycle related gene CDKN2A, tumor suppressor genes like TET2, NOTCH1 and NRG1 were all found to have a high mutation percentage in the cohort. Of note, all patients with NOTCH1 mutations were in the relapsed group ( sFigure 1 ). Co-occurring mutations in SETD2/NRG1, SETD2/ ankyrin repeat domain containing 11 (ANKRD11), SOX17/ kelch like ECH associated protein 1 (KEAP1) and spectrin alpha, erythrocytic 1 (SPTA1)/spen family transcriptional repressor (SPEN) were also found (p<0.05) ( Figure 3A ). Overall, a total of five major mutation types were detected: 273 mis-sense mutations (58.2%), 110 copy number variations (23.4%), 43 frame-shift mutations (9.2%), 35 nonsense mutations (7.5%) and 8 other mutation types (1.7%). Separate group analyses showed no significant difference in the proportion of mutation types between the two defined groups ( Figure 3B ).

TMB is thought to be related to the number of neoantigens in tumors and plays an important role in predicting the efficacy of immunotherapy (8, 9). TMB-H has been shown to be a poor prognostic factor for recurrence in patients with resected NSCLC. Therefore, we investigated whether TMB-H is associated with decreased RFS in the laryngeal cancer cohort. Interestingly, the TMB values in the relapsed group were higher than those in the non-relapsed group (7.03 vs. 4.69, p=0.045, Figure 4A ), suggesting that there were more somatic gene mutations in the relapsed group. In line with the findings reported for resected NSCLC, a higher TMB (top 25%) was associated with poorer RFS outcomes in the resected laryngeal cancer cohort (median RFS: not reached vs. 11.5 months in TMB-L vs. TMB-H; Hazard ratio (HR): 2.4, 95% CI: 1.05-5.48, p=0.036; Figure 4B ). Next, we further examined the clinical impact of the presence of multiple somatic mutations (>1 mutation) and only one somatic mutation in this cohort. In line with the findings for TMB-H, patients with multiple somatic gene mutations (>1 mutation) (n=37) showed worse RFS than those (n=5) with only one somatic gene mutation (p=0.033) ( sFigure 1B ). These results indicate that TMB-H, which is characterized by a high number of somatic gene mutations, may be related to the degree of malignancy of laryngeal cancer, leading early recurrence in laryngeal cancer.

We further explored the rate of individual somatic mutated genes between the two groups. Among all the candidate genes (mutated in more than five cases, >10% patients), LRB1P mutation (30.4% vs. 5.3%, p=0.039, Figure 4C ) and NOTCH1 mutation (21.7% vs. 0%, p=0.030, Figure 4C ) were occurred more frequently in the relapsed group, with significant differences observed. Specifically, all NOTCH1 mutations occurred in the relapsed group. However, Kaplan–Meier survival curves showed that LRB1P mutation was not associated with RFS (median RFS: not reached vs. 25.0 months in LRB1P-WT vs. LRB1P-Mu; HR: 1.53, 95% CI: 0.62-3.75, p=0.335; Figure 4D ) or overall survival (OS) (median OS: not reached vs. 58.0 months in LRB1P-WT vs. LRB1P-Mu; HR: 1.14, 95% CI: 0.23-5.65, p=0.877; Figure 4E ) in our laryngeal cancer cohort. Of note, patients with NOTCH1 mutation presented significantly shorter RFS (median RFS: 40.0 vs. 8.0 months in NOTCH1-WT vs. NOTCH1-Mu; HR: 4.73, 95% CI: 1.68-13.27, p=0.001; Figure 4F ) and OS (median OS: not reached vs. 25 months in NOTCH1-WT vs. NOTCH1-Mu; HR: 5.37, 95% CI: 0.97-29.61, p=0.031; Figure 4G ). In addition, TMB values in the NOTCH1-mutant group were higher than those in the wild-type group, although this difference was not significant, possibly due to the small number of samples (9.5 vs. 5.2, Mu vs. WT, p=0.098). Consistent with the findings in our cohort, the TCGA laryngeal cancer dataset confirmed that patients with NOTCH1 mutation across different tumor stages, mainly advanced stages, had significantly shorter progression-free interval (PFI) (median PFI: 68.7 months vs. 14.5 months in NOTCH1-WT vs. NOTCH1-Mu; HR: 2.08, 95% CI: 0.99-4.37, p=0.049; Figure 4I ), disease-specific survival (DSS) (median DSS: 213.9 months vs. 32.7 months in NOTCH1-WT vs. NOTCH1-Mu; HR: 3.08, 95% CI: 1.34-7.04, p=0.005; Figure 4J ) and OS (median OS: 70.7 months vs. 15.7 months in NOTCH1-WT vs. NOTCH1-Mu; HR: 2.61, 95% CI: 1.35-5.05, p=0.003; Figure 4K ). These results suggest that NOTCH1 mutation is a poor prognostic factor for laryngeal cancer and that patients with resected laryngeal cancer whose tumors carried NOTCH1 mutation have a higher risk of relapse.

Clinical factors and genetic alterations are potential predictors of prognosis; therefore, we used univariate and multivariate Cox regression models to examine the correlation between genetic changes and patients’ RFS and OS. Baseline clinical characteristics including sex, age, smoking history, drinking history, clinical stage and anterior commissure involvement status were also examined. The results of univariate analysis showed that NOTCH1 mutation (HR: 4.44, 95% CI: 1.59-12.4, p=0.004) and TMB-H (HR: 2.40, 95% CI: 1.05-5.48, p=0.038) were adverse prognostic factors for RFS ( Table 1 ). To determine whether NOTCH1 mutation and TMB-H were independent predictors of survival outcomes, variables with p<0.05 in the univariate analysis were then included in multivariate analysis. The results of multivariate analysis showed that both NOTCH1 mutation (HR: 4.18, 95% CI: 1.20-14.6, p=0.025) and TMB-H (HR: 2.85, 95% CI: 1.16-7.01, p=0.023) were independent genetic factors found to be significantly associated with shorter RFS ( Tables 1 , 2 ). In addition, univariate and multivariate analyses revealed that only NOTCH1 mutation was identified as an independent prognostic factor for OS (univariate, HR: 5.37, 95% CI: 0.97-23.6, p=0.054; multivariate, HR: 10.2, 95% CI: 1.49-69.7, p=0.018) ( sTable 5 ).

Previous studies suggested that a weakened immune phenotype was present in relapsed stage IA NSCLC tumors (30, 31). Thus, to elucidate the role of NOTCH1 mutation on recurrence in resected laryngeal cancers, we further investigated the variation in the anti-tumor immunity of the TiME in relapsed tumors with NOTCH1 mutation and NOTCH1 wild-type tumors by Nano-String gene expression assay. After excluding specimens exhibiting RNA degradation or unqualified quality control, tumors from 4 NOTCH1 mutation patients with relapse and 17 NOTCH1 wild-type relapsed patients were included in the analysis. We first explored the DEGs between the two groups. In total, seven upregulated DEGs were identified in NOTCH1 mutation group, which are shown as a heatmap and volcano plot ( Figures 5A, B ). In DEGs, most of these genes were related to specific growth-stimulating molecules, such as MYC (MYC proto-oncogene), ADM (adrenomedullin), TGFβ1 (transforming growth factor beta 1) and TAP1 (transporter associated with antigen processing 1). Besides, the expression of CD44, a cell-surface glycoprotein gene involved in cell-cell interactions, cell adhesion and migration, was higher in patients with NOTCH1 mutation. The expression of TNFSF10 (TNF superfamily member 10) and PSMB9 (proteasome 20S subunit beta 9) were also higher in patients with NOTCH1 mutation. In particular, these upregulated genes were determined to have the potential to be small-molecule targets for precision therapy in patients who had recurrent laryngeal cancer with NOTCH1 mutations.

Furthermore, we explored the function of DEGs by GO and KEGG analysis and found that the enrichment of DEGs by GO was mainly concentrated in receptor ligand activity, signaling receptor activator activity pathways in molecular function ( Supplementary Figure 2A ). KEGG analysis showed that DEGs enrichment was concentrated in Natural killer cell-mediated cytotoxicity and cell adhension molecules pathways ( Supplementary Figure 2B ). In addition, the GSEA results demonstrated that cell cycle pathway were enriched in the NOTCH1 mutation group (NES: 1.865, p.adjust=0.004), indicating that NOTCH1 mutation might be positively correlated with the cell cycle process leading to carcinogenesis, which is associated with poorer survival outcomes ( Supplementary Figure 2C ).

Next, we measured the intratumoral abundance of various immune cell populations with differential gene expression based on the nCounter immune profile panel. At the metagene level, the 14 immune cells composition were assessed, and the scores of some immune cell subsets, such as B cells (p=0.031), CD45⁺ cells (p=0.035), dendritic cells (p=0.024), macrophages (p=0.006), mast cells (p=0.009), neutrophils (p=0.031) and T cells (p=0.018), were significantly decreased in relapsed samples with NOTCH1-mutant compared to NOTCH1 wild-type ( Figures 5C, D ), indicating a weakened immune phenotype in NOTCH1-mutant relapsed patients.

Several biomarkers have been developed to evaluate different aspects of the adaptive immune response, such as the downstream consequences of immune activation, as measured by the presence of immune cells, and the immunoscore or gene expression signature related to the immune environment (32). We next evaluated 10 predefined immune gene signatures included in the immuno-oncology panel to assess the biological differences in the tumor immune environment. The results showed that the expression of related gene signatures, such as the T-cells marker score (CD2, CD3D, CD3E, HLA-E, IL2RG, NKG7) (WT vs. Mu: 8.6 vs. 6.7, p=0.018, Figures 6A, B ), B-cells score (CXCL13, CXCR6, IL18, IL2RG, LCK, PSMB10, TNFRSF4, TNFRSF14) (WT vs. Mu: 7.3 vs. 6.4, p=0.024, Figures 6C, D ), and immune signature score (CD2, CD247, CD3E, GZMH, GZMK, NKG7, PRF1) (WT vs. Mu: 7.0 vs. 4.1, p=0.024, Figures 6E, F ), were significantly reduced in recurrent patients with NOTCH1 mutation. In addition, the TILs score, which incorporates different immune cell populations (B cells, CD8 T cells, cytotoxic cells, exhausted CD8 cells, macrophages, NK CD56dim cells, total T cells), was further analyzed and found to be significantly decreased in the NOTCH1-mutant group (WT vs. Mu: 6.6 vs. 4.5, p=0.006, Figures 6G, H ). However, no significant differences were found in terms of the cytotoxic T lymphocyte activity (CTL) score (WT vs. Mu: 7.5 vs. 6.43, p=0.2), cytolytic activity (CYT) score (WT vs. Mu: 7.9 vs. 6.8, p=0.2), chemokine score (WT vs. Mu: 7.8 vs. 6.8, p=0.14), angiogenesis score (WT vs. Mu: 11.6 vs. 11.9, p=0.57), IFN-γ signature score (WT vs. Mu: 2.6 vs. 2.4, p=0.76), or T-cell-inflammatory gene expression profile (GEP) score (WT vs. Mu: -0.4 vs. -0.5, p=0.46) ( Supplementary Figures 2D–I ). These results highlighted that the attenuation of these immunophenotypes may lead to immunosurveillance escape and postoperative recurrence in patients with NOTCH1 mutation.

To further confirm the alterations observed in the various immune cell populations, we used the TIMER algorithm to estimate the TCGA immune cell infiltration status in laryngeal cancer datasets. In line with the findings of the NanoString gene expression assay, the proportion of activated B cells in NOTCH1-mutant tumors was lower than that in wild-type tumors ( Figures 7A ). To further explore the potential mechanism behind NOTCH1 mutation and the adaptive immune response in laryngeal cancer, immune-related signatures and signaling genes were analyzed based on RNA data from the TCGA database. The GSEA results revealed prominent enrichment in signatures related to the activation of the immune response pathway (NES: -1.57, adjusted p<0.001), adaptive immune response pathway (NES: -1.49, adjusted p<0.001), B-cell activation pathway (NES: -1.570, adjusted p=0.009) and immune response-activating cell surface receptor signaling pathway (NES: -1.58, adjusted p<0.001), all of which were downregulated in the NOTCH1-mutant group ( Figure 7B ), further indicating a weakened adaptive immune response. Moreover, cell-cell adhesion via the plasma membrane adhesion pathway (NES: -1.91, adjusted p<0.001) was also downregulated in the NOTCH1-mutant group ( Figure 7B ), highlighting the impaired cell-cell adhesion functions in these patients. These results suggest that preexisting tumor immunity characteristics are a determinant of recurrence, and the weakened adaptive immune response and impaired adhesion functions in the tumor microenvironment may contribute to recurrence.