This article is mentioned in:

Introduction

Oral squamous cell carcinoma (OSCC) is one of the most common cancers worldwide. Annually, over 350,000 new patients are diagnosed with oral cancer worldwide, and over 170,000 patients die from this disease (1). Tongue squamous cell carcinoma (TSCC) is the most common type of OSCC (2,3). For example, in Japan, TSCC accounted for 30.4–43.8% of OSCC between 2000 and 2014 (4). TSCC presents with early metastasis and poor prognosis due to the abundance of blood and lymphatic vessels in the tongue, as well as the continuous movement of the tongue (5). Although surgery, radiotherapy and chemotherapy have demonstrated great progress in cancer treatment in recent decades, the 5-year survival rate (40–50%) remains low, and the prognosis of most patients is poor (6–10). Among the pathogenic factors of several TSCCs, reports have demonstrated an association between areca nut chewing and TSCC (11,12).

The areca nut is classified as a group 1 carcinogen, which can induce oral precancerous lesions and oral cancer (13). Oral submucous fibrosis (OSF) is a precancerous disorder caused by areca nut chewing (14). Areca nut extract can stimulate oral keratinocyte to secrete fibrosis-related cytokines, promote fibroblast proliferation and collagen synthesis, cause oral mucosal fibrosis, and ultimately induce OSF (15). Previous studies have demonstrated that OSF transforms into a malignant tumor in 1.5–15% of all cases (14,16) and increases the risk of oral cancer by 19.1-fold (17,18). In India, Pakistan and other countries or regions, oral cancer is at the forefront of malignant tumors due to the high prevalence of areca nut chewing, and consequently, the high morbidity of OSF (19,20). For example, 36.5% of men and 8.4% of women in India have areca nut chewing habits between the ages of 15 and 49 years in 2006, which led to a prevalence of OSF increased from 0.03% in 1968 to 6.42% in 2006, and >100,000 cases of OSCC are registered each year in India (20–24). Approximately 13.3% of residents in the Hunan province of China chew areca nuts, which has significantly increased the local oral cancer mortality (age-standardized mortality rates in 2013: 5.50 per 100,000) (25–27). However, the molecular mechanisms whereby areca nut chewing induces OSF and TSCC remain unclear.

A large number of genes are abnormally expressed during carcinogenesis, which promotes tumor initiation and progression, including p53 and MYC (28–30). Increasing evidence indicates that long non-coding (lnc)RNAs play important roles in carcinogenesis (31–40). However, the differential lncRNA expression profiles and related gene-regulatory networks associated with carcinogenesis of areca nut chewing-induced TSCC have not yet been reported. Thus, the present study collected five pairs of TSCC and adjacent normal tissues and monitored the resultant lncRNA- and mRNA-expression profiles using a lncRNA microarray. In addition, lncRNAs associated with areca nut chewing-induced TSCC and their potential gene expression-regulatory networks were screened, and the functions of these lncRNAs in TSCC were investigated.

Materials and methods

Sample collection

A total of five pairs of TSCC and adjacent normal tissues were randomly collected for a chip study at The Second Xiangya Hospital of Central South University (Changsha, China) between August 2015 and October 2015. All patients were male, with a mean age of 44.6 years (age range, 35–57 years). The inclusion criteria were as follows: Patients who reported a history of areca nut chewing, newly diagnosed with TSCC, who had not undergone preoperative radiotherapy or chemotherapy, and who did not have any other systemic disease. The exclusion criteria were as follows: Patients without a history of areca nut chewing and patients who had undergone preoperative radiotherapy or chemotherapy, or who had any other systemic diseases. In addition, the functions of the vital organs of all patients were in the normal range. Based on the International Union Against Cancer Eighth Edition TNM classification system for malignant tumors (41), 1 case was identified as T1N1M0, 1 case as T1N2M0, 2 cases as T2N0M0 and 1 case as T2N1M0 (Table SI). All tissue specimens were immediately frozen and stored in liquid nitrogen until subsequent experiments. The present study was approved by the Ethics Committee of Cancer Research Institute of Central South University (approval no. 2015-068; Changsha, China) and written informed consent was provided by all patients.

Gene chip

