Daniel Hu,Diana V. Messadi
Immune-Related Long Non-Coding RNA Signatures for Tongue Squamous Cell Carcinoma
tongue squamous cell carcinoma,immune-related genes,lncRNA,survival analysis,prognosis
Background: Tongue squamous cell carcinoma (TSCC) represents one of the major subsets of head and neck cancer, which is characterized by unfavorable prognosis, frequent lymph node metastasis, and high mortality rate. The molecular events regulating tongue tumorigenesis remain elusive. In this study, we aimed to identify and evaluate immune-related long non-coding RNAs (lncRNAs) as prognostic biomarkers in TSCC. Methods: The lncRNA expression data for TSCC were obtained from The Cancer Genome Atlas (TCGA) and the immune-related genes were downloaded from the Immunology Database and Analysis Portal (ImmPort). Pearson correlation analysis was performed to identify immune-related lncRNAs. The TCGA TSCC patient cohort was randomly divided into training and testing cohorts. In the training cohort, univariate and multivariate Cox regression analyses were used to determining key immune-related lncRNAs, which were then validated through Cox regression analysis, principal component analysis (PCA), and receiver operating characteristic (ROC) analysis in the testing cohort. Results: Six immune-related signature lncRNAs (MIR4713HG, AC104088.1, LINC00534, NAALADL2-AS2, AC083967.1, FNDC1-IT1) were found to have prognostic value in TSCC. Multivariate and univariate cox regression analyses showed that the risk score based on our six-lncRNA model, when compared to other clinicopathological factors (age, gender, stage, N, T), was an important indicator of survival rate. In addition, Kaplan–Meier survival analysis demonstrated significantly higher overall survival in the low-risk patient group than the high-risk patient group within both training and testing cohorts. The ROC analysis indicated that the AUCs for 5-year overall survival were 0.790, 0.691, and 0.721, respectively, for training, testing, and entire cohorts. Finally, PCA analysis demonstrated that the high-risk and low-risk patient groups presented significant deviation regarding their immune status. Conclusions: A prognostic model based on six immune-related signature lncRNAs was established. This six-lncRNA prognostic model has clinical significance and may be helpful in the development of personalized immunotherapy strategies.
Immune-Related Long Non-Coding RNA Signatures for Tongue Squamous Cell Carcinoma

Background: Tongue squamous cell carcinoma (TSCC) represents one of the major subsets of head and neck cancer, which is characterized by unfavorable prognosis, frequent lymph node metastasis, and high mortality rate. The molecular events regulating tongue tumorigenesis remain elusive. In this study, we aimed to identify and evaluate immune-related long non-coding RNAs (lncRNAs) as prognostic biomarkers in TSCC. Methods: The lncRNA expression data for TSCC were obtained from The Cancer Genome Atlas (TCGA) and the immune-related genes were downloaded from the Immunology Database and Analysis Portal (ImmPort). Pearson correlation analysis was performed to identify immune-related lncRNAs. The TCGA TSCC patient cohort was randomly divided into training and testing cohorts. In the training cohort, univariate and multivariate Cox regression analyses were used to determining key immune-related lncRNAs, which were then validated through Cox regression analysis, principal component analysis (PCA), and receiver operating characteristic (ROC) analysis in the testing cohort. Results: Six immune-related signature lncRNAs (MIR4713HG, AC104088.1, LINC00534, NAALADL2-AS2, AC083967.1, FNDC1-IT1) were found to have prognostic value in TSCC. Multivariate and univariate cox regression analyses showed that the risk score based on our six-lncRNA model, when compared to other clinicopathological factors (age, gender, stage, N, T), was an important indicator of survival rate. In addition, Kaplan–Meier survival analysis demonstrated significantly higher overall survival in the low-risk patient group than the high-risk patient group within both training and testing cohorts. The ROC analysis indicated that the AUCs for 5-year overall survival were 0.790, 0.691, and 0.721, respectively, for training, testing, and entire cohorts. Finally, PCA analysis demonstrated that the high-risk and low-risk patient groups presented significant deviation regarding their immune status. Conclusions: A prognostic model based on six immune-related signature lncRNAs was established. This six-lncRNA prognostic model has clinical significance and may be helpful in the development of personalized immunotherapy strategies.

