Screening of tumor grade-related mRNAs and lncRNAs for Esophagus Squamous Cell Carcinoma.

BACKGROUND: The goal of our study was to screen tumor grade-related lncRNAs and mRNAs to reveal the underlying molecular mechanism of esophagus squamous cell carcinoma (ESCC).
METHODS: The lncRNA and mRNA sequencing data were obtained from The Cancer Genome Atlas (TCGA). Tumor grade correlation analysis of lncRNAs and mRNAs was executed, followed by the functional enrichment analysis of all tumor grade-related mRNAs. The differentially expression mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) were obtained. PPI network and DEmRNA-DElncRNA interaction analysis were constructed. The functional annotation of the DEmRNAs co-expressed with DElncRNAs was performed. The expression levels of the candidate genes were validated using qRT-PCR.
RESULTS: A total of 1864 tumor grade-related mRNAs (846 positively related and 1018 negatively related) and 552 tumor grade-related lncRNAs (331 positively related and 221 negatively related) were obtained. The top 10 significantly grade-related mRNAs and lncRNAs included CA12, FABP4, DECR1, BAIAP2, IL1RAPL2, PPARD, LAD1, TSPAN10, LDOC1, ZNF853, RP11-25G10.2, RP11-557H15.3, RP11-521D12.5, CHKB-AS1, RP11-219B4.3, CH17-335B8.4, RP11-99 J16-A.2, CTB-111H14.1, ADNP-AS1, and JHDM1D-AS1. SFN, IL1RAPL2, and RP11-25G10.2 were overlapped from grade 1, grade 2, and grade 3. PPI network showed that top 10 proteins with higher degrees, including GNAI1, RAP2B, GNAZ, SHH, ADCY1, PRKAR2B, SH3GL1, GNA15, and ARRB1. A DElncRNAs-nearby DEmRNAs network was constructed to obtain hub lncRNAs including ADAMTS9-AS2, RP11-210 M15.2, RP11-13 K12.1, ZBED3-AS1, and RP11-25G10.2. Except for RP11-25G10.2, ADAMTS9-AS1, ZBED3-AS1, SFN, ATP1A2, and GNA15 were consistent with our TCGA analysis.
CONCLUSIONS: Alterations of DEmRNAs and DElncRNAs may provide key insights into the molecular mechanisms of ESCC.

INTRODUCTION

Esophageal cancer is one of the most common malignancies worldwide. The incidence of esophageal cancer is threefold higher in men than in women. The predominant histological types of esophageal cancer are adenocarcinoma and squamous cell carcinoma. Esophageal squamous cell carcinoma (ESCC) is one of aggressive squamous cell carcinomas with high incidence. The high mortality rate has been attributed to late diagnosis and poor treatment response, among others. The World Health Organization ESCC histological tumor classification is as follows: grade 1 (well differentiation), grade 2 (moderately differentiated), and grade 3 (lowly differentiated). Therefore, screening of tumor grade‐related mRNAs and lncRNAs in ESCC is key issue to the treatment of ESCC patients.

Long non‐coding RNA (lncRNA), of more than 200 nucleotides in length, can regulate gene expression during several biological processes. Cumulative reports of abnormal lncRNA expression have shown that lncRNA may be used as a new independent biomarker for the early diagnosis, prognosis, and metastasis prediction of various cancer types. , , HOTAIR is a well‐known lncRNA whose expression level strongly predicts the metastatic and survival of different cancer types, including ESCC. , , However, the potential molecular mechanisms of ESCC remain unclear.

Here, we aimed to investigate the differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) associated with the different tumor grades including grade 1, grade 2, and grade 3. We first performed a tumor grade correlation analysis of lncRNAs and mRNAs expression data in patients with ESCC using the Cancer Genome Atlas (TCGA). Then, the screening of DEmRNAs and DElncRNAs between different tumor grades was performed, and this was followed by the DElncRNAs‐target DEmRNAs interactions network. The qRT‐PCR was applied to validate the expression of candidate gene. Our study identified potentially key mRNAs and lncRNAs in ESCC and provided a basis for further understanding the mechanisms of ESCC and the predictive performance of lncRNAs.

METHODS

Integrated profiles in TCGA

The TCGA data portal included the data of 185 patients with Esophagus Squamous Carcinoma, including clinical data of 185 patients, RNA sequencing data of 164 patients up to 24 May 2018. In this study, ESCC patients were allocated to grade 1, grade 2, and grade 3 groups in accordance with the recorded tumor grade. The lncRNA and mRNA gene expression profiles and clinical data of ESCC were downloaded from TCGA.

