Abstract

The study discusses the challenges associated with the HNSCC metastasis and aims to identify hub genes related to the Cancer Stem Cells and metastasis in HNSCC through bioinformatics analysis. HNSCC is a prevalent type of cancer having a low five-year rate of survival, largely due to advanced stages at diagnosis and distant metastasis. The study hypothesizes that CTCs with CSC characteristics are responsible for metastasis. To investigate, we utilized mRNA-based stemness index (mRNAsi) and WGCNA on transcriptomic data by TCGA HNSCC cohort. The study involved multiple analytical steps, including DGE analysis, co-expression network, functional enrichment, hub gene identification, PPI network, and survival evaluation. The results revealed 18 hub genes associated with CSCs and metastasis in HNSCC. Further screening narrowed the hub gene list to nine genes, CDC45, MCM5, ASF1B, RFC4, E2F1, TK1, CHTF18, CENPM, and CDCA3.These genes have been discovered to be up-regulated in the samples of HNSCC compared to normal tissues. Highlighted for potential prognostic significance in HNSCC patients, CHTF18 and CDC45, with CDC45 playing a role in DNA replication, are associated with metastasis and CSC-like behavior. The study offers insights into potential hub genes associated with CSCs and metastasis in HNSCC. However, further experimental validation and clinical data integration are necessary to confirm these findings and establish clinical relevance. This study utilizes bioinformatics tools to identify genes influencing metastatic behavior in HNSCC, revealing potential therapeutic targets and underscoring the intricacies of cancer metastasis. Other hub genes are linked to various cancers, promoting invasion, proliferation, and migration.

Keywords: Head and neck squamous cell carcinoma, Metastasis, Hub genes, Cancer stem cells, mRNA-based stemness index, Weighted gene Co-expression network analysis

Introduction

Metastasis remains a major challenge in most solid tumors due to the lack of adequate treatment regimens and subsequent poor prognosis of the patients (Rus et al., 2020). HNSCC, the sixth most common cancer in the world, has a 50 % 5-year rate of survival, one of the major reasons being the advanced stage at presentation (Dahiya & Dhankhar, 2016). Patients who present having advanced disease are found to develop distant metastasis at a rate of twenty percent, which results in a poor prognosis for the patient (Schaeffer et al., 2014); the 5-year survival drops from 80% for the stage-I-II disease to about 20% in stage-III-IV. Distant metastasis rates in HNSCC range between 25-50%, the common sites being the brain, lungs, and bones, the rates of metastasis also differ between the various sub-sites of HNSCC (Nicoară et al., 2023; Radu et al., 2023). The current mode of detection of metastasis is imaging-based, PET scan be the best with a sensitivity of detecting 1-1.5 cm of tumor at the distant site, however, this does not provide for appropriate treatment management. Detection of patients with metastatic potential at presentation is a way forward to increase the rate of survival and patient’s treatment outcomes with advanced disease. Investigations down the decades have delineated the major processes underlying distant metastasis; dissemination of epithelial tumor cells into the blood, epithelial-mesenchymal transition (EMT), migration/invasion, and secondary tumor initiation. Circulating tumor cells (CTCs) in the bloodstream, considered to be the aggressive, migratory cells that disseminate from the primary tumor, are said to be the potential seeds that lead to distant metastasis. CTCs were majorly registered in the patient's peripheral blood with cancers of the breast, colorectal, lung, prostate, and oral cancer (Valiyaveedan et al., 2015; Kulsum et al., 2017) with associated prognostic efficacy, possibility of treatment surveillance, and disease progression (Maxime et al., 2016; Rutkowski et al., 2017; Lop et al., 2018; Matos et al., 2023). Cancer stem cells (CSCs) are known to attain EMT-like characteristics, with higher migratory and invasive capacity and retaining their tumorigenic property (Dahiya & Dhankhar, 2016). Further, reports suggesting the presence of non-migratory and migratory (Schaeffer et al., 2014) subgroups within CSC emphasize the need for a detailed investigation. Studies from the different groups have also shown that CSCs enriched by exposure to drugs do attain high migratory and invasive capacity, EMT marker profile, and increased properties of tumor initiation (Valiyaveedan et al., 2015; Kulsum et al., 2017). With the background of this evidence, we hypothesize that the CTCs with CSC characteristics are the ‘metastatic CTCs’ and the use of CSC-specific hub genes might provide insights into their accurate detection/targeting in HNSCC.

