This article is mentioned in:

Introduction

Retinoblastoma (RB) typically occurs in pediatric patients under the age of 6 years and is the most common type of childhood intraocular malignancy (1). RB is an uncommon pediatric cancer, with a global incidence of 1 in every 15,000-18,000 live births (2,3). In recent years, with earlier diagnosis and the development of multimodal treatments, the survival rates of patients with RB have markedly improved (4). However, numerous survivors of RB experience blindness or eye loss. Furthermore, individuals with hereditary RB continue to have an increased risk of developing metastases, trilateral RB and secondary malignant neoplasms, which leads to high mortality rates (2,5). Therefore, it is essential to identify novel specific markers for RB, elucidate its molecular mechanisms and develop better, targeted therapeutic strategies for RB.

Long noncoding RNAs (lncRNAs) are a class of mRNA-like transcripts that are longer than 200 nucleotides and lack protein-coding capability (6). In the past decade, lncRNAs have been implicated in the pathogenesis of a variety of diseases and biological processes, such as metabolic disorders, tumorigenesis, cellular differentiation, X-chromosome imprinting and even metastasis (7-9). Furthermore, lncRNAs are widely expressed in various cancer tissues and may be exploited as diagnostic or prognostic indicators for various cancer types, including RB (10-12). However, the full lncRNA expression profile in RB has not been determined and the detailed lncRNA regulatory mechanisms in RB have remained elusive. The present study aimed to explore the lncRNA expression profile in RB and identify mechanistically relevant lncRNAs by bioinformatics analysis of the Gene Expression Omnibus (GEO) dataset GSE111168, which may allow for the identification of novel biomarkers for the diagnosis and treatment of RB.

Materials and methods

Microarray data

The gene expression dataset GSE111168 was downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo) (13). The dataset includes the mRNA expression profiles of three RB samples and three adjacent normal tissue samples determined using the GeneChip Human Genome U133 Plus 2.0 platform (Affymetrix; Thermo Fisher Scientific, Inc.). The annotation files were downloaded from the platform.

Data preprocessing and determination of differential expression

Based on the annotation information, the probe levels were converted into lncRNA or mRNA expression values and preprocessed using log2 transformation and Z-score normalization in R software 3.6 (R Core Team) for each study. Thereafter, the differentially expressed RNAs in the RB and healthy retina samples were selected using the linear models for microarray analysis package in R (Affymetrix; Thermo Fisher Scientific, Inc.; http://www.affymetrix.com/analysis/) (14). The cutoff values for the differentially expressed lncRNAs (DELs) and differentially expressed mRNAs (DEMs) were P<0.01 and |log2(fold-change)|≥2.

Functional enrichment analysis

To predict potential biological functions of DEL-targeted and other DEMs in RB, Gene Ontology (GO) analysis (www.geneontology.org) was performed using GO Association tools (www.github.com/tanghaibao/GOatools). Functional annotations were performed and terms in the categories biological process, cellular component or molecular function were determined. To predict the signaling pathways in which DEL-targeted and other DEMs participated, the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology-Based Annotation System was used and KEGG (www.genome.ad.jp/kegg) pathway enrichment analysis was performed. The GO and KEGG items with P<0.05 were considered statistically significant.

DEL-DEM coexpression network analysis

Weighted correlation network analysis (www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA) was used to construct the coexpression network of DELs and DEMs. In brief, Pearson's correlation coefficients (PCCs) for comparisons between random single DELs and DEMs were calculated. The DEL-DEM pairs were filtered for network construction using PCCs; PCCs ≥0.99 were considered meaningful values. The b parameter was simulated to weigh the network and a meaningful biological coexpression network was established. Cytoscape (http://cytoscape.org/) software was used to visualize the network.

Cell culture

The human retinal pigment epithelial cell line ARPE19 (The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences) was cultured in Dulbecco's modified Eagle's medium/F12 medium (Gibco; Thermo Fisher Scientific, Inc.). The RB cell lines Y79, SO-RB50 (The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences) and WERI-RB1 (Zhejiang Ruyao Biotechnology Co., Ltd.) were cultured in RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc.). The media were supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin/streptomycin and cells were incubated at 37˚C in a humidified atmosphere with 5% CO2.

Reverse transcription-quantitative (RT-q)PCR

