Using RNA sequencing to identify a putative lncRNA-associated ceRNA network in laryngeal squamous cell carcinoma

Abstract

Accumulating evidence indicates that lncRNAs can interact with miRNAs to regulate target mRNAs through competitive interactions. However, this mechanism remains largely unexplored in laryngeal squamous cell carcinoma (LSCC). In this study, transcriptome-wide RNA sequencing was performed on 3 pairs of LSCC tissues and adjacent normal tissues to investigate the expression profiles of lncRNAs, miRNAs and mRNAs, with differential expression of 171 lncRNAs, 36 miRNAs and 1709 mRNAs detected. Seven lncRNAs, eight mRNAs and three miRNAs were identified to be dysregulated in patients' tissues by using qRT-PCR. GO and KEGG pathway enrichment analyses were performed to elucidate the potential functions of these differentially expressed genes in LSCC. Subsequently, a ceRNA (lncRNA-miRNA-mRNA) network including 4631 ceRNA pairs was constructed based on predicted miRNAs shared by lncRNAs and mRNAs. Cis- and transregulatory lncRNAs were analysed by bioinformatics-based methods. Importantly, mRNA-related ceRNA networks (mRCNs) were further obtained based on potential cancer-related coding genes. Coexpression between lncRNAs and downstream mRNAs was used as a criterion for the validation of mRCNs, with the ZNF561-AS1-miR217-WNT5A and SATB1-AS1-miR1299-SAV1/CCNG2/SH3 KBP1/JADE1/HIPK2 ceRNA regulatory interactions determined, followed by experimental validation after siRNA transfection. Moreover, ceRNA activity analysis revealed that different activities of ceRNA modules existing in specific pathological environments may contribute to the tumorigenesis of LSCC. Consistently, both downregulated SATB1-AS1 and ZNF561-AS1 significantly promoted laryngeal cancer cell migration and invasion, indicating their important roles in LSCC via a ceRNA regulatory mechanism. Taken together, the results of this investigation uncovered and systemically characterized a lncRNA-related ceRNA regulatory network that may be valuable for the diagnosis and treatment of LSCC.

Keywords: CeRNA; RNA sequencing; laryngeal carcinoma; lncRNA; noncoding RNA.

Figure 1.

Heat maps of differentially expressed lncRNA (A), mRNA (B) and miRNA (C) profiles. The heat maps describe the most significantly changed lncRNAs (top 40), mRNAs (top 40) and miRNAs (top 40) in 3 pairs of LSCC tissues compared with adjacent nontumor tissues, with fold change ≥1.5 for lncRNAs, fold change ≥2.0 for mRNAs and miRNAs, p < 0.05 and FDR < 0.05. The expression values are depicted in accordance with the colour scale. The intensity increases from green to red. LSCC: laryngeal squamous cell carcinoma; T: tumour tissues. N: adjacent nontumor tissues.

Figure 2.

Cis-regulation by lncRNAs of nearby coding genes. The distances between lncRNAs and their cis-regulated coding genes within the regions 10 kb upstream and downstream are described. The left vertical axis shows the sequence names of the coding genes, and the right vertical axis exhibits the gene symbols of the sequences. The middle vertical axis presents the cis-regulatory lncRNAs.

Figure 3.

Validation of differentially expressed lncRNAs (A), mRNAs (B), miRNAs (C) by quantitative RT-PCR. The relative expression levels of seven lncRNAs (A), eight mRNAs (B) and three miRNAs (C) are shown by comparing tumour tissues with normal tissues via qRT-PCR. The data are presented as the mean ±SEM, n = 36. *p < 0.05, **p < 0.01, ***p < 0.001..

Figure 4.

ceRNA networks in laryngeal squamous cell carcinoma. The ceRNA network was constructed based on the shared MREs of lncRNAs and mRNAs predicted by the Mireap, miRanda and TargetScan software programs, followed by determination using a hypergeometric cumulative distribution function test. This figure shows three networks containing miR-1299 ceRNA networks (A), miR-592 ceRNA networks (B) and miR-5187-5p ceRNA networks (C), which are hub miRNAs in the global ceRNA networks.

Figure 5.

KEGG pathway analysis of ceRNA networks. The KEGG pathway analysis results corresponding to miR-1299 ceRNA networks (A), miR-592 ceRNA networks (B) and miR-5187-5p ceRNA networks (C), based on coding genes in the networks, are presented.

Figure 6.

Experimentally validated mRNA-related ceRNA triplets in LSCC. Two sets of mRNA-related ceRNA networks were experimentally validated by determining the coexpression coefficients of lncRNAs and mRNAs using qRT-PCR, including lncRNA SATB1-AS1-related ceRNA clusters, which target SAV1, CCNG2, SH3KBP1, JADE1 and HIPK2 by interacting with miR-1299 (A) and the ZNF561-AS1-miR217-WNT5A triplet (B).

Figure 7.

Coexpression validation of ceRNA triplets after siRNA transfection in cell lines. SATB1-AS1-miR-1299-SAV1-JADE1/HIPK2/SH3KBP1/SAV1/CCNG2 was validated by qRT-PCR measurement after siRNA- SATB1-AS1 transfection in the SNU899 and SNU1066 LSCC cell lines (A-B). ZNF561-AS1-miR217-WNT5A was validated by qRT-PCR after siRNA-ZNF561-AS1 transfection in the SNU899 and SNU1066 LSCC cell lines (C-D).

Figure 8.

Characterization of different ceRNA patterns in LSCC tissue and nontumor tissue. Putative mRNA-related ceRNA triplets were further validated by calculating the correlation coefficients of lncRNAs and their target coding genes. SATB1-AS1-miR1299-SAV, SATB1-AS1-miR1299-SH3KBP1, SATB1-AS1-miR1299-JADE1, SATB1-AS1-miR1299-CCNG2, SATB1-AS1-miR1299-HIPK2 and ZNF561-AS1-miR217-WNT5A were observed to be ceRNAs in LSCC tissue (A-F). The SATB1-AS1-miR1299-SAV and SATB1-AS1-miR1299-SH3KBP1 triplets were also observed to have ceRNA abilities in nontumor cancer (G-H). The SATB1-AS1-miR1299-CCNG2 triplet had a relatively low level of competitive activity in noncancer tissues (I); SATB1-AS1-miR1299-JADE1, SATB1-AS1-miR1299-HIPK2 and ZNF561-AS1-miR217-WNT5A had no competitive activity in noncancer tissues (J-L). The ceRNA pair (lncRNA-miRNA-mRNA) with a correlation coefficient ≥ 0.5 determined by qRT-PCR was defined as the final ceRNA triplet.

Figure 9.

Downregulated SATB1-AS1 or ZNF561-AS1 contributes to the EMT process in LSCC.Wound healing assay was performed in SNU1066 cells transfected with siRNA-NC (siRNA-negative control), siRNA-SATB1-AS1 or siRNA-ZNF561-AS1 (A). The migration rate in the wound healing assay (B). Transwell migration assay (C) and the corresponding migrated cell numbers (E) after siRNA transfection. Transwell invasion assay (D) and the corresponding invaded cell numbers (F) after SNU1066 cells were transfected with siRNA-NC, siRNA-SATB1-AS1 or siRNA-ZNF561-AS1.