Single-Cell Transcriptome Analysis of Stepwise Progression of Laryngeal Squamous Cell Carcinoma
- 1. The Affiliated Lihuili Hospital, Ningbo University, China
- 2. School of Medical, Ningbo University, China
- 3. Institute of Drug Discovery Technology, Ningbo University, China
Abstract
Laryngeal squamous cell carcinoma (LSCC) evolves from normal tissues through precancerous conditions, ultimately culminating in primary laryngeal squamous cell carcinoma. To elucidate the heterogeneity of tumor cells and their interactions during this progression, we performed RNASeq analysis on normal and primary tumor samples. This technique allowed us to distinguish the diverse expression profiles of malignant epithelial cells, highlighting shared patterns and abnormal activation pathways during tumorigenesis. We discovered in situ cancer cells in lesions that were not previously detected by pathological examinations, revealing DNA copy number variations of ALG1L and the associated transcriptional regulation that promotes tumor progression. Furthermore, we uncovered the transcriptional dynamics in T cell exhaustion, identifying key transcription factors related to T cell depletion and highlighting new therapeutic targets to enhance anti-tumor immunity. Our analysis of cellular communication has illuminated the complex interactions between immune and cancer cells across the cancer continuum, providing valuable insights into the biological underpinnings of laryngeal cancer and suggesting that various cell lineages contribute to potential carcinogenic pathways. In essence, our study offers insights into the pathobiology of intercellular interactions during the development of laryngeal squamous cell carcinoma, proposing potential targets for the diagnosis and treatment of LSCC.
Keywords
• Laryngeal carcinoma
• Cancer
• Single-cell transcriptomics
• RNA sequencing
• Tumor heterogeneity
CITATION
Yidian C, Chen L, Yanguo L, Shanshan G, Zhisen S (2024) Single-Cell Transcriptome Analysis of Stepwise Progression of Laryngeal Squamous Cell Carcinoma. J Cancer Biol Res 11(1): 1145.
INTRODUCTION
Laryngeal squamous cell carcinoma (LSCC) ranks as the 11th most common malignant tumor worldwide and the second most common in the head and neck region [1]. Accounting for over 90% of laryngeal cancer cases are of the LSCC type. The five-year survival rate for laryngeal cancer in China stands at approximately 77% [2]. Hence, early diagnosis, treatment, disease monitoring, and prognosis evaluation are imperative for improving survival rates. The development of LSCC involves a multistep process originating from laryngeal epithelial precursor lesions, such as chronic hypertrophic laryngitis, laryngeal leukoplakia, and vocal cord leukoplakia, which exhibit similar clinical symptoms, pathological features, and physical signs [3]. Specifically, vocal cord leukoplakia emerges as the most common precancerous lesion of laryngeal cancer, with an incidence ranging from 6% to 7%. The early detection and prevention of malignant transformation in precancerous lesions are effective strategies against LSCC (Molecular Biomarkers of Malignant Transformation in Head and Neck Dysplasia PubMed, n.d.). However, the molecular mechanisms underlying the transformation from precancerous leukoplakia to malignancy remain unknown.
Single-cell technology surpasses the limitations of traditional sequencing methods by providing accurate and detailed cell state information at the individual cell level [4]. This technology helps us understand the immune microenvironment composition in LSCC comprehensively and guides individualized immunotherapy [5]. Additionally, this approach offers distinctive opportunities to evaluate the regulation, evolution, and interaction of individual cells. [6,7]. There have been a study that utilized single-cell RNA sequencing to construct a gastric cell atlas, identifying cellular and molecular characteristics of premalignant and early- malignant lesions, including the acquisition of an intestinal-like stem cell phenotype during metaplasia, marker identification for unique endocrine cells and pre-goblet cell clusters, and the discovery of EGC-specific signatures for precise diagnosis [8]. These studies indicate that cells in precancerous lesions exhibit distinct transcriptome characteristics compared to normal cells. Identifying changes in the tumor microenvironment and single- cell transcriptome early on is a scientific challenge that needs to be addressed.
Although scRNA seq has gained traction in the study of head and neck tumors, its application in LSCC had been mainly focused on revealing heterogeneity among tumor cells and the complexity of the TME [9-12]. To investigate the molecular evolution mechanism of LSCC and precancerous lesions and address gaps in previous studies, we conducted single-cell transcriptome sequencing using three laryngeal squamous cell carcinoma precancerous lesions, three laryngeal squamous carcinoma tumor tis- sues, and two vocal cord polyps’ tissues obtained from clinical LSCC patients. Cluster analysis of 39,733 cells yielded 24 cell subsets, which were then annotated for their functions. Our comprehensive analysis included characterizing the functional attributes of each cell subset, delineating cell development and differentiation trajectories, examining cellular interactions, investigating gene regulatory networks, and identifying key genes associated with cancer progression. This methodology enabled us to provide a meticulous description of cell types in LSCC and precancerous tissues, showcasing the gradual changes in cell composition during LSCC progression. We pinpointed epithelial cells, stromal cells, and immune cells, unraveling their interactions within the tumor microenvironment and offering crucial mechanistic and clinical insights into LSCC.
In conclusion, our findings have elucidated the significant molecular alterations that occur during the progressive development of LSCC and have mapped the transcriptome changes from normal cells to tumor cells for the inaugural time.
MATERIALS AND METHODS
Specimen collection and processing
The ethical consent form for this study received approval from the Ethics Committee of Ningbo Lihuili Hospital and written informed consent from the participants. The samples originated from patients at Ningbo Lihuili with pathologically verified laryngeal cancer, vocal cord leukoplakia, and vocal cord polyps. The procedure for gathering and processing tumor samples is as follows: once the fresh samples are cleaned, part of them is immersed in liquid nitrogen and moved to our hospital’s Human Resource Specimen Bank for cryopreservation at −80?C; another part of the tissue undergoes fixation for 24 hours prior to paraffin embedding, and subsequently, a single-cell nucleus suspension is prepared.
scRNA-seq library preparation and sequencing
The library for single-cell nuclear sequencing was prepared using the Chromium Single Cell 3’ v3 (10x Genomics) for a tenfold Chromium Single Cell 3’ library construction. The library for single cell nuclear sequencing was prepared using the Chromium Single Cell 3’ v3 (10x Genomics) for a tenfold Chromium Single Cell 3’ library construction. Following the preparation and quality assessment of the library, bidirectional sequencing was conducted on the Illumina (Nova 6000) platform using a 2 ∗ 150bp format.
scRNA-seq data preprocessing and quality control
The original BCL files are assembled into raw fastQ file data using Illumina’s bcl2fastq converter. After primary quality control and data filtering, high-quality Clean Reads are obtained. The filtering criteria are: 1) removal of reads containing poly A; 2) removal of reads with more than 3Ns; 3) removal of low- quality reads (where bases with a quality value ≤ 5 make up more than 20% of the entire Read). Subsequently, Cell Ranger (version 3.1.0), the official software of 10x Genomics, is used for cell quality control, detection, and comparison with the reference genome (GRCh38) to obtain the original expression matrix.
We then perform downstream analysis. Using the Seurat R package (version 4.3.0) for further quality control, we apply the following filtering criteria: 1) removal of cells detected with < 200 genes; 2) removal of genes expressed in ≤ 3 cells; 3) removal of cells with > 20⌊% mitochondrial gene content. Additionally, we use the R package Doublet Finder to remove potential doublet cells. After quality control, a total of 48,694 high-quality cells were obtained from 8 patient samples Table S1. Batch effects are reduced using the Seurat R package with default settings, by selecting genes for integration, preparing for integration, finding anchor genes, and integrating samples in four steps. The patient characteristics of the study samples are summarized in [Table S2].
Non-linear dimension reduction of unsupervised cell type annotation and subtype (tSNE/UMAP) to identify clusters and cell type annotation
In this study, the integrated data were dimensionally reduced using PCA, and the Elbow plot function of the Seurat package was used for visualizing the reduction of PC variance. The top 30 principal components were selected for downstream analysis. Then, the Find Neighbors and Find Clusters functions (resolution= 0.15) were used for clustering all cells. Uniform Manifold Approximation and Projection (UMAP) was applied to the PCA results for a two-dimensional visual display, resulting in the identification of 24 cell subgroups.
Combining cell-specific marker genes reported in the literature and specific marker genes of different cell types from the Cell Marker database (http://biocc.hrbmu.edu.cn/ CellMarker/), eight cell groups were identified: epithelial cells (KRT14, KRT16, S100A8); myeloid cells (LYZ, APOE); T cells (GNLY, IL32); B cells (MS4A1, CD79A); fibroblasts (COL3A1); endothelial cells (PECAM1); mast cells (GATA2); dendritic cells (JCHAIN, GZMB). Additionally, epithelial cells, immune cells, and stromal cells were extracted for further subgroup analysis to explore their internal consistency. The steps for the subgroup analysis were similar to the aforementioned analysis.
Differential analysis and functional enrichment
We used the Find Markers function with default parameter Settings to identify differentially expressed genes between subgroups. The cluster Profiler package (version 4.10.0) was used to perform gene ontology analysis and precision gene and Genome Encyclopedia signal pathway enrichment analysis for differential genes, so that gene set variation analysis could reveal the subclub specific biological processes.
Analysis of transcription factor expression
To identify the key regulators that determine the fate of each epithelial subpopulation, regulatory networks and regulatory activities were analyzed using SCENIC (v 1.1.0).
Cell transition trajectory and diffusion map analysis
For the analysis of cell trajectory in high-quality epithelial cells (UMI ¿ 5500), we utilized Monocle 2 (version 2.10.1) [13]. We employed the top 100 markers from each cluster to sequence the cells. The process involved reducing dimensionality and constructing trajectories using the default methods and parameters on the chosen genes. To calculate diffusion components, we applied the Run Diffusion function from the Seurat package, adhering to standard parameters. We used the first three dimensions for creating diffusion maps. Each cluster’s mean cell coordinates were designated as the cluster’s centroid. The starting point was determined as the cell furthest from the combined distance of other stages in the NOR stage.
Analysis of interaction between cell types
To uncover interactions among different cell types, we employed Cell chat (version 1.6.1), utilizing its default settings. Due to the smaller number of cells in some epithelial clusters, we downsampled each cluster to a consistent count of 100 cells. This approach enabled us to effectively demonstrate the interactions between epithelial and immune cells, as well as the reciprocal interactions from immune cells towards epithelial cells.
Immunohistochemical
Immunohistochemistry Procedure Overview: Tissue sections undergo dewaxing followed by antigen retrieval and endogenous peroxidase blocking. Subsequently, serum blocking is applied to tissue surfaces, followed by incubation with primary and secondary antibodies for specific staining. Positive signals are stained using DAB, then counterstaining of cell nuclei highlights morphology. Finally, tissue sections are dehydrated and mounted to achieve transparency and fixation. The immunohistochemically stained samples were independently assessed by two pathologists employing a double-blind evaluation method.
Cell culture
The AMC-HN-8 cell line, obtained from the Cell Bank of the Chinese Academy of Sciences, was cultured under specific conditions. For the culture medium, we utilized Dulbecco’s Modified Eagle Medium provided by Gibco Thermo Fisher Scientific, Inc., USA, and supplemented it with 10% fetal calf serum from the same supplier. Additional supplements to the medium included 100U/ml penicillin and 100µg/ml streptomycin, sourced from Solarbio, Beijing, China. Cells were cultured in a humidified controlled incubator at 37?C and 5%CO2.
Gene overexpression
To induce overexpression of ALG1L in AMC cells, we designed and synthesized small interfering RNAs (siRNAs) targeting ALG1L (si ALG1L), along with a negative control RNA (si NC), through GenePharma Co., Ltd. For the transfection process, cells were plated in 6-well plates at a density of 1 ×10∧5 cells/well and incubated at 37?C until they reached about 70% confluence. Subsequent to this, the cells were transfected with either the siRNAs or si NC using Lipofectamine®2000 (Invitrogen; Thermo Fisher Scientific, Inc.), adhering to the guidelines provided by the manufacturer. The antisense sequences for the three siRNAs were as follows:
si ALG1L: 5’- CUCAAACUUUCCUGAUCCUTT -3’.
si NC: 5’-TTCTCCGAACGTGTCACGTTT-3’.
Knockdown efficiency evaluation by RT-qPCR
The relative expression of ALG1L in AMC cells was assessed using RT-qPCR. Total RNA was extracted from the cells using TRIzol reagent (Invitrogen, USA). Subsequently, cDNA was synthesized employing the RevertAid First Strand cDNA Synthesis kit (Promega Corporation, USA), adhering to the manufacturer’s protocol. The RT- qPCR analysis was carried out with the iQ? SYBR Green Supermix (Bio-Rad, USA). The following primers were used for qPCR:
ALG1L-F: 5′-TCGGCGGGAAGGACAATC-3′
ALG1L-R: 5′-TAGACTCCACGCGGGATGAA-3′
GAPDH-F: 5′-AAGGTGAAGGTCGGAGTCAA-3′
GAPDH-R: 5′-AATGAAGGGGTCATTGATGG -3′
Cell proliferation assay
The proliferation of AMC cells was evaluated using the Cell Counting Kit-8 (CCK-8) assay from Dojindo Molecular Technologies, Kumamoto, Japan. Following 24 hours post- transfection with ALG1L mimics, we seeded 4000 cells in a volume of 100µL into each well of a 96-well plate. We established three replicates for each treatment group. The cells were then incubated at 37?C in a 5 ? % CO2 environment for 24 hours, after which 10µL of CCK-8 solution was added to each well. Two hours later, the optical density (OD) at 450 nm for each well was measured using a microplate reader (Spectra Max M5, Molecular Devices, CA, USA). This experiment was conducted a minimum of three times for consistency.
Edu incorporation assay
The 5-ethynyl-2’-deoxyuridine (EdU) proliferation assay is an effective method for assessing cell proliferation. Following the manufacturer’s guidelines, we prepared a 50µM solution of EdU (Beyotime Biotechnology, Ningbo, China). We added 100µL of this solution to each well and incubated at room temperature for 2 hours. Subsequently, cells were fixed with 4% paraformaldehyde for 30 minutes at room temperature. They were then incubated with Apollo staining solution and Hoechst 33342 for an additional 30 minutes, followed by 2-3 washes with PBS, each lasting 5 minutes. Cell proliferation was observed under an Olympus fluorescence microscope in five randomly selected fields. Post-staining, nuclei appeared blue due to Hoechst 33342, and proliferating cells appeared red from EdU Apollo staining, with the blue and red-stained cells overlapping in the images. The total number of cells and the number of proliferating cells were counted using Image J software.
Cell Migration Assay
We conducted a cell migration assay using a 24-well Transwell chamber, without the inclusion of Matrigel. Following 48 hours of transfection with si ALG1L and si NC, the cells were suspended in serum-free medium and seeded into the upper chamber of the Transwell at a density of 1 × 10∧5 cells per well. The lower chamber was filled with culture medium supplemented with 10% FBS. After a 24-hour incubation period at 37?C, cells that had migrated to the lower chamber were fixed with 4 | % paraformaldehyde and stained with 0.5% crystal violet. Any cells remaining on the upper surface of the filter membrane were carefully removed using a cotton swab. We used a microscope to capture images of the cells on the lower surface, and the counting process was performed three times to ensure accuracy.
Statistical analysis
To visualize the outcomes of statistical analyses, we employed R software (version 4.0.3) along with GraphPad Prism 9.0. Additionally, various statistical tests including the Chi-squared test, Student? s t-test, Wilcoxon test, and Kruskal-Wallis test were appropriately used to determine the significance of differences between groups in our data. A p-value of less than 0.05 was regarded as statistically significant, except in cases where a different threshold was specifically mentioned.
RESULTS
The integration of scRNA-seq profiles from multiple patients maps the cellular diversity in tumors, precancerous lesions and normal laryngeal mucosa at high resolution.
We obtained scRNA-seq profiles from single-cell sequencing of eight samples, each from a unique untreated patient. These included Normal mucosa next to vocal cord polyps (NM, n = 2), Laryngeal Papillomatosis (LP, n = 3 ), and Laryngeal Cancer ( LC, n = 3 ) [Figure 1A].
Figure 1: Establishing of the Cell Atlas for LSCC Based on Single-Cell RNA-Seq. (A) An overview of the study design and workflow. (B) Uniform Manifold Approximation and Projection (UMAP) plot of the analyzed single cells. Each color represents one cluster. Annotated cell types are listed below. (C) The distribution of cells derived from different sample origins. (D) Pie chart of the proportion of cells. (E) Scale plots show the proportion of cells from different sample sources. (F) The heat map shows differentially expressed genes for eight cell types in LSCC. The color keys from purple to yellow indicate low to high expression levels. Typical cell marker gene expression that defines each cluster is shown at the right.
Post mass filtration, a total of 48,694 cells were isolated, comprising 13,438 (27.61%) from LC, 22,522 (46.3)%) from LP, and 12,734 (26.1%) from NM [Figure 1B]. We performed Canonical Correlation Analysis (CCA) on all the high- quality single cells to detect anchors or mutual nearest neighbors (MNN) [14]. After integrating all the units, we conducted an unsupervised clustering analysis. Using Uniform Manifold Approximation and Projection (UMAP) with a resolution setting of 1.1, we successfully classified eight primary cell lineages: epithelial cells, myeloid cells, endothelial cells, fibroblasts, mast cells, B cells, T cells, and dendritic cells, as shown in Figure 1C. The identification of these cell types was primarily based on canonical cell markers and their functional categories, as determined by the genes significantly expressed in different clusters (refer to Supplementary Figure 1A-1D). These cell types were derived from normal, precancerous, and laryngeal cancer tumor tissues, and we noted varied proportions in each cell type. Notably, significant differences were evident between the tumors [Figure 1D,E]. The highest differentially expressed genes in each group were shown in the gene expression heat map [Figure 1F]. Epithelial cells represented the most abundant cell type. Their proportion notably increased in precancerous lesions compared to normal tissues. However, in tumor tis- sues, epithelial cell representation diminished, while immune cell populations increased significantly, aligning with existing studies [15]. Since proliferation is one of the main characteristics of tumor cells, we performed cell cycle analysis using a cell cycle gene score approach (Supplementary Figure 1E).
Characterization of transcriptional heterogeneity of epithelial cells in LSCC
Our study elucidates the transition from normal laryngeal mucosal epithelium to invasive cancer by examining the changes in expression and functional alterations in laryngeal squamous carcinoma epithelial cells across normal, precancerous, and cancerous stages. We analyzed 23,283 epithelial cells from these three stages, integrating and reaggregating them to identify distinct subpopulations based on marker gene expression levels. These were categorized into groups EP-C1 through EP-C12 [Figure 2A], with EP-C1 to EP-C10 present in all samples.
Figure 2: LSCC Epithelial Cells’ Diverse Expression Profiles. (A) UMAP Plot: Distinct subsets of epithelial cells are visualized, with color coding indicating different cell types. (B) Tissue Distribution: Depicts the expression of epithelial cell subsets across various sample tissues. (C) Active Factor Heat Map: Details active biological factors within each epithelial cell subset. (D) Oncogenic Pathway Heat Map: Illustrates the signaling in different epithelial cell subsets, focusing on oncogenic path- ways. (E) Gene Expression Heat Map: Displays characteristic genes in the epithelial gene co-expression module. (F) Co-expression Module Enrichment: Bubble plots show functional enrichment in various epithelial gene co-expression modules. (G) Cell Differentiation Trajectories: Maps the developmental pathways of epithelial cells based on differing expression programs.
Notably, EP-C12 was exclusive to tumor samples, while EP-C11 was found only in normal samples [Figure 2B]. Using the DoRothEA package, we analyzed transcription factors (TFs) in different epithelial cell clusters. The C10 subgroup expressed several tumor- related highly active TFs, including BACH1, associated with enhanced cancer cell invasion, metastasis, and poor prognosis in various cancers. Another key TF, NFE2L2 (NRF2), regulates antioxidant enzymes and plays a crucial role in oxidative stress and inflammatory response management. Clusters EP-C7, EP-C9, and EP-C10 highly expressed SREBF1, crucial for SCC survival and migration, with over- expression correlating with poor SCC patient survival. The Progeny package revealed strong signals in EGFR, p53, TNF, and other pathways in EP-C10, indicating regulation of cellular proliferation, differentiation, division, survival, and cancer development. EP-C12, derived from tumor samples, showed strong tumor signaling activation, while EP-C11, derived from normal samples, exhibited weak tumor signaling. GO and KEGG enrichment analyses revealed overlapping biological functions in some epithelial cell subsets. Building on Song L et al.’s categorization of malignant cells in laryngeal squamous cell carcinoma (LSCC), we conducted a secondary grouping of epithelial cells using differentially highly expressed functional genes as markers. We identified expression programs encompassing common patterns and abnormal activation pathways in epithelial cells during tumorigenesis [Figure 2E,F].
Further co-expression pattern analysis led to the classification of epithelial cells into Ciliated EPs, Metastatic EPs, Keratinocyte EPs, Immortal EPs, Proliferative EPs, and Immune EPs. Notably, Metastatic EPs expressed high levels of TSPAN18 and TM4SF1, while Keratinocyte EPs expressed SPRR2A, KRT6A, and KRT6B, potentially regulating proliferation and epithelial-mesenchymal transition. The Immortal program was characterized by the expression of MTRNR2L1, MTRNR2L8, and MTRNR2L12, and the Metastasis program by SERPINB5, GJA1, and TM4SF1. Functional enrichment analysis revealed associations with receptor antagonist activity and membrane permeability regulation [Figure 2F]. Finally, differentiation trajectory analysis showed that Immortal EPs appeared at the terminal differentiation stage [Figure 2G], Metastasis EPs were distributed throughout the trajectory, and Immune EPs, although scarce, expressed genes indicative of lymphocyte activation and cell adhesion regulation.
DNA copy number variation at each stage of LSCC development
We utilized the infer CNV method to identify recurring DNA copy number variations (CNVs) in epithelial cells. Profound CNVs were noted in premalignant lesion (PL) and laryngeal carcinoma (LC) tissues, contrasting with normal mucosa (NM) tissues. In PL tissues, 50 epithelial cells exhibited CNVs akin to malignant cells, termed carcinoma in situ (CIS) cells. CIS cells displayed significant DNA copy number gains at chromosomes 3, 18, 20, and 22, and losses at chromosomes 9,6, and 13, correlating with prior studies [Figure 3A].
Figure 3: Copy-number variation analysis of epithelial cells. (A) Copy number variation (CNV) in epithelial cells using normal tissue as a template, with CNV frequency displayed in the order of chromosomes by tissue type. (B) A heatmap illustrating the enrichment pathways of hallmark genes across different tissue types of epithelial cells in NM, PL, CIS, and LC. (C) Dot plot representation of differential expression of hallmark genes in epithelial cells of tissue types NM, PL, CIS, and LC. (D) Differential expression of CAN (Copy Number Alterations) frequency in the tissue types NM, PL, CIS, and LC.(E) Expression levels of the ALG1L gene across different tissue types. (F) Top ten pathways in the GSVA (Gene Set Variation Analysis) enrichment score for the ALG1L gene. (G) Validation of ALG1L expression levels in adjacent normal and head and neck tumor tissues from the TCGA-HNSCC (The Cancer Genome Atlas - Head and Neck Squamous Cell Carcinoma) database. (H) Immunohistochemistry (IHC) staining results in NM, PL, and LC tissues.(I) Quantification of relative proliferative activity in cells post-transfection. The bar graph shows mean values ± SEM of three independent experiments. The asterisks indicate a statistically significant difference between NC and Si ALG1L transfected cells (p < 0.01)
These cells initiated expression of genes linked to TNF- α/NF − kB, P53, estrogen response, and PI3 K/AKT/mTOR pathways, a trend intensified in malignant cells [Figure 3B]. Gene-level analysis showed CIS cells overexpressed known oncogenic markers (e.g., CD247, BARX1, CBX6, ARTN), affirming their malignant nature [Figure 3c]. CNA- driven transcriptional dysregulation is implicated in cancer progression. We identified specific CNA-dependent genes in CIS cells, showing altered DNA copy numbers and transcriptional changes across tissue types [Figure 3D]. Notably, ALG1L was overexpressed in CIS cells [Figure 3E], a finding corroborated by TCGA database analysis of head and neck tumors, revealing elevated ALG1L transcription compared to normal tissue [Figure 3G]. GSVA single gene functional enrichment analysis suggested ALG1L’s role in DNA replication, oxidative phosphorylation, and chromosome segregation. It may modulate DNA replication processes, contributing to chromosomal instability and mutation accumulation, hence advancing precancerous lesion development. ALG1L might also regulate cellular energy metabolism and oxidative stress responses, impacting cell growth and survival, thereby affecting lesion progression [Figure 3E]. siRNA-mediated overexpression experiments confirmed this, as laryngeal squamous cancer cells treated with siRNA targeting ALG1L (siALG1L) showed increased proliferation and migration [Figures 3J,I,K]. Immunohistochemical analysis verified ALG1L protein expression in cancerous epithelial cells of laryngeal squamous carcinoma and leukoplakia tissues, with higher expression than adjacent normal cells, localized in the cytoplasm and cell membrane [Figure 3H]. These findings collectively underscore ALG1L’s pivotal role in CIS cell mediated LSCC progression.
Tumor microenvironment plays an important role in the multi-stage development of LSCC
In our comprehensive analysis of the immune landscape within the tumor microenvironment, we utilized a series of hallmark marker genes and distinct gene expression profiles as a foundational framework to reclassify immune cells, successfully delineating 14 unique subsets [Figure 4A].
Figure 4: Tumor microenvironment plays an important role in the multi-stage development of LSCC (A) Label each of the 14 unique immune cell subsets. Highlight the gene expression profiles and hallmark marker genes used for this classification. (B) Annotate the prevalence of NK cells in pre-cancerous lesions. Emphasize their role in innate immunity and their capability to detect and destroy emerging cancer cells. (C) Show the increase in CD8+ISG+ T cells within the tumor microenvironment. Describe their activation state and role in the anti-tumor response. (D) Illustrate the use of DoRothEA in analyzing transcriptional dynamics. Focus on transcription factors (TFDP1, E2F2, FOXM1) and their roles in T cell exhaustion. (E) Distinguish between the three fibroblast phenotypes: CAFs, myofibroblasts, and NFs. Highlight key markers like FAP, α-SMA, and PDGFRβ for CAFs. (F) Use scatter plots to show the expression patterns of CXCL14, IFI27, and IL6. Include density plots for pseudotemporal ordering of fibroblast evolution. (G) Focus on the differential expression of JUNB, MGP, and NFKBIA.
This granular subcategorization allowed us to chart the complex roles played by a diverse array of immune cells during the multi- faceted stages of tumor development. During the formative stages of tumorigenesis, we noted a pronounced prevalence of Natural Killer (NK) cells within precancerous lesions [Figure 4B]. These cells, instrumental in innate immunity, possess the intrinsic capability to detect and destroy emerging cancer cells without the prerequisite of prior sensitization-a critical component of the body’s initial immune surveillance against malignancy [16]. Their targeting mechanisms are specialized to identify cells exhibiting abnormal expressions of self-molecules, a hallmark of early malignant transformation. As the tumor matures, our observations revealed a conspicuous increase in the presence of CD8+ISG+ T cells within the tumor microenvironment [Figure 4C]. Characterized by the expression of Interferon-Stimulated Genes (ISGs), these cells suggest an activated state, potentially instigated by interferon signaling a key defensive strategy against both viral pathogens and neoplastic cells [16]. The cytotoxic nature of CD8+ T cells serves to potentiate the anti- tumor response, addressing the heightened antigenic burden imposed by the progression of tumor cells [17]. This dynamic immunological interplay evolves as the malignancy advances; NK cells, initially serving as a front- line defense either through direct cytolytic activity against nascent tumor cells or by shaping the activity of other immune cell types, gradually cede their predominant role to CD8+ISG+ T cells [18]. These latter cells emerge as central figures in the adaptive immune response, engaging more sophisticated and targeted responses against advanced tumor cells that may exhibit immunoreactive traits. This transition from an innate to an adaptive immune response encapsulates the intricate, stage-dependent immunological shifts occurring during tumor progression. Our deep dive into the transcriptional dynamics of the tumor microenvironment was significantly propelled forward by leveraging DoRothEA, an expansive gene regulatory network (GRN), as shown in Figure 4D. Through this analytical lens, we could unravel the transcriptional underpinnings of T cell exhaustion - a state marked by diminished effector function and the upregulation of inhibitory receptors. Our data revealed a substantial overexpression of transcription factors such as TFDP1, E2F2, and FOXM1 within the exhausted T cell cohort. TFDP1, operating in concert with the E2F family of transcription factors, appears to orchestrate the modulation of cell cycle-related gene expression, potentially preserving a state of T cell quiescence amidst an immunosuppressive environment [19]. E2F2, a regulator of cell cycle checkpoints, may impose a senescent phenotype upon T cells, diminishing their responsiveness [20]. Meanwhile, the link between FOXM1 and cellular proliferation suggests its contributory role in sustaining the viability of functionally impaired T cells [21]. Together, these transcription factors sketch a regulatory landscape that perpetuates T cell exhaustion, underlining them as putative therapeutic targets. Modulating these factors offers a strategy to rejuvenate antitumor immunity, presenting a promising avenue for cancer immunotherapy. The implication of our findings paves the way for the development of novel therapeutics aimed at re-energizing the exhausted T cells, potentially restoring their capability to combat tumor growth and progression.
In our meticulous analysis, fibroblasts within the laryngeal carcinoma milieu were classified into three distinct phenotypes: cancer-associated fibroblasts (CAFs), myofibroblasts, and normal fibroblasts (NFs) [Figure 4E]. CAFs were distinguished by their robust expression of fibroblast activation protein (FAP), α-smooth muscle actin (α-SMA), and platelet-derived growth factor receptor beta (PDGFR β), markers emblematic of their heightened pro-tumorigenic activity [22]. Myofibroblasts, on the other hand, were marked by their expression of α-SMA, indicative of their contractile function; however, they lacked the full spectrum of CAF-specific markers [Supplementary Figure 2,3]. Leveraging Monocle for our single-cell trajectory analysis, we meticulously tracked the fibroblast state transitions through the cancer continuum [Figure 4F]. Scatter plots underscored the expression patterns of pivotal genes such as CXCL14, IFI27, and IL6 across stromal subpopulations, implicating their roles in inflammation, immune modulation, and cellular response. Density plots paired with pseudo-temporal ordering enriched our understanding of the fibroblast evolution within the tumor landscape. Notably, CAFs demonstrated a progressive activation within the tumor microenvironment, as opposed to NFs, which maintained a relatively subdued profile suggestive of a chronic inflammatory backdrop. Myofibroblasts appeared to occupy an intermediary niche, potentially orchestrating interactions between CAFs and NFs. Gene expression shifts in JUNB, MGP, and NFKBIA were particularly salient [Figure 4G]. JUNB is associated with cellular stress responses and proliferation, while MGP and NFKBIA play roles in immune and inflammatory path- ways [23,24]. Differential expression across fibroblast subsets implies varied contributions to tumor progression or suppression. The pronounced expression signature of CAFs may correlate with their roles in structural support, immune regulation, and facilitation of tumor cell survival within the TME. The expression profiles of myofibroblasts and NFs suggest their involvement in early tumor development stages, particularly in immune surveillance and the initial inflammatory response, underscoring the intricate interplay of stromal components in cancer pathophysiology.
Complex intercellular and molecular interaction networks in LSCC
Our study investigates the evolving landscape of cell communication networks through- out laryngeal cancer progression. Changes in signaling patterns underscore the dynamic crosstalk between immune cells and tumor cells within the tumor microenvironment (TME) [Figure 5A].
Figure 5: Analysis of intercellular communication networks (A) The network diagram displays the interaction outcomes between various cell subpopulations based on ligand- receptor paracrine signaling, with the thickness of the connecting lines representing the quantity of ligand-receptor interactions. The colored bar graph at the top represents the sum of column values shown in the heatmap (incoming signals). The colored bar graph on the right indicates the sum of row values (outgoing signals). In these color bars, red or blue denotes increased or decreased signals in the second dataset compared to the first. (B) The total number of interactions between NM, PL, LC in a bar graph. (C) Ranking important signaling pathways based on the differences in overall information flow in the network inferred between NM and PL. The top path ways in red are enriched in NM, while those in green are enriched in PL. (D) (E) The outgoing (or incoming) signal patterns between the NM and PL datasets, thereby identifying signaling pathways/ligand-receptor pairs that exhibit distinct signal pat- terns. (F) Identifying up-regulated (increased) and down-regulated (decreased) signal ligand-receptor pairs in a PL compared to NM. (G) Signal gene expression distribution.
Compared to normal mucosa (NM), premalignant lesion (PL) tissues exhibit augmented information flow across several key signaling pathways, notably those governing immune regulation and inflammation, such as IL- 16, EGF, and WNT [Figure 5B,C]. T cells in NM tissue predominantly engage the PI3K/AKT pathway, while PL tissues demonstrate a more diverse signaling repertoire, with pronounced T cell activity within the IFN- γ pathway. This distinction may reflect the immune cells’ strategic signaling adaptations at various tumor development stages [Figure 5D,E]. The amplified IFN- γ signaling in PL tissues likely represents an adaptive response to the increasing antigenic challenge posed by cancer cells [Figure 5G]. As the condition progresses to LC, signaling pathway activities, particularly NAMPT-INR and MIF-ACKR3, become more pronounced, potentially indicating an evolving tumor immune microenvironment and immune activation [18] [Supplementary Figure 4]. Macrophage migration inhibitory factor (MIF), a non-cellular TME constituent, is significantly elevated in several cancers, including liver, colorectal, and pancreatic cancers [25,26]. By binding with receptors like CD74, CXCR3, CXCR4, and CXCR7, MIF facilitates tumor proliferation, migration, invasion, angiogenesis, and chemo resistance. It contributes to anti-inflammatory responses, immune escape, and immune tolerance phenotypes in both innate and adaptive immune cells, ostensibly shielding cancer cells from immunogenic death and attenuating tumor immune responses [27]. MIF has been shown to activate the PI3K/AKT pathway via CD74 and CRCX7 interaction, suppress NR3C2, p53, and p27 expression, and elevate downstream tar- gets such as cyclins, cyclin-dependent kinases, MMP-7, and VEGF, thus fostering tumor progression and dampening apoptosis [28] [Figure 5F]. The AREG- EGFR pathway’s upregulation may correspond with tumor proliferation and survival, echoing the cancer cells’ tactic of exploiting growth signals for dissemination and apoptosis resistance [29]. Nonetheless, pathways like LGALS9-HAVCR2 and CXCL8-ACKR1 show diminished activity in LC tissues, suggesting possible evasion tactics by cancer cells against immune surveillance. Moreover, an overarching analysis of signaling patterns indicates a diminished information flow in certain cell types within the LC stage, possibly reflecting a redistribution of cellular populations in the TME or alterations in communication patterns triggered by tumor adaptation. Comparative analysis of signaling between PL and LC tissues unveils distinct information flow characteristics [Supplementary Figure 5]. For instance, while NAMPT and MIF pathways’ escalation in LC tissues may correlate with immune escape and inflammation during cancer progression, the AREG-EGFR pathway’s enhancement is linked to cell proliferation and survival. Concurrently, the down-trending interaction between LGALS9 and CD44 could signify a decline in immune surveillance efficacy. In summary, the progression of laryngeal cancer appears intimately linked with the dynamics of intercellular communication networks. The intricate interplay between immune and cancer cells throughout the cancer continuum offers valuable perspectives on laryngeal cancer’s biological underpinnings and could illuminate paths to novel therapeutic approaches.
DISCUSSION
This study, leveraging single-cell RNA sequencing data from eight untreated patients, has uncovered significant cellular diversity within normal laryngeal mucosa, laryngeal papillomatosis, and laryngeal carcinoma. Our analysis of 48,694 cells from diverse laryngeal tissues identified 24 cellular subgroups and eight primary lineages, underscoring the marked heterogeneity in both precancerous lesions and cancerous tissues. This variance was particularly evident in the shifting proportions of epithelial and immune cells across different tissue types. To further understand the impact of this diversity, we delved into the transcriptional heterogeneity of these cells and its implications for the progression of laryngeal squamous cell carcinoma.
Further exploration into the transcriptional heterogeneity of epithelial cells revealed their intricate relationship with the progression of laryngeal squamous cell carcinoma. By classifying 23,283 epithelial cells into distinct subgroups based on the expression levels of marker genes, we observed that certain subgroups, such as EP- C12, were exclusively present in tumor samples, whereas others, like EP-C11, were found only in normal samples. The analysis of transcription factors within these sub- groups highlighted the expression of key transcription factors- BACH1, NFE2L2, and SREBF1-suggesting these subgroups play pivotal roles in the survival, migration, and prognosis of laryngeal squamous cell carcinoma. These findings underscore the complex interplay between cellular diversity and cancer progression, leading us to investigate the molecular mechanisms further.
The study also delved into the impact of CNA-dependent transcriptional dysregulation on cancer progression, particularly emphasizing the role of the ALG1L gene. CNAs in epithelial cells were found to be significantly increased in precancerous lesions and laryngeal carcinoma tissues compared to normal mucosa, with specific CNAs possessing malignant features. This underlines the importance of CNAs in cancer development, as highlighted by the functional enrichment analysis which revealed these CNA cells’ genes were enriched in critical signaling pathways like TNF- α/ NF-kB, P53, estrogen response, and PI3K/AKT/mTOR.
Immunohistochemical analysis of 15 cases of vocal cord leukoplakia and 20 cases of laryngeal squamous cell carcinoma further validated the overexpression of ALG1L protein in cancerous epithelial cells, indicating its significant role in promoting proliferation and migration in laryngeal squamous cell carcinoma cells. By comprehensively analyzing the immune landscape within the tumor microenvironment of laryngeal squamous cell carcinoma, we classified immune cells into 14 unique subgroups, revealing the intricate roles these cells play across different stages of tumor progression. This in-depth analysis of both the cellular and molecular landscape provides critical insights into the nature of laryngeal cancer.
Despite the groundbreaking findings, this research is not without its limitations. The small sample size necessitates validation with a larger cohort to mitigate potential selection bias. Additionally, the lack of matched normal tissue, with only non-neoplastic samples from patients with vocal cord polyps, is a constraint. Paired normal tissue would serve as an ideal control to minimize individual differences. The ethical and practical challenges in obtaining normal tissue samples during biopsy procedures also add to this limitation. Furthermore, while we have identified key genes in laryngeal cancer progression through bioinformatics and cell function experiments, additional animal studies are needed to elucidate the underlying biological mechanisms more robustly. Ultimately, by combining various single-cell multi-omics approaches, including single-cell DNA sequencing, single-cell proteomics, and spatial transcriptomic, we can deepen our insight into the underlying mechanisms that influence the development and progression of laryngeal cancer.
CONCLUSION
In conclusion, our study illuminates the complex nature of LSCC, laying a foundation for future research and potential shifts in clinical strategies. By exploring tumor cell heterogeneity, progression models, and molecular alterations, we open new avenues for understanding and combating this intricate disease.
COMPETING INTERESTS
The authors announce that they have no known competing personal relationships or financial interests may influence this paper.
FUNDING
This research project is supported by Ningbo Top Medical and Health Research Pro- gram (No.2023030514) ; Ningbo medical and health brand discipline (No.PPXK2018- 02); Ningbo Clinical Research Center for Otolaryngology Head and Neck Dis- ease (No.2022L005); Ningbo ”Technology Innovation 2025” Major Special Project (2020Z097).
DATA AVAILABILITY STATEMENT
Any further inquiries about the data and materials used in this study should be directed to yidianchu@gmail.com.
ETHICS APPROVAL
For detailed information on relevant ethical standards and criteria, please refer to the sections on “Specimen collection and processing”.
REFERENCES
- Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin. 2019; 69: 7-34.
- Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68: 394-424.
- Rodrigo JP, Garc´?a-Pedrero JM, Su´arez C, Takes RP, Thompson LDR, Slootweg PJ, et al. Biomarkers predicting malignant progression of laryngeal epithelial precursor lesions: a systematic review. Eur Arch Otorhinolaryngol. 2012; 269: 1073-1083.
- Navin NE. The first five years of single cell cancer genomics and beyond. Genome Res. 2015; 25: 1499-1507.
- Campbell JD, Yau C, Bowlby R, Liu Y, Brennan K, Fan H, et al. Genomic, Pathway Network, and Immunologic Features Distinguishing Squamous Carcinomas. Cell Rep. 2018; 23: 194-212.e6.
- Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016; 352: 189-196.
- Kumar MP, Du J, Lagoudas G, Jiao Y, Sawyer A, Drummond DC, et al. Analysis of Single-Cell RNA-Seq Identifies Cell-Cell Communication Associated with Tumor Characteristics. Cell Rep. 2018; 25: 1458- 1468.e4.
- Zhang P, Yang M, Zhang Y, Xiao S, Lai X, Tan A, et al. Dissecting the Single-Cell Transcriptome Network Underlying Gastric Premalignant Lesions and Early Gastric Cancer. Cell Rep. 2019; 27: 1934-1947.e5.
- Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell. 2017; 171: 1611-1624. e24.
- Cavaliere M, Bisogno A, Scarpa A, D’Urso A, Marra P, Colacurcio V, et al. Biomarkers of laryngeal squamous cell carcinoma: a review. Ann Diagn Pathol. 2021; 54: 151787.
- Janiszewska J, Bodnar M, Paczkowska J, Ustaszewski A, Smialek MJ, Szylberg L, et al. Loss of the MAF Transcription Factor in Laryngeal Squamous Cell Carcinoma. Biomolecules. 2021; 11: 1035.
- Mo B-Y, Li G-S, Huang S-N, He W-Y, Xie L-Y, Wei Z-X, et al. The underlying molecular mechanism and identification of transcription factor markers for laryngeal squamous cell carcinoma. Bioengineered. 2021; 12, 208-224.
- Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014; 32: 381-386.
- Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, et al. Comprehensive Integration of Single-Cell Data. Cell. 2019; 177: 1888-1902.e21.
- Song L, Zhang S, Yu S, Ma F, Wang B, Zhang C, et al. Cellular heterogeneity landscape in laryngeal squamous cell carcinoma. Int J Cancer. 2020; 147: 2879-2890.
- Kucuksezer UC, Aktas Cetin E, Esen F, Tahrali I, Akdeniz N, Gelmez MY, et al. The Role of Natural Killer Cells in Autoimmune Diseases. Front Immunol. 2021; 12: 622306.
- Giles JR, Globig A-M, Kaech SM, Wherry EJ. CD8+ T cells in the cancerimmunity cycle. Immunity. 2023; 56: 2231-2253.
- Gonzalez H, Hagerling C, Werb Z. Roles of the immune system in cancer: from tumor initiation to metastatic progression. Genes Dev. 2018; 32: 1267-1284.
- Ostroumov D, Duong S, Wingerath J, Woller N, Manns MP, Timrott K, et al. Transcriptome Profiling Identifies TIGIT as a Marker of T-Cell Exhaustion in Liver Cancer. Hepatology. 2021; 73: 1399-1418.
- Ghram M, Bonnet-Magnaval F, Hotea DI, Doran B, Ly S, Des- Groseillers L. Staufen1 is Essential for Cell-Cycle Transitions and Cell Proliferation Via the Control of E2F1 Expression. J Mol Biol. 2020; 432: 3881-3897.
- Yu C, Sheng Y, Yu F, Ni H, Qiu A, Huang Y, et al. Foxm1 haploin- sufficiency drives clonal hematopoiesis and promotes a stress- related transition to hematologic malignancy in mice. J Clin Invest. 2023; 133: e163911.
- Nurmik M, Ullmann P, Rodriguez F, Haan S, Letellier E. In search of definitions: Cancer associated fibroblasts and their markers. Int J Cancer. 2020; 146: 895-905.
- P´erez-Benavente B, Fathinajafabadi A, de la Fuente L, Gand´?a C, Mart´?nez- F´erriz A, Pardo-S´anchez JM, et al. New roles for AP-1/ JUNB in cell cycle control and tumorigenic cell invasion via regulation of cyclin E1 and TGF-β2. Genome Biol. 2022; 23: 252.
- Rong D, Sun G, Zheng Z, Liu L, Chen X, Wu F, et al. MGP promotes CD8+ T cell exhaustion by activating the NF-κB pathway leading to liver metastasis of colorectal cancer. Int J Biol Sci. 2022; 18: 2345-2361.
- Balogh KN, Templeton DJ, Cross JV. Macrophage Migration Inhibitory Factor protects cancer cells from immunogenic cell death and impairs anti-tumor immune responses. PLoS One. 2018; 13: e0197702.
- Osipyan A, Chen D, Dekker FJ. Epigenetic regulation in macrophage migration inhibitory factor (MIF) mediated signaling in cancer and inflammation. Drug Discov Today. 2021; 26: 1728-1734.
- Noe JT, Mitchell RA. MIF-Dependent Control of Tumor Immunity. Front Immunol. 2020; 11: 609948.
- Liu R-M, Sun D-N, Jiao Y-L, Wang P, Zhang J, Wang M, et al. Macrophage migration inhibitory factor promotes tumor aggressiveness of esophageal squamous cell carcinoma via activation of Akt and inactivation of GSK3β. Cancer Lett. 2018; 412: 289-296.
- Chen Z, Chen J, Gu Y, Hu C, Li J-L, Lin S, et al. Aberrantly activated AREG-EGFR signaling is required for the growth and survival of CRTC1-MAML2 fusion-positive mucoepidermoid carcinoma cells. Oncogene.2014; 33, 3869-3877.