Over the past decade, the advancement of high-throughput technology has led to a substantial accumulation of data within public databases, furnishing a wealth of high-quality information for more in-depth exploration (Hajjo et al., 2021). The realm of oncological research has greatly benefited from the successful application of bioinformatics based machine-learning techniques. Among these, the One-Class Logistic Regression (OCLR) (Hemanth et al., 2020) stands out as a prominent method for encapsulating the prospects of stem cells. OCLR introduced the concept of a stemness index known as mRNAsi. In the context of this investigation, we delved into the clinical implications of mRNAsi and identified genes associated with CSCs (Cancer Stem Cells) or metastasis in HNSCC samples sourced from the TCGA dataset. This exploration was undertaken through the adept utilization of mRNAsi and WGCNA, harnessing the power of bioinformatics techniques. To execute this study, a variety of R-programming packages were employed, enabling a comprehensive analysis.

Materials and Methods

Software and R Packages

In this project, we utilized R version 4.3.1 running on the Windows platform. All the given R-packages were sourced from Bioconductor (Gentleman et al., 2004; Otasek et al., 2019), which provides a collection of open-source packages for bioinformatics research. The software and platform we employed are both open source and readily accessible.

Data Sources

The RNA-Seq data utilized in this research were sourced through TCGA, HNSC cohort, available at the UCSC database. The data set consisted of 500 primary site tumor samples, 2 metastatic samples, and 44 normal samples of tissue. Concurrently, relevant clinical data, including age, sex, life status, and stages, have been extracted through the same TCGA HNSC cohort. Subsequently, the expression values underwent log2 transformation (Ritchie et al., 2015). In instances where multiple expression values were available for a particular gene, we averaged the values to obtain are present active expression level. Additionally, the corresponding mRNAsi data, along with EREG-mRNAsi, were obtained from a referenced journal (Hemanth et al., 2020). To maintain data integrity, samples lacking clinical information and stemness index data were excluded from the subsequent analyses.

Correlation of the mRNAsi and Clinical Characteristics

The mRNAsi scores, based on the OCLR algorithm, were determined for stemless and stemness aspects, referred to as EREG-mRNAsi, and both scores ranged from 0 to 1. Our investigation involved a comparative analysis of mRNAsi across various aspects: normal vs. tumor group and tumor vs. metastasis group, conducted using an unpaired t-test. To assess significant changes in mRNAsi among different groups of variables, including treatment effects and different TNM (Tumor Node and Metastasis) stages, we utilized one-way ANOVA. All the above evaluations have been conducted by utilizing the GraphPad Prism version 8 (64-bit). For further analysis, the mRNAsi values (stemness index) were categorized into two subgroups as per the median cutoff value. We then conducted a Kaplan-Meier (K-M) evaluation for comparing the OS among the higher and lower subgroups, employing the log-rank test with the "survival" and "survminer" packages in R software. A p<0.05 was considered indicative of a statistically major difference.

Differential Gene Expression Analysis Statistically Significant Difference Was Presumed to Exist When the p< 0.05.

In this research, we employed the R edgeR (Robinson et al., 2010) package to conduct a comprehensive evaluation of DEGs (Mirza et al., 2019) among Tumor + Metastasis Vs Normal (T+M Vs N), Metastasis Vs Normal (M Vs N), Metastasis Vs Tumor (M Vs T). The identification of DEGs was carried out using the Wilcoxon signed-rank test. Our selection criteria for DEGs encompassed three key factors: a minimum absolute log2 fold change (FC) greater than 0.1, a p<0.05, and FDR (False Discovery Rate) below 0.05 for combinations M Vs N and T+M Vs N. For M vs. T a minimum absolute log2 fold change (FC) >F 0.1, a p-value less than 0.05, and an FDR below 0.34. To enhance the scope of our investigation and capture as many potential DEGs as possible between the sample groups, we adopted a criterion focused on low log fold change. Subsequently, we utilized the ensemble database R package "org.Hs.eg.db" to convert the Ensembl IDs of the differentially expressed genes into their corresponding gene names. To present our findings in a visually informative manner, we created volcano plots using the R-packages "ggplot" and "ggplot2”. These plots offered a clear representation of the DEGs identified in each of the examined groups.