Total RNA was extracted from cells using the TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.). Complementary DNA was generated by RT of total RNA using the PrimeScript RT reagent kit (Takara Biotechnology, Co. Ltd) according to the manufacturer's protocol. qPCR was performed using SYBR Premix Ex Taq II (Takara Biotechnology, Co. Ltd) on a 7500 real-time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). The thermocycling conditions were as follows: 45 cycles of 95˚C for 7 sec, 57˚C for 10 sec and 72˚C for 15 sec. Relative expression levels of DELs were calculated using the 2-ΔΔCq method (15) and GAPDH was used as an internal control. The PCR primers used in the present study are listed in Table SI.

Statistical analysis

RT-qPCR experiments were performed in triplicate and the data are presented as the mean ± standard deviation. The differences between two groups were analyzed using one-way ANOVAs and Scheffe tests were performed for the post-hoc tests. P<0.05 was considered to indicate statistical significance.

Results

DELs and DEMs in RB

Simultaneous genome-wide analyses of lncRNA and mRNA expression profiles were performed and 3,915 DELs and 3,715 DEMs were identified between the RB and adjacent normal tissue samples. Among the DELs, 1,774 were upregulated and 2,141 were downregulated in RB tissues compared to those in healthy tissues (Fig. 1A and B). Among the DEMs, 1,492 were upregulated and 2,223 were downregulated in RB tissues compared to those in healthy tissues (Fig. 2A and B). The top 10 up- and downregulated DELs are listed in Table I. TCONS_00008015 was the most significantly upregulated DEL and lnc-LRRC59-3:3 was the most significantly downregulated DEL.

Figure 1 Expression of DELs in retinoblastoma samples and adjacent normal tissue samples. (A) Clustered DEL heatmap. Rows represent specific DELs and columns represent individual samples (red, upregulated; blue, downregulated). (B) Volcano plot of DELs. The x-axis presents log2(fold-change) and the y-axis -log10(P-value) (red, upregulated; green, downregulated). DEL, differentially expressed long noncoding RNA.

Figure 2 DEM expression in retinoblastoma and adjacent normal tissue samples. (A) Clustered DEM heatmap. Rows represent specific DEMs and columns represent individual samples (red, upregulated; blue, downregulated). (B) Volcano plot of DEMs. The x-axis presents log2(fold-change) and the y-axis -log10(P-value) (red, upregulated; green, downregulated). DEM, differentially expressed mRNA.

Table I Top 10 most up and downregulated mapped differentially expressed lncRNAs and mRNAs.

Table I

Top 10 most up and downregulated mapped differentially expressed lncRNAs and mRNAs.

ID	Chromosome	Position	log2FC	P-value
TCONS_00008015	12	116405401-116463531	Inf	3.61x10-27
lnc-DAZ1-161:1	Y	56675831:56678602	Inf	3.07x10-23
lnc-BICRA-12:1	19	47840709:47843109	Inf	1.44x10-21
miR124-2HG:15	8	64379148:64382668	Inf	5.20x10-19
lnc-NEPRO-36:2	3	114314500:114329714	Inf	8.59x10-18
TCONS_00024597	3	32716474-32835950	Inf	2.75x10-17
lnc-OTUD7A-24:1	15	32765545:32767369	Inf	2.41x10-16
lnc-HDAC7-21:1	12	48688463:48692507	Inf	1.48x10-15
lnc-CCDC39-11:1	3	179900672:180037035	Inf	2.15x10-15
Lnc-CCDC85C-26:2	14	100361702:100375473	Inf	2.33x10-15

[i] FC, fold-change; lnc/lncRNA, long noncoding RNA; Inf, infinite.

Functional and pathway enrichment analyses

To explore potential biological functions and signaling pathways associated with the DEL-targeted DEMs, GO and KEGG enrichment analyses were performed. The results indicated that the DEL-targeted DEMs were highly enriched by genes associated with hexose transport (GO:0008645, biological process), muscle tissue morphogenesis (GO:0060415, biological process), the stereocilium membrane (GO:0060171, cellular component), endothelin B receptor binding (GO:0031708, molecular function) and γ-filamin/ABP-L, α-actinin and telethonin binding protein of the Z-disc binding (GO:0051373, molecular function) (Fig. 3). The KEGG pathway analysis suggested that the DEL-targeted DEMs were most enriched in pathways such as the PI3K/AKT, Hippo and cancer pathways, as well as extracellular matrix (ECM)-receptor interaction and neuroactive ligand-receptor interaction pathways (Fig. 4).

