The workflow used in this study is shown in Figure 1 . The search words “laryngeal carcinoma” OR “LSCC” OR “laryngeal cancer” AND “human” AND “Expression profiling by array” were applied for dataset retrieval. Datasets that included non-tumor samples and the total sample number of over 100 could be selected. Therefore, GSE142083 was selected to perform the following bioinformatics analysis. The gene expression matrix of 106 LSCC samples (53 LSCC and paired adjacent normal mucosa tissues) was downloaded from the GSE142083 dataset obtained from the GEO database (21). The GEO expression matrix was annotated with gene symbols using information from the GPL20301 Illumina HiSeq 4000 (Homo sapiens) platform file (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL20301) as transcripts per million (TPM), and this information was log₂ transformed in R (version 4.0.4) and RStudio (version 1.2.5033) when necessary. Principal component analysis (PCA) using the “ggfortify” package was performed, and the outliers were excluded (GSM4219698 and GSM4219730) ( Supplementary Figure 1 ). The remaining data from the 52 LSCC and 52 paired normal tissue samples were included in the subsequent analysis. GSE27020 (22), GSE39366 (23), and GSE127165 (24) were also downloaded from GEO for validation.

The gene expression RNA-Seq data and the corresponding clinical data from the TCGA-HNSC dataset were downloaded from the University of California, Santa Cruz (UCSC) Xena (http://xena.ucsc.edu/).

scRNA-Seq data from the GSE150321 dataset were obtained from the GEO database (25). The “Seurat” R package (version 4.0.2) was used to process the data (26). The scRNA-Seq data of GSM4546857 and GSM4546858 from two LSCC patients were included and processed as described previously. Cells with less than 200 genes detected genes were excluded. Data were normalized by the method “LogNormalize” with a scale factor equal to 10,000. After determining the statistically significant principal components, cell clusters were annotated using the information from previous literatures. Next, we implemented the t-distributed stochastic neighbor embedding (t-SNE) algorithm to explore and visualize the cluster classification across cell samples. The trajectory analysis and pseudotime analysis of tumor cells were performed using the “monocle” R package (version 2.18.0) (27).

After removing the outliers, the expression data were standardized using R software. The “limma” (version 3.46.0) software package (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (28) was used to perform the DEG analysis between the LSCC and adjacent normal control tissues from the GSE142083 dataset. Genes with an expression-adjusted p value < 0.05 and log₂ (fold change) > 1.5 were regarded as significant DEGs. The volcano and heatmap plots were generated using “ComplexHeatmap” and “ggplot2” package. Extracellular matrix-associated genes were annotated using the Matrisome Project (http://matrisomeproject.mit.edu) (29).

GO and KEGG enrichment analyses of DEGs and hub modules were performed using the “clusterProfiler” (version 3.18.1) R package (http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (30). The three main processes in GO analysis were as follows: biological process (BP), molecular function (MF), and cellular component (CC). Statistical significance was set at adjusted-p < 0.05. Plots were generated using the “GOplot” (version 1.0.2) package.

The “WGCNA” package (version 1.70-3) in R (14) was used to construct co-expression networks of LSCC and adjacent normal tissue samples. Genes with a mean fragments per kilobase per million mapped reads (FPKM) > 0.5 were selected. The adjacency matrices that stored information of the entire co-expression network were created based on Pearson’s correlation matrices. A topological overlap matrix (TOM) was created from the adjacency matrix to estimate the connectivity of the network. Average linkage hierarchical clustering was used to construct a clustering dendrogram of the TOM matrix with a minimum module size of 30. Finally, similar gene modules were merged with a threshold of 0.25.

The hub genes were identified as intersecting between the midnight-blue, blue, and green-yellow modules from the WGCNA and DEGs in the GSE142083 dataset. These hub genes were validated for their prognostic signature.

First, we identified the LSCC patient cohort through the clinical data from the TCGA-HNSC dataset (TCGA-HNSC-Larynx). Univariate Cox regression analysis of these hub genes was performed to screen the significant genes associated with overall survival (OS) in the training cohort. Genes with p < 0.1 were included in the subsequent analysis. Finally, multivariate Cox regression model analysis was performed to establish a prognostic model. The coefficients of multivariate Cox regression analysis were used to calculate the risk score (RS) of each sample using the following formula:

According to the RS, samples were divided into a high-risk group and a low-risk group with a cutoff value of 50% in the TCGA-HNSC-Larynx cohort. Receiver operating characteristic (ROC) curve analysis and Kaplan–Meier analysis were conducted between the high-risk and low-risk groups. For ROC analysis, “survivalROC” package (version 1.0.3) was performed to draw 5-year OS ROC curves with the survival status (dead or alive) as outcome variable. Verification was also performed using the GSE27020, GSE39366, and GSE127165 datasets.

The R package “maftools” (version 2.6.05) was used to analyze the mutation data of TCGA-HNSC-Larynx dataset. The tumor mutation burden (TMB) score for each sample from the high-risk and the low-risk groups was calculated.

The CIBERSORT (https://cibersort.stanford.edu/) algorithm was used to assess the proportions of 22 types of infiltrating immune cells based on the TCGA-HNSC-Larynx dataset and GSE142083 data following a previously reported procedure (31). A bar plot was drawn to show the differences in the composition of 22 kinds of tumor-infiltrating immune cells between the high-risk group and low-risk group.

KEGG analysis was conducted on the high- and low-risk groups via Gene Set Variation Analysis (GSVA). The reference information was downloaded from the Molecular Signature Database v7.4 (MSigDB v7.4, http://software.broadinstitute.org/gsea/msigdb/index.jsp) (32). Enriched pathways with p-value < 0.05 were considered to be statistically significant.

We conducted single-gene gene set enrichment analysis (GSEA) (33) of the prognostic genes using data from TCGA. All samples were divided into high-expression and low-expression groups based on the median of the prognostic genes, and GSEA was conducted to explore the enrichment of GO biological process (BP) and KEGG pathways in different groups. Statistical significance was set at p < 0.05.

The laryngeal cancer cell line AMC-HN-8 and the dysplastic oral keratinocyte cell line DOK were acquired from the Department of Otolaryngology-Head and Neck Surgery, Xiangya Hospital Central South University. AMC-HN-8 cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM) with high glucose (HyClone, Logan, UT, USA) and 10% fetal bovine serum (FBS) (Biological Industries, CT, USA). DOK cells were cultured in RPMI-1640 medium (HyClone, Logan, Utah, USA) and 10% FBS (Biological Industries, CT, USA). Cells were maintained at 37°C in a humidified incubator with 5% CO₂.

The TEDC2 siRNA (siTEDC2 #1 and siTEDC2 #2) was produced by RiboBio Inc. (Guangzhou, China). SiRNA was transfected into cells using Lipofectamine 3000 Reagent (Invitrogen, USA) with Opti-MEM medium (Gibco, USA).

The total RNA was extracted using TRIzol reagent (Solarbio, Beijing, China) and subjected to reverse transcription with random primers using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, USA). The expression level of targeted genes was measured with Maxima SYBR Green/ROX qPCR Mix (Thermo Scientific, USA) using real-time PCR system (Roche, Basel, Switzerland). The relative RNA expression levels were calculated using the 2^{(−△△CT)} method. The 18 s rRNA was used as internal controls. The sequences of primers will be provided upon request.

For CCK-8 assay, 1,000 cells per well were seeded into 96-well plates. For each well, 10 μl CCK-8 (Solarbio, Beijing, China) was added into the medium. Then, incubated at 37°C for 2 h, the OD value at 450 nm was measured at different time.

For cell migration assay, 2.5 × 10⁵ cells per well that were suspended in DMEM were seeded into a 24-well 8.0-μm transwell top chamber (Jet Biofil, Guangzhou, China). DMEM medium supplemented with 12% FBS was added to the bottom chambers. After incubating at 37°C for 12 h, cells at top chambers were fixed with 4% paraformaldehyde for 30 min, followed by permeabilization with methyl alcohol for 20 min. Then, cells were stained with 0.1% crystal violet (Solarbio, Beijing, China) for 15 min. Cells that did not migrate through the pores were removed by a cotton swab. Cells on the bottom of the chamber were counted using an inverted phase-contrast microscope at low magnifications (×5) (at least three randomly selected fields were quantified).

Statistical analysis was conducted using R (version 4.0.4) or GraphPad Prism (version 8.0). Student-t test was used to compare two normal-distribution groups. NS: not statistically significant; *p < 0.05; **p < 0.01; ***p < 0.001.

We analyzed the DEGs in the GSE142083 dataset using the limma package with a threshold of |log₂(fold-change)| > 1.2 and an adjusted-p < 0.05. A total of 701 DEGs, including 329 upregulated genes and 372 downregulated genes, were screened between the LSCC and adjacent normal samples. The heatmap and volcano plots show the expression patterns of these DEGs ( Figures 2A, B ).

The clusterProfiler package was then used to determine the role of DEGs in the pathogenesis of LSCC using a cutoff criterion of adjusted-p < 0.05. For the biological process group, skin development, epidermis development, cornification, keratinization, and keratinocyte differentiation were significantly enriched. The genes in the cellular component group were significantly enriched in cornified envelope, collagen-containing extracellular matrix (ECM), and ECM component. In the molecular function group, the DEGs were primarily enriched in ECM structural constituent and peptidase regulator activity. KEGG analysis revealed that the IL-17 signaling pathway, salivary secretion, ECM–receptor interaction, cell cycle, and protein digestion and absorption were significantly enriched ( Figure 2C ). We also noticed that ECM-associated genes were enriched in the DEGs, and most of them were highly expressed in the LSCC samples ( Supplementary Figure 2 ). Moreover, hsa04970:Salivary secretion was one of the few terms that had decreasing Z-scores. This may indicate a significant decrease in salivary secretion in patients with LSCC.

To detect the functional module in LSCC, we applied the WGCNA package based on the GSE142083 dataset to establish the gene co-expression networks. To ensure that the network was scale-free, an empirical analysis was conducted to choose an optimal parameter β. Both the scale-free topology model fit index (R²) and mean connectivity reached a steady state when β = 9 ( Figures 3A, B ).

A total of 28 modules were identified via average linkage hierarchical clustering, and each module is represented using a different color in Figure 3C . Among the modules, the modules midnight-blue, blue, and green-yellow exhibited a strong correlation with cancer traits (cor ≥ 0.8, p-value < 0.001) ( Figures 3D – G ). Therefore, these modules were selected as clinically significant modules for further analysis. A set of 400 selected genes was identified for the network heatmap ( Supplementary Figure 3 ).

GO and KEGG analyses were used to assess the biological processes, molecular functions, cellular components, and KEGG pathways of genes in the modules. The midnight-blue module was primarily related to glycosylation events and glycosphingolipid biosynthesis, whereas the blue module was primarily associated with the DNA replication and the cell cycle pathway. The green-yellow module was associated with RNA splicing ( Supplementary Figure 4 ).

Hub genes were screened using a |module membership (MM) score| > 0.8 and |gene significance (GS) score| > 0.2 as the cutoff criteria (34–36). Based on this criterion, we subsequently sorted the genes according to their connectivity to select candidate genes ( Figures 3E – G ). A total of 41 genes in the midnight-blue module, 120 genes in the blue module, and 99 genes in the green-yellow module were screened out.

To identify the “real” key genes, we compared these candidate genes obtained from WGCNA with the 701 DEGs. A total of 85 key genes from the midnight-blue module, blue module, and DEGs were identified ( Figure 3H ).

To validate and explore the prognostic values of these “real” key genes in LSCC, we used the TCGA-HNSC-Larynx dataset (N = 111) as training cohort. Univariate Cox regression analysis of the 85 key genes was conducted with regard to the OS of samples from the training cohort ( Supplementary Table 1 ). Six genes with p < 0.1 were then included for multivariate Cox regression model analysis to establish a prognostic model. Finally, we identified three genes—SLC35C1, HOXB7, and TEDC2—as potential independent prognostic markers for LSCC ( Supplementary Figure 5 ). The risk score for each patient was calculated using the following formula: risk score = 0.703437*SLC35C1 + 0.833199*HOXB7 + (-0.891338)*TEDC2.

According to the level of risk score, samples from the two cohorts were divided into a high-risk group and a low-risk group based on a cutoff value of 50% in both training and testing cohorts, as shown in Figures 4A – G . The risk score was significantly associated with the survival time of patients with LSCC in the training cohort (TCGA-HNSC-Larynx) ( Figure 4A ) and testing cohort (GSE27020) ( Figure 4F ). The ROC curve was then used to evaluate the accuracy of the survival analysis. In the training and testing cohorts, the area under the curve (AUC) was 0.726 and 0.645, respectively, indicating that the prediction effect was good ( Figures 4B, G ). Besides, patients with continuous risk scores harbored various clinical outcomes in different groups ( Figures 4C, D, H , I ). Additionally, all three genes (SLC35C1, HOXB7, and TEDC2) were significantly associated with poor prognosis and unhealthy living habits, in either the training cohort (TCGA-HNSC-Larynx), the testing cohort (GSE27020), or other LSCC datasets (GSE39366 and GSE127165) ( Figures 4E, J and Supplementary Figure 6 ).

Analysis of differences in gene mutations between high-risk and low-risk groups showed that the proportion of patients with gene mutations was 92.73% (51 in 55) and 96.30% (52 in 54) in the high-risk and low-risk score groups, respectively ( Figures 5A, B ). The mutation frequency of TP53 was observably higher in the high-risk group (82%) than in the low-risk group (69%) ( Figures 5A, B ).

The distribution of tumor-infiltrating immune cells is an important indicator of a patient’s lymph node status and prognosis. In order to explore the relationship of the prognostic model and the tumor immune microenvironment, we analyzed the tumor-infiltrating immune cell abundance in the TCGA-HNSC-Larynx dataset using the CIBERSORT algorithm ( Figure 5C ). The results showed that the infiltration levels of the following cells: “T cells CD4 memory activated”, and “T cells regulatory” were significantly higher in the low-risk group than in the high-risk group. In contrast, the infiltration levels of “Macrophages M2” and “Mast cells resting” were significantly lower in the low-risk group. A similar abundance was also estimated in the GSE27020 dataset ( Supplementary Figure 7 ).

GSVA and GSEA were conducted to search for GO and KEGG pathways in which the prognostic genes or risk scores were enriched in samples with high-risk levels from the TCGA-HNSC-Larynx dataset. Metabolism and cancer-associated pathways including “Glyoxylate and dicarboxylate metabolism”, “Protein export”, “Sphingolipid metabolism”, “Lysosome”, and “Bladder cancer” were significantly enriched with low expression of most prognostic genes and high-risk score ( Figure 6A ). Besides, GSEA based on DEGs between the high- and low-risk score groups also revealed that “Apical part of cell”, “Anion transmembrane transporter activity”, and “Tissue homeostasis” were significantly enriched with GO BP ( Figure 6B ) and “Metabolism of xenobiotics by cytochrome p450”, “Cardiac muscle contraction”, and “Systemic lupus erythematosus” were significantly enriched with KEGG analysis ( Figure 6C ).

To further verify the relationship between the prognostic genes and risk score in LSCC, we employed single-cell RNA-Seq data from the GSE150321 dataset (25), which comprised data from two LSCC samples. We identified five cell groups, namely tumor cells, immune cells, epithelial cells, mesenchymal cells and endothelial cells ( Figure 7A ). We then calculated the risk score for each cell and plotted them in a t-SNE plot and the violin plots ( Figure 7B and Supplementary Figure 8 ). Tumor cells and some immune cells showed higher risk scores than other cell types.

Next, we profiled the tumor cells and arranged them into five clusters ( Figure 7C and Supplementary Figures 9A, B ). Pseudotime and trajectory analyses of the tumor cells suggested a continuous cell fate, starting at tumor 3 and progressing toward tumor 1 and tumor 4, with tumor 2 and tumor 5 being a transitioning state ( Figures 7D, E ). Diffusion mapping placed tumor 4 and tumor 1 within populations of high-risk score cells, indicating a differentiation trajectory from low-risk tumor cells to high-risk tumor cells ( Figure 7F ). Together with previous results, these results further supported that the identified prognostic genes influence the LSCC progression.

Based on the results of the previous analysis, we selected TEDC2, one of the three prognostic genes with less research involving LSCC, for further analysis. We first verified the expression level of TEDC2 in different datasets, and they both indicated that TEDC2 was higher in LSCC than in paired normal tissues ( Figure 8A ). Similar results were observed in the LSCC cell line and a precancerous cell line ( Figure 8A ). We therefore used siRNA to inhibit the expression of TEDC2. Results from RT-PCR indicated that both of the two TEDC2 siRNAs in AMC-HN-8 showed a good silenced efficiency ( Figure 8B ).

CCK-8 assays showed that silencing the expression of TEDC2 promoted AMC-HN-8 cell proliferation ( Figure 8C ). Cell-cycle genes, including PCNA, CNNB1, and CNND1, were both downregulated when TEDC2 was knocked down ( Figure 8D ). Results of the transwell assays showed that knockdown of TEDC2 promoted migration ( Figure 8E ). Based on these validations, we concluded that silenced TEDC2, to a certain extent, promoted LSCC cell proliferation and migration.