WGCN Analysis

Co-Expression Network and Identification of Significant Modules

In this study, we constructed co-expression networks by utilizing the WGCNA package (Langfelder & Horvath, 2008) in the R-software (v4.3.1), based on the DEGs attained in the former step. The goodSamplesGenes function was employed to verify the completeness of the sample data and identify any outliers through sample clustering, which were subsequently removed. Pearson correlation coefficients were calculated between groups of genes and their absolute values have been utilized to create a gene expression similarity matrix by utilizing the mentioned formula: aij=|cor(xi,xj)|^β, xi&xj denotes the nodes i&j of the network, correspondingly. The parameter 'β' was carefully chosen to ensure the gene distribution adhered to a scale-free network based on connectivity. With the chosen 'β' value, we made the adjacency matrix and the TOM (Topological Overlap Matrix) (Pei et al., 2020). The TOM has been then subjected to gene dissimilarity-based clustering, and the resulting trees have been divided into distinct modules utilizing the dynamic shear approach. Each module was required to contain at least 50 genes. To enhance module identification, we combined some modules based on correlation coefficients. From the pool of modules, we have chosen the hierarchical clustering module which exhibited the closest connection to mRNAsi & EREG-mRNAsi for additional investigation. To gauge the importance of genes within each module, we utilized Module Eigengene (ME), Module Significance (MS), and Genetic Significance (GS). GS was determined by correlating gene expression with mRNAsi and EREG-mRNAsi, employing a log10 to transform the p-value in linear regression. This approach sheds light on the significance of every gene to specific traits. MS represented the average significance of all genes within a module, providing an overall measure of the module's importance. By merging similar modules by utilizing a cut-off value of 0.25, we identified modules having a great MS, deeming them most pertinent to sample traits. The first principal component discovered by evaluating the principal component on the gene expression matrix for every module is known as ME. The module with the greatest MS was chosen for further study because it was thought to be connected with mRNAsi and EREG-mRNAsi. Following the identification of major modules, we computed the GS and module membership (MM), which represents the association among the genes that make up a module also the gene expression patterns. To determine hub genes related to stemness indices, we set the inclusion criteria as follows: MM > 0.6 and correlation between gene GS > 0.4 and 0.1.

Functional Enrichment Analysis

Functional enrichment for GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) evaluates the module that is selected using the clusterProfiler package (Yu et al., 2012; Aschenbrenner et al., 2021). The GO evaluation involved three terms: molecular function (MF), cellular component (CC), and biological process (BP). To filter the results, the thresholds were set as follows: a modified P< 0.05 and a q<0.05.

Relationships and Interactions Among Hub Genes

The network of protein interactions may be studied and visualized using the online database STRING. Genes with a score of less than 0.4 in the STRING database were selected to create a network model that Cytoscape display (v3.10.0). CytoHubba, a plugin for Cytoscape, determined the MNC, MCC, Degree, and EPC of each node. The genes with the highest MNC, MCC, Degree, and EPC values in this analysis were regarded as hub genes.

Correlation of Hub Genes with Survival Data

The TCGA database's tumors and normal tissues analyzed for differential expression profiling, patient survival, and pathological staging using the web tool GEPIA2. To make up for the dearth of normal tissue samples in the TCGA database, GEPIA data are additionally taken from the sequencing data of normal tissues in the Genotype-Tissue Expression (GTEx) project 6 (Tang et al., 2017). To determine if the hub gene was differentially expressed in various TNM stages and whether its expression was greater in tumor tissue than in normal tissue, this database was employed. The hub genes were also analyzed using HPA (Human Protein Atlas).

Verification of the Hub Genes Using Human Protein ATLAS (HPA)