Tongue squamous cell carcinoma (TSCC) is the most common malignancy of the oral cavity [1,2]. A recent study suggested that about 45% of all oral cavity cancers were TSCCs [3]. TSCC generally affects the older age group of patients (e.g., 60 to 80 years of age) who have been exposed to tobacco and/or drinking for certain period of time. However, recent studies indicate that there has been an increasing incidence of TSCC in the younger age group [4,5,6]. TSCC in young age patients may not necessarily be due to direct exposure to smoking and drinking alcohol but might be caused by human papilloma virus infection [7]. Etiologically, it represents a distinct disease entity from TSCC in older patients. According to the National Cancer Institute Cancer Statistics, an estimated 17,860 new cases of TSCC and 2790 related deaths are projected to occur in the United States in 2022 [8]. The rate of new cases of tongue cancer is 3.6 per 100,000 men and women annually, with a death rate of 0.7 per 100,000 men and women per year. These rates are based on age-adjusted data from 2015–2019 cases and 2016–2020 deaths. In 2017–2019, approximately 0.4% of men and women were diagnosed with tongue cancer at some point in their lifetime. TSCC is generally an aggressive cancer with a relatively poor prognosis [2,9], and with a 5-year relative survival rate of 68.8% in the United States between 2012–2018 [8]. Despite advancements in treatment options, the survival rates for TSCC patients have remained relatively unchanged over the past few decades. Factors that contribute to the poor prognosis of TSCC include late diagnosis, high rates of metastasis, and the tendency for tumors to be aggressive and resistant to treatment. Early detection plays a critical role in the survival of TSCC patients, as the prognosis and five-year survival rate improve significantly with early diagnosis. Among TSCC patients, 29% are diagnosed at the local stage where the tumor is confined to the primary site, 53% are diagnosed at the regional stage where the tumor has spread to regional lymph nodes, 15% are diagnosed at the distant stage where the tumor has metastasized, and the remaining 4% are diagnosed with an unknown stage [8]. Based on data from the National Cancer Institute Cancer Statistics, the 5-year relative survival rates for localized, regional, distant, and unknown stage tongue cancer are 84.2%, 69.8%, 40.7%, and 53.0%, respectively [8]. These rates highlight the importance of early detection and timely treatment, as patients with localized tumors have significantly better survival rates compared to those with advanced disease. Therefore, regular screening and monitoring of high-risk individuals, such as those with a history of smoking and heavy alcohol consumption, can help with early detection and improve patient outcomes. The clinicopathologic diagnosis of TSCC does not always provide accurate predictions for clinical outcomes, highlighting the need for reliable prognostic biomarkers to guide treatment strategies. Recent advancements in molecular pathology have introduced thousands of tumor biomarkers associated with the progression and/or prognosis of various cancers, many of which have been evaluated for their prognostic value in OTSCC. However, despite the presentation of several promising biomarkers that could provide additional value to traditional factors such as stage, tumor grade, and depth of invasion, none have been approved for daily clinical use [10]. The aim of this study was to investigate the expression of immune-related long non-coding RNAs (lncRNAs) in patients with TSCC and identify a panel of lncRNAs that could serve as prognostic biomarkers for TSCC. LncRNAs are a type of RNA that are longer than 200 nucleotides and do not encode proteins [11]. Over the past decade, studies have revealed that lncRNAs have diverse roles in gene regulation, including modulating chromatin function, regulating nuclear bodies, and altering mRNA stability and translation. These functions can ultimately affect gene expression in various signaling pathways, such as those involved in immune responses and cancer [12]. While lncRNAs have been linked to head and neck cancer, including TSCC [13,14], little is known about their clinical applications or functional roles of immune-related lncRNAs in TSCC.

