Four Autophagy-Related Long Noncoding RNAs Provide Coexpression and ceRNA Mechanisms in Retinoblastoma through Bioinformatics and Experimental Evidence.

Retinoblastoma (RB) is the most common type of intraocular malignant tumor that lowers the quality of life among children worldwide. Long noncoding RNAs (lncRNAs) are reported to play a dual role in tumorigenesis and development of RB. Autophagy is also reported to be involved in RB occurrence. Although several studies of autophagy-related lncRNAs in RB have been explored before, there are still unknown potential mechanisms in RB. In the present study, we mined dataset GSE110811 from the Gene Expression Omnibus database and downloaded autophagy-related genes from the Human Autophagy Database for further bioinformatic analysis. By implementing the differential expression analysis and Pearson correlation analysis on the lncRNA expression matrix and autophagy-related genes expression matrix, we identified four autophagy-related lncRNAs (namely, N4BP2L2-IT2, SH3BP5-AS1, CDKN2B-AS1, and LINC-PINT) associated with RB. We then performed differential expression analysis on microRNA (miRNA) from dataset GSE39105 for further analyses of lncRNA-miRNA-mRNA regulatory mechanisms. With the miRNA-lncRNA module on the StarBase 3.0 website, we predicted the differentially expressed miRNAs that could target the autophagy-related lncRNAs and constructed a potential lncRNA-miRNA-mRNA regulatory network. Furthermore, the functional annotations of these target genes in regulatory networks were presented using the Cytoscape and the Metascape annotation tool. Finally, the expression pattern of the four autophagy-related lncRNAs was evaluated via qRT-PCR. In conclusion, our findings suggest that the four autophagy-related lncRNAs could be critical molecules associated with the development of RB and affect the occurrence and development of RB through the lncRNA-miRNA-mRNA regulatory network. Genes (GRP13B, IFT88, EPHA3, GABARAPL1, and EIF4EBP1) may serve as potential novel therapeutic targets and biomarkers in RB.

Introduction

Retinoblastoma (RB) is the most common primary intraocular malignancy in childhood, although it comprises a rare pediatric cancer (worldwide incidence of 1:15000 to 1:20000 live births, approximately 8000 cases per year).[1,2] Patient survival is approximately 30% globally; the low survival rate is related to a lack of early detection and effective therapeutic strategies, particularly in low-income countries.[3] The early detection and effective therapeutic strategies are critical for improving the survival rates.[4] As is known, the biallelic mutation and inactivation located in the chromosome 13q14 tumor suppressor gene Rb1 lead to the occurrence of hereditary RB; approximately 1% of RB (usually unilateral) is caused by somatic MYCN amplification, rather than the loss of RB1 function.[5,6] However, with further research, diverse genomic alterations and epigenetic modifications like microRNAs, DNA methylation, and lncRNAs are also involved in tumorigenesis, many of which mechanisms are still unclear.[7,8] Due to the biggest success of the targeted therapy in the treatment of cancer in the past few decades, various specific drugs developed for molecular targets have played a corresponding role in cancer treatment.[9] Therefore, to provide more clinical treatments, it is imperative to study new molecular targets in RB.

Long noncoding RNAs (lncRNAs) are molecules longer than 200 nucleotides, which are transcribed by RNA polymerase II.[10] LncRNAs affect the pathophysiological development of diverse systems (e.g., nervous, muscular, cardiovascular, adipose, hematopoietic, and immune systems); they also play essential roles in cellular processes such as epigenetic regulation, cell cycle regulation, cell differentiation regulation, and post-transcriptional regulation.[11,12] According to the previous studies, various lncRNAs contribute to the proliferation, apoptosis, migration, and invasion of RB.[13] For example, the lncRNA SNHG16 exerts an oncogenic role in RB by sponging miR-140-5p, thus suppressing retinoblastoma progression.[14] Additionally, the knockdown of LINC00324 decreased RB cell proliferation, colony formation, migration, and invasion while promoting apoptosis and cell cycle arrest in vitro.[15] In our research, we are committed to finding some lncRNAs unstudied in RB.