Immunohistochemistry (IHC) by the HPA database has been utilized for comparing the survival-related gene protein expression between HNSCC and normal tissues. For researchers, HPA is a useful database that offers a significant quantity of proteomics and transcriptomic data in certain human tissues and cells. Additionally, the most used method of using immune staining for determining the relative amount as well as the protein location is the IHC based on a pattern of protein expression.

Results and Discussion

mRNAsi Index is Significantly Associated with the HNSCC

We obtained the clinical data and transcriptome profiling for 546 Samples. The dataset included 44 samples of normal tissue, 2 samples from metastatic disease, and 500 samples of tumors. Age, gender, status of life, length of survival, illness stage, and classification of TNM (Tumor-Node-Metastasis) stage are all included in the clinical information. The mRNAsi index of each instance was obtained from the appendix of Malta's paper (Huang et al., 2019) and was then combined with the dataset from the TCGA. A total of six TCGA samples without mRNAsi information were excluded from the study.

Scatter plots were created, as shown in Figure 1, to show how the global mRNAsi profiles and clinical characteristics correlate. The difference in the mRNAsi between normal and HNSCC tumors is 44 and 494 samples respectively (Figure 1a). The difference in the mRNAsi among metastasis is 2 samples and for normal it is 44 samples (Figure 1b). The difference in the mRNAsi between metastasis is 2 samples and the tumor is 494 samples (Figure 1c). A comparison of mRNAsi between Tumor samples (494), normal samples (44), and metastasis samples (2) (Figure 1d). Variations in the mRNAsi according to the treatment impact; complete responder 217 samples, Partial responder 6 samples, stable disease 11 samples, progressive disease 10 samples, and potentially resettable disease 68 samples (Figure 1e). Variations in the mRNAsi among the male 394 samples and Female146 samples (Figure 1f). Comparison of the mRNAsi among 4 various T-stages (T0, 1 sample; T1, 48 samples; T2, 150 samples; T3, 103 samples and T4, 183 samples; Tx, 33 samples (Figure 1g). By comparing the mRNAsi among the 4 various N-stages (N0, 185 samples; N1, 70 samples; N2, 173 samples; and N3, 7 samples; Nx, 80 samples) (Figure 1h). Comparison of the mRNAsi among various M stages (M0, 187; M1, 1; Mx, 64) (Figure 1i). Comparison of The mRNAsi among 4 various stages of Tumor, stage I 27 samples; stage II 85 samples; stage III 84 samples and stage IV 276 samples (Figure 1j). Comparison of the mRNAsi between the vital status (Alive 292 samples; Dead 248) (Figure 1k). Kaplan–Meier survival curves depicted the low mRNAsi group having a good rate of OS in comparison to the high mRNAsi group (Figure 1l).

Differential Gene Expression Analysis

Differential genes observed with Tumor, metastasis Vs Normal 2532 genes (up-regulated 1244, down-regulated 1288) (Figure 2a). Differential genes observed between Metastasis Vs Normal 994 (up-regulated 618, down-regulated 376) (Figure 2b). Differential genes observed between Metastasis Vs Tumor 2823 (up-regulated 1397, down-regulated 1426) (Figure 2c).

WGCNA

Identification of modules associated with clinical information from the TCGA HNSCC patient sample data. After eliminating 25 outlier samples, the co-expression network included all DEGs (Figure 3a). β=5 met the soft threshold parameter of the construction requirements for scale-free distribution, and the curve reached R² = 0.85 (Figure 3b). To combine the comparable modules, the MEDiss Threshold was adjusted to 0.45, and 9 modules were eventually acquired (Figure 3c). The relationships between the modules and the characteristics of mRNAsi and metastasis were then assessed. By the mRNAsi and metastasis correlation coefficient R², we selected the positively correlated mRNAsi Brown (R²=0.57; P-Value = 2e-41) module and metastasis (R² =0.15; P-Value = 0.001) module (Figure 3d). Scatter plot of module eigengenes in the brown co-expression modules for mRNAsi and metastasis (Figures 3e and 3f).

For identifying hub gens the threshold for mRNAsi used is (GS>0.4 and MM>0.6) and for metastasis (GS>0.1 and MM>0.6). To identify the overlapping hub genes between mRNAsi Brown module, Metastasis Brown module, DEGs of M vs. T, and DEGs of M vs N, we found 18 common hub genes (Figure 4a).

Functional Enrichment Analysis

GO and KEGG pathway enrichment studies have been conducted for the modules that have been closely associated with the metastasis and mRNAsi, and the top 5 enriched outcomes are shown in (Figure 4b). The brown modules, that show the strongest positive correlation with the metastasis and mRNAsi, are greatly enriched in nuclear division, DNA geometric change, DNA duplex unwinding, BP (Biological Process), chromosomal region, kinetochore, condensed chromosome centromeric region, CMG complex and DNA replication pre-initiation complex under cellular composition (CC) and ATP-dependent activity, DNA helicase activity, acting on DNA, single-standard DNA binding and DNA replication origin binding Molecular Functions (MF). According to the outcome of the KEGG evaluation (Figure 4c), the DEGs have been mostly concentrated in the replication of DNA, oocyte meiosis, cell cycle, and cellular senescence.

PPI Network Construction and Screening of the Hub Genes

To further screen the hub genes, we have constructed a PPI network by utilizing the STRING (Search Tool for the Retrieval of Interacting Genes). For the common hub genes (Figure 4a) 18 nodes, and 59 edges have been created in the network, and the PPI improved the P-Value:<1.0e-16 (Figure 5a).

To create a network module that Cytoscape could see, the genes with a score of less than 0.4 were selected (v3.10.0). A plug-in called CytoHubba is used to compute the nodes. We utilized all four algorithms: Degree, MNC, MCC, and EPC. when comparing the top 10 genes between these algorithms 9 genes were found to be common in all and were considered as the top hub genes CDC45, MCM5, ASF1B, RFC4, E2F1, TK1, CHTF18, CENPM, and CDCA3 (Figure 5b).

Survival Analysis

Following the CytoHubba plugin's screening of the nine hub genes (CDC45, MCM5, ASF1B, RFC4, E2F1, TK1, CHTF18, CENPM, and CDCA3). We checked the hub genes levels of expression in the TCGA database's patient database. When compared to normal tissues, all nine of the hub genes are highly up-regulated in HNSCC carcinoma. To examine the predictive significance of the genes in the patients of HNSCC, the Kaplan-Meier plotter has been utilized to conduct OS, DFS, and pathological stage evaluation of the 9 hub genes by utilizing the R-survival package as well as GEPIA2 database (Figure 6).

Of the 9 hub genes, the Kaplan–Meier evaluation proposed that the greater expression level of CHTF18 has been majorly linked with the poor OS of the patients of the HNSCC (P<0.044) (Figure 6d), and for CDC45 a significant difference was observed in the DFS (Figure 7b) in HNSCC patients. Furthermore, the protein expression levels of MCM5, RFC4, CDC45, ASF1B, E2F1, and TK1 (Figures 7 J,H) show moderate to high expression of HNSCC and other cancers.

HNSCC is a group of malignant tumors that arise in the mucosal linings of the oral cavity, pharynx, and larynx. Despite advancements in treatment modalities, HNSCC continues to pose a significant clinical challenge due to its high potential for metastasis. In recent years, research has focused on understanding the underlying metastasis mechanisms as well as the role of CSCs in driving this process. In this research, we utilized a bioinformatics method to identify the hub genes that have CSC sand involved metastasis using mRNA stemness index data. When performing differential analysis, we found 18 hub genes and when further evaluated with PPI and Cytohub we identified 9 hub genes (CDC45, MCM5, ASF1B, RFC4, E2F1, TK1, CHTF18, CENPM, and CDCA3) which are all up-regulated in HNSCC TCGA data. Among these nine hub genes, CHTF18 and CDC45 have shown high significance in OS and DFS in HNSCC patients. Where the expression of MCM5, RFC4, CDC45, ASF1B, E2F1, and TK1 in the HPA data base is moderate to high in Neck and head cancer and other cancers (Seo & Kang, 2018). The CDC45 gene, known for its role in DNA replication, is now implicated in cancer metastasis and stem cell-like behavior. Increased CDC45 expression in hepatocellular carcinoma (HCC) has been associated with larger tumor size and advanced stage, potentially facilitating metastasis by promoting epithelial-mesenchymal transition (EMT) (Yang et al., 2020; Yang et al., 2021). Additionally, CDC45 appears to influence cancer stem cell (CSC)-like behavior, as observed in gastric cancer cells where CDC45 .knockdown reduced the CSC-like population (Huang et al., 2019; Rommasi, 2023). CDC45 and MCM5 were reported for their interaction in the human replicative helicase, where CDC45 is termed as the guard for the MCM5 (Petojevic et al., 2015; Seo & Kang, 2018; Nowinska et al., 2019; Schmit & Bielinsky, 2021). ASF1B enhances migration, invasion, and cell proliferation by influencing the P53-mediated EMT signaling pathway in cancer of the lung, male thyroid cancer, renal cancer, and pancreatic adenocarcinoma. This highlights the potential of targeting ASF1B for effective cancer therapy (Jiangqiao et al., 2019; Qiu et al., 2022; Wang et al., 2022; Huang et al., 2023). RFC4 is a DNA replication factor in NSCLC it promotes metastasis and stemness properties when it couples with NICD1 (Liu et al., 2021). Transcriptomic analysis of tumors lacking the E2F1 gene revealed its significance as a central controller of several genes promoting metastasis in human breast cancer. Moreover, studies unveiled certain genes influenced by E2F1 that have uncharted functions in pulmonary metastasis (Hollern et al., 2019). EF1 also has been reported for its progression and metastasis in melanoma and prostate cancer (Davis et al., 2006; Alla et al., 2010; Putzer et al., 2010; Marquardt et al., 2018; Jamroze et al., 2021; Singh et al., 2023). TK1 appears to play a role in the intricate procedure for melanoma metastasis and breast cancer metastasis by facilitating its spread through the deep lymphatic system and contributing to its progression. This suggests that TK1 might be involved in promoting the movement of tumor cells (CTC) through the primary tumor site to distant lymph nodes and potentially other parts of the body. This mechanism could be a critical factor in the aggressive nature of melanoma and breast cancer and its ability to establish secondary growths in distant locations (Wu et al., 2013; Del Re et al., 2019; Qian et al., 2023). CHTF18, an evolutionarily conserved protein involved in DNA replication, assumes vital functions in both mammalian fertility and the process of meiosis (Mukerji et al., 2016). Increased expression of CENPM promotes the metastatic spread of tumors through the mTOR/p70S6K signaling pathway in the context cancer of pancreatic (Zheng et al., 2020). Cell division cycle-associated protein 3 (CDCA3) plays a crucial role in bladder cancer (BLCA) as MYBL2, a critical transcription factor, binds to its promoter, transactivates it, and promotes BLCA cell proliferation and migration. Elevated MYBL2 levels in BLCA tissues are associated with clinic pathological parameters and poor cancer-specific survival. This study establishes CDCA3 as a downstream target of MYBL2 and highlights their collaboration with FOXM1 to activate Wnt/β-catenin signaling, ultimately driving BLCA malignancy. CDCA3 is also involved in the promotion of tumorigenicity and metastasis in many cancers (Dou et al., 2020; Yu et al., 2020; Liu et al., 2022).

Conclusion

This study has successfully employed mRNAsi and WGCNA techniques to identify pivotal hub genes within HNSCC, offering the potential to predict metastatic tendencies. These hub genes are acknowledged for their significance in diverse cancer prognoses and metastasis scenarios. However, it's important to acknowledge that this investigation primarily constitutes a bioinformatics analysis, and the absence of clinical trials is notable. Future research endeavors will be directed towards the accumulation of an expanded clinical dataset, Accompanied by the integration of pertinent clinical insights. This strategic approach aims to bolster the credibility of the hub gene’s association with HNSCC metastasis and CSC behavior. This will be achieved through a series of targeted cellular and molecular experiments, ultimately a factor in a more comprehensive understanding of the complex dynamics underlying HNSCC metastasis.

Acknowledgments: None

Conflict of interest: None

Financial support: None

Ethics statement: None