The workflow of this study is illustrated in Figure 1. Gene expression RNA-Seq data, including count and TPM data, were obtained from the University of California, Santa Cruz (UCSC) browser Xena (http://xena.ucsc.edu/) and were associated with TCGA-HNSC clinical information. Gene expression microarray datasets (log₂-transformed Robust Multi-Array normalized data) and their associated prognostic datasets, GSE27020 (Fountzilas et al., 2013) and GSE42743 (Lohavanichbutr et al., 2013), were selected and downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). These GEO datasets were annotated with gene symbols using information from the GPL96 Affymetrix Human Genome U133 A Array and GPL570 Affymetrix Human Genome U133 Plus 2.0 Array platform files, respectively. All the data were processed using R (version 4.0.4) and RStudio (version 1.2.5033). TCGA-HNSC contained 502 HNSCC tumor 44 normal tissues, GSE27020 contained 109 LSCC samples with survival data, and GSE42743 contained 74 OSCC and 29 normal samples with survival data (Supplementary Table S1).

ScRNA-Seq data from the GSE150321 (Song et al., 2020) and GSE172577 (Peng et al., 2021) datasets were obtained from the GEO database. The R packages “Seurat” (version 4.0.2) and “scCancer” (Guo et al., 2021) were used to process the data. The uniform manifold approximation and projection (UMAP) algorithm was used to explore and visualize cluster classifications across the cell samples.

A list of 274 core matrisome genes was obtained from the Matrisome Project (http://matrisomeproject.mit.edu) (Shao et al., 2020). Genes with expression in TCGA-HNSC, GSE27020, and GSE42743 datasets were included for further analysis. A total of 190 common core matrisome genes (Supplementary Table S2), including 33 collagens, 127 ECM glycoproteins, and 30 proteoglycan-coding genes, were selected.

Differential expression analysis was performed using edgeR (Robinson et al., 2010) in the OmicShare tool, a free online platform for data analysis (www.omicshare.com/tools). Protein-coding genes with mean counts of >5 were selected. The default parameters of edgeR were used, and DEGs were selected according to log₂ |fold–change (FC)| ≥ 2.0 and Q value >0.05. As there were no replicates in this study, the biological coefficient of variation (BCV), which is the square root of dispersion, was set to 0.01 following the suggestion from the edgeR official manual.

The Metascape tool (Zhou et al., 2019b) (https://metascape.org) was used to perform Gene Ontology (GO) analysis for biological processes (BP), along with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Statistical significance was set at p < 0.01.

Univariate Cox regression analysis was performed to screen the differential matrisome genes that were significantly associated with overall survival (OS) in the TCGA-HNSC dataset. Genes with p < 0.1 were included for subsequent research.

LASSO and multivariate Cox regression model was widely used for high-dimensional predictor identification. LASSO regression analysis was performed using the “glmnet” R package (version 4.1-2). For multivariate Cox regression model analysis, OS-associated matrisome genes were used to select the significant genes associated with the OS of patients with HNSCC according to the coefficient value. These factors were incorporated in the multivariate Cox regression model to construct the HNSCC prognostic signature. The risk score was calculated using the following formula: Risk Score = ∑ i = 1 n β i × E x p i

In this formula, β denotes the regression coefficient, and Exp denotes the expression levels of each core matrisome gene, i (Hou et al., 2020; Huang et al., 2021b). Samples in the TCGA-HNSC cohort were divided into high- or low-risk groups depending on their median risk scores. Receiver operating characteristic (ROC) and Kaplan–Meier analyses were conducted between the high- and low-risk groups. The R packages “pROC” (version 1.18.0), “timeROC” (version 0.4), “survival” (version 3.2-13) and “survminer” (version 0.4.9) were used for visualization.

Single-sample gene set enrichment analysis (ssGSEA) was used to evaluate the relative proportion of 23 infiltrating immune cell types in TCGA-HNSC, including adaptive immune cells (activated B, activated CD4⁺ T, activated CD8⁺ T, gamma delta T, immature B, regulatory T, T follicular helper, type 1 T helper [Th1], Th17, and Th2 cells), and innate immune cells (activated dendritic cells [DCs], CD56bright natural killer cells, CD56dim natural killer cells, eosinophils, immature DCs, macrophages, mast cells, myeloid-derived suppressor cells (MDSCs), monocytes, natural killer cells, natural killer T cells, neutrophils, and plasmacytoid DCs). A bar plot was drawn to evaluate the differences in the composition of these 23 types of tumor-infiltrating immune cells between the high and low-risk groups. The correlation analysis of the relationship between risk score and immune cells was visualized by the “corrplot” R package (version 0.92).

The R package “maftools” (version 2.6.05) was used to calculate the tumor mutation burden score for each sample from the high- and the low-risk groups of the TCGA-HNSC dataset and generate the oncoplot waterfall plot.

KEGG analysis was conducted on the high- and low-risk groups using GSVA. Reference information was downloaded from the Molecular Signature Database v7.4 (MSigDB v7.4, http://software.broadinstitute.org/gsea/msigdb/index.jsp). Statistical significance was set at adjusted-p < 0.05.

The HPA is a database that maps all human proteins in cells, tissues, and organs using an integration of various omics technologies (https://www.proteinatlas.org/). We verified the protein expression levels of survival-related matrisome genes by IHC using the HPA database.

The laryngeal cancer cell line AMC-HN-8 and OSCC cell line JHU011 were acquired from the Department of Otolaryngology-Head and Neck Surgery, Xiangya Hospital, Central South University. The cell lines were maintained in Dulbecco’s Modified Eagle’s Medium (DMEM) with high glucose (Procell Life Science&Technology Co., Ltd., Wuhan, China) and 10% fetal bovine serum (FBS) (Procell Life Science&Technology Co., Ltd.). Cells were maintained at 37°C in a humidified incubator with 5% CO₂.

The LAMB4 siRNA was produced by GenePharma Inc. (Suzhou, China). The LAMB4 siRNA sequence was 5′-GCC​UUC​GAU​GUU​UGC​ACA​ATT-3′, and the control siRNA sequence was 5′-UUC​UCC​GAA​CGU​GUC​ACG​UTT-3′. The siRNA was transfected into AMC-HN-8 or JHU011 using Lipofectamine 3,000 Reagent (Invitrogen, United States) with the Opti-MEM medium (Gibco, Waltham, MA, United States).

The cell sample total RNA was extracted using the TRIzol reagent (Solarbio, Beijing, China) and subjected to reverse transcription with random primers using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, United States). The expression level of the targeted genes was measured with the Maxima SYBR Green/ROX qPCR Mix (Thermo Fisher Scientific) using a 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 an internal control. The primer sequences of LAMB4 were 5′-AAA​GAG​AAC​GTG​GAA​GGA​GC-3′ (forward) and 5′-TCA​CAG​GTC​AAG​AAT​GGC​AG-3′ (reverse), and the primer sequences of 18 s rRNA were 5′-CCT​GGA​TAC​CGC​AGC​TAG​GA-3′ (forward) and 5′-GCG​GCG​CAA​TAC​GAA​TGC​C-3′ (reverse).

A total of 2,000 cells per well were seeded into 96-well plates. For each well, 10 μL CCK-8 reagent (Solarbio, Beijing, China) was added into the medium. After incubation at 37°C for 2 h, the optical density (OD) was measured at different time points at 450 nm.

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

The statistical analyses were performed using R software (version 4.0.4) or GraphPad Prism (version 8.0). Gene expression levels were compared using Student’s t-test. P < 0.05 was considered statistically significant. *p < 0.05; **p < 0.01; ***p < 0.001.

DEG analysis was conducted by genetic screening of TCGA-HNSC microarray data with a threshold of |log₂(FC)| > 2.0 and an adjusted-p < 0.05. A total of 1,779 DEGs between tumor and normal samples were identified, including 939 upregulated genes and 840 downregulated genes (Figure 2A).

The DEGs were mostly enriched in “muscle structure development,” “supramolecular fiber organization,” “skeletal system development,” “cilium movement,” “pattern specification process,” and “extracellular matrix organization” in GO BP analysis (Figure 2B and Supplementary Figure S1). In the KEGG analysis, DEGs were enriched in “protein digestion and absorption,” “salivary secretion,” “ECM-receptor interaction,” “calcium signaling pathway,” “focal adhesion,” and “hypertrophic cardiomyopathy” (Supplementary Figures S2A,B).

Intersection analysis of the DEGs in TCGA-HNSC and 190 common core matrisome genes identified 52 differentially expressed matrisome-related genes (Figures 2C,D). These genes were used for subsequent analysis.

To establish a prognostic assessment model for patients with HNSCC, we selected the TCGA-HNSC dataset (only tumor samples with survival data, N = 502) as the training cohort and GSE27020 (LSCC samples, N = 109) and GSE42743 (OSCC samples, N = 74) datasets as testing cohorts. Univariate Cox regression analysis of the 52 differentially expressed matrisome genes was performed to evaluate the OS of the samples from the training cohort (Supplementary Table S3). Twelve genes (MATN3, LAMC2, FN1, SPP1, LAMB4, LAMB3, DMP1, IBSP, DMBT1, TGFBI, CHAD, and MMRN1) with p < 0.1 were then included for LASSO analysis. All 12 matrisome genes achieved minimum partial expression (Supplementary Figure S2A,B) and were further incorporated into multivariate Cox regression analysis with their expression levels and prognostic data to identify genes involved in signature construction. Multivariate Cox regression analysis established a prognostic model consisting of a risk signature comprising six genes (FN1, LAMB4, LAMB3, DMP1, CHAD, and MMRN1). The formula for the risk score calculation was as follows: risk score = 0.075123 × FN1 + (−0.348,219) × LAMB4 + 0.090359 × LAMB3 + 0.363,187 × DMP1 + (−0.314,325) × CHAD + (−0.125,206) × MMRN1 (Supplementary Figure S3).

According to the risk score, samples from the training and testing cohorts were divided into two groups, a high-risk and low-risk group, based on a cutoff value of 50% (Figures 3A−O). Kaplan–Meier analysis showed that patients with high-risk scores had significantly shorter survival times than those in the low-risk group, both in the training (TCGA-HNSC) and testing cohorts (GSE27020 and GSE42743) (Figures 3A,F,K). In addition, the area under the curve (AUC) values of the ROC curves for the 1-, 2-, and 5-years survival rates were 0.636, 0.635, and 0.571, respectively, indicating that the risk score can be used to predict prognosis in the training cohort (TCGA-HNSC) (Figure 3B). In the testing cohorts, the AUC values for the 1-, 2-, and 5-years survival rates were 0.724, 0.691, and 0.650 in GSE27020 (Figure 3G) and 0.707, 0.602 and 0.859 in GSE42743 (Figure 3L), respectively. Additionally, all six genes (FN1, LAMB4, LAMB3, DMP1, CHAD, and MMRN1) were significantly associated with poor prognosis and unhealthy living habits in both the training (TCGA-HNSC) and the testing cohorts (GSE27020 and GSE42743) (Figures 3E,J,O). Therefore, a matrisome-associated prognostic model for HNSCC was established and verified.

To evaluate the landscape of the 23 infiltrating immune cell types in TCGA-HNSC, ssGSEA was used. Based on the clinical data of TCGA-HNSC and the risk score of each sample, a heatmap of the 23 infiltrating immune cell types were drawn (Figure 4A). Next, correlation analysis between tumor-infiltrating immune cells was performed. As shown in Figure 4B, the highest significantly positive correlation was between activated CD8 T cells and neutrophils, whereas the highest significantly negative correlation was between CD56bright natural killer cells and monocytes. The comparison analysis of tumor-infiltrating immune cell abundance between the different risk groups showed that the presence of activated and immature B cells, eosinophils, macrophages, monocytes, neutrophils, and Th2 cells were significantly reduced in the high-risk group (Figure 4C). Finally, we analyzed the correlation between 23 infiltrating immune cell types and the risk scores. The results showed that Th17 cells (p = 0.015) were the only significantly positively correlated cell type with the risk score, whereas the risk score had a significantly negative correlation with immature B cells, monocytes, macrophages, eosinophils, activated B cells, Th2 cells, gamma delta T cells, MDSCs, and neutrophils (Figure 4D). Collectively, our results show that the 7 cell types (activated B cells, eosinophils, immature B cells, macrophages, monocytes, neutrophils, and Th2 cells) may play an important role in the matrisome-related HNSCC microenvironment.

Furthermore, the somatic mutation profiles of TCGA-HNSC patients were used to explore common somatic mutations in the high- and low-risk groups. Among these patients, 239 (97.15%) belonged to the high-risk group, and 215 (87.40%) belonged to the low-risk group. The frequency of gene mutations was generally higher in the high-risk group compared to the low-risk group. Alterations in the mutation landscape in high- or low-risk group were as follows: eight genes were mutated in >15% of tissues with high-risk score: TP53 (72%), TTN (35%), FAT1 (26%), CDKN2A (19%), MUC16 (17%), CSMD3 (16%), NOTCH1 (16%), and LRP1B (17%), while eight genes were mutated in >15% of tissues with low-risk score: TP53 (52%), TTN (35%), FAT1 (16%), CDKN2A (17%), MUC16 (16%), PIK3CA (18%), CSMD3 (16%), and SYNE1 (16%). Notably, TP53 was one of the most commonly mutated genes in cancer, occurring more frequently in the high- (72%) than in the low-risk group (52%) (Figures 5A,B).

GSVA was performed to identify prognosis-related KEGG pathways. As shown in the heatmap, the high-risk group was significantly enriched in the P53 signaling pathway, thyroid cancer, colorectal cancer, endometrial cancer, prostate cancer, basal cell carcinoma, and melanoma KEGG pathways, while glycosaminoglycan biosynthesis keratan sulfate, glycosylphosphatidylinositol GPI anchor biosynthesis, and folate biosynthesis KEGG pathways were significantly enriched in the low-risk group (Figure 5C). Therefore, gene expression in the high-risk group was significantly enriched in tumor-related pathways.

Using HPA online datasets, we then verified the protein expression of these prognostic genes through IHC. The FN1, LAMB3, and DMP1 proteins were found to be upregulated, while the LAMB4, CHAD, and MMRN1 proteins were found to be downregulated in HNSCC samples compared with the normal controls (Figure 6). These findings were consistent with our results obtained using the TCGA dataset.

Next, we used single-cell RNA-Seq data from the GSE150321 and GSE172577 datasets to further verify the relationship between the prognostic model and genes in HNSCC. For the GSE150321 dataset, which comprised data from two LSCC samples, a total of 5 cell clusters (tumor, immune, epithelial, mesenchymal, and endothelial cells) were identified based on previous literature (Song et al., 2020) (Figure 7A). We then calculated the risk score for each cell and plotted it in a UMAP plot and violin plots (Figure 7B and Supplementary Figure S4). LSCC tumor cells had higher risk scores compared to the non-tumor cells.

For the GSE172577 dataset, four OSCC samples, GSM5258385, GSM5258386, GSM5258387, and GSM5258388, were combined. The cell clusters were annotated manually based on the “scCancer” R package (Figures 7C,D). Notably, the expression profile of most prognostic genes was similar to that in LSCC, except that in GSE172577 CHAD was expressed at lower levels (Figure 7E and Supplementary Figure S5). Collectively, these results further supported that this signature of matrisome-related prognostic genes influences HNSCC progression.

Based on the previous analysis results, we selected LAMB4 for further analysis since it was one of the six prognostic genes that is less researched in the context of HNSCC. LAMB4 siRNA was used to silence LAMB4 expression in the AMC-HN-8 and JHU011 cell lines (Figure 8A). The CCK-8 assay demonstrated that silenced LAMB4 levels promoted both AMC-HN-8 and JHU011 cell proliferation (Figure 8B). We also detected the migration change after LAMB4 silencing. Transwell assays showed that LAMB4 knockdown promoted migration in both LSCC and OSCC cell lines (Figures 8C,D). Therefore, as a potential tumor suppressor gene, silencing LAMB4 promoted the proliferation and migration of HNSCC cell lines.