Autophagy is the process by which cellular components are degraded in the lysosome; it includes macro-autophagy, micro-autophagy, and chaperone-mediated autophagy.[16] Mutation or loss of function of key autophagy genes is associated with the incidence and progression of neuropathies, cardiovascular diseases, autoimmune diseases, and malignant tumors.[17] Although the roles of autophagy in tumors have not been fully elucidated, autophagy is demonstrated to be involved in the metabolism of various cancers and tumors, affecting the tumor microenvironment and significantly controlling their immune responses.[18,19] Some studies have shown that autophagy can affect retinoblastoma development through distinct pathways. The protein E2F1 in RB, a transcriptional regulator of autophagy, can upregulate autophagy-related genes and proteins in DNA damage-induced autophagy, thus reducing DNA damage.[20] MicroRNA MIR34A-dependent high mobility group box 1 (HMGB1) downregulation promotes oxidative damage and DNA damage by inhibiting autophagy and finally inducing apoptosis in the retinoblastoma cell.[21] Although the inhibition of autophagy can promote the apoptosis and proliferation of RB,[22,23] there are still many other unsolved mechanisms in RB.

More and more studies have confirmed that the interaction between lncRNA and autophagy affects mechanisms of tumors. LINRIS knockdown prevented the degradation of IGF2BP2 through the autophagy–lysosome pathway and attenuated the downstream effects of IGF2BP2 in colorectal cancer cells.[24] The lncRNA HAGLROS activated the mTORC1 signaling pathway, which inhibited the expression of autophagy-related genes ATG9A and ATG9B, thereby promoting excessive proliferation and maintaining the malignant phenotype of GC cells.[25] However, there are a few studies on the role and mechanism of autophagy-related lncRNAs in retinoblastoma. In this study, we identified four autophagy-related lncRNAs in RB based on data from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/)[26] and the Human Autophagy (HADb) Database (http://www.autophagy.lu/).[27] Subsequently, through the further functional enrichment analyses and the construction of lncRNA–miRNA–mRNA regulatory networks, we concluded four regulatory networks associated with RB cells, which may provide a certain theoretical basis for fundamental experimental research studies and clinical target therapies.

Materials and Methods

Data Acquisition and Processing

An independent RB gene expression profile (GSE110811) containing 28 RB samples and 3 retinal samples from healthy controls (HCs) was obtained from the GEO database using a GEO query package in R software. First, we normalized gene expression profiles using the limma package in R software to obtain normalized boxplots for subsequent research and analysis. Then, we annotated and classified those gene expression profiles to acquire a lncRNA expression matrix and an mRNA expression matrix by processing data using Bombyx_mori.GCA_000151625.1.23.gtf files. A list of autophagy-related genes was obtained by means of HADB and GSEA analyses. By integration of autophagy-related genes with the mRNA expression matrix, an expression matrix of autophagy-related genes was produced. Next, the lncRNA expression matrix and the autophagy-related gene expression matrix were separately analyzed using the limma package in R software. Differentially expressed autophagy-related genes and lncRNAs were identified using the following criteria: |log2 fold change (FC)| > 1 and P < 0.05.

GO and KEGG Analyses of Autophagy-Related Genes and Autophagy-Related lncRNAs

The cluster profiler package in R software was used to perform Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses of autophagy-related genes. A P value of <0.05 was considered statistically significant; the overall results of enrichment analyses were selected for the KEGG visualization analysis of autophagy-related genes to identify molecular pathways related to RB. Subsequently, Pearson correlation analysis was carried out using R software to analyze correlations between abnormally expressed lncRNAs and autophagy-related genes in retinoblastoma using |r| > 0.5 and P < 0.05 as the thresholds to screen out autophagy-related lncRNAs in RB. Autophagy-related lncRNA and mRNA regulatory networks were constructed, and Metascape annotation tools were used to analyze the molecular pathways related to this regulatory network in RB.

Construction of the lncRNA–miRNA–mRNA Regulatory Network and Enrichment Analysis

Based on the gene–ceRNA theory, lncRNAs may be targeted to bind miRNAs to release specific miRNAs from their target mRNAs, thereby promoting mRNA expression.[28] Therefore, we selected lncRNAs and mRNAs with positive regulatory relationships in the autophagy-related lncRNA–mRNA regulatory network. First, we used the GSE39105 dataset to analyze differences in miRNA expression between RB patients and HCs. Using |logFC| > 1 and P < 0.05 as thresholds, we screened out abnormally expressed miRNAs, visualized by heatmaps and volcano plots. Then, we predicted miRNAs that could target autophagy-related lncRNAs through the miRNA–lncRNA module on the StarBase 3.0 website (http://starbase.sysu.edu.cn/) in the Encyclopedia of RNA Interactomes.[29] We used the intersection of miRNAs targeting lncRNAs and their related mRNAs, thereby constructing a potential lncRNA–miRNA–mRNA regulatory network.

Functional Analysis

Based on the constructed ceRNA regulatory network, we performed Pearson correlation analysis using R software to analyze correlations between the regulatory network and autophagy-related genes. We identified autophagy-related coexpressed genes related to this regulatory network with the standards of |r| > 0.5 and P < 0.05. Then, we performed functional annotations of these genes using the Metascape annotation tool.

Cell Culture and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Due to the universal usability and cell line stability of the retinal pigment epithelium cell line (ARPE-19) and the human retinoblastoma cell lines (Y79), we will use these two types of cells for experiments. ARPE-19 and Y79 were purchased from ATCC (Rockville, MD, USA). The ARPE-19 cell line was maintained in DMEM (Gibco, Grand Island, NY, USA) with 10% fetal bovine serum (Gibco, Grand Island, NY, USA) and antibiotics (100 U/mL penicillin and 100 μg/mL streptomycin, Gibco, Grand Island, NY, and Scotland, UK). The Y79 cell lines were cultured using RPMI 1640 medium (Gibco, NY, USA) containing 100 U/L penicillin/streptomycin and 10% fetal bovine serum (FBS). All cells were cultured in an incubator at 37 °C in an atmosphere containing 5% CO2.

Total RNA was extracted using TRIzol (Invitrogen). PrimeScript RT Master Mix (Yeasen, Shanghai, China) was used for cDNA synthesis. The PCR primers are listed as follows: H-N4BP2L2-IT2, forward: 5′-TTG AAT GCC TTC ACC TGT GC-3′; H-N4BP2L2-IT2, reverse: 5′-CAG ACT CAG CAA AGA AGG CG-3′; H-SH3BP5-AS1, forward: 5′-ATC AGG CTC AGG TTT GCT CC-3′; H-SH3BP5-AS1, reverse: 5′-AGG CTA GCA GGG TAG TCT TCA-3′; H-CDKN2B-AS1, forward: 5′-GCA GAA ACC ACA TCC CTT GG-3′; reverse: 5′-TAG TGC GTT AGG CAT CTG TGT-3′; H-LINC-PINT, forward: 5′-ATG AGG TAG GAG GCT CAG CA-3′; reverse: 5′-CAA GAG GTA GCT GGC GGA AA-3′; GAPDH, forward: 5′-GAG AAG GCT GGG GCT CAT TT-3′; reverse: 5′-TAA GCA GTT GGT GGT GCA GG-3′. The 2–ΔΔ method was conducted to calculate the lncRNA expression. The Student t-test was used to compare the expression level of each lncRNA between different groups.

Results

Differentially Expressed Autophagy-Related Genes and lncRNAs in RB

The results of standardized preprocessing of gene expression profiles are presented in the form of normalized boxplots (Figure A,B: before and after standardization, respectively). Based on these data, 17 significantly differentially expressed lncRNAs (Figure C,D) and 10 significantly differentially expressed autophagy-related genes (Figure E,F) were identified between RB and HC groups by cluster analysis, visualized by heatmaps and bar graphs, respectively. The figure shows 7 significantly upregulated lncRNAs and 10 significantly downregulated lncRNAs, as well as 1 significantly upregulated autophagy-related gene and 9 significantly downregulated autophagy-related genes.

Figure 1

Differentially expressed autophagy-related genes and lncRNAs. (A) Boxplot of RB gene expression profiles (GSE110811) before standardization. (B) Boxplot of RB gene expression profiles (GSE110811) after standardization. (C, E) Heatmaps of lncRNAs and autophagy-related genes (|log2 fold change| > 1; P < 0.05). (D, F) Bar graphs of lncRNAs and autophagy-related genes (|log2 fold change| > 1; P < 0.05).

Functional Enrichment Analysis

Ten autophagy-related differential genes were assessed by cluster analysis. The GO analysis of the 10 autophagy-related genes showed biological functions (BPs) related to response to starvation, positive regulation of endothelial migration, endothelial autophagy, etc. In addition, the molecular functions (MFs) of those genes included integrin binding, transmembrane receptor protein tyrosine kinase activity, and heat shock protein binding (Figure A,B). The KEGG enrichment analysis of autophagy-related genes showed enrichment in EGFR tyrosine kinase inhibitor resistance, PI3K-Akt signaling pathway, MAPK signaling pathway, etc. (Figure C,D).

Figure 2

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms of differentially expressed genes (DEGs). (A, B) GO analysis based on DEGs. BP, biological process; MF, molecular function. (C, D) KEGG analysis based on DEGs.

Identification Autophagy-Related lncRNAs and Functional Analysis

To identify autophagy-related lncRNAs, Pearson correlation analysis was conducted based on autophagy-related genes and lncRNAs. A total of 17 autophagy-related lncRNAs were identified. The autophagy-related lncRNA–mRNA regulatory network is shown in Figure A. The Metascape annotation tool analysis of the regulatory network function revealed an association with response to starvation, autophagy, and positive regulation of endothelial cell migration (Figure B).

Figure 3

(A) lncRNA–target network showing autophagy-related genes targeted by lncRNAs. (B) Functional diagram of the lncRNA–mRNA regulatory network. Red: response to starvation; blue: autophagy; green: positive regulation of endothelial cell migration; purple: regulation of plasma membrane-bound cell project.

Construction of the lncRNA–miRNA–mRNA Regulatory Network

In the StarBase 3.0 website, differentially expressed miRNAs (DEMs) were identified based on GSE39105 between RB and HC samples, including 121 upregulated and 67 downregulated DEMs (Figure A,B). The upregulated miRNAs intersected with miRNAs targeting downregulated lncRNAs and miRNAs targeting downregulated mRNAs; the downregulated miRNAs intersected with miRNAs targeting upregulated lncRNAs and miRNAs targeting upregulated mRNAs. Seven significantly upregulated miRNAs were identified (Figure C): hsa-miR-125a-5p, hsa-miR-125b-5p, hsa-miR-181a-5p, hsa-miR-181c-5p, hsa-miR-194-5p, hsa-miR-26b-5p, and hsa-miR-485-3p. Three significantly downregulated miRNAs were identified (Figure D): hsa-miR-15b-5p, hsa-miR-16b-5p, and hsa-miR-195-5p. According to their targeting mechanisms, a lncRNA–miRNA–mRNA regulatory network related to autophagy in RB is constructed (Figure E).

Figure 4

Differentially expressed miRNAs (DEMs) and construction of (A) a heatmap of DEMs (|log2 fold change| > 1; P < 0.05) and (B) a volcano plot of DEMs in GSE39105. (C, D) Bar graphs of DEMs (|log2 fold change| > 1; P < 0.05). (E) ceRNA network of the lncRNA–miRNA–mRNA network.

Functional Analyses of the ceRNA Network

Following completion of the ceRNA network of the lncRNA–miRNA–mRNA network, Pearson correlation analysis was performed between autophagy-related genes and ceRNAs, with the following criteria: |log2 FC| > 1 and P < 0.05. Five coexpressed genes were explored: GPR137B, IFT88, EPHA3, GABARAPL1, and EIF4EBP1. Functional annotation of GPR137B and its coexpressed genes suggested that LINC-PINT in RB may control the expression of GPR137B through hsa-miR-26b-5p, subsequently regulating the MAPK signaling pathway, invasion, autophagy, etc. (Figure A). Functional annotation of IFT88 and its coexpressed genes suggested that N4BP2L2-IT2 in RB may control the expression of IFT88, which in turn regulates membrane transport, Ras protein signal transduction, and cell polarity establishment, by targeting hsa-miR-125a-5p, hsa-miR-125b-5p, hsa-miR-181a-5p, and hsa-miR-181c-5p (Figure B).

Figure 5

Genes of coexpression network and their enrichment analyses. (A) Functional annotation of GPR137B and its coexpressed genes; (B) functional annotation of IFT88 and its coexpressed genes; (C) functional annotation of EPHA3 and its coexpressed genes; (D) functional annotation of GABARAPL1 and its coexpressed genes; (E) functional annotation of EIF4EBP1 and its coexpressed genes.

Functional annotation of EPHA3 and its coexpressed genes suggested that SH3BP5-AS1 in RB may control the expression of EPHA3 through hsa-miR-194-5p, thereby regulating the cell cycle, cell differentiation, and retinoblastoma-related genes (Figure C). Functional annotation of GABARAPL1 and its coexpressed genes suggested that SH3BP5-AS1 in retinoblastoma may control the expression of GABARAPL1 through hsa-miR-485-3p, thereby regulating cell differentiation, cell cycle, and protein autophosphorylation (Figure D). Functional annotation of EIF4EBP1 and its coexpressed genes suggested that CDKN2B-AS1 in RB may control the expression of EIF4EBP1 through hsa-miR-15b-5p, hsa-miR-195-5p, and hsa-miR-16-5p, thereby regulating cell cycle, cell differentiation, DNA repair, etc. (Figure E).

Four Autophagy-Related lncRNAs in Y79 Cell Lines

We evaluated the relative expressions of four lncRNAs (N4BP2L2-IT2, SH3BP5-AS1, CDKN2B-AS1, and LINC-PINT) in Y79 cell lines and ARPE-19 cell lines. The results revealed that the mRNA level of N4BP2L2-IT2, SH3BP5-AS1, and CDKN2B-AS1 are significantly up-regulated in Y79 cell lines than ARPE-19 cell lines (Figure A–C). LINC-PINT has a lower expression level in Y79 cell lines than in AREP-19 cell lines (Figure D).

Figure 6

qPCR result of four lncRNAs in Y79 and ARPE-19 cell lines: (A) qPCR result of N4BP2L2-IT2; (B) qPCR result of SH3BP5-AS1; (C) qPCR result of CDKN2B-AS1; (D) qPCR result of LINC-PINT.

Discussion

RB is the prototype genetic cancer, which is due to germline mutations in the RB1 gene.[30] But with changes in complex epigenetic and genetic events, it is necessary to explore other regulatory mechanisms in RB. Autophagy is a complex process of capturing and degrading damaged proteins and senescent or malfunctioning organelles.[31] Autophagy in tumor cells and the host could sustain essential metabolic functions of tumor cells by impairing DNA, activating transcription programs, and influencing protein synthesis.[42] Substantial evidence suggests that the autophagy process could be regulated by RB/E2F1-dependent transcriptional activation, strengthening the interaction between autophagy, apoptosis, and aging in RB.[32] LncRNAs can influence chromatin modification and regulate gene and genome activities of distinct tumors at various levels.[33] They can regulate many proteins required for autophagy. Here, we systematically studied the associations between autophagy-related lncRNAs and RB through bioinformatic analysis. We aimed to identify significant differentially expressed biomarkers and analyzed the network regulation they are involved in, which could provide a certain degree of diagnostic value and guide clinical treatment.

In this study, we first mined GEO and HADb for data concerning RB samples, revealing 17 differentially expressed lncRNAs and 10 differentially expressed autophagy-related genes (NFE2L2, GABARAPL1, FGE1, EPHA1, EDIL3, KDR, GPR137B, IFT88, DAPL1, and EIF4EBP1). Subsequent GO and KEGG analyses of the 10 autophagy-related differential genes showed their biological processes, molecular functions, and pathways. Analysis of these genes will facilitate future identification of pathways between lncRNAs and target genes. Subsequently, we analyzed differentially expressed lncRNAs and autophagy-related genes through Pearson correlation analysis, which revealed four autophagy-related lncRNAs (LINC-PINT, N4BP2L2-IT2, SH3BP5-AS1, and CDKN2B-AS1) that participated in an autophagy-related lncRNA–miRNA–mRNA coexpression network. Metascape annotations were used to analyze the pathways and mechanisms by which each coexpression network pathway contributed to the occurrence and development of RB. To our knowledge, most of the network nodes were identified for the first time in this study.

The four autophagy-related lncRNAs regulate different network mechanisms. LINC-PINT was found to target gene GRP137BT through has-miR-26b-5p. It affects the occurrence and development of different cancers, including glioblastoma, lung cancer, thyroid cancer, esophageal cancer, and gastric cancer.[34−38] In addition, LINC-PINT inhibits miR-523-3p to reverse malignant phenotypes, thus inhibiting retinoblastoma progression by upregulation of Dickkopf-1 (DKK1).[39] Similarly, CDKN2B-AS1 is found to target gene EIF4EBP1 through has-miR-15b-5p, has-miR-195-5p, and has-miR-16-5p. CDKN2B-AS1 sponges let-7c-5p, thus promoting nucleosome assembly protein 1-like 1 (NAP1L1) expression and activating PI3K/AKT/mTOR signaling in HCC cells, which inhibits the growth and metastasis of human hepatocellular carcinoma.[40] CDKN2B-AS1 regulates tumor progression and metastasis of renal cell carcinoma by direct interaction with miR-141.[41] Although CDKN2B-AS1 has been studied in other tumors, it has not yet been investigated in RB. The functions of N4BP2L2-IT2 and SH3BP5-AS1 are yet unclear. In this study, we found that N4BP2L2-IT2 could target gene IFT88 by regulating has-miR-125a-5p, has-miR-181a-5p, has-miR-181c-5p, and has-miR-125b-5p; SH3BP5-AS1 could target two genes (EPHA3 and GABARAPL1) by regulating has-miR-194-5p and has-miR-485-5p, respectively. These coexpression network pathways provide crucial points for follow-up research studies. A series of studies based on the obtained genes and RNAs (lncRNAs and microRNAs), such as targeted gene therapy for tumor epiregulation and control of protein expression from the research of validated RNA modification, can be conducted.[43] Overall, we identified four autophagy-related lncRNAs, which may provide new insights into finding tumor biomarkers and therapeutic targets. They will also be useful in future in vitro experiments (e.g., western blotting, polymerase chain reaction analyses, etc.).

Several limitations still exist in our study. First, our study only verified the expression patterns of four autophagy-related lncRNAs in cell experiments and lacked cell or animal experiments to verify biological functions and the whole regulatory networks. Second, due to the lack of follow-up information of RB patients, prognostic analysis of autophagy-related lncRNAs was not conducted in our study. Thus, we cannot evaluate the clinical value and significance of these lncRNAs in RB objectively. Third, the disparity in sample size between RB and the normal control may bring potential deviations to the results of analyses.