The present study used a 4×180K lncRNA array chip (Agilent Technologies, Inc.), which contains 78,243 lncRNA microarray probes and 32,776 mRNA microarray probes, and every microarray probe consists of 60 oligonucleotides. Data were uploaded to the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/; accession no. GSE139869).

Obtaining lncRNA expression profile data

Tissue RNA was extracted using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.; cat. no. 15596026), and all materials used were processed using RNAase-free diethyl pyrocarbonat water (Sigma-Aldrich; Merck KGaA; cat. no. D5758). Briefly, 50–100 mg of human tissue samples were ground into tissue homogenate in liquid nitrogen. After volatilizing the liquid nitrogen, 1 ml TRIzol reagent was immediately added. RNA was subsequently extracted according to the manufacturer's instructions (Invitrogen; Thermo Fisher Scientific, Inc.).

Sample labeling, microarray hybridization, and washing were performed according to the manufacturer's protocols (42). Briefly, total RNA was reverse transcribed into double-strand cDNA by using Low Input Quick Amp Labeling Kit (Agilent Technologies, Inc.; cat. no. 5190-2305). The reverse transcription procedure was as follows: 65°C for 10 min, 40°C for 2 h, 70°C for 15 min, at 4°C. Next complementary RNA (cRNA) was synthesized and labeled with cyanine-3-CTP. The labeled cRNAs were hybridized onto the microarray. After washing by using Gene Expression Wash Buffer Kit (Agilent Technologies, Inc.; cat. no. 5188-5327) twice, the arrays were scanned using a G2505C Agilent Scanner (Agilent Technologies, Inc.) and raw data was obtained using Feature Extraction software (version 10.7.1.1; Agilent Technologies, Inc.).

Microarray data were filtered to remove ‘background noise’ for non-expressed or very poorly expressed targets. After filtering the data, 64,556 probes remained, of which 24,149 corresponded to mRNAs and 40,407 corresponded to lncRNAs.

Principal component analysis (PCA) and differential expression analysis