Figure 1 illustrates our workflow in developing an immune-related signature lncRNA model for TSCC prognosis. First, we extracted the lncRNA profiling data of TSCC patients from the TCGA, a landmark cancer genomics program with >20,000 primary cancer and matched normal samples spanning a large number of cancer types, and identified the lncRNAs at significantly differential levels in TSCCs. Subsequently, the Pearson correlation analysis was performed to identify differentially expressed immune-related lncRNAs. Next, the TSCC patients who have meaningful survival data (n = 145, excluding 4 TSCC patients with missing survival data) were randomly separated into training and testing cohorts (training to testing cohort ratio ~2:1). Both univariate and multivariate cox regression analyses were performed to build the signature lncRNA model. Finally, the testing cohort was used to validate the efficacy of our developed lncRNA model, and subsequent prognosis analysis, receiver operating characteristic (ROC) analysis, and PCA were performed to evaluate the signature lncRNA model.

As previously stated, our data were downloaded from the TCGA’s GDC portal where we extracted the transcriptome profiling data for 149 tongue cancer subjects from the TCGA head and neck squamous carcinoma project. The lncRNA expression levels were extracted and the list of immune genes was imported from the Immunology Database and Analysis Portal.

The downloaded gene expression and clinical data from the TCGA were organized and merged so that each patient’s transcriptome data corresponded to their clinicopathological data using Perl. Since this study was based on the premise of discovering key immune-related lncRNAs with prognostic value, patients with either unknown survival time or survival times of less than 30 days were removed. We then used the edgeR package in the R program (version 4.2.0) to filter out significantly differentially expressed lncRNAs between the cancer and normal samples, with a threshold of FDR < 0.05 and |log2FC| ≥ 1 [15]. With the list of both immune genes from the ImmPort and significantly differentially expressed lncRNAs, we conducted the Pearson correlation coefficient analysis, with a threshold of p < 0.001 and |R| > 0.4, to extract differentially expressed immune-related lncRNAs.

To identify potentially valuable lncRNAs for TSCC prognosis, we utilized the “survival” package in the R 4.2.3 to conduct univariate Cox regression analysis on the immune-related lncRNAs obtained. This analysis was performed to identify key lncRNAs (p < 0.01), which were then further screened using a stepwise multivariate Cox regression model. Based on the expression level of the lncRNAs and the regression coefficient β of the weighted linear combination in the multivariate analysis, the risk score of each patient was calculated using the following formula: Risk score = βgene1 × exprgene1 + βgene2 × exprgene2 + … + βgenen × exprgenen (βgenen: coefficient; exprgenen: expression level). Using the median risk score, the patients were divided into low-risk and high-risk groups.

The Kaplan–Meier curves were created to estimate the survival capability for the low- and high-risk groups being tested. Independent prognostic analysis was then carried out through univariate and multivariate cox proportional hazards regression analysis. Subsequently, the receiver operating characteristics (ROC) analyses were performed to compute the area under curve (AUC), verify our model, and to evaluate the accuracy of our model in comparison with other clinicopathological factors.

With the R packages “limma” and “scatterplot3d”, the PCA analysis was performed to visualize and compare the predictive power of all immune-related lncRNAs or the six immune-related signature lncRNAs.

The clinicopathological and transcriptomic data of 149 TSCC subjects were downloaded from the TCGA-HNSC project initially. A total of 16,902 lncRNAs were extracted from the complete transcriptome profiles, and through differential expression analysis with the edgeR package, 1108 lncRNAs were found to be significantly differentially expressed between the TSCC and normal tissues (adjusted p < 0.05 and |log2FC| ≥ 1). The heat map and volcano plot of significantly changed lncRNAs in TSCCs are shown in Figure 2. Using the list of immune related genes (n = 2483) from the ImmPort, the expression data of 1502 immune-related genes was extracted from the TSCC dataset, and then, through Pearson correlation analysis of the expression between 1108 lncRNAs and 1502 immune-related genes (mRNAs), 695 differentially expressed, immune-related lncRNAs were obtained (p < 0.001 and |R| > 0.4). These immune-related lncRNAs had either positive or negative modulation on the expression of the immune-related genes.

Patients with missing survival data or with survival times of less than 30 days were excluded from the model building. Therefore, 145 TSCC patients were included for subsequent analysis in this study. The demographic data of the 145 TSCC patients, including age, gender, and TNM stage, are shown in Table 1. Using the R software, the TSCC subjects were randomly divided into training and testing cohorts by a ratio of 2:1. In the training cohort, univariate analysis was performed on the expression profiles of the 695 differentially expressed, immune-related lncRNAs, resulting in 10 immune-related signature lncRNAs that had potential prognostic value (Table 2). For the ten immune-related signature lncRNAs, stepwise multivariate Cox regression analysis was performed to further screen the immune-related lncRNAs, resulting in six immune-related signature lncRNAs (Table 3). These six signature lncRNAs were then used to construct the prognostic model for TSCC. Using the following formula described in Section 2.4, the risk score was calculated for each patient. The TSCC patients were then divided into low-risk and high-risk groups based on the median risk score. The p value for each of the six lncRNA genes is <0.05, indicating that these genes are statistically significant.

Once the risk scores for the training cohort patients were all calculated using the formula above, the median risk score was then calculated and used to separate patients into high-risk and low-risk groups (Figure 3A). The distribution of the survival status between high-risk and low-risk patient groups is visualized in Figure 3B. The Kaplan–Meier survival curve demonstrated a relatively clear difference in patient survival probability between the high-risk and low-risk groups (Figure 3C). We then validated the model with the testing cohort by performing the same type of analysis, and the data analysis showed similar results to those of the training cohort (Figure 4).

The 3-year and 5-year ROC curves for the six immune-related signature lncRNAs based on the data of the training cohort were computed to evaluate the performance of the six-lncRNA model, resulting in AUC values of 0.773 and 0.795, respectively (Figure 5A). The same analysis was carried out for both the testing cohort (Figure 5B), which resulted in the 3-year and 5-year AUC values of 0.677 and 0.754, respectively, and the entire cohort, which resulted in the AUC values of 0.677 and 0.754, respectively (Figure 5C).

We performed both univariate and multivariate Cox regression analyses to determine if our model can be independently used for prognosis in TSCC. We included other clinicopathological characteristics such as age, gender, and stage (Figure 6) for comparison. Both univariate and multivariate Cox regression analyses indicated higher median hazard ratios with less variation of our six-lncRNA model when compared to the clinicopathological factors (age, gender, and stage). The results suggested that our six-lncRNA risk score model was a prognostic predictor independent of age, gender, and stage.

The ROC analysis was performed to determine the clinical significance of the six lncRNA signatures in comparison with other clinicopathological factors such as age, gender, stage, T, and N (Figure 7). The ROC curves indicate that our six-lncRNA model had the highest AUC value when compared to the other factors for both 3-year (Figure 7A) and 5-year overall survival (Figure 7B).

PCA analysis was performed to show the difference between the high- and low-risk groups based on the expression profiles of six immune-related signature lncRNAs or all immune-related lncRNAs (Figure 8). Clear separation between the two patient groups was demonstrated based on the expression profile of the six immune-related signature lncRNAs. The two patient groups tended to be distributed in different directions, indicating that there was a significant difference in the immune status between high-risk and low-risk patient groups.

lncRNAs represent a class of RNA molecules that are longer than 200 nucleotides and lack the ability to encode proteins. Once thought to be “junk” RNA, lncRNAs have been found to play important roles in a number of biological processes, including transcriptional regulation, RNA processing, post-transcriptional regulation, and protein synthesis, through various mechanisms, such as chromatin remodeling, transcriptional activation or repression, mRNA splicing, and protein interaction [16]. Moreover, lncRNAs have been implicated in diverse biological functions, such as cell proliferation, differentiation, apoptosis, and migration, as well as in development and aging [16,17]. One of the key functions of lncRNAs is their ability to interact with chromatin and regulate gene expression through chromatin remodeling and/or epigenetic modifications. For example, the lncRNA X-inactive specific transcript (XIST) plays a crucial role in X-chromosome inactivation in females by recruiting chromatin modifiers to silence gene expression on one of the two X chromosomes [18]. LncRNAs also have the ability to modulate the activity of other RNA species, including mRNAs and miRNAs, through various mechanisms such as competing for miRNA binding or serving as decoys to sequester proteins and prevent their interaction with other RNAs [19]. The dysregulation of lncRNA expression has been implicated in the pathogenesis of various diseases, including cancer. Many lncRNAs have been found to be differentially expressed in various types of cancer and are thought to be involved in tumorigenesis and cancer progression [20]. For example, HOTAIR (HOX transcript antisense RNA), a lncRNA that is overexpressed in breast, liver, and colorectal cancers, has been shown to promote tumor growth, metastasis, and drug resistance through various mechanisms, including the epigenetic silencing of tumor suppressor genes [21]. Numerous studies have provided insights into the mechanisms underlying the functions of lncRNAs in cancer. LncRNAs have been found to regulate gene expression at various levels, including chromatin remodeling, transcriptional regulation, post-transcriptional regulation, and epigenetic modification [17]. LncRNAs can interact with DNA, RNA, and proteins to exert their regulatory functions. For example, HOTAIR has been shown to interact with polycomb repressive complex 2 (PRC2) and lysine-specific demethylase 1 (LSD1) to regulate gene expression by modulating histone methylation and acetylation [22]. LncRNAs can regulate the expression and activity of oncogenes and tumor suppressor genes in cancer cells. For example, the lncRNA MALAT1 (metastasis-associated lung adenocarcinoma transcript 1) is overexpressed in various types of cancer and involved in cancer cell proliferation, migration, and invasion via promoting the expression of oncogenes such as TGF-β1 [23]. Similarly, the lncRNA HULC (highly upregulated in liver cancer) is overexpressed in hepatocellular carcinoma and promotes liver cancer growth through increasing the expression of the HMGA2 oncogene via sequestration of the microRNA-186 [24]. In addition, lncRNA CCAT1-L, which is transcribed specifically in human colorectal cancers, plays an important role in MYC transcriptional regulation and promotes long-range chromatin looping. Knockdown of CCAT1-L reduces long-range interactions between the MYC promoter and its enhancers. These results suggested an important role of CCAT1-L in gene regulation at the MYC locus in colon cancer [25]. On the other hand, the lncRNA PANDA (p21-associated noncoding RNA DNA damage-activated) inhibits the expression of the oncogene MYC and promotes the expression of the tumor suppressor p21, thereby inhibiting cancer cell proliferation [26]. In addition, lncRNAs may regulate the activity of signaling pathways involved in cancer pathogenesis, such as the Wnt/β-catenin, PI3K/Akt/mTOR, and MAPK/ERK pathways. For instance, the lncRNA H19 has been shown to promote tumor growth and metastasis by activating the Wnt/β-catenin signaling pathway in several types of cancer [27]. Similarly, the lncRNA MALAT1 has been shown to promote proliferation and metastasis of hepatocellular carcinoma cells by activating the MAPK/ERK pathway [28]. lncRNAs can serve as potential biomarkers for cancer diagnosis, prognosis, and therapy. Several lncRNAs have been shown to be associated with clinical outcomes, such as tumor stage, metastasis, recurrence, and survival, in different cancer types [29]. For instance, the lncRNA prostate cancer-associated transcript 1 (PCAT-1) is upregulated in prostate cancer and associated with poor prognosis and tumor progression [30]. LncRNAs have also emerged as promising therapeutic targets for cancer treatment. LncRNA-based therapeutics can be designed to target specific lncRNAs involved in cancer pathogenesis or to modulate their expression levels or activity. For instance, a peptide nucleic acid- based approach has been used to block the ability of lncRNA HOTAIR to interact with EZH2 and subsequently inhibit HOTAIR-EZH2 activity and resensitize resistant ovarian tumors to chemotherapy [31]. Furthermore, recent studies have shown that immune-related lncRNAs may be used to predict cancer patients’ survival and response to immunotherapy [32,33]. This highlights the potential of immune-related lncRNAs as prognostic biomarkers or therapeutic targets for clinical applications in cancer. Researchers have identified immune-related lncRNAs in various cancer types, including breast cancer [34], hepatocellular carcinoma [35], colorectal carcinoma [36], and renal cancer [37]. However, few studies have explored the clinical significance of immune-related lncRNAs in TSCC. As tongue cancer is one of the deadliest forms of oral cavity cancer, it is crucial to develop effective prognostic biomarkers for predicting malignancy behavior and response to treatments such as immunotherapy. Therefore, this study aimed to identify an effective immune-related lncRNA signature model for TSCC prognosis by performing a bioinformatic analysis of lncRNA gene expression data from TSCC patients. Through the use of univariate Cox regression analysis and multivariate Cox regression analysis of the lncRNA gene expression data of TSCC patients, we were able to identify six immune-related lncRNAs—MIR4713HG, AC104088.1, LINC00534, NAALADL2-AS2, AC083967.1, and FNDC1-IT1—as signature biomarkers for TSCC. Among the six lncRNA biomarkers, AC083967.1 and MIR4713HG were found to be of significant prognostic value in oral squamous cell carcinoma in previous studies, corroborating our results [38,39]. In fact, lncRNA MIR4713HG acted as a pro-tumor factor facilitating cell proliferation and metastasis of TSCC via the let-7c-5p/TMC7 signaling pathway, suggesting MIR4713HG might serve as a promising therapeutic target in TSCC [39]. To the best of our knowledge, the other four lncRNA biomarkers have not been studied in TSCC. However, lncRNA AC104088.1 was significantly overexpressed in non-small cell lung cancer (NSCLC) and suggested as a prognostic biomarker for NSCLC [40]. The expression level of lncRNA LINC00534 was found to be significantly increased in the serum samples of breast cancer patients [41] and also in lung squamous cell carcinoma tissues [42]. LncRNA FNDC1-IT1 was significantly overexpressed in breast squamous cell carcinoma and lung adenocarcinoma tissues and showed prognostic values predicting overall survival [43,44]. Studies also revealed significant upregulation of LncRNA NAALADL2-AS2 in metastatic castration-resistant prostate cancer (mCRPC) and demonstrated NAALADL2-AS2 as a biomarkers of response measure for survival in patients with mCRPC [45]. Androgen regulated lncRNA NAALADL2-AS2 promoted the survival of prostate cancer cells, and transcriptome and pathway analyses revealed that NAALADL2-AS2 modulated the expression of genes involved with cell cycle control and glycogen metabolism [46]. In addition, upregulation of NAALADL2-AS2 was found to be highly dependent on Epstein–Barr virus nuclear antigen 2 (EBNA2) [47]. Based on the six immune-related signature lncRNAs, we constructed a risk score model for the prognosis of TSCC s. To verify our model, the Kaplan Meier survival curves were utilized to illustrate the survival capability of the high- and low-risk patient groups. The p-value of the Kaplan Meier survival curve was 1.393 × 10−5 indicating a strong correlation between our risk score model and the survival outcome of TSCC patients. We further verified our lncRNA prognostic model through the testing cohort. In our ROC analysis, the AUC value for 3-year ROC of the training and testing cohorts were 0.773 and 0.677, respectively, demonstrating the proficient accuracy of the model. Univariate and multivariate Cox regression analysis on the training cohort, testing cohort, and entire cohort, demonstrated that our risk score model based on the six lncRNA biomarkers can be used as an independent prognostic factor. Further ROC analysis demonstrated that our six-lncRNA model had a significantly greater AUC in comparison to other clinicopathological characteristics such as age, gender, stage, N, and T. Lastly, PCA analysis demonstrated that our six-lncRNA predictive model could classify TSCC patients into different risk groups based on their immune status. Our study had one limitation, mainly, that the six immune-related signature lncRNAs were only tested on the patient samples from the TCGA database. Further validation studies are warranted to confirm the clinical utility of our developed lncRNA model. Moreover, in vitro and in vivo experiments would be needed to elucidate the molecular mechanisms of the six immune-related signature lncRNAs in TSCC.

In this study, we have demonstrated and verified six immune-related lncRNAs that may collectively serve as an independent prediction model for the prognosis of TSCC. In addition, the model can be used to distinguish TSCC patients between high- and low-risk groups in TSCC patient population, which are known to have significantly different survival outcomes. To the best of our knowledge, this is the first study on the identification and construction of an immune-related lncRNAs prediction model for TSCC. The discovered six immune-related signature lncRNAs may provide further insight regarding both TSCC prognosis and development of immunotherapy strategies.