Laryngeal Squamous Cell Carcinoma: Potential Molecular Mechanism and Prognostic Signature Based on Immune-Related Genes

Background

Laryngeal squamous cell carcinoma (LSCC) is currently the most common pathological classification of laryngeal carcinoma (LC), accounting for more than 90% of cases [1–6]. Globally, more than 177 000 cases of LC are diagnosed every year, and more than 94 000 deaths from the disease [7]. Currently, treatments such as surgery, radiation therapy, and chemotherapy are used for early-stage LC [8,9], with favorable effects on survival [10–12]. However, more than two-thirds of patients are diagnosed with LC when the disease is already at an advanced stage [13–16]. This delay facilitates the growth and spread of LC cells [17], posing challenges for successful treatment of the disease.

Several pathogenic factors and mechanisms have been proposed for LC. Previous research has found that the incidence and development of it may be related to multiple risk factors, such as drinking [18–20], smoking [13], asbestos exposure [21], and human papillomavirus infection [22]. The mechanism may involve multiple molecules (such as MYC Target 1 [23], p16 [24], and NLK [25]), and multiple pathways (such as Wnt/β-catenin [26], ATR/p53/caspase-3 [27], and phosphatidylinositol-3-kinase/Akt/mTOR [28]); clearly, the onset and development of LC are complex. In summary, it is necessary to further explore the molecular mechanism of the pathogenesis and development of LC and identify novel and effective prognostic markers to assist in clinical screening and treatment of LC.

The immune system has been considered to affect the development, growth, and metastasis of tumors. Immune-related genes (IRGs) are considered to be relevant to tumor prognosis. For example, high expression of FGF19 may promote breast cancer progression by activating Akt signaling in cells, resulting in a discouraging prognosis for individuals with the disease [29]. Similarly, in esophageal squamous cell carcinoma, patients with TUBB3-negative carcinomas had higher rates of disease-free survival (DFS) and overall survival (OS) rates than those with TUBB3-positive disease [30]. In glioblastomas, another IRG, RBP1, is considered to be associated with poor prognosis [31]. Studies of LSCC and IRGs have found that programmed death-ligand 1 (PD-L1) expression is associated with favorable prognosis [32], and STC2 [33] overexpression is associated with poor prognosis. Clearly, IRGs are likely to contribute to the incidence, growth, and spread of LSCC. Therefore, exploring the relationship between IRGs and LSCC prognosis is crucial to understanding the pathogenesis of the disease and facilitating clinical screening and treatment.

High-throughput sequencing technologies are widely used in cancer multi-omics. In this study, we combined high-throughput data from public databases and published studies to explore the potential molecular mechanism of development of LSCC and to identify a prognostic signature based on IRGs to assist in clinical screening and treatment of LSCC.

Material and Methods

COLLECTION OF SAMPLE DATA AND IRGS:

The RNA sequencing (RNA-Seq) dataset from The Cancer Genome Atlas (TCGA) was obtained through the Genomic Data Commons Data Portal database (https://portal.gdc.cancer.gov/repository). The dataset contained 111 LSCC samples and 12 non-LSCC control samples. Clinical information from the samples in the RNA-Seq dataset was obtained from the University of California Santa Cruz Xena database (http://xena.ucsc.edu/). In addition, a list of 1811 IRGs was obtained from the ImmPort database (https://www.immport.org/home).

DIFFERENTIAL EXPRESSION ANALYSIS:

RNA-Seq count data were input into R (v3.6.1), and through the edgeR package [34], differentially expressed genes (DEGs) that had a log2-fold change of 2 or more and an adjusted P<0.05 were selected for further study. Genes were identified as upregulated (up-DEGs) if they showed a log2-fold change of 2 or more, while genes with a log2 fold change less than −2 were considered downregulated (down-DEGs). The same methods were used to screen immune-related differentially expressed genes (IRDEGs), immune-related upregulated DEGs (up-IRDEGs), and immune-related downregulated DEGs (down-IRDEGs).

ANALYSIS OF THE POTENTIAL MOLECULAR MECHANISM OF IRDEGS:

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of up-IRDEGs and down-IRDEGs were performed using the clusterProfiler package [35] in R (v3.6.1). GO terms and KEGG signaling pathways were then filtered by identifying those that simultaneously satisfied the following 3 conditions: (1) adjusted P<0.05; (2) q<0.05; and (3) number of enriched genes ≥5. The STRING database (https://string-db.org/) was used for analysis of protein-protein interactions (PPIs) of up-IRDEGs and down-IRDEGs. The minimum required interaction score for a PPI was set at 0.7 in the STRING database, which is considered by the database as high confidence. The CytoHub plug-in in Cytoscape software (v3.7.2) was applied to screen for hub genes, and the hub genes of up-IRDEGs and down-IRDEGs were identified according to the degree algorithm.

TRANSCRIPTIONAL REGULATOR (TR) PREDICTION FOR IRDEGS AND CONSTRUCTION OF TR-MEDIATED REGULATORY NETWORK:

To further understand the potential molecular mechanisms of up-IRDEGs and down-IRDEGs in LSCC, we attempted to explore the upstream regulation mechanism of these genes through the correlation between the genes and transcriptional regulators (TRs). TRs for genes were predicted using Epigenetic Landscape In Silico (Lisa, http://lisa.cistrome.org/), a bioinformatics analysis tool [36] that contains a large amount of H3K27ac ChIP-seq data from Homo sapiens. We predicted TRs related to up-IRDEGs and down-IRDEGs with the Lisa tool and selected TRs with P<0.01. Pearson coefficients of expression levels of TRs and IRDEGs were calculated in R, while the TR-IRDEG regulatory network was constructed using Cytoscape.

DEVELOPMENT AND ASSESSMENT OF IRDEG SIGNATURE:

Univariate Cox analysis and least absolute shrinkage and selection operator (LASSO) were used to select the candidates for signature from IRDEGs, and a multivariate Cox proportional hazards model was applied to construct the IRDEG signature of LSCC based on candidate IRDEGs (Figure 1). Using univariate Cox analysis, IRDEGs related to the survival of patients with LSCC (P<0.05) were identified and named as preliminary candidate IRDEGs. Then, based on the LASSO regression analysis, which can quickly and effectively extract important variables from many variables, secondary candidate IRDEGs were obtained from the preliminary candidate IRDEGs. Eventually, through the multivariate Cox proportional hazard model, some of the secondary candidate IRDEGs were used to construct the IRDEG signature of LSCC.

Curves for receiver operating characteristics (ROCs) and survival analysis (Kaplan-Meier) were used to evaluate the relationship between the IRDEG signature and the effects on screening and prognosis, respectively. Univariate and multivariate Cox regression analyses were used to identify factors related to independent prognoses in the IRDEG signature and clinical parameters. Development and assessment of the IRDEG signature were completed in R with the help of the glmnet [37], survival, and survivalROC package [38].

RELATIONSHIP BETWEEN RISK SCORE OF IRDEG SIGNATURE AND IMMUNE INFILTRATION LEVEL:

The Tumor IMmune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) is an online tool for comprehensive analysis of tumor infiltration of immune cells. The table matrix of TCGA estimation in TIMER was used to analyze the correlation between the IRDEG signature infiltration level risk score and 6 types of immune cells: B cells, CD4 T cells, CD8 T cells, macrophages, neutrophils, and dendritic cells.

APPLICATION OF IRDEG SIGNATURE:

Combining the IRDEG signature and the clinical parameters related to the independent prognosis of LSCC, a nomogram was constructed to quantify the risk of LSCC in individuals in a clinical environment. A calibration curve was used to detect the difference between the predicted and actual survival rates in the nomogram. The concordance index (C-index; range 0–1) was used to objectively evaluate the predictive power of the signature; the larger the C-index, the better the predictive power of the signature. Construction and verification of the nomogram were performed in R.

STATISTICAL ANALYSIS:

Except for differential expression analysis, for which original counts should be used, the RNA-Seq data used in the rest of the analysis were log₂ (counts+1) converted data. All statistical analyses were performed in R.

The main processes in the present study are shown in Figure 1. Unless otherwise specified, P<0.05 for a difference was considered statistically significant.

Results

IDENTIFICATION OF DEGS AND IRDEGS IN LSCC:

After filtering, 4274 DEGs were included: 2990 up-DEGs and 1284 down-DEGs. Among the 1811 IRGs, there were 302 IRDEGs, including 224 up-IRDEGs and 78 down-IRDEGs. Supplementary Figure 1 shows the expression of DEGs and IRDEGs in LSCC.

GO, KEGG, AND PPI ANALYSES OF IRDEGS:

To analyze the potential molecular mechanism of LSCC, we conducted GO, KEGG, and PPI analyses of IRDEGs.

GO analysis was performed on 224 up-IRDEGs and 78 down-IRDEGs. Among up-IRDEGs, the most common biological terms were “complement activation, classical pathway” (“biological process”), “blood microparticle” (“cell component”), and “antigen binding” (“molecular function”) (Figure 2A–2C). For down-IRDEGs, the most common biological process, cell component, and molecular function terms were “antimicrobial humoral response,” “azurophil granule lumen,” and “receptor ligand activity,” respectively (Supplementary Figure 2A–2C). Table 1 shows the 4 most enriched terms among the 3 types of terms.

KEGG analysis results showed that the most enriched KEGG pathway for both up-IRDEGs (Figure 3) and down-IRDEGs (Supplementary Figure 3) was the “cytokine-cytokine receptor interaction.” The 4 most enriched pathways for up-IRDEGs and down-IRDEGs are shown in Figure 3, Supplementary Figure 3, and Table 1. The GO and KEGG analyses show that the biological terms enriched in IRDEGs and the KEGG signaling pathway are closely related to immunity. In fact, these results are not surprising, given that the IRDEGs used for KEGG analysis were originally genes closely related to immunity.

The results of the PPI analysis and the degree algorithm showed that in up-IRDEGs, Kininogen 1 (KNG1) and C-X-C motif chemokine ligand 10 (CXCL10) were hub genes in the development of LSCC, while in down-IRDEGs, the hub genes were elastase, neutrophil expressed (ELANE) and lysozyme (LYZ). Figure 4A and 4B, respectively, show the interaction between the first 30 up-IRDEGs and down-IRDEGs obtained with the Degree algorithm.

CONSTRUCTION OF A TR-MEDIATED REGULATORY NETWORK:

Based on the prediction results from the Lisa tool, we screened 103 up-TRs related to up-IRDEGs with a threshold of P<0.01. To reduce false positives in TR prediction results, we conducted a correlation analysis on the 96 up-TRs included in up-IRDEGs and RNA-Seq and finally selected 15 real up-TRs and 50 up-IRDEGs with Pearson correlation coefficients >0.5 and P<0.05 for constructing up-TR-IRDEGs. The results showed that SPI1, SP140, and signal transducer and activator of transcription 4 (STAT4) were likely to be the hub TRs regulating up-IRDEGs (Figure 5A).

In the same way, we screened 159 down-TRs, 147 of which were included in the RNA-Seq dataset used in this study. Finally, we selected 17 real down-TRs and 24 down-IRDEGs with correlation coefficients >0.5 to construct down-TRs-IRDEGs. The results indicated that zinc finger E-box binding homeobox (ZEB1) and Ikaros family zinc finger 2 (IKZF2) may be important TRs in regulating down-IRDEGs (Figure 5B).

CONSTRUCTION OF IRDEGS SIGNATURE:

Through univariate Cox analysis, 30 preliminary candidate IRDEGs related to the survival of patients with LSCC were obtained from 302 IRDEGs, as shown in Supplementary Figure 4. Based on LASSO regression analysis and selecting the largest lambda with an average error within 1 standard deviation, 14 secondary candidate IRDEGs for construction of an IRDEG signature were selected from as many as 30 preliminary IRDEGs (Supplementary Figure 5). In total, 14 secondary candidate IRDEGs – beta cellulin (BTC), EPO, fibroblast growth factor (FGF)19, FGF5, gastrin (GAST), immunoglobulin heavy chain variable (IGHV)3–7, IGHV6–1, interleukin 13 receptor subunit alpha 2 (IL13RA2, retinol binding protein 1 (RBP1), RAR related orphan receptor C (RORC), semaphoring 6C (SEMA6C), TNR receptor superfamily member 4 (TNFRSF4), TUBB3, and urocortin (UCN) – were used to construct the IRDEG signature. Finally, we used the multivariate Cox proportional hazards model to construct the IRDEG signature of LSCC. The final IRDEG signature contained 9 secondary candidate IRDEGs for LSCC. As shown in Table 2, of the 9 IRDEGs, 7 (FGF19, GAST, IGHV3–7, IL13RA2, RBP1, TUBB3, UCN) were highly expressed in LSCC, while IL13RA2, RBP1, and TUBB3 are prognostic risk factors for patients with LSCC (hazard ratio [HR] >1 and 95% confidence interval [CI] do not contain 1). Of the 9 IRDEGs, BTC and RORC are expressed at low levels in LSCC and BTC is a risk factor for patients with LSCC, while RORC is a prognostic protective factor (HR <1 and 95% CI do not contain 1). The 9 genes included in the IRDEG signature as well as the coefficients of each gene are listed in Table 2. The product of the expression value of each gene and the gene’s corresponding coefficient is the contribution of the gene to the risk score, and the sum of the contribution of all genes is the risk score, and we performed a risk score calculation for 111 LSCC samples.

VALIDATION OF IRDEG SIGNATURE:

To evaluate the screening and prognostic effects of the signature in LSCC prognosis, patients with LSCC were divided into high- and low-risk groups according to scores higher or lower than the median risk score (0.859).

In terms of screening, the ROC curve showed that the risk score calculated according to the signature was the best measure for determining the survival of patients with LSCC, with a higher area under the curve than any clinical parameter (age, sex, neoplasm histologic grade, stage, tumor stage, node stage; Figure 6A) and any single gene included in the IRDEG signature (Figure 6B). Moreover, we also compared the performance of the signature to potential markers related to the prognosis of LSCC identified in research recently added to PubMed (published between January 1, 2020 and May 15, 2020), while 9 of the potential markers – activated leukocyte cell adhesion molecular (ALCAM) [39], ATM [40], BAF chromatin remodeling complex subunit (BCL11A) [41], B-cell lymphoma 2 (BCL-2) [42], desmoglein 2 (DSG2) [43], epidermal growth factor receptor (EGFR) [44], FGFR1 [45], homeobox A13 (HOXA13) [46], insulin like growth factor 1 receptor (IGF-1R) [47] – were included in the RNA-Seq data used in this study. Compared with any gene from these 9 potential prognostic markers, the risk score based on the IRDEG signature showed the best ability to screen for prognosis in patients with LSCC (Figure 6C).

In terms of prognosis, the OS rate for the high-risk group was significantly lower than that for the low-risk group (Figure 6D, P=1.234e-12). Moreover, the risk curve showed that the majority of patients with LSCC who died were those with a higher risk score, whereas there were significantly fewer deaths in the group with low-risk scores (Figure 7A). Figure 7B shows the expression of the 9 IRDEGs in the signature in the high- and low-risk groups.

In addition, we used Cox regression analysis to explore prognostic factors in LSCC. The results of univariate and multivariate Cox regression analysis demonstrated that the risk score based on the IRDEG signature can be considered a risk factor for the independent prognosis of LSCC (Figure 8A, 8B). In terms of risk of death, the risk for men with LSCC was 0.432 times that for women, whereas the node stage of the tumor was a risk factor related to the prognosis of LSCC (Figure 8B).

RELATIONSHIP BETWEEN IRDEG SIGNATURE AND CLINICAL PARAMETERS IN PATIENTS WITH LSCC:

To explore the clinical significance of IRDEGs in the signature of LSCC, we studied the expression of each IRDEG in the signature and its relationship with clinical parameters. Among the 9 IRDEGs included in the signature, BTC and RORC had low levels of expression in LSCC and the remaining 7 genes were highly expressed in LSCC. A Wilcoxon test showed that there were statistically significant differences in expression of 6 of the remaining 7 genes in the LSCC and non-LSCC samples. The exception was IL13RA2 (Figure 9A). Subsequently, based on data from 102 patients with LSCC for whom we had complete clinical information (Supplementary Table 1), we explored the relationship between risk score and expression level of the 9 genes and the clinical parameters in the patients (age, sex, neoplasm histologic grade, stage, and tumor and node stage), and found that the risk score and GAST were related to the neoplasm histologic grade of LSCC (Figure 9B, 9C), while the levels of expression of BTC, RORC, and TUBB3 were related to the age of patients with LSCC (Figure 9D–9F).

In addition, we analyzed the relationship between the IRDEG signature and the infiltration levels of B and CD4 and CD8 T cells; macrophages; neutrophils; and dendritic cells and found that the risk score based on the IRDEG signature was negatively correlated with the level of B-cell immune infiltration (Supplementary Figure 6).

CONSTRUCTION AND VERIFICATION OF NOMOGRAM BASED ON IRDEG SIGNATURE:

To investigate a preliminary application of the IRDEG signature, we constructed a nomogram of it and the clinical parameters related to independent prognosis of LSCC and predicted 1-, 3-, and 5-year survival rates in patients with LSCC (Figure 10A).

The calibration curve revealed that although the 3- and 5-year survival rate predictions by the nomogram were not in good agreement with actual patient survival (data not shown), the 1-year survival rate predictions were close to actual survival (Figure 10B). Moreover, the more objective C-index was 0.777 (standard deviation=0.031), indicating that the nomogram could predict the survival rate with moderate accuracy.

Discussion

According to current research, a number of IRGs are closely related to the incidence and progression of tumors, potentially indicating that IRGs play an important role in cancer [32]. Therefore, new markers for screening and treatment of tumors could be discovered through exploration of and research about IRGs.

Using the RNA-Seq dataset containing 123 samples, we first explored the potential molecular mechanism of LSCC from the perspective of IRGs, and for the first time, constructed an IRDEG signature related to the prognosis of LSCC based on IRGs. Through GO and KEGG analysis, we found, as expected, that IRDEGs may participate in the progression of LSCC through immune-related reactions. PPI analysis identified hub genes among IRDEGs in LSCC. We then constructed a regulatory network of TRs and IRDEGs to explore the underlying molecular mechanism of LSCC at the upstream level. In addition, the risk score based on the IRDEG signature that we constructed did well in distinguishing the prognosis of patients with LSCC. The risk score based on the IRDEG signature was an independent prognostic risk factor for LSCC, indicating that the IRDEG signature can be used as a novel marker to predict prognosis in patients with LSCC.

Initially, we explored the underlying molecular mechanism of LSCC from the perspective of IRGs. To this end, we conducted enrichment and PPI analyses on IRDEGs. In up-IRDEGs, the most common GO terms in biological processes, cellular components, and molecular functions were “complement activation, classical pathway,” “blood microparticle,” and “antigen binding,” respectively, while in down-IRDEGs, the terms were “antimicrobial humoral response,” “azurophil granule lumen,” and “receptor ligand activity,” respectively. The KEGG analysis revealed that the most enriched pathways for both up-IRDEGs and down-IRDEGs was the “cytokine-cytokine receptor interaction.” It is worth mentioning that quite a few GO terms and KEGG signaling pathways enriched in IRDEGs were closely related to tumor progression and immune response [48,49]; this, to a certain extent, indicates that IRDEGs may be involved in immune response, and thus, they participate in the occurrence and progression of LSCC. Through PPI analysis, we found that in up-IRDEGs, KNG1 and CXCL10 were likely the hub genes in the development of LSCC, and in down-IRDEGs, ELANE and LYZ were the hub genes. Interestingly, each of these 4 hub genes were considered a potential marker for screening or treating at least 1 cancer [50–53], suggesting that they are essential genes in LSCC. In addition, to explore the upstream regulatory mechanism of IRDEGs, we predicted TRs related to up-REDEGs and down-IRDEGs and plotted a TR-mediated regulatory network according to the co-expression relationship between TRs and IRDEGs. SPI1, SP140, and STAT4 were likely to be key TRs regulating up-IRDEGs; ZEB1 and IKZF2 may be the important TRs regulating down-IRDEGs. Although in previous studies, SPI1 [54] and STAT4 [55] have been identified as influencing tumors, insight into the regulatory mechanisms of these 5 TRs in LSCC is limited in the literature and more research is needed to verify their role.

To explore new markers that can be used for the screening and treatment of LSCC, we constructed a signature from 9 IRDEGs (BTC, FGF19, GAST, IGHV3-7, IL13RA2, RBP1, RORC, TUBB3, and UCN), and the screening and prognostic effects of it in LSCC were evaluated. In terms of screening, according to the ROC curves, the risk score calculated for the signature proved better at predicting the prognosis of patients with LSCC than any single parameter including age, sex, neoplasm histologic grade, stage, or tumor or node stage, as well as any single IRDEG used to create the signature. Moreover, the risk score also was better for screening and predicting prognosis than any of the recently reported prognostic markers for LSCC, namely, ALCAM [39], ATM [40], BCL11A [41], BCL-2 [42], DSG2 [43], EGFR [44], FGFR1 [45], HOXA13 [46], and IGF-1R [47]. In terms of prognosis, not only was the OS rate in the high-risk score group significantly lower than in low-risk score group, but the risk curve showed that there were significantly fewer deaths in patients with lower risk scores. The risk score also proved to be an independent prognostic factor for LSCC in multivariate Cox regression analysis. In summary, the IRDEG signature may be a novel marker for the screening and treatment of LSCC.

All the IRDEGs included in the signature other than IGHV3-7 have been identified as playing a role in tumors. Among them, BTC is considered to promote the development of head and neck squamous cell carcinoma (HNSCC) and has been shown to be a prognostic risk factor in patients with the disease [56]; this was confirmed by our research. We also observed a high level of expression of FGF19, which also has been identified in research about HNSCC [57] and which activates the ERK signaling cascade to stimulate the progression of nasopharyngeal carcinoma [58]. Another IRDEG signature gene, GAST, may contribute to the progression of gastric cancer at the hormone level [59] through participation in neuroendocrine function. In papillary thyroid carcinoma, with which IL13RA2 has been associated, the gene may be involved in cell migration by enhancing the epithelial-mesenchymal transition. RBP1 may be a factor associated with poor prognosis in glioblastoma [31]. RORC is thought to trigger proliferation of bladder cancer cells, and this process likely results from its ability to inhibit the PD-L1/integrin beta-6/signal transducer and activator of transcription 3 signal axis [60]. In patients with esophageal squamous cell carcinoma, those who were TUBB3-negative had considerably higher DFS and OS rates than those who were TUBB3-positive [30]. UCN has a dual effect on liver cancer cell metastasis due to its potential involvement in the regulation of cytosolic phospholipase A2 (cPLA2) and calcium-independent phospholipase A2 (iPLA2), that is, contributing to cell migration via promoting cPLA2 expression while preventing cell migration through reducing iPLA2 expression [61]. It is worth mentioning that, with the exception of TUBB3 [62] (which contributes to shorter progression-free survival and cancer-specific survival in patients), none of the remaining 8 IRDEGs have been reported in LSCC. This demonstrates the novelty of our work and reveals that further research into these 9 IRDEGs in LSCC is still needed.

To further explore the clinical significance of the IRDEG signature in LSCC, we analyzed the relationship between the signature, clinical parameters, and immune cell infiltration level in patients with LSCC. The results demonstrated that the risk score and the levels of expression of GAST were related to the neoplasm histologic grade of LSCC, while the levels of expression of BTC, RORC, and TUBB3 were related to the age of patients with LSCC; furthermore, neither the IRDEG signature nor the 9 IRDEGs were relevant to sex, stage, or tumor or nodes stage, but larger samples would be required to verify these observations. At the same time, the risk score based on the IRDEG signature was negatively correlated with the level of B-cell immune infiltration, suggesting that the signature may affect the occurrence and progression of LSCC by affecting the B-cell infiltration level. In short, we found that the clinical significance of the IRDEG signature in LSCC was mainly reflected in the neoplasm histologic grade and the B-cell immune infiltration level.

In this study, although we identified potential molecules involved in LSCC from the perspective of IRGs and established the first IRDEG signature related to the prognosis of the disease, there were several limitations. For example, the sample size was relatively small (n=123), and a study with a larger sample would be needed to verify the results. Moreover, because of the limited sample data, we could not externally verify the IRDEG signature, nor could we identify the relationship between the IRDEG signature and specific treatment methods, such as chemotherapy. Further experimental research is needed to verify the clinical significance and molecular mechanisms of the IRDEG signature in LSCC.

Conclusions

The IRDEG signature constructed in this study had a strong effect on the screening and prognosis of patients with LSCC and it may be considered an independent prognostic factor for use as a new marker in early screening and treatment of LSCC.