Figure 3 Functional enrichment analysis of DEL-targeted DEMs. Enriched GO terms for DEL-targeted DEMs in retinoblastoma. The x-axis presents the term ID and the y-axis the negative logarithm of the P-value. DEL, differentially expressed long noncoding RNA; DEM, differentially expressed mRNA; GO, Gene Ontology.

Figure 4 Pathway enrichment analysis of DEL-targeted DEMs. Enriched KEGG signaling pathways of DEL-targeted DEMs in retinoblastoma. DEL, differentially expressed long noncoding RNA; DEM, differentially expressed mRNA; KEGG, Kyoto Encyclopedia of Genes and Genomes; hsa, Homo sapiens.

DEL-DEM coexpression network analysis

Next, a DEL-DEM coexpression network was constructed to investigate the potential roles of lncRNAs in RB. The coexpression network consisted of 152 DELs and 285 DEMs and comprised 671 nodes (Table SII). The top three DEMs, namely, EIF1AY, GSTM1 and NLRP11, formed coexpression modules with 33, 50 and 41 DELs, respectively (Fig. 5A-C). The top three DELs, namely, lnc-DAZ1-161, lnc-HDAC7-21 and lnc-OR52A1-55, formed coexpression modules with 181, 156 and 210 DEMs, respectively (Fig. 5D-F). The expression of each DEL correlated with that of several DEMs and the expression of each DEM correlated with that of multiple DELs.

Figure 5 DEL-DEM coexpression network. Rectangular and elliptical nodes represent DELs and DEMs, respectively. (A-C) Coexpression network of the top three DEMs: EIF1AY, GSTM1 and NLRP11. (D-F) Coexpression network of the top three DELs: lnc-DAZ1-161, lnc-HDAC7-21 and lnc-OR52A1-55. DEL, differentially expressed long noncoding RNA; DEM, differentially expressed mRNA.

RT-qPCR validation of DELs

To validate the accuracy of the microarray data from the GEO database, RT-qPCR was used to evaluate the expression levels of lnc-DAZ1-161, lnc-HDAC7-21 and lnc-OR52A1-55 in human retinal pigment epithelial cells and RB cells. Among these lncRNAs, lnc-DAZ1-161 and lnc-HDAC7-21 were the top upregulated DELs and lnc-OR52A1-55 was a downregulated DEL, as obtained through bioinformatics analysis. Both lnc-DAZ1-161 and lnc-HDAC7-21 were significantly overexpressed in RB cells, which was consistent with their expression pattern in the database (Fig. 6A and B). However, lnc-OR52A1-55 was also significantly overexpressed in RB cells, which was inconsistent with its expression pattern in the database (Fig. 6C).

Figure 6 Validation of differentially expressed lncRNAs by reverse transcription-quantitative PCR in tumor cell lines compared with normal cell lines. Relative expression of (A) lnc-DAZ1-161, (B) lnc-HDAC7-21 and (C) lnc-OR52A1-55 in human retinal pigment epithelial cells and retinoblastoma cells. *P<0.05; **P<0.01; ***P<0.001 vs. APRE19. lnc/lncRNA, long noncoding RNA.

Discussion

RB is the leading global pediatrics eye malignancy based on the mortality rate of patients (16); therefore, its underlying molecular mechanisms have been widely investigated and the treatment of RB has evolved significantly in the past half-century (4). However, as a rare childhood cancer, the pathogenesis of RB has remained largely elusive. Recent studies have demonstrated that numerous lncRNAs are dysregulated in RB and drive its development and progression. Hu et al (7) reported that the lncRNA XIST promoted RB progression, in part by modulating the miR-124/STAT3 axis. Shang (11) indicated that the lncRNA THOR was upregulated in RB, acting as an oncogene by enhancing the expression of the MYC mRNA and IGF2BP1 protein. Yang and Peng (17) demonstrated that ANRIL depletion suppressed RB progression by activating the ATM/E2F1 signaling pathway. However, comprehensive analyses of lncRNA and related mRNA expression profiles, as well as potential biological functions of these RNAs in RB, have not been previously reported, to the best of our knowledge. In the present study, a comprehensive analysis of GEO datasets was performed to identify lncRNAs and mRNAs that are potentially involved in the pathogenesis of RB. A total of 1,774 upregulated DELs, 2,141 downregulated DELs, 1,492 upregulated DEMs and 2,223 downregulated DEMs were identified between the RB and para-tumor samples, which comprise the lncRNA/mRNA profile of RB.

