Abstract Optic nerve invasion (ONI) is an important high-risk feature and prognostic indicator of retinoblastoma (RB). Emerging evidence has revealed that non-coding RNAs (ncRNAs) play important roles in tumor perineural invasion (PNI). Nevertheless, the regulatory role of ncRNAs in the ONI of RB is poorly understood. In the current study, whole-transcriptome sequencing was performed to assess the expression profiles of ncRNAs and mRNAs in RB tissues, with or without ONI. Based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, we predicted the biological functions of differentially expressed (DE) mRNAs. We then constructed competing endogenous RNA (ceRNA) regulatory networks based on bioinformatics analysis. The hsa_circ_0015965/lncRNA MEG3-hsa-miR-378a-5p-NOTCH1 pathway was selected and validated by real-time qPCR, western blotting, and dual luciferase reporter assays. Moreover, we demonstrated that NOTCH1 promotes the malignant progression of RB. Taken together, our results provide novel insights into the mechanism underlying optic nerve invasion in RB. Keywords: ceRNA, ncRNA, Optic nerve invasion, Perineural invasion, Retinoblastoma 1. Introduction Retinoblastoma (RB) is the most common intraocular tumor in children [[33]1], and its prognosis has improved over the past 50 years [[34][2], [35][3], [36][4]]. Advances have been made based on the combination of surgery, chemotherapy, and radiotherapy, and the five-year survival rate of patients with stage III RB with orbital extension is more than 50%; however, the prognosis of stage IV patients with systemic metastasis or central nervous system (CNS) involvement remains unfavorable [[37]5]. Optic nerve invasion (ONI) has been confirmed as a significant risk factor for extraocular metastasis and recurrence of RB [[38][6], [39][7], [40][8]] and can cause intracranial metastasis via contiguous spread or seeding tumor cells into the leptomeninges and cerebrospinal fluid [[41]9]. However, the mechanisms underlying ONI in RB remain unclear. Therefore, identifying additional potential biomarkers is crucial for devising more effective strategies for preventing and treating RB with ONI. Non-coding RNAs (ncRNAs) were initially thought to be the products of random transcription with no function. However, recent emerging evidence has revealed the functional importance of ncRNAs in cellular physiology and pathology processes [[42]10,[43]11]. Generally, ncRNAs can be divided into several types, including microRNA (miRNA), long ncRNA (lncRNA), circular RNA (circRNA), endogenous small interfering RNA, Piwi-associated small RNA, small nucleolar RNA (snoRNA), sno-derived RNA, transcription-initiation RNA, and miRNA-offset RNA [[44]12]. Among them, miRNA, lncRNA and circRNA have attracted considerable interest and have been proven to be closely related to tumor progression [[45][13], [46][14], [47][15], [48][16]]. Perineural invasion (PNI) is an aggressive form of cancer cell invasion along nerves that occurs in many malignant tumors [[49]17]. The function of ncRNAs in PNI has been previously documented. Zhang et al. [[50]18] revealed that in salivary adenoid cystic carcinoma, CXCR5 facilitated PNI through the miR-187/S100A4 axis. Moreover, in cholangiocarcinoma, upregulated miR-21 expression promotes PNI by inhibiting RECK expression [[51]19]. However, the function of ncRNAs in the ONI of RB has seldom been reported. According to the competing endogenous RNA (ceRNA) hypothesis, mRNAs, lncRNAs and circRNAs affect target mRNAs by competing for the same miRNA [[52]20]. Accumulating evidence has demonstrated a critical role for ceRNAs in RB progression. Huang et al. [[53]21] reported that Circ-E2F3 promotes the progression of RB via the miR-204-5p/ROCK1 pathway. Li et al. [[54]22] indicated that TRPM2-AS acts as an oncogenic lncRNA in RB via the miR-497/WEE1 axis. Wang et al. [[55]23] revealed that DANCR promotes MMP-9-mediated progression and metastasis by working as a ceRNA that targets miR-34c and miR-613 in RB. However, to date, a comprehensive expression profile has not been developed for the miRNAs, lncRNAs, circRNAs and mRNA involved in the ONI of RB or their underlying regulatory networks. In this study, we performed whole-transcriptome sequencing in snap-frozen RB tissues with or without ONI to identify differentially expressed (DE) lncRNA, circRNA, miRNA, and mRNA. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on DE mRNA. Subsequently, we constructed ceRNA networks and a ceRNA pathway was selected and verified using real-time qPCR, western blotting, and dual-luciferase reporter gene assays. This is the first study using whole-transcriptome sequencing to depict ncRNAs and related ceRNA networks in the ONI of RB, and the findings may provide novel perspectives and strategies towards the development of targeted RB treatments. 2. Material and methods 2.1. Human tissue samples Tumor tissue samples from 22 RB patients who underwent enucleation surgery were included in this study. Six tissues were selected for whole-transcriptome sequencing analysis, and 16 samples were used to verify and confirm the accuracy of the microarray analysis. Normal retinal tissues were collected from eight patients who received enucleation due to ophthalmorrhexis. Informed consent for diagnostic and research studies was obtained from all subjects in accordance with the Declaration of Helsinki protocols, and the research protocol was approved by the Institutional Review Board of Eye & ENT Hospital of Fudan University. Prior to the surgery, none of the patients had received any other treatment, such as radiotherapy or chemotherapy. The diagnosis of the patients was histologically confirmed by two independent pathologists, and the tissues were stored in liquid nitrogen until use. The clinical and pathological features of the patients in this study are shown in [56]Table S1. 2.2. Cell culture Human retinoblastoma cell line (WERI-Rb-1) and human embryonic kidney cells (293 T) were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). For cancer cells, RPMI-1640 medium (Gibco, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (FBS; Invitrogen, Carlsbad, CA, USA) was used for cell culture. 293 T cells were cultured in DMEM (Gibco) supplemented with 10% FBS. All cells were grown at 37 °C in 5% CO[2] incubator. 2.3. RNA extraction and quality control Total RNA isolation from RB cells and tissues was performed using TRIzol reagent (Invitrogen). The concentration and integrity of the total RNA were assessed using a NanoDrop ND-1000 (Thermo Fisher, Wilmington, DE, USA) and Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). 2.4. Library construction, RNA sequencing and data analysis For long RNA library construction, total RNA was subjected to rRNA removal, fragmentation, first- and second-strand cDNA synthesis, end repair, poly A tail addition, ligation and enrichment. Small RNA sequencing used next-generation sequencing technology to sequence miRNAs (18–30 nt). The 18-30 nt small RNAs were ligated to a 5’- and 3’-adaptor. The adapter-ligated small RNAs were then transcribed into cDNA. RNA-seq sequencing was performed using DNBSEQ (BGI, Shenzhen, China). Long RNA reads were filtered using SOAPnuke (v1.5.2) and mapped to the reference genome using HISAT2 (v2.0.4). Bowtie2 (v2.2.9) was used to align the clean reads to known mature human miRNA sequences. Novel miRNAs was predicted using mirDeep2 (v2.0.0.8) software. The discrepancy between different groups was analyzed using DEGseq software. Detailed information was uploaded to the SRA database under accession number PRJNA861968 ([57]https://www.ncbi.nlm.nih.gov/sra/PRJNA861968). 2.5. Construction of competing endogenous RNA (ceRNA) network Three search platforms, miRanda ([58]http://www.microrna.org/microrna/home.do), RNAhybrid ([59]https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/submission.htm l/) and TargetScan ([60]http://www.targetscan.org/), were used to predict mRNA-miRNA, lncRNA-miRNA, and circRNA-miRNA interactions. mRNA was selected for the miRNA-target pairs if a negative correlation occurred between the expression of miRNA and mRNA. LncRNAs and circRNAs were analyzed in the same manner. Cytoscape (v3.6.1) software was used to visualize the ceRNA network. 2.6. Functional enrichment analysis The Database for Annotation, Visualization and Integrated Discovery (DAVID), an online functional annotation tool ([61]https://david.ncifcrf.gov/), was used to conduct Gene Ontology (GO) and Kyoto Genome Encyclopedia (KEGG) pathway analyses. The GO terms and KEGG pathways with p values of <0.05 were regarded as significantly enriched. 2.7. Real-time qPCR For cDNA synthesis, total RNA was reverse-transcribed using the PrimeScript RT Reagent Kit (Takara, Beijing, China). GAPDH served as the internal reference gene. miRNAs were detected using miDETECT A Track™ miRNA qRT-PCR Kit (RiboBio, Guangzhou, China) normalized by U6. The sequences are covered by a patent. All values were calculated by 2^−ΔΔCt method. The sequences of primers were shown in [62]Table S2. 2.8. Western blot Total proteins were extracted from cells and tissues were subjected to SDS/PAGE and transferred onto the PVDF membrane. Membranes were blocked at room temperature for 1 h with 5% skim milk, followed by overnight incubation at 4 °C with primary antibody (NOTCH1, 20687-1-AP, Proteintech, Wuhan, China; c-Myc, 10828-1-AP, Proteintech, Wuhan). After incubated with HRP-conjugated secondary antibody, immunoreactive bands were detected using the ECL western blotting system. Anti-GAPDH (5174, Cell Signaling Technology, Danvers, MA, USA) was used as an internal loading control. 2.9. Cell transfection Short interfering RNAs (siRNAs) targeting lncRNA MEG3, hsa_circ_0015965 overexpression plasmid and NOTCH1 overexpression plasmid, as well as hsa-miR-378a-5p mimics and inhibitor, were designed and synthesized by RiboBio (Guangzhou, China). Cell Transfection was performed using Lipofectamine 3000 (Invitrogen). RNA sequences of siRNAs are listed in Table S3. 2.10. Luciferase gene reporter assay Negative control mimics or hsa-miR-378a-5p mimics were co-transfected with the dual-luciferase reporter plasmid into 293 T cells using Lipofectamine 3000. Forty-eight hours later, cells were collected and lysed. Luciferase activity was detected using a dual-luciferase reporter assay kit (Promega, Madison, WI, USA). 2.11. Cell proliferation For the Cell Counting Kit-8 (CCK-8) assay, cells were seeded into 96-well plates with 10,000 cells/well (in 100 μl medium) and cultured for 24, 48, 72, 96, and 120 h, respectively. 10 μl CCK-8 (Selleck, shanghai, China) was then added to each well. After 2 h incubation in cell incubator, the absorbance at 450 nm was measured. EdU Cell Proliferation Kit (Beyotime, shanghai, China) was utilized to determine the ability of cell proliferation according to the manufacturer’s manual. Images were taken using fluorescence microscope. 2.12. Transwell migration and invasion assays Cells were plated into the top chambers with serum-free medium and migration-inducing medium containing 10% FBS was added into the bottom chambers. For Matrigel invasion assay, Matrigel was used to coat the top chamber before cells were seeded. After 24 h of cell culture, the filters were fixed using 4% paraformaldehyde. Cells on the upper side of the membrane were scraped and the bottom side stained using crystal violet. 2.13. Xenograft tumor model Six-week-old BALB/c nude male mice were subcutaneously injected with 1 ⨯ 10^6 RB cells. Tumors were visible at seven days after inoculation. Tumor size was measured weekly using caliper. Around 4 weeks after injection, all mice survived and were later euthanized, the tumors were dissected and weighted. No tumor abscess was observed in xenograft tumors. Immunohistochemistry staining was used to evaluate NOTCH1 expression levels in xenograft tumors. Animal experiments were approved by the Institutional Research Ethics Committee of Eye & ENT Hospital of Fudan University (Date April 2021/No 2021056). 2.14 Immunohistochemistry. Immunohistochemistry was performed as previously described [[63]24]. 2.14. Statistical analysis Statistical analyses were performed using SPSS 13.0 software (SPSS Inc., Chicago, IL, USA), and graphics were generated using GraphPad Prism 6.0 (GraphPad Software, Inc., San Diego, CA, USA). Data are presented as means ± standard deviation (SD). Student’s t-test was used to analyze the statistical differences between the two groups. The significance was defined as the p-value of <0.05. 3. Results 3.1. Expression profiles of circRNAs, lncRNAs, miRNAs and mRNAs The whole-transcriptome sequencing data of the two groups, including circRNAs, lncRNAs, miRNAs, mRNAs, were obtained using the DNBSEQ platform. Six RB patients who underwent primary enucleation were included in this study. Of these six patients, three patients were pathologically diagnosed with ONI and designated “RB with ONI” and the other three patients without ONI were designated “RB without ONI.” We analyzed DE ncRNAs and mRNAs using significance analysis of sequencing technique based on the cutoff fold changes ≥2 or ≤0.5 along with q value < 0.05. The analytical design and workflow of this study are shown in [64]Fig. 1. Fig. 1. [65]Fig. 1 [66]Open in a new tab Flow Chart of This Study Design. DE, differentially expressed. Volcano plots ([67]Fig. 2A) and heatmaps ([68]Fig. 2B) were used to demonstrate the DE ncRNAs and mRNAs in the RB samples. High-throughput sequencing revealed that 94 circRNAs, 64 lncRNAs, 67 miRNAs, and 141 mRNAs were differentially expressed between RB with ONI and without ONI group. As shown in [69]Table 1, the most upregulated RNAs were circRNA hsa_circ_0002493, lncRNA LOC101928517, miRNA hsa-miR-544a, and mRNA DLK1. The most downregulated RNAs were circRNA hsa_circ_0091543, lncRNA LOC105375301, miRNA hsa-miR-4707-3p and mRNA NLRP2. The top five upregulated and downregulated ncRNAs and mRNAs are listed in [70]Table S4-7. Fig. 2. [71]Fig. 2 [72]Open in a new tab Differentially expressed ncRNAs and mRNAs in RB tissues with or without ONI. Volcano plot (A) and heatmap (B) displaying ncRNA and mRNA expression of differentially expressed genes between the two groups. ONI, optic nerve invasion. (C) The gene ontology enrichment analysis covered three domains: biological process, cellular component and molecular function. (D) KEGG pathway enrichment analysis of DE mRNAs. Enrichment score values were calculated as − log [10] (P values). Table 1. Statistical Analysis of All of the Differentially Expressed ncRNAs and mRNAs. DE RNAs Total No. No. of Upregulated No. of Downregulated The Most Upregulated (Fold Change) The Most Downregulated (Fold Change) circRNA 94 47 47 hsa_circ_0002493 (21.61) hsa_circ_0091543 (21.53) lncRNA 64 33 31 LOC101928517 (22.04) LOC105375301 (21.09) miRNA 67 59 8 hsa-miR-544a (11.32) hsa-miR-4707-3p (11.02) mRNA 141 73 68 DLK1 (10.92) NLRP2 (9.37) [73]Open in a new tab 3.2. GO and KEGG analyses of DE mRNAs The circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA networks were constructed based on the analysis design and workflow. GO and KEGG pathways were used to elucidate the potential functions of the DE mRNAs involved in the ceRNA network. The results of the GO enrichment analysis showed that in the biological process (BP) category, the DE mRNAs were mainly related to nervous system development, forebrain development and regulation of membrane potential; in the cell component (CC) category, they were mainly enriched in MHC class II protein complex, plasma membrane and apical plasma membrane; and in the molecular function (MF) category, they were enriched in ionotropic glutamate receptor activity, MHC class II receptor activity and MHC class II protein complex binding ([74]Fig. 2C). In addition, the KEGG pathway analysis revealed that there were 13 pathways enriched in these DE mRNAs ([75]Fig. 2D). Among them, the antigen processing and presentation, Th1 and Th2 cell differentiation and asthma pathway were the most abundant. 3.3. ceRNA network construction CircRNA and lncRNA can generally regulate the activity of mRNAs by acting as a miRNA sponge. The miRanda, RNAhybrid, and TargetScan platforms were used to identify potential RNAs (including circRNAs, lncRNAs and mRNAs) that share the miRNA-recognition elements (MREs) with miRNAs. Based on the sequencing data, we built circRNA–miRNA–mRNA and lncRNA-miRNA-mRNA interaction networks. A total of 306 circRNA-miRNA-mRNA networks were established, including 23 circRNAs, 10 miRNAs, and 46 mRNAs ([76]Fig. 3A). Additionally, 515 lncRNA-miRNA-mRNA networks were constructed, including 20 lncRNAs, 17 miRNAs, and 49 mRNAs ([77]Fig. 3B). Fig. 3. [78]Fig. 3 [79]Open in a new tab ceRNA networks construction. (A)cirRNA-miRNA-mRNA network. (B) lncRNA-miRNA-mRNA network. (C) lncRNA/circRNA-miRNA-mRNA integrated network. Analysis of the lncRNA and circRNA ceRNA networks revealed that nine miRNAs were included in both; therefore, an lncRNA/circRNA-miRNA-mRNA integrated network was also established ([80]Fig. 3C). To identify the key nodes associated with ONI in RB, the degree, betweenness centrality, and closeness centrality of all circRNAs, lncRNAs, and miRNAs were calculated using cytoHubba plugin based on Cytoscape. The results for the top there circRNAs, lncRNAs, and miRNAs are listed in [81]Table 2. Table 2. The node scores of circRNAs, lncRNAs, and miRNAs calculated by cytoHubba plugin based on Cytoscape. ID Degree Betweenness Closeness hsa_circ_0015965 3 30.95 27.33 hsa_circ_0139601 2 17.23 26.42 hsa_circ_0113239 1 0 24.55 lncRNA MEG3 4 87.26 30.00 LINC02593 4 87.26 30.00 DIO3OS 4 87.26 30.00 hsa-miR-6724-5p 39 1957.92 45.17 hsa-miR-1275 23 950.18 34.50 hsa-miR-378a-5p 20 669.71 32.50 [82]Open in a new tab 3.4. Verification of key DE ncRNAs and mRNAs in RB tissues To confirm and verify the sequencing data, we first performed real-time qPCR to detect the expression levels of the top there key circRNAs, lncRNAs, and miRNAs in eight RB tissues with ONI and 8 RB tissues without ONI. The results revealed the upregulation of hsa_circ_0015965, hsa_circ_0139601, hsa_circ_0113239, lncRNA MEG3, LINC02593, and DIO3OS and the downregulation of hsa-miR-378a-5p, hsa-miR-1275 and hsa-miR-6724-5p expression in RB tissues with ONI compared to the control ([83]Fig. 4A-I). These were consistent with the RNA sequencing results. We also examined the expression of the top there key circRNAs, lncRNAs, and miRNAs in normal retinal tissues. Results showed that hsa_circ_0015965, hsa_circ_0139601, hsa_circ_0113239, lncRNA MEG3, LINC02593, and DIO3OS expression was enhanced and hsa-miR-378a-5p, hsa-miR-1275 and hsa-miR-6724-5p expression were downregulated in RB tissues compared with normal retinal tissues ([84]Fig. 4A-I). These results indicated that these ncRNAs may contribute to the progression and ONI of RB. Fig. 4. [85]Fig. 4 [86]Open in a new tab Verification of the differentially expressed key transcripts. The relative expression levels of 3 circRNAs (A–C), 3 lncRNAs (D–F), 3 miRNAs (G–I), and 2 mRNAs (J–K) in normal retinal tissues and RB tissues were detected by qRT-PCR. The protein level of NOTCH1 and c-Myc in normal retinal tissues and RB tissues were detected by Western Blot (L). ONI, optic nerve invasion. Data are presented as the mean ± SD. Data were analyzed with Student’s t-test. * indicates p < 0.05, ** indicates p < 0.01 and *** indicates p < 0.005. Of these verified ncRNAs, hsa_circ_0015965 was the most significantly upregulated circRNA, and lncRNA MEG3 was the most significantly upregulated lncRNA. The most significantly down-regulated miRNA was hsa-miR-378a-5p. Hsa_circ_0015965 and lncRNA MEG3 could bind to hsa-miR-378a-5p in the ceRNA network. In addition, NOTCH1 binds to hsa-miR-378a-5p. NOTCH1 and c-Myc, a direct target of NOTCH1 [[87]25], were overexpressed in RB tissues. Compared with RB tissues without ONI, RB tissues with ONI revealed much higher NOTCH1 and c-Myc levels ([88]Fig. 4J-L). NOTCH1 is a member of the Notch family, and deregulation of the NOTCH1 pathway occurs in multiple cancers [[89]26]. Strong NOTCH1 expression is significantly related to the presence of perineural invasion in poorly differentiated oral squamous cell carcinoma (OSCC) and colorectal cancer (CRC) [[90]27,[91]28]. Recent studies showed that Myc is implicated in subsets of high-risk RB [[92]29]. Consequently, the hsa_circ_0015965/lncRNA MEG3-hsa-miR-378a-5p-NOTCH1 pathway was selected for validation. 3.5. Verification of the circRNA/lncRNA-miRNA-mRNA pathway through dual luciferase reporter gene system Transfection of the hsa_circ_0015965 overexpression vector or the control vector was conducted in the WERI-Rb-1 cell line, and the overexpression effect was confirmed by qPCR ([93]Fig. 5A). The effects of the miR-378a-5p mimics or miR-378a-5p inhibitor were also examined ([94]Fig. 5B). We found that compared with the control vector, hsa-miR-378a-5p expression was significantly decreased after transfection with hsa_circ_0015965 overexpression vector ([95]Fig. 5C) while NOTCH1 and c-Myc mRNA and protein expression were increased ([96]Fig. 5D-F). Additionally, the upregulation of NOTCH1 and c-Myc mediated by hsa_circ_0015965 was suppressed by hsa-miR-378a-5p mimics ([97]Fig. 5D-F). Collectively, hsa_circ_0015965 activated NOTCH1 by targeting hsa-miR-378a-5p. Fig. 5. [98]Fig. 5 [99]Open in a new tab The hsa_circ_0015965/lncRNA MEG3-hsa-miR-378a-5p-NOTCH1 pathway was verified by dual luciferase reporter gene system. (A) qRT-PCR analysis of the effect on overexpression of hsa_circ_0015965 by vector transfection in WERI-Rb-1 cell line. (B) qRT-PCR analysis of the effect on mimics or inhibitor of miR-378a-5p in WERI-Rb-1 cell line. (C) The expression level of hsa-miR-378a-5p in hsa_circ_0015965 overexpressed RB cells determined with qRT-PCR. (D–E) The mRNA expression level of NOTCH1 (D) and c-Myc (E) in RB cells after hsa_circ_0015965 overexpressing with or without hsa-miR-378a-5p mimics or inhibitor determined by qRT-PCR. (F) The protein expression level of NOTCH1 and c-Myc in RB cells after hsa_circ_0015965 overexpressing with or without hsa-miR-378a-5p mimics or inhibitor determined using Western Blot. (G) qRT-PCR analysis of the effect on knockdown of MEG3 expression by siRNA in WERI-Rb-1 cell line. (H) The expression level of hsa-miR-378a-5p in MEG3 silenced RB cells determined with qRT-PCR. (I–J) The mRNA expression level of NOTCH1(I) and c-Myc (J) in RB cells after MEG3 silencing with or without hsa-miR-378a-5p mimics or inhibitor determined using qRT-PCR. (K) The protein expression level of NOTCH1 and c-Myc in RB cells after MEG3 silencing with or without hsa-miR-378a-5p mimics or inhibitor determined using Western Blot. (L) Potential binding sites of hsa-miR-378a-5p to hsa_circ_0015965, lncRNA MEG3 and NOTCH1 with mutation sites for specific binding assay. (M) The relative luciferase activities were determined in 293 T cells after co-transfection with WT or MUT vectors with or without hsa-miR-378a-5p mimics. NC, negative control; OE, overexpression. Data are presented as the mean ± SD. Data were analyzed with Student’s t-test. * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.005 and **** indicates p < 0.001. Then, siRNA were designed to knock down lncRNA MEG3 expression and siMEG3-1 (called siMEG3) with highest inhibition efficacy was chosen for further analysis ([100]Fig. 5G). The qPCR results showed that lncRNA MEG3 knockdown increased hsa-miR-378a-5p levels ([101]Fig. 5H). Downregulation of lncRNA MEG3 decreased NOTCH1 and c-Myc expression, which could be rescued by hsa-miR-378a-5p inhibitor ([102]Fig. 5I–K). The results indicated that lncRNA MEG3 could also regulate NOTCH1 expression via hsa-miR-378a-5p. The bioinformatics prediction tools were next utilized to confirm that hsa-miR-378a-5p can target hsa_circ_0015965, lncRNA MEG3 and NOTCH1 mRNA with complementary binding sites ([103]Fig. 5L). We performed dual luciferase assays to investigate whether hsa-miR-378a-5p could interact with hsa_circ_0015965, lncRNA MEG3 and NOTCH1. The wild-type and mutant 3′-UTR sequences of hsa_circ_0015965, lncRNA MEG3 and NOTCH1 were cloned to construct wild-type (WT) and mutant vectors (MUT) luciferase reporter vectors ([104]Fig. 5L). Co-transfection with WT or MUT luciferase reporter plasmids and hsa-miR-378a-5p mimics were then performed in 293 T cells. Results showed that the wild-type luciferase reporter activity of hsa_circ_0015965, lncRNA MEG3 and NOTCH1 was significantly decreased by hsa-miR-378a-5p mimics, while the MUT luciferase reporter activity showed no obvious change ([105]Fig. 5M). These results suggest that hsa_circ_0015965, lncRNA MEG3 and NOTCH1 are the direct targets of hsa-miR-378a-5p. Taken together, these findings confirmed the hsa_circ_0015965/lncRNA MEG3-hsa-miR-378a-5p-NOTCH1 pathway. 3.6. NOTCH1 promotes RB cell proliferation, migration and invasion in vitro and RB growth in vivo To explore the biological functions of NOTCH1 in RB, we overexpressed NOTCH1 in RB cells. The overexpression efficiency was confirmed by qPCR and Western blot ([106]Fig. 6A–B). CCK8 and EdU assays showed that NOTCH1 overexpression significantly promoted cell proliferation ([107]Fig. 6C–D). Moreover, transwell migration and invasion assays revealed that the migration and invasion abilities of RB cells were elevated by NOTCH1 overexpression ([108]Fig. 6E). Next, NOTCH1 overexpressing and control RB cells were subcutaneously injected into nude mice to evaluate the effects of NOTCH1 overexpression in vivo. The results showed that NOTCH1 overexpression increased the rate of tumor growth ([109]Fig. 6F). The overexpression efficiency of NOTCH1 in xenografted tumors was confirmed by immunohistochemistry ([110]Fig. 6G). Taken together, these results demonstrate that NOTCH1 overexpression potentiates the malignant progression of RB. Fig. 6. [111]Fig. 6 [112]Open in a new tab NOTCH1 promotes RB cell proliferation, migration and invasion in vitro and RB growth in vivo. (A) qRT-PCR analysis of NOTCH1 mRNA expression in RB cells. (B) Western blots of NOTCH1 protein in RB cells. (C) Cell proliferation was detected by CCK8 assay. (D) Representative images (left) and quantification (right) of EdU incorporation assay. Scale bar, 100 μm. (E) Representative images (left) and quantification (right) of the migration and invasion assay. Scale bar, 100 μm. (F) Representative tumor images (left), tumor volume curves and tumor weight quantification (right) were shown. (G) Representative IHC staining for NOTCH1 expression in xenografts (left). Scale bars, 50 μm. Quantitative analysis of NOTCH1 protein expression based on staining index in xenografts (right). OE, overexpression. Data are presented as the mean ± SD. Data were analyzed with Student’s t-test. * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.005 and **** indicates p < 0.001. 4. Discussion RB is the most common eye cancer in children [[113]30]. CNS metastasis is the most important cause of death in RB and occurs owing to ONI [[114][31], [115][32], [116][33]]. In the past few decades, numerous studies have attempted to diagnose ONI in RB using imaging to select better treatments [[117][34], [118][35], [119][36]]. However, only a few studies have explored the underlying molecular mechanism of ONI in RB. Nikia Laurie et al. verified that changes in cadherin-mediated cell adhesion are associated with ONI of RB prior to metastasis [[120]37]. Laura Asnaghi et al. founded that ACVR1C mRNA levels were significantly increased in RB tissues which invaded the optic nerve [[121]38]. PNI is an important adverse pathological feature of many malignancies. Recent works have revealed that PNI of cancer is a continuous multistep process among tumor cells, nerves and stromal cells in perineural niche. In the tumor background environment of hypoxia, high glucose level, inflammatory reaction, and sympathetic system activation, tumor cells, nerve cells, supporting cells, recruited inflammatory cells, blood vessels, immune cells, ncRNAs and various soluble signaling molecules and their receptors in a perineural niche comprise a complex signal network, which achieves the dynamic interaction between nerve and tumor [[122]17]. Among these molecules, ncRNAs can regulate PNI by orchestrating oncogene and anti-oncogene, and they are hot topics of current studies [[123]39]. Downregulated hsa_circ_0065149 was significantly correlated with PNI in gastric cancer tissues [[124]40]. In colorectal cancer, lncRNA-CYTOR may enhance PNI via CYTOR/miR-3679-5p/MACC1 axis [[125]41]. However, no comprehensive analysis of expression profiles of mRNAs and ncRNAs in RB optic nerve infiltration has been reported to date. In the present study, whole-transcriptome sequencing techniques were used to identify the mRNAs and ncRNAs expression profiles in RB with or without ONI. The data revealed that 94 circRNAs, 64 lncRNAs, 67 miRNAs and 141 mRNAs were significantly differentially expressed. Subsequently, bioinformatics methods, including GO, KEGG pathway and ceRNA network analyses were carried out to explore the potential functions of the DE transcripts. GO and KEGG pathway analyses showed that these DE transcripts were mainly enriched in nervous system development, MHC class II protein complex, ionotropic glutamate receptor activity, antigen processing and presentation, Th1 and Th2 cell differentiation, and asthma pathway, suggesting that tumor neurogenesis, neurotransmitters-glutamate and immune inflammatory response are involved in the process of RB optic nerve infiltration. In PNI, reciprocal communication occurs between tumors and the nervous system. On the one hand, nerve fibers stimulate tumor progression by releasing neurotransmitters [[126]42]. Glutamate is the major excitatory neurotransmitter in the CNS. Many studies have proved the role of glutamate and its receptors in tumor growth, metastasis and PNI. Glutamate-related receptors and their signaling pathways may provide novel treatment strategies for malignant tumors [[127]43]. On the other hand, cancer stem cells differentiate to autonomic neuron-like cells generating tumor-derived neurogenesis [[128]44] and neural progenitor cell migration and differentiation of neurons into tumors producing tumor-induced neurogenesis [[129]45]. Tumor neurogenesis is a potential target for cancer therapy. Immune cells can also communicate with tumor cells and nervous systems leading to poor prognosis [[130]46]. Thus, it would be important to explore the cellular and molecular basis of communication among tumors, nervous and immune cells, which may help us gain an in-depth understanding of RB optic nerve infiltration. To explore the functions of ncRNAs in ONI of RB, we constructed an integrated lncRNA/circRNA-miRNA-mRNA network. CytoHubba plugin based on Cytoscape was used to identify key ncRNAs. The results showed that the three top-ranked circRNAs, lncRNAs, and miRNAs were hsa_circ_0015965, hsa_circ_0139601, and hsa_circ_0113239; lncRNA MEG3, LINC02593, and DIO3OS; and hsa-miR-378a-5p, hsa-miR-1275 and hsa-miR-6724-5p. Next, the hsa_circ_0015965/lncRNA MEG3-hsa-miR-378a-5p-NOTCH1 pathway was selected and verified to demonstrate the credibility of our results. Several reports have indicated that NOTCH1 is associated with PNI in cancer [[131]27,[132]28]. However, its role in RB ONI has not been demonstrated. In the present study, we found that NOTCH1 was substantially upregulated in RB tissues with optic nerve infiltration and NOTCH1 overexpression contributes to the malignant progression of RB. However, our study had some limitations. Up date, in vitro or in vivo models of optic nerve infiltration of RB have not been constituted, which makes it impossible for us to verify the function of ceRNA network in ONI; thus, further investigation is required. In conclusion, the present study is the first to identify whole transcriptome expression profiles in RB with ONI. Furthermore, GO, KEGG pathway and ceRNA network analyses were conducted to explore the potential mechanisms of ONI. Additionally, the hsa_circ_0015965/lncRNA MEG3-hsa-miR-378a-5p-NOTCH1 pathway was selected and verified. Our findings preliminarily revealed the ceRNA regulatory networks involved in ONI, which could provide new novel perspectives and strategies towards the target of RB treatment. 5. Ethics statement All experiments have been carried out adhering to the tenets of the Helsinki Declaration and were approved by the Investigational Review Board of the Eye and ENT Hospital of Fudan University. Informed consent was obtained from all subjects involved in the study. Declarations Author contribution statement Jie Sun: Performed the experiments; Wrote the paper. Lu Gan; Jie Ding: Analyzed and interpreted the data. Ruiqi Ma: Contributed reagents, materials, analysis tools or data. Jiang Qian; Kang Xue: Conceived and designed the experiments. Funding statement This work was supported by the Science and Technology Commission of Shanghai Municipality [20Y11911200]; the National Natural Science Foundation of China [82171099, 81970835]; and the Natural Science Foundation of Shanghai [20ZR1409500]. Data availability statement Data associated with this study has been deposited at SRA database under the accession number PRJNA861968. Declaration of interest’s statement The authors declare no conflict of interest. Footnotes ^Appendix A Supplementary data to this article can be found online at [133]https://doi.org/10.1016/j.heliyon.2023.e13813. Contributor Information Jiang Qian, Email: qianjiang@fudan.edu.cn. Kang Xue, Email: xuekang@foxmail.com. Abbreviations RB Retinoblastoma ONI Optic nerve invasion PNI Perineural invasion ncRNA non-coding RNA GO Gene ontology KEGG Kyoto Encyclopedia of Genes and Genomes DE Differentially expressed ceRNA Competing endogenous RNA CNS Central nervous system Appendix A. Supplementary data The following are the Supplementary data to this article. Multimedia component 1 [134]mmc1.docx^ (19.5KB, docx) Multimedia component 2 [135]mmc2.docx^ (15KB, docx) Multimedia component 3 [136]mmc3.docx^ (14KB, docx) Multimedia component 4 [137]mmc4.docx^ (14.8KB, docx) Multimedia component 5 [138]mmc5.docx^ (14.9KB, docx) Multimedia component 6 [139]mmc6.docx^ (14.8KB, docx) Multimedia component 7 [140]mmc7.docx^ (14.9KB, docx) Multimedia component 8 [141]mmc8.docx^ (671.6KB, docx) References