High-throughput sequencing profile of laryngeal cancers: analysis of co-expression and competing endogenous RNA networks of circular RNAs, long non-coding RNAs, and messenger RNAs

Abstract

Background: Circular RNAs (circRNAs) and long non-coding RNAs (lncRNAs) have been recently identified as new classes of non-coding RNAs which participate in carcinogenesis and tumor progression. However, the functions of these non-coding RNAs and gene expression patterns are largely unknown.

Methods: We carried out high-throughput sequencing to analyze the differential expression of RNAs in 5 coupled laryngeal cancer (LC) and corresponding adjacent noncancerous tissues. Bioinformatics analyses were performed to predict the functions of these non-coding RNAs via co-expression, competing endogenous RNA networks and pathway enrichment analysis. The differential expression of the selected RNAs were confirmed using RT-qPCR. The CCK8, EDU, Transwell, and wound healing assays were conducted to validate the biological functions of SNHG29 in LC. Western blot assay was performed to identify the effects of SNHG29 having on the epithelial to mesenchymal transition process. Kaplan-Meier analysis was used to investigate whether the expression level of SNHG29 correlated with survival in LC patients. One-way ANOVA was used to analyze the correlation between the expression of SNHG29 and clinicopathological parameters of the included patients.

Results: Compared to normal laryngeal tissues, 31,763 non-coding RNAs were upregulated and 11,557 non-coding RNAs were downregulated in cancer tissues. SNHG29 expression was low in the LC cell lines and tissues predicting a better clinical prognosis. SNHG29 was also found to inhibit the proliferation, migration, and invasion ability of LC, exerting a suppressive role in the epithelial to mesenchymal transition process as well. SNHG29 downregulation was significantly correlated with differentiation (P=0.026), T-stage (P=0.041), lymphatic metastasis (P=0.044), and clinical stage (P=0.037). We found that the biological functions of differentially expressed transcripts included cell adhesion, biological adhesion, and migration and invasion related to adherens junction pathways.

Conclusions: Our study was the first to describe the non-coding RNA profile of LC, and suggested that dysregulated non-coding RNAs could be involved in LC tumorigenesis. SNHG29 was demonstrated to play crucial roles in inhibiting the pathogenesis and progression of LC. Our findings provide a new approach for further analyses of pathogenetic mechanisms, the detection of novel transcripts, and the identification of valuable biomarkers for this tumor.

Keywords: Differential expression analysis; biomarker identification; novel transcript detection.

Figure 1

Detailed characteristics of transcripts in the sequencing data. (A) Bar chart showing types and counts of all kinds of transcripts in the sequencing data. (B) Box plot of the transcript FPKM values between tumor and normal tissues. Ordinate represents log10 FPKM. Maximum, upper quartile, median, lower quartile, and minimum are displayed from top to bottom. (C) Heatmap showing inter-sample correlation. Correlation was assessed by Pearson’s correlation coefficient of all transcript expression levels. Blue color indicates significant differences among different samples, red color indicates slight differences.

Figure 2

Differentially expressed transcripts. (A) Bar chart showing counts of differentially expressed transcripts. Red colors, upregulated transcripts; blue colors, downregulated transcripts. (B) Volcano plots showing significantly differentially expressed transcripts between the tumor and normal groups. Vertical lines represent 2-fold changes for upregulation and downregulation. Horizontal line represents statistical significance, P<0.05. Red colors indicate differentially expressed transcripts with statistical significance. (C) Clustering analysis showing all of the differentially expressed transcripts. Heatmap is based on transcript expression values with the criteria of fold-change>2.0 and P<0.05. T: tumors; N: normal.

Figure 3

Co-expression network of coding and non-coding RNAs. The absolute value of Pearson’s correlation coefficient was limited to >0.99, the P value was limited to <0.05. Red nodes indicate upregulation; gray nodes indicate downregulation. All included transcripts were positively correlated.

Figure 4

Competing endogenous RNA network of non-coding and coding RNAs. This network was based on the correlations among circRNA/miRNAs, lncRNA/miRNAs, and miRNA/mRNAs. Red color indicates upregulation; gray color indicates downregulation.

Figure 5

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed genes. (A) GO analysis of biological process, cellular components, and molecular functions (q<0.0001). Left Y-axis indicates counts of genes; right Y-axis indicates P values. (B) and (C) Enriched bubble diagrams of the top 20 differentially expressed genes. Enrichment factor indicates the ratio between the differentially expressed genes enriched in the pathway. The scale of the bubble indicates counts of genes; the depth of bubble color indicates the P value.

Figure 6

The expression trends of selected genes in sequencing results and corresponding verification by quantitative real-time PCR (qRT-PCR). (A,B,C,D,E) The expression levels of selected RNAs were verified by qRT-PCR using 47 patient samples. T, laryngeal cancer tissues; N, normal adjacent tissues. (F) The comparison of selected genes between sequencing data and qRT-PCR. The vertical axis represents the mean value of fold change (log2 scale) of the selected RNAs. (G) The expression levels of SHKBP1, CLDN1, LINC00519, and SNHG29 were detected by qRT-PCR in TU212, LCC, and HaCat cells. Data are shown as mean ± SD, n=3.

Figure 7

The characteristics and anticancer role of SNHG29. (A) Fluorescence in situ hybridization (FISH) analysis for SNHG29 was performed in TU212 and LCC cells. Scale bar: 100 µm. (B) RT-qPCR was performed to evaluate the expression level of SNHG29 in LC cells transfected with si-RNA, si-NC, mock, or SNHG29. (C) The CCK-8 assay was conducted to measure the viability of LC cells transfected with si-RNA or SNHG29 vector. (D) The EDU assay was performed to measure the cell proliferation ability. Scale bar: 50 µm. (E) Cell invasion or migration assays were performed to validate cell motility with the indicated vectors using a transwell chamber with or without matrigel. Cells were dyed by 0.1% crystal violet. Scale bar: 100 µm. (F) The wound healing assay was performed to validate the cell motility of LC cells transfected with si-SNHG29 or SNHG29. Scale bar: 100 µm. (G) The effects of SNHG29 on the epithelial to mesenchymal transition process was measured in LC cells transfected with SNHG29 overexpressed vector using the western blot assay. (H) The expression levels of DDIT4L and PFKFB3 were evaluated using RT-qPCR in LC cells transfected with si-RNA, si-NC, mock, or overexpression vectors. (I) Kaplan-Meier curves indicated that the upregulation of SNHG29 was related to the higher survival rates. Data indicate mean ± SD, n=3; *P<0.05, **P<0.01, ***P<0.001.