Identification of tumor grade‐related mRNAs and lncRNAs

The TCGA datasets were downloaded and transformed from Fragments PerKilobase Million (FPKM) data into Transcripts PerKilobase Million (TPM) data. Here, we used log2 of the TPM value to measure of mRNA and lncRNA expression level. The undetectable lncRNAs and mRNAs (with read count value =0 in more than 20% ESCC case) were filtered and deleted. The linear by linear association test was applied to analyze the correlation of the expression of mRNAs and lncRNAs with tumor grade by using the lbl.test function of the coin package in R‐3.3.3. p<0.05 was defined as the criteria of statistical significance.

Screening of differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) related to tumor grade

The difference in the expression of grading degree related lncRNAs and mRNAs for the different clinical tumor grades (grade 1 vs. grade 2, grade 1 vs. grade 3, grade 2 vs. grade 3) was performed. The Tukey's Honest Significant Difference test was used for multiple comparisons between different clinical tumor grades. The p < 0.05 was considered as statistical significance. The overlapped DElncRNA and DEmRNA for all clinical tumor grades were also analyzed.

Functional annotation

Gene Ontology (GO) classification and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using GeneCoDis3 (http://genecodis.cnb.csic.es/analysis). FDR <0.05 was considered to show statistically significant differences.

Protein‐protein interaction (PPI) network construction

A total of 332 DEmRNAs (after overlapping for grade 1 vs. grade 2 and grade 1 vs. grade 3) in ESCC were used to construct the PPI network by using String (https://string‐db.org/) and Cytoscape 3.5.0 (http://www.cytoscape.org/). The nodes and edges represent the proteins and interactions between the two proteins, respectively.

DEmRNA‐DElncRNA interaction analysis

The correlation between DElncRNAs and DEmRNAs was analyzed by the Spearman correlation coefficient. DElncRNA‐DEmRNA pairs with p<0.05 and | r |>0.6 were served as significant mRNA‐lncRNA co‐expression pairs. The Cytoscape 3.5.0 (http://www.cytoscape.org/) was used to build the DElncRNA‐DEmRNA co‐expression network.

qRT‐PCR Confirmation

Six genes (RP11‐25G10.2, ADAMTS9‐AS1, ZBED3‐AS1, SFN, ATP1A2, and GNA15) were selected as candidate genes. Eighteen tissue samples were obtained, including nine samples of ESCC para‐cancer, two tumor samples of grade 3 ESCC, and seven tumor samples of grade 2 ESCC. This study was approved by the ethics institute of our hospital. The signed informed consent was obtained from all the participants.

Total RNA was isolated using a TRizol kit (Tiangen, China). A Fast Quant RT Kit (Tiangen, China) was utilized to obtain the complementary DNA. Using Super Real PreMix Plus SYBR Green (Tiangen, China), quantitative real‐time PCR was generated using the LightCycler 96. Relative gene expression was analyzed by the 2−ΔΔCt method.

Validation in a published GEO dataset (GSE53625)

The expression patterns of the expression of top 10 upregulated and top 10 downregulated mRNAs and lncRNAs were validated using the GSE53625 dataset. The Gene Expression Omnibus (GEO) dataset GSE53625 dataset was were downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/). The GSE53625 dataset was published and the tissue sample consisted of 179 normal controls, 32 grade 1 ESCC patients, 98 grade 2 patients, and 49 grade 3 patients.

RESULTS

Tumor grade‐related mRNAs and lncRNAs

A total of 1864 tumor grade‐related mRNAs (846 positively related and 1018 negatively related) and 552 tumor grade‐related lncRNAs (331 positively related and 221 negatively related) were identified. The top 10 significantly grade‐related mRNAs, including CA12, FABP4, DECR1, BAIAP2, IL1RAPL2, PPARD, LAD1, TSPAN10, LDOC1, and ZNF853, were listed in Figure 1 and Table 1. The top 10 significantly grade‐related lncRNAs, including RP11‐25G10.2, RP11‐557H15.3, RP11‐521D12.5, CHKB‐AS1, RP11‐219B4.3, CH17‐335B8.4, RP11‐99 J16‐A.2, CTB‐111H14.1, ADNP‐AS1, and JHDM1D‐AS1, were shown in Figure 2 and Table 2. Besides, the heat map of top 100 tumor grade‐related mRNAs and lncRNAs was shown in Figure 3.

FIGURE 1

The box plots of the top 10 tumor grade‐related mRNAs. The x‐axis and y‐axis represent the clinical tumor grade and expression levels, respectively

TABLE 1

The top 10 significantly related mRNAs

mRNA	Mean_G1	Mean_G2	Mean_G3	PValue	Association
CA12	15.0268	14.1534	13.1209	3.69E−05	Negative
FABP4	9.6833	8.9306	6.1734	3.69E−05	Negative
DECR1	11.8339	11.5243	11.0303	7.08E−05	Negative
BAIAP2	13.7962	13.0492	12.5766	7.08E−05	Negative
IL1RAPL2	1.4241	3.0979	4.3787	9.89E−05	Positive
PPARD	13.7714	12.8994	12.6284	1.33E−04	Negative
LAD1	15.4383	14.8062	14.1925	1.33E−04	Negative
TSPAN10	8.7461	8.0707	7.1247	1.33E−04	Negative
LDOC1	8.1131	10.4949	10.9634	1.33E−04	Positive
ZNF853	7.0001	8.6198	9.3348	1.33E−04	Positive

G1 (grade 1); G2 (grade 2); G3 (grade 3).

FIGURE 2

The box plots of the top 10 tumor grade‐related lncRNAs. The x‐axis and y‐axis represent the tumor grade and expression levels, respectively

TABLE 2

The top 10 significantly related lncRNAs

lncRNA	Mean_G1	Mean_G2	Mean_G3	PValue	Association
RP11‐25G10.2	0.7837	2.0551	2.9416	1.39E−05	Positive
RP11‐557H15.3	8.3860	6.5275	5.0147	1.88E−05	Negative
RP11‐521D12.5	7.2330	6.3871	5.2905	7.08E−05	Negative
CHKB‐AS1	6.8058	6.4156	5.8739	1.33E−04	Negative
RP11‐219B4.3	1.6527	2.4740	3.4196	2.47E−04	Positive
CH17‐335B8.4	4.3742	4.0931	2.7438	3.52E−04	Negative
RP11‐99J16__A.2	0.7723	1.4853	2.0785	6.05E−04	Positive
CTB−111H14.1	1.9260	2.8209	3.0849	6.93E−04	Positive
ADNP‐AS1	6.5575	6.9455	7.1321	7.49E−04	Positive
JHDM1D‐AS1	10.1829	9.6409	8.7570	7.73E−04	Negative

G1 (grade 1); G2 (grade 2); G3 (grade 3).

FIGURE 3

The heat map of the top 100 tumor grade‐related mRNAs and lncRNAs. The diagram presents the results of a two‐way hierarchical clustering of 100 tumor grade‐related mRNAs/lncRNAs and samples. The clustering is constructed using the complete‐linkage method together with the Euclidean distance. Each row represents a tumor grade‐related mRNAs/lncRNAs and each column, a sample. The tumor grades related mRNAs/lncRNAs clustering tree is shown on the right. The color scale illustrates the relative expression level of tumor grade‐related mRNAs/lncRNAs: red, below the reference channel; green, higher than the reference

Screening of DEmRNAs and DElncRNAs related to grading degree

Based on the expression difference analysis of the tumor grade‐related lncRNAs and mRNAs in different tumor grades, there was an overlap of 455 DEmRNAs and 118 DElncRNAs between grade 1 and grade 2; 928 DEmRNAs and 254 DElncRNAs between grade 1 and grade 3; 186 DEmRNAs and 46 DElncRNAs between grade 2 and grade 3. It is worth mentioning that two overlapped DEmRNAs including SFN and IL1RAPL2 were found between grade 1, grade 2, and grade 3. It is noted that one overlapped DElncRNAs including RP11‐25G10.2 was identified between grade 1, grade 2, and grade 3. The overlap of DEmRNAs and DElncRNAs correlated with the tumor grade of ESCC between the grade 1 and grade 2, between the grade 1 and grade 3, and between the grade 2 and grade 3 were displayed in Figure 4.

FIGURE 4

Venn diagram showing the overlap of DEmRNAs and DElncRNAs for the tumor grade of ESCC between the grade 1 and grade 2, between the grade 1 and grade 3 and between the grade 2 and grade 3. Numbers represent the number of mRNAs and lncRNAs

Functional annotation of tumor grade degree related mRNAs

A total of 332 DEmRNAs (after overlapping for grade 1 vs. grade 2 and grade 1 vs. grade 3) in ESCC were used to perform the GO and KEGG enrichment analysis. According to GO enrichment analysis, signal transduction (FDR=1.57E‐09), protein phosphorylation (FDR=6.20E‐05), plasma membrane (FDR=1.75E‐13), and nucleotide binding (FDR=4.54E‐05) were significantly enriched GO terms. The KEGG pathway enrichment analysis displayed that Inositol phosphate metabolism (FDR=4.39E‐06), Gastric acid secretion (FDR=0.000134781), Phosphatidylinositol signaling system (FDR=0.00134861), Pathways in cancer (FDR=0.00852066), and Salivary secretion (FDR=0.00999689) were four significantly enriched pathways. The top 15 most significantly enriched GO and KEGG pathways of mRNAs were shown in Figure 5.

FIGURE 5

The top 15 significant enrichment GO terms and KEGG pathways of DEmRNAs. The x‐axis shows ‐log FDR and the y‐axis shows GO terms or KEGG pathways. (A) Biological processes. (B) Cellular components. (C)Molecular functions. (D) KEGG pathways

PPI network construction

A total of 332 DEmRNAs (after overlapping for grade 1 vs. grade 2 and grade 1 vs. grade 3) in ESCC were used to construct the PPI network construction. The PPI network was consisted of 153 nodes and 274 edges (Figure 6). The top 10 proteins with higher degrees, including GNAI1 (degree =18), RAP2B (degree =14), GNAZ (degree =13), SHH (degree =12), ADCY1 (degree =12), PRKAR2B (degree =12), SH3GL1 (degree =12), GNA15 (degree =11), and ARRB1 (degree =11).

FIGURE 6

The PPI network. Ellipses represent nodes and lines represent edges. Red and green represent positive correlation and negative correlation with tumor grade, respectively

The 87 DElncRNAs (after overlapping for grade 1 vs. grade 2 and grade 1 vs. grade 3) and 332 DEmRNAs (after overlapped of grade 1 vs. grade 2 and grade 1 vs. grade 3) were used to build the DEmRNA‐DElncRNA interaction analysis. A total of 456 DElncRNA‐DEmRNA co‐expression pairs including 55 lncRNAs and 175 mRNAs were identified with absolute value of the Spearman correlation coefficient |r|>0.6 and p‐value<0.05. The co‐expressed DElncRNA‐DEmRNA network was shown in the Figure 7. Among which, ADAMTS9‐AS2 (degree=95), RP11‐210 M15.2 (degree=61), RP11‐13 K12.1 (degree=39), ZBED3‐AS1 (degree=30), and RP11‐25G10.2 (degree=17) were lncRNAs with higher degrees.

FIGURE 7

The co‐expressed DElncRNAs‐DEmRNAs network. The ellipses and inverted triangles represent the DEmRNAs and DElncRNAs, respectively. Red and green represent positive correlation and negative correlation with tumor grade, respectively

Validation of qRT‐PCR.

As shown in Figure 8, RP11‐25G10.2, SFN, and GNA15 were downregulated and ADAMTS9‐AS1, ZBED3‐AS1, and ATP1A2 were upregulated in the grade 2 ESCC compared with grade 3 ESCC. Except for RP11‐25G10.2, the other five genes were consistent with our TCGA analysis.

FIGURE 8

The validation of the expression levels of selected genes. The x‐axis shows DEGs and the y‐axis shows log2 (fold change) between AD and normal controls. *p<0.05, **p<0.01

Validation in GSE53625

As presented in Figure 9 and Figure 10, the expression of eight DEmRNAs (CA12, FABP4, BAIAP2, PPARD, LAD1, TSPAN10, LDOC1, and ZNF853) and three DElncRNAs (CHKB‐AS1, RP11‐219B4.3, and JHDM1D‐AS1) was consistent with our TCGA results, generally. This validation may enhance the credibility of our results in integrated analysis. Further confirmation and experiments with larger sample sizes will be conducted in our following study.

FIGURE 9

Validation in GSE53625 the expression levels of selected mRNAs. (A) CA12. (B) FABP4. (C) BAIAP2. (D) PPARD. (E) LAD1. (F) TSPAN10. (G) LDOC1. (H) ZNF853. The x‐axis and y‐axis represent the clinical tumor grade and expression levels, respectively. *p<0.05, **p<0.01

FIGURE 10

Validation in GSE53625 the expression levels of selected lncRNAs. (A) CHKB‐AS1. (B) RP11‐219B4.3. (C) JHDM1D‐AS1. The x‐axis and y‐axis represent the clinical tumor grades and expression levels, respectively. *p<0.05, **p<0.01

DISCUSSION

ESCC has a higher incidence among blacks and Asians, accounting for 70% of esophageal cancers. Although its diagnosis and treatment have advanced recent years, ESCC ranks among the fourth leading cause of cancer‐related death. In this study, two differentially expressed mRNAs (SFN and IL1RAPL2) and one common differentially expressed lncRNAs (RP11‐25G10.2) overlapped for ESCC in grade 1, grade 2, and grade 3.

In the PPI network, GNAI1, RAP2B, GNAZ, SHH, ADCY1, PRKAR2B, SH3GL1, GNA15, and ARRB1 were top 10 proteins with higher degrees. GNAI1, G protein subunit alpha i1, is differentially expressed mRNA in pancreatic ductal adenocarcinoma, and plays a key role in the pathogenesis of pancreatic ductal adenocarcinoma. GNAI1 was associated with the tumor grade of ESCC and significantly enriched in pathways of cancer. ADCY1 is highly expressed in multidrug‐resistant in ESCC cell lines and facilitated ESCC cells resistance to most chemotherapy drugs. In this study, GNAI1 and ADCY1 were DElncRNAs related to ESCC grading. GNAI1 and ADCY1 were significantly enriched pathway of Gastric acid secretion. Therefore, we hypothesized that GNAI1 and ADCY1 may play pivotal roles in ESCC by regulating signaling pathway of Gastric acid secretion.

GNA15, a G protein subunit alpha 15, is differentially expressed mRNA in human squamous cell carcinoma and pancreatic ductal adenocarcinoma. , It has been reported that GNA15 is negatively correlated with the tumor grade of ESCC patients. In the current study, GNA15 was negatively correlated with the tumor grade of ESCC. The interaction network analysis results showed that GNA15 was co‐expressed with RP11‐210 M15.2 and ZBED3‐AS. Hence, we further hypothesized that P11‐210 M15.2 and ZBED3‐AS regulate GNA15, which may be crucial to ESCC tumor grading.

Phospholipase C beta 4 (PLCB4) is overexpressed in primary gastrointestinal stromal tumors, and it is a novel overexpressed enzyme that regulates lipid catabolism, promotes cell proliferation, and independently confers a worse prognosis. PLCB4 hotspot mutation is similar to a gain‐of‐function mutation leading to the activation of the same signaling pathway, promoting uveal melanoma tumorigenesis. In the present study, PLCB4 was positively correlated with tumor grade in patients with ESCC. The KEGG results showed that PLCB4 was a significantly enriched pathway of salivary secretion. The interaction network analysis results showed that PLCB4 was co‐expressed with RP11‐210 M15.2. Therefore, we presumed that RP11‐210 M15.2 may be involved in ESCC tumor grade by regulating the salivary secretion pathway.

Two DEmRNAs (SFN and IL1RAPL2) and one overlapped DElncRNAs (RP11‐25G10.2) overlapped for grade 1, grade 2, and grade 3. Stratifin (SFN) is originally identified as a p53‐inducible gene that is responsive to DNA‐damaging agents. SFN is a differentially expressed gene related to human lung adenocarcinoma cell proliferation, and SFN may facilitate lung tumor development and progression. SFN promoter methylation may be associated with the carcinogenesis of breast cancer and that the use of SFN promoter methylation may represent a useful blood‐based biomarker for the clinical diagnosis of breast cancer. Hepatocellular carcinoma‐secreted SFN promotes the expression of matrix metalloproteinases in cancerous surrounding cells via an aminopeptidase N dependent mechanism, and SFN has a paracrine effect on educating stromal cells in the tumor‐associated microenvironment. Here, SFN was co‐expressed with RP11‐25G10.2 and ZBED3‐AS1. Therefore, we hypothesized that both RP11‐210 M15.2 and ZBED3‐AS1 may be involved in regulation of SFN and may underlie the tumor grade.

In summary, we identified several differentially expressed mRNAs and lncRNAs associated with different tumor grades of ESCC. Interestingly, the overlapped DEmRNAs (SFN and IL1RAPL2) and DElncRNA (RP11‐25G10.2) overlapping for grade 1, grade 2, and grade 3 may play a crucial role in the regulation of tumor grade in ESCC. However, some limitations should be addressed. The number of samples used qRT‐PCR confirmation was small. More studies are needed to uncover the precise mechanisms of ESCC.

CONSENT FOR PUBLICATION

The subjects gave written informed consent for the publication of any associated data and accompanying images.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

AUTHOR CONTRIBUTIONS

Yufeng Cheng and Xin Gao contributed to the conception of the study. Qian Liu and Xue Chen contributed the materials and performed the experiment. Shaoping Chen, Jianmei Yang, and Qiang Liu performed the data analyses. Yufeng Cheng and Xin Gao contributed significantly in writing the manuscript. All authors read and approved the final manuscript.

ETHICAL APPROVAL

This study was approved by the ethical committee of the Qilu Hospital, Cheeloo College of Medicine, Shandong University.