lncRNAs mainly participate in various biological processes by regulating their target mRNAs. To gain insight into potential functions of the DELs in RB, GO and KEGG analyses of the pathways that were most highly represented among the differentially regulated RNAs were performed. The results indicated that the highly differentially expressed RNAs were commonly associated with PI3K/AKT signaling, Hippo signaling, cancer pathways, ECM-receptor interactions and neuroactive ligand-receptor interactions. The PI3K/AKT signaling pathway has a central role in regulating apoptosis, cell proliferation, metabolism and migration, as well as in maintaining the biological characteristics of malignant cells, affecting carcinogenesis and tumor progression (18). The PI3K/AKT signaling pathway also has an important role in RB (19,20). The Hippo signaling pathway is an evolutionarily conserved network that promotes uncontrolled cell proliferation, impairs differentiation and is associated with cancer (14). Furthermore, biologically active ECM fragments are able to regulate tissue injury, remodeling and cell growth and induce inflammatory responses, which are characteristic of carcinogenesis; thus ECM-receptor interactions may have a role in cancer progression (21-23). The present results suggested that the potential functions of the identified DELs were closely related to the development and pathogenesis of RB.

The lncRNA and mRNA coexpression network consisted of 671 nodes among 152 DELs and 285 DEMs. Of note, the top three DELs, lnc-DAZ1-161, lnc-HDAC7-21 and lnc-OR52A1-55, have not been previously studied. The network analysis indicated that lnc-DAZ1-161 formed a coexpression module with 181 DEMs, one of which, EIF1AY, was previously reported to be upregulated in dilated ischemic cardiomyopathy and during neural differentiation (24,25). However, the association of EIF1AY with RB has not been previously reported. Thus, the possibility that lnc-DAZ1-161 has a role in RB development by regulating EIF1AY expression requires further investigation. Furthermore, lnc-HDAC7-21 formed a coexpression module with 156 DEMs, including NLRP11. NLRP11, a primate-specific member of the NOD-like receptor family, is highly expressed in the testis, ovary, liver and immune cells (26). Upregulation of NLRP11 has pivotal roles in attenuating Toll-like receptor signaling and preventing the dysregulation of inflammatory responses (27). The possibility that lnc-HDAC7-21 regulates NLRP11 expression during RB development requires further investigation. Furthermore, lnc-OR52A1-55 formed a coexpression module with 210 DEMs. The null genotype of GSTM1, one of the coexpressed mRNAs, alters enzymatic activity and exerts indirect effects on the development of various cancer types (28,29). In addition, the GSTM1 null genotype appears to be associated with the risk of glaucoma (30). Further research is required to elucidate the potential relationship between lnc-OR52A1-55 and GSTM1 in RB.

Furthermore, RT-qPCR analysis was performed to compare expression levels of the three representative lncRNAs in human retinal pigment epithelial cells and RB cells. Among these lncRNAs, lnc-DAZ1-161 and lnc-HDAC7-21 were the top upregulated DELs in RB and were indicated to be significantly overexpressed in RB cells. However, the expression levels of lnc-OR52A1-55 in the cells were inconsistent with those in the tissue microarray data, which may be due to an insignificant log2(fold-change) of lnc-OR52A1-55 in the GEO database that was used in the present study compared with other top deregulated genes. Thus, the three lncRNAs require to be further studied in clinical RB tissues.

In conclusion, the present study used bioinformatics analyses to determine the lncRNA and mRNA expression profiles associated with RB. The potential pathological roles of the aberrantly expressed lncRNAs and mRNAs were predicted using GO and KEGG analyses. Furthermore, RT-qPCR validation confirmed that the expression of certain representative lncRNAs in RB cells was consistent with that in the tissue microarray database. These results may contribute to the elucidation of the molecular mechanisms of RB and provide novel targets for the development of improved clinical therapies. However, further studies are required to determine the biological functions of the dysregulated lncRNAs in RB.

Supplementary Material

Primer sequences used for PCR.

Supplementary Data

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available from the GEO repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111168). The experimental data (PCR) used and/or analyzed during the current study are available from the corresponding author on reasonable request..

Authors' contributions

XF and FC designed the study and contributed to the editing and proofreading of the manuscript. JG, QL, CX and JP retrieved and processed the data. JG, RZ and LZ analyzed the data. All authors read and approved the final manuscript. XF and FC checked and approved the authenticity of the raw data in the present study.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References