To understand whether there are different characteristics between the areca nut chewing and the non-areca nut chewing TSCC expression profiles, PCA was performed on the two groups of data. PCA is commonly used for data dimensionality reduction, which can easily reveal the characteristics of different types of samples (43,44). The non-areca nut chewing TSCC lncRNA expression profile GSE34105 dataset was downloaded from the GEO database (45). These data are from Swedish samples, where there is no habit of chewing areca nut. PCA was performed using the programming language R (version 3.6.0, http://www.r-project.org).

Differentially expressed mRNAs and lncRNAs were screened using the significance analysis of microarrays (SAM) technique. The screening criteria were fold-change of ≥2.0 and P<0.05. Heat maps were constructed to identify differentially expressed mRNAs and lncRNAs, using Genesis software.

Chromosomal localization of differentially expressed mRNAs and lncRNAs by gene set-enrichment analysis (GSEA)

GSEA is a method for analyzing gene expression information, which helps determine the effects of gene expression changes on several biological functions and pathways, including T cell receptor signaling pathway (46,47). The present study used TSCC and adjacent normal tissues as group variables and used GSEA to perform enrichment analysis of positional gene sets (c1) with differentially expressed genes between TSCC and adjacent normal tissues, to predict the regional chromosomal changes (chromosomal deletions or amplifications). ‘c1: Positional gene sets (c1.all.v6.2. symbols)’ were used as the given sets of genes from the Molecular Signatures Database (MSigDB; http://www.gsea-msigdb.org/gsea/msigdb), and the number of permutations was set to repeat 1,000 times.

Establishment of mRNA and lncRNA co-expression networks using weighted gene co-expression network analysis (WGCNA)

In addition to chromosomal locations of mRNAs and lncRNAs, the associations between mRNAs and lncRNAs can be analyzed via their expression trends (48,49). WGCNA is a commonly used algorithm for constructing gene co-expression networks by gene expression trends (50). First, soft threshold was set to 24 to make network construction obeyed a scale-free network distribution, and a heat map of the topological overlap matrix (TOM) between genes was constructed. The thresholds were set and highly correlated co-expression networks were selected and visualized using Cytoscape software (version 3.7.0; http://www.cytoscape.org).

Predicting upstream regulatory factors and selecting related lncRNAs that may be regulated by these transcription factors

GSEA was performed to determine whether some of the differentially expressed lncRNAs and mRNAs are regulated by the same transcription factors and to predict the upstream regulatory factors of these differentially expressed genes. GSEA was also performed to determine whether the related lncRNAs are regulated by these transcription factors. ‘c3: Motif gene sets (c3. tft. v6.2. symbols. gmt)’ was used as the reference gene set (from the GSEA MSigDB database; http://www.gsea-msigdb.org/gsea/msigdb) for enrichment analysis to predict the upstream regulatory factors of differentially expressed genes, as previously described (46). The number of permutations was set to repeat 1,000 times. Genes regulated by these transcription factors were identified in the GSEA reference gene set and selected lncRNAs that are highly related to these genes (the topological-overlap values between these lncRNAs and more than half of mRNAs of the GSEA gene set exceeded the threshold of 0.4) in the co-expression network constructed by WGCNA. The related regulatory networks were constructed by Cytoscape software (version 3.7.0, http://www.cytoscape.org).

Analyzing canonical pathways of differentially expressed genes via ingenuity pathway analysis (IPA) and selecting pathway-related lncRNAs in the WGCNA constructed co-expression network

To identify functional lncRNAs in complex co-expression networks, the IPA tool was used to perform canonical pathway analysis on differentially expressed genes, and lncRNAs associated with significantly enriched pathways were screened in WGCNA-constructed co-expression networks. Based on the Ingenuity Knowledge Base (51), IPA can be used to identify biological networks and signaling pathways for specific datasets (51,52). By importing gene-expression profiles into the IPA program, signaling pathways that were significantly enriched for differentially expressed genes, as well as the genes involved in these pathways were identified. lncRNAs that were highly correlated with these genes (topological-overlap values between these lncRNAs and more than half of mRNAs of the IPA gene sets exceeded the threshold of 0.4) were selected in the WGCNA-constructed gene co-expression network.

Results

Differential expression of lncRNAs and mRNAs

Gene expression profiles were obtained from five pairs of TSCC and adjacent normal tissues, via mRNA- and lncRNA-gene chip analysis, and we eventually detected 78,243 lncRNAs and 32,276 mRNAs. Following preliminary screening to remove data for lowly expressed lncRNAs and mRNAs, 40,407 lncRNAs and 24,149 mRNAs remained. PCA analysis was subsequently performed and the results demonstrated that the areca nut chewing expression and non-areca nut chewing expression profiles obtained from the database (GSE34105) had different characteristics (Fig. S1). In addition, differentially expressed mRNAs and lncRNAs were screened via SAM, with the screening criteria set as follows: Fold-change ≥2.0 and P<0.05. A total of 3,860 differentially expressed RNAs were screened (Fig. 1A), including 1,667 mRNAs (866 upregulated and 801 downregulated; Fig. 1B) and 2,193 lncRNAs (1,057 upregulated and 1,136 downregulated; Fig. 1C).

Figure 1. Differentially expressed molecular heat maps of TSCC and adjacent normal tissues. (A) Heat map of differentially expressed mRNAs and lncRNAs. (B) Heat map of differentially expressed mRNAs. (C) Heat map of differentially expressed lncRNAs. Fold-change ≥2 and P<0.05. T represents TSCC tissue and N represents adjacent normal tissue. Red indicates upregulated genes and green indicates downregulated genes. TSCC, tongue squamous cell carcinoma; lncRNA, long non-coding RNA.

Chromosomal localization of differentially expressed mRNAs and lncRNAs via GSEA

GSEA positional gene sets can be used to identify regional chromosomal changes (chromosomal deletions or amplifications) that affect gene expression (53). The results of the present study demonstrated significant changes in three chromosomal segments via GSEA, including 22q13, 8p21 and 3p21. The upregulated genes were significantly enriched in the 22q13 segment (Fig. 2A), whereas the downregulated genes were significantly enriched in the 8p21 and 3p21 segments (Figs. S2 and S3). Taken together, these results suggest that chromosome amplification may occur in the 22q13 region, while chromosome deletion may occur in the 8p21 and 3p21 regions. In addition, according to the annotation information of the gene chip, 21 differentially expressed lncRNAs were in segment 22q13 (Fig. 2B). It is speculated that the mRNAs and lncRNAs expressed from this region were upregulated due to chromosomal amplification of the segment. Similarly, 13 downregulated lncRNAs from the 8p21 region and 51 downregulated lncRNAs from the 3p21 region were identified (Figs. S2 and S3), further indicating that chromosomal deletion occurs within these regions.

Figure 2. GSEA revealed that differentially expressed genes in TSCC and adjacent normal tissues are enriched in chromosome 22q13. (A) GSEA gene sets demonstrated that several genes located in the 22q13 segment were upregulated in TSCC. (B) Heat map of differentially expressed long non-coding RNAs (asterisks) and mRNAs from the 22q13 segment. GSEA, gene set-enrichment analysis; TSCC, tongue squamous cell carcinoma.

Construction of a co-expression network of differentially expressed lncRNAs and mRNAs

In addition to chromosomal co-localization of mRNAs and lncRNAs, a co-expression network based on similar expression trends between differentially expressed mRNAs and lncRNAs was constructed using WGCNA. First, the soft threshold for network construction was determined (Fig. S4), which made the constructed network conform to a scale-free distribution. Networks constructed in such a manner are similar to biological networks found in nature (50). Subsequently, the topological-overlap values among the molecules were calculated, and a hierarchical-clustering tree was constructed based on these values (Fig. 3A). Based on the topological-overlap values, genes with similar expression patterns were grouped into the same module. A TOM of the differentially expressed lncRNAs and mRNAs was constructed using the associated topological-overlap values, and the TOM was visualized in the form of a heat map (Fig. 3A). The expression network was divided into seven color modules in the present study. By setting the threshold, the highly correlated co-expressed molecules were screened (topological-overlap values >0.40) and a visual co-expression network was constructed using Cytoscape software (Fig. 3B). Due to the high threshold, four modules with fewer molecules and lower topological-overlap values between the molecules were almost filtered out. Thus, only three modules remained in the network, which consisted of 913 nodes, including 538 lncRNAs and 375 mRNAs, and 60,223 linkages (mRNA-lncRNA associations).

Figure 3. WGCNA-constructed mRNA-lncRNA co-expression network in TSCC. (A) Topological overlap matrix heat map of differentially expressed lncRNAs and mRNAs, which demonstrated that the network exhibits a scale-free topology. Red represents a lower overlap, while yellow represents a higher overlap. The top and left sides of the heat map are hierarchical cluster trees, with different branches of the cluster tree representing different gene modules, while the corresponding colors represent different modules. (B) Co-expression network of differentially expressed mRNAs and lncRNAs in TSCC and adjacent normal tissues, obtained via WGCNA and visualized using Cytoscape software. The network consisted of 913 nodes, including 538 lncRNAs and 375 mRNAs, and 60,223 linkages (mRNA-lncRNA associations). The topological-overlap value was higher than the threshold of 0.40. WGCNA, weighted gene co-expression network analysis; lncRNA, long non-coding RNA; TSCC, tongue squamous cell carcinoma.

Prediction of upstream regulatory factors and related lncRNA-mRNA network that these factors may drive

To determine which transcription factors regulate the differentially expressed lncRNAs and mRNAs in TSCC, the c3: Motif gene sets (c3.tft.v6.2.symbols.gmt; GSEA) was used as a reference gene set for enrichment analysis to predict the upstream regulatory factors (Fig. 4A). The most significantly enriched molecules were interferon regulatory factors (IRFs; P=0.003) and nuclear factor-kappa B (NF-κB; P=0.006; Fig. 4A). The downstream mRNA-lncRNA networks regulated by IRFs and NF-κB were constructed (Fig. 4B). The results demonstrated that the transcription of four differentially expressed genes (CXCL10, MSC, UBD and HOXB3) was simultaneously regulated by IRFs and NF-κB (Fig. 4B).

Figure 4. Prediction of IRF- and NF-κB-regulated mRNAs via GSEA, screening related lncRNAs and construction of regulatory networks. (A) The upstream regulators of differentially expressed mRNAs were predicted via GSEA. IRF- and NF-κB-regulated mRNAs were the most significantly enriched. (B) mRNA regulatory networks driven by IRFs and NF-κB. (C) The co-expression heat map of IRF-regulated mRNAs and lncRNAs (asterisks). In the heatmap, some mRNAs were derived from the reference gene set that was regulated by IRFs in the GSEA database, whereas others were significantly associated with IRFs in the WGCNA-constructed co-expression network. The lncRNAs that were significantly associated with IRFs in the WGCNA-constructed co-expression network are presented. (D) The co-expression heat map of NF-κB-regulated mRNAs and lncRNAs (asterisks). In the heatmap, some mRNAs were derived from the reference gene set that was regulated by NF-κB in the GSEA database, whereas others were significantly associated with NF-κB in the WGCNA-constructed co-expression network. The lncRNAs that were significantly associated with NF-κB in the WGCNA-constructed co-expression network are presented. (E) mRNA-lncRNA regulatory networks driven by IRFs and NF-κB. Yellow represents IRFs or NF-κB; green represents mRNAs derived from the reference gene set that was regulated by IRFs or NF-κB in the GSEA database; blue represents lncRNAs significantly associated with IRFs or NF-κB in the WGCNA-constructed co-expression network and red represents mRNAs significantly associated with IRFs or NF-κB in the WGCNA-constructed co-expression network. IRF, interferon regulatory factor; NF-κB, nuclear factor kappa B; GSEA, gene set-enrichment analysis; lncRNA, long non-coding RNA; WGCNA, weighted gene co-expression network analysis.

In addition to these downstream mRNAs, the present study investigated whether IRFs and NF-κB drive the expression of some lncRNAs in TSCC. The correlation between mRNAs and lncRNAs can be analyzed via their expression. Highly related molecules may participate in the same regulatory network and be regulated by the same regulatory factors. Thus, the present study screened for lncRNAs (Tables SII and SIII) whose expression trends were highly consistent with those of IRF-regulated or NF-κB-regulated mRNAs in the WGCNA-constructed co-expression network. As presented in Fig. 4C, the expression trends of the lncRNAs (Tables SII) and IRF-regulated mRNAs were similar, suggesting that the lncRNAs may also be regulated by IRFs. Similarly, as presented in Fig. 4D, the mRNAs and lncRNAs (Tables SIII) that were regulated by NF-κB in TSCC have similar expression trends. Based on the IRF- and NF-κB-regulated mRNAs enriched by GSEA and the mRNA-lncRNA co-expression network, the downstream mRNA-lncRNA networks regulated by IRFs and NF-κB in TSCC were constructed (Fig. 4E).

Canonical pathway analysis of differentially expressed genes via IPA and construction of related regulatory-network pathways

To further understand what functions and pathways the differentially expressed lncRNAs participate in, the gene expression profile dataset was imported into IPA to identify significantly enriched biological pathways in areca nut chewing-induced TSCC (Fig. S5A). The majority of the top five pathways (Th cell activation pathway, antigen presentation pathway, role of NFAT in regulation of the immune response, GP6 signaling pathway and T cell exhaustion signaling pathway) were associated with inflammatory and immunity. Collectively, these results suggest that inflammatory immune responses play important roles in the carcinogenesis of areca nut chewing induced TSCC. Given that the Th cell activation pathway was the most significantly enrichment pathway, lncRNAs that were associated with the Th cell activation pathway genes in the WGCNA-constructed co-expression network were screened (Table SIV). The expression trends of these lncRNAs were similar to mRNAs in the Th cell activation pathway (Fig. S5B). A network of lncRNAs that potentially regulates the Th cell activation pathway was derived (Fig. 5) by integrating the lncRNAs into the Th cell activation pathway (Fig. S6).

Figure 5. Network of lncRNAs involved in regulating the Th cell activation pathway in tongue squamous cell carcinoma. The lncRNAs were derived from the co-expression network and had similar expression trends compared with the mRNAs in the Th cell activation pathway. The triangles represent lncRNAs. Round and oval represent mRNAs. Red and pink represent upregulation, green represents downregulation and white represents no difference in expression. lncRNA, long non-coding RNA.

Discussion

Areca nut has been identified as a group 1 carcinogen that is closely associated with the development of OSCC and its precancerous lesions (11). Epidemiological evidence has demonstrated that in the Indian mainland, Southeast Asia countries and Hunan province of China, extensive areca nut chewing has caused a high incidence of OSCC, including TSCC (10,24,54). In certain countries, such as Sri Lanka, India and Pakistan, OSCC is the most common cancer, accounting for 25% of all new cancer cases (24). For example, more than 100,000 cases of OSCC are registered each year in India. The age-adjusted cancer incidence rates of OSCC in India were between 3.4 per 100,000 and 6.0 per 100,000 in 2001–2002 (24). However, the molecular mechanisms underlying areca nut chewing induced TSCC remain unclear. lncRNAs are non-coding RNA transcripts that exceed 200 nucleotides, and regulate gene expression at the epigenetic, transcriptional and post-transcriptional levels (55). Increasing evidence has demonstrated that aberrant lncRNA expression and dysfunctional lncRNA activity are associated with tumor initiation and progression (56–64). However, a full transcriptome analysis of the interactions between mRNAs and lncRNAs, and related regulatory networks in areca nut chewing induced TSCC have not yet been reported.

The lncRNA and mRNA profiles from TSCC were obtained by lncRNA microarray technology to determine the potential molecular mechanisms of areca nut chewing-induced TSCC. To determine whether there are different characteristics between the areca nut chewing and the non-areca nut chewing TSCC expression profiles, PCA analysis was performed and the results demonstrated that the components of the two data groups have obvious differences in clustering. This indicates that the expression profiles of patients with areca nut chewing and patients with non-areca nut chewing have different characteristics, and there may be specific differential expression genes involved in areca nut chewing induced TSCC. A total of 3,860 differentially expressed genes were identified, including 1,667 mRNAs and 2,193 1ncRNAs via SAM analysis.

OSF is a potentially malignant disorder associated with areca nut chewing and many patients with OSCC who habitually chew areca nuts develop OSCC from OSF malignant transformation (19). Through microarray analysis, Li et al (65) discovered several differentially expressed genes between OSF and normal tissues, were as Loricrin, COMP, CXCL9, KRT19 and CYP3A5 were the most differentially expressed genes and may play important roles in OSF. In the present study, CXCL9 and CYP3A5 were also differentially expressed and the expression trend was consistent with that of the report by Li et al indicating that CXCL9 and CYP3A5 may be involved in the entire process leading from normal mucosa to OSF and TSCC. Qiu et al (66) also constructed a lncRNA-expression profile for OSCC tissues obtained from individuals who did not habitually chew areca nuts. A total of 2,294 differentially expressed lncRNAs were obtained and the results demonstrated that lnc-TPP2-7:2, lnc-MANSC4-8:1, TMPRSS11BNL, GAS5 and DANCR are significantly downregulated in OSCC. The lncRNAs, GAS5 and DANCR were also downregulated in the present study, suggesting that they all play causative roles in OSCC induced by areca nut chewing or other reasons. However, most lncRNAs in the present differentially expressed lncRNA profile have not been previously reported. Some of these lncRNAs may be specifically associated with areca nut chewing induced TSCC and warrant further study. Regarding the lncRNAs that have not been previously studied, it is vital to determine their functional roles.

Genome instability is one of the common characteristics of cancer (67). During carcinogenesis, exogenous stimulation may induce chromosomal rearrangement, amplification, deletion, ectopic recombination or inversion (68). GSEA can be used to identify regional chromosomal deletions or amplifications. The present study performed chromosomal-localization analysis via GSEA and demonstrated that the differentially expressed genes were predominantly enriched in the chromosomal regions 22q13, 8p21 and 3p21. The upregulated genes were significantly enriched in the 22q13 region, whereas the downregulated genes were significantly enriched in the 8p21 and 3p21 regions. Taken together, these results indicate that chromosome amplification may occur in the 22q13 region, while chromosome deletion may occur in the 8p21 and 3p21 regions. If chromosome amplification or deletion occur in these regions, the lncRNAs located in such regions should also be abnormally expressed. The annotated positional information for the lncRNA chip was analyzed and 21 additional lncRNAs that were significantly upregulated were identified in the 22q13 region. Thus, it is speculated that chromosome amplification occurred in 22q13 and several mRNAs and lncRNAs were upregulated simultaneously in this region. Consistent with this result, Janjetovic et al (69) identified chromosomal amplification of 22q13 in patients with OSCC using comparative expressed sequence hybridization. In addition, Sakai et al (70) identified chromosomal amplification of 22q13 in patients with esophageal cancer by performing comparative genomic hybridization (CGH). The results of the present study demonstrated that multiple mRNAs and lncRNAs were significantly downregulated in the 8p21 and 3p21 regions, such as CCDC25 and FAM3D, which suggests that these 2 chromosomal segments may have been deleted in the areca nut chewing induced TSCC specimens, particularly considering that these deletions have been reported in different types of cancers, including OSCC (71–75). This finding needs to be further verified in more areca nut chewing induced TSCC tissues using CGH, fluorescence in situ hybridization or immunohistochemistry.

In addition to chromosomal colocalization analysis, the association between mRNAs and lncRNAs can be studied by constructing a co-expression network based on their expression trends. WGCNA is commonly used for constructing gene co-expression networks (50,76). When constructing co-expression networks, thresholds usually need to be set to determine whether two genes have similar expression patterns (50). Unlike the traditional ‘hard-threshold’ algorithm, WGCNA generates connections between genes in the network that obey scale-free networks by setting soft thresholds, which is more in line with biological characteristics (50). Thus, the present study constructed a co-expression network of differentially expressed lncRNAs and mRNAs using WGCNA. However, this co-expression network only provided lncRNA-mRNA associations at the transcriptome level, and so the functions of these lncRNAs remained unclear. Thus, the set of differentially expressed genes were imported into data analysis software programs based on known gene functions, such as IPA (Ingenuity Systems, www.ingenuity.com) and GSEA (version 3.0; http://software.broadinstitute.org/gsea/index.jsp), to find common transcriptional regulators and cell-signaling pathways.

GSEA for known upstream transcription factor-regulated gene sets demonstrated that IRFs and NF-κB that transcriptionally regulated gene sets were the most significantly enriched. IRFs and NF-κB are immunoregulatory factors that regulate the expression of several genes, particularly those involved in inflammation and immune responses, such as IFN-α and TNF-α (77,78). In addition, they are closely associated with tumor initiation and progression (79–81). For example, the chemokine, CCL5 promotes the invasion and metastasis of oral cancer cells by activating NF-κB (82). IRFs are also associated with the expression of immune evasion-related molecules, such as CTLA-4 and IDO (83,84). Jiang et al (84) reported that IRF-1 is an important factor for IFN-γ-induced IDO expression. These previous reports (82–84) are consistent with the results of the regulatory networks constructed in the present study, which proves that IRFs or NF-κB regulate areca induced TSCC through these downstream mRNAs derived from the reference gene set in the GSEA database. We constructed an mRNA-lncRNA co-expression network to establish the relationship between mRNAs and lncRNAs. Thus, lncRNAs that are associated with the downstream mRNAs regulated by IRFs or NF-κB in the co-expression network, may also be regulated by IRFs or NF-κB. A lncRNA and mRNA co-expression network was constructed to determine which lncRNAs are regulated by IRFs or NF-κB and the results identified 30 lncRNAs that were associated with IRFs or NF-κB, including NONHSAT085150. However, the functions of these lncRNAs have not yet been reported, and how these lncRNAs participate in IRF- and NF-κB-driven downstream mRNA-lncRNA-signaling pathways or promote TSCC merits further investigation.

IPA pathway enrichment analysis demonstrated that four of the top five pathways were associated with inflammation and immunity, including the Th cell activation pathway, antigen presentation pathway, role of NFAT in regulation of the immune response pathway and T cell exhaustion signaling pathway. These pathways differ from those of patients with OSCC that do not habitually chew areca nuts. Qiu et al (66) reported that the differentially expressed genes in patients with OSCC who do not chew areca nut are predominantly enriched in metabolic pathways, such as steroid biosynthesis, glycosaminoglycan biosynthesis and the synthesis and degradation of ketone bodies. These findings suggest that immune inflammatory responses may play more important carcinogenic roles in patients with oral cancer who chew areca nuts compared with patients who develop OSCC due to other causes. Autoimmune reaction is one of the pathogenic factors of OSF which is the main precancerous lesion of areca nut chewing induced OSCC (14). Long-term stimulation of exogenous antigens, and being repeatedly damaged and repaired places the oral mucosa in a state of chronic inflammation in patients who have chewed areca nuts for a long period of time (14,85). Long-term unresolved inflammation triggers cancer progression (86,87), thus, inflammation-mediated carcinogenesis may be one of the pathogenic factors of areca nut chewing induced TSCC. In addition to chronic inflammation, immune escape may also be one of the molecular mechanisms underlying areca nut chewing induced TSCC (88,89). Immune evasion plays important roles in oncogenesis and tumor development by protecting tumors from attack by the immune system (90,91). Th cell activation, antigen extraction and T cell exhaustion play important roles in immune escape, particularly for exhausted T cells (90). Exhausted T cells participate in the negative regulation of tumor immunity through inhibitory molecules and cytokines on the cell surface, thus causing tumor immune escape (90–101). In addition, some T cell exhaustion-related immune checkpoints were upregulated in the present study, such as CTLA4, TIGIT and Tim3. The differential expression of these immune checkpoints further suggests that immune escape occurs during areca nut chewing induced carcinogenesis. The present study aimed to investigate whether lncRNAs play important roles in these pathways via construction of a WGCNA co-expression network, and the results demonstrated that 74 lncRNAs were significantly associated with Th cell-associated pathways, including NR_026691.1. These 74 lncRNAs have similar expression trends to those of genes involved in the Th cell pathway, indicating that they may have the same function (48,49). Thus, these 74 lncRNAs were selected to study the regulatory networks that they may be involved in. However, the functions of these lncRNAs have not yet been reported, and their underlying molecular mechanisms require further investigation.

To the best of our knowledge, the present study was the first to investigate the differential lncRNA-expression profiles and the lncRNAs involved in signaling pathways and regulatory networks in areca nut chewing induced TSCC. The present study identified 1,667 differentially expressed mRNAs and 2,193 differentially expressed lncRNAs in five pairs of areca nut chewing induced TSCC and adjacent normal tissues via lncRNA gene chip detection and analysis. The majority of the differentially expressed lncRNAs have not yet been studied. Construction and analysis of a lncRNA-mRNA co-expression network generated a preliminary set of lncRNAs involved in biological signaling pathways and regulatory networks. The results of the present study provided insight for further investigation to determine the molecular mechanisms underlying areca nut chewing induced TSCC, provided candidates for selecting potential molecular markers and provided potential targets for clinical interventions.

A potential limitation of this study is the small sample size. Some genes that are less different between cancer and adjacent normal tissues but are important for the induction of cancer by areca nut may not have been screened. In addition, due to the small sample size, it is easy to screen some genes with false positive differences. With this regard, larger sample study will be carried out in the future to verify our results.

Supplementary Material

Supporting Data

Acknowledgements

The authors would like to thank Mr. Yanshan Ge (NHC Key Laboratory of Carcinogenesis, Cancer Research Institute and School of Basic Medical Sciences, Central South University, Changsha, China) and Mr. Zhihao Qu (Key Laboratory of Tea Science of Ministry of Education, Hunan Agricultural University, Changsha, China) for their technical support.

Funding

The present study was supported by the National Natural Science Foundation of China (grant nos. 81803025 and 81772901) and the “225” Project of Hunan Province High Level Health Talents.

Availability of data and materials

The datasets used and analyzed during the present study are available from the corresponding author upon reasonable request.

Authors' contributions

ZG conceived and designed the experiments. JL, YX, PL, SZ and YM performed the experiments. PL, SZ, YM, LZ, YW, FX and SZ analyzed the data. PL, JL, YX and ZZ drafted the initial manuscript. ZZ, WX, YL and ZG contributed significantly to data analyses and manuscript revision. All authors read and approved the final manuscript and agree to be accountable for all aspects of the research in ensuring that the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Ethics approval and consent to participate

The present study was approved by the Ethics Committee of Cancer Research Institute of Central South University (Changsha, China) and written informed consent was provided by all patients.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References