Comprehensive Analysis of Serum Small Extracellular Vesicles-Derived Coding and Non-Coding RNAs from Retinoblastoma Patients for Identifying Regulatory Interactions

Simple Summary The diagnosis of retinoblastoma (RB) is usually made by clinical examination and imaging modalities. Routine tissue biopsy is not done due to the risk of extraocular spread. Blood-based RNA cargoes could be promising surrogate markers for RB diagnosis and prognostication. Our data indicated that the size, morphology, and zeta potential (ZP) of RB and non-RB serum extracellular vesicles (EVs) met standard exosome properties with similar concentrations. MALTA1, AFAP1-AS1, miR-145, and miR-101 were identified as hub non-coding RNAs that promote RB progression by targeting cyclins, cyclin-dependent kinases, c-MYC, EZH2, ZEB1, TP53, and BCL2. Along with these, the aberrantly expressed miRNAs, lncRNAs, and their target mRNAs of RB EVs were implicated in cell cycle, metabolism, and tumor-associated signaling pathways. The differential expression of EV RNAs in RB compared to controls may aid in the identification of possible serum prognostic biomarkers for RB. Abstract The present study employed nanoparticle tracking analysis, transmission electron microscopy, immunoblotting, RNA sequencing, and quantitative real-time PCR validation to characterize serum-derived small extracellular vesicles (sEVs) from RB patients and age-matched controls. Bioinformatics methods were used to analyze functions, and regulatory interactions between coding and non-coding (nc) sEVs RNAs. The results revealed that the isolated sEVs are round-shaped with a size < 150 nm, 5.3 × 1011 ± 8.1 particles/mL, and zeta potential of 11.1 to −15.8 mV, and expressed exosome markers CD9, CD81, and TSG101. A total of 6514 differentially expressed (DE) mRNAs, 123 DE miRNAs, and 3634 DE lncRNAs were detected. Both miRNA-mRNA and lncRNA-miRNA-mRNA network analysis revealed that the cell cycle-specific genes including CDKNI1A, CCND1, c-MYC, and HIF1A are regulated by hub ncRNAs MALAT1, AFAP1-AS1, miR145, 101, and 16-5p. Protein-protein interaction network analysis showed that eye-related DE mRNAs are involved in rod cell differentiation, cone cell development, and retinol metabolism. In conclusion, our study provides a comprehensive overview of the RB sEV RNAs and regulatory interactions between them.

The collected blood samples were centrifuged at 2000× g for 15 min to remove any cellular debris. The supernatants containing the cell-free serum samples were stored in aliquots at −80 • C. Fresh RB tumor tissues (n = 5) were collected following enucleation of the eye as part of the treatment protocol for advanced intraocular tumor. An area of maximum tumor volume based on orbital imaging was identified, a 5 mm sclero-choroidal incision was given with a blade, the tumor was identified, and 2-3 mm of fresh tumor was obtained and transferred to an aliquot containing 400 µL of RNA later. The entire procedure was performed in the operation theater immediately after enucleation, and under sterile conditions. The enucleated globe was then submitted for routine histopathological examination. The control retinas (n = 2) were acquired from human cadaveric eyeballs, which were collected within 6 h of death and preserved in a sterile moist chamber in Ramayamma International Eye Bank at LVPEI, Hyderabad. The globe was bisected with blade adjacent to the optic nerve, the retina was identified, carefully dissected, and 2 mm of it was transferred to an aliquot containing 400 µL of RNA later and stored at −80 • C until RNA isolation. The remaining retina was sectioned and stained to rule out any evidence of retinal disease.

Small Extracellular Vesicles Isolation from Serum
sEVs were recovered from serum samples using the commercial kit Total Exosome Isolation™ from serum (Invitrogen by Thermo Fisher Scientific, Vilnius, Lithuania) as per manufacturer's instructions. Briefly, 100 µL of reagent was added to 500 µL of serum sample, and the solution was incubated overnight at 4 • C. Then precipitated vesicles were recovered by centrifugation at 10,000× g for 10 min. The pellet was dissolved in 100 µL PBS, aliquoted and stored at −80 • C for further experiments. The detailed methodology for comprehensive analysis of serum sEVs is illustrated in Figure S1. From all the 9 RB patients, sEVs were isolated and characterized for physical properties, but for RNA sequencing R4, R5 and R6 aged 5, 2, and 4 were selected. All three of them were male unilateral RB cases (Table S2). Among the five controls, we randomly picked first three controls aged 5, 5 and 4.

Transmission Electron Microscopy
The morphology of sEVs was examined by TEM according to the technique described by Ahmed et al. [24]. Briefly, 20 µL of EVs PBS solution (1:100 dilutions) drop was loaded onto carbon coated copper grids and permitted to stand overnight for air drying. The absorbed sEVs were negatively stained with 2% uranyl acetate for 10 min. Finally, the images of sEVs were captured under TEM (JEOL JEM-1400Flash, Road Peabody, MA, USA) at 80 kV after the grids were dried.

Nanoparticle Tracking Analysis by Zeta View
All the purified single freeze-thawed sEV fractions were analyzed for particle size, concentration, and zeta potential using a ZetaView device (Particle Metrix GmbH, Meerbusch, Germany) and its corresponding software (version 8.05.12 SP1) [25]. Each sEV sample (1 µL in triplicates) was diluted (1:4000) in PBS, and 2 µL of this solution was injected into the cell. The instrument measured each sample at 11 different positions throughout the cell, with two cycles of readings at each position. The detection threshold of the zetaView software was set to 5, and the maximum jump distance, and the minimum track segment length were both set to auto. After automated analysis of all 11 positions, the mean, median and mode (indicated as diameter) sizes, ZP, as well as the concentration of the sample were calculated by the inbuilt software. As described earlier [26], we selected the mode as the measurement for size in our analysis. Final concentration was calculated by multiplying the observed concentration with the dilution factor. The concentration of EVs present in each sample was expressed as (particles/mL).

Immunoblotting for Exosome Specific Proteins
To demonstrate the presence of exosomal protein markers in serum small EV samples, immunoblotting was performed. Total exosomal protein concentration was estimated by a commercially available BCA kit (Thermo Fisher Scientific. Waltham, MA, USA) as per manufacture instructions. Then denatured proteins were separated on 10% SDS PAGE and transferred onto a PVDF membrane using a Trans-Blot ® SD Semi-Dry Transfer Cell (Bio Rad, Hercules, CA, USA). Membranes were blocked with 5% nonfat milk TBST for 1 h at room temperature and then incubated with primary antibodies anti-CD9 (D8O1A, Cell Signaling Technology, Danvers, MA, USA), anti-CD81 (EPR4244, Abcam, Cambridge, MA, USA), and anti TSG101 (EPR7130(B)), Abcam, Cambridge, MA, USA), overnight at 4 • C. The membranes were washed three times with 1 × TBST for 10 min and incubated with an HRP-conjugated secondary antibody (7074, Cell Signaling Technology, Danvers, MA, USA) for 2 h at room temperature, and washed in TBST. Signals were visualized after incubation with enhanced chemiluminescence kit by ChemiDoc (Bio-Rad, Hercules, CA, USA).

RNA Isolation, Library Preparation, Sequencing, and Data Processing
Large RNA (>200 nt) and small RNA (<200 nt) fractions were extracted from pooled sEV preparations of RB patients (n = 3) and age-matched controls (n = 3) using Total Exosome RNA & Protein Isolation Kit (Invitrogen™) following the manufacturer's protocol. One mL PBS was added to the enriched sEV pellet. To this, 2% denaturation solution was added and mixed thoroughly. The mixture was incubated for 5 min at 4 • C. Next, one volume of acid-phenol-chloroform solution was added, the samples mixed by vortexing for 30-60 S and then centrifuged for 5 min at 10,000× g at room temperature. The aqueous fraction was collected, and 1/3 volume of 100% ethanol was added. The lysate was then passed through the filter cartridge and centrifuged for 30 S. At this stage, large RNA and small RNA were collected: the flow through contained small RNA and filter contained large RNA. Then, the large RNA filter was centrifuged for an additional 1 min at 10,000× g, transferred into a fresh collection tube, and 50 µL of nuclease-free water was added to the center of the filter. The samples were centrifuged for 30 S at 10,000× g. The elute containing the large RNA was collected and stored at −80 • C.
For small RNA recovery, a 2/3rd volume of ethanol was added to the flow-through containing the small RNA and mixed thoroughly. Ethanol mixture containing the small RNA was dropped onto the second filter cartridge, centrifuged for 30 S, and the same step repeated twice. The flow-through was discarded. The small-RNA bound cartridge was washed with wash buffer I, II and III, centrifuged, and the flow-through discarded. A volume of 50 µL of nuclease free water was added to the small RNA-bound cartridge and centrifuged for 30 S to recover the RNA. The eluate, which contained small RNAs was stored at −80 • C. Quality and quantity of the sEV RNAs were analyzed by NanoDrop TM 2000 (Thermo Scientific). The large RNA fraction was subjected to whole transcriptome analysis (WTA) and the small RNA portion to small RNA sequencing.

Differential Expression Analysis
Differential expression analysis of mRNAs, lncRNAs and miRNAs was performed using EdgeR package (v 3.28.1) that counted data using a negative binomial distribution, and implemented statistical methods proposed by Robinson and Smyth [29]. Individual gene expression was calculated as the mean expression of each gene averaged over all samples of each group and presented as the logarithm of counts per million reads. A p-value cutoff of 0.05 and log2foldchange of (+/−) 2 were used for identifying significantly differentially regulated transcripts. The Benjamini Hochberg procedure, proposed by Benjamini and Hochberg was used to control the false discovery rate (FDR) [30].

Functional Enrichment Analysis
The ClusterProfiler package (http://bioconductor.org/packages/release/bioc/html/ clusterProfile-r.html) [36] was used to identify and visualize the top GO terms and KEGG pathways enriched by DE mRNAs and DE lncRNAs. The target genes of miRNAs were selected for GO analysis using GOnet tool http://tools.dice-database.org/GOnet/ accessed on 16 February 2022) [37]. Gene set enrichment analysis (GSEA) was performed using gene pathways extracted from GO and KEGG. Values were derived for the set of genes by permuting the gene sets for 'n' number of times within the available datasets. Cohesive rankings and differences in pathways were assigned by the GSEA algorithm. Validation of rankings and differences was done by calculating the statistical significance of the normalized enrichment score using nominal p-value and FDR q-value. The cutoffs for p-value and FDR q-value were set to 0.05 and 0.25.

Construction of RNA Interaction Networks
The interaction network of DE miRNAs and their target genes identified in RB serum small EVs was generated by miRTargetLink 2.0 (https://www.ccb.uni-saarland.de/ mirtargetlink2 accessed on 3 March 2022) with parameter setting to strong interaction and minimum 5 shared targets options [38]. The associated biological processes were extracted by GeneTrail3.0 (http://genetrail.bioinf.uni-sb.de accessed on 3 March 2022) [39]. The miRNA-target network analysis was carried out using the tool cytoscape with 10 selected miRNAs that target the gene of interest RB1. An LncRNA-miRNA-mRNA interaction network was constructed on the basis that 23 selected lncRNAs bind to their target genes (either to mRNAs or miRNAs or to both). The network was visualized using cytoscape software (v 3.8.2; accessed on 15 March 2022) [40]. In this network, each node represented a biological molecule, and the edges were defined as interactions between nodes [41]. LncRNAs, mRNAs, and miRNAs in the network were presented as green triangles, blue rectangles, and pink ovals, respectively. Hub RNAs was selected based on the topological features of the network such as betweenness centrality, closeness centrality network and degree layout calculated by a built-in NetworkAnalyzer tool in Cytoscape software [42]. The functionally enriched KEGG pathway and GO terms for the network were visualized using the ClueGO/CluePedia plugin from Cytoscape [43,44].

Construction of Protein-Protein Interaction Network
A total of 39 DE mRNAs (20 up and 19 down regulated coding mRNAs) related to ocular development were filtered from core DE mRNAs of RB sEVs. We used a string database (version: 11.5 available at: https://string-db.org/ accessed on 20 March 2022) for construction of the protein interaction (PPI) network [45]. The protein names for the selected genes were given as input data by selecting multiple proteins by names/identifiers and interactions pertaining to Homo sapiens. The topological and functional properties of the PPI network were analyzed by K-means clustering, an unsupervised learning algorithm [46]. The resulting clusters were separated manually for better visual representation and comprehension of the interaction network. The functional significance and statistical analysis of the network were investigated automatically by inbuilt String software against the statistical background of the whole genome.

Quantitative Reverse Transcriptase-Polymerase Chain Reaction
Among the total DE mRNAs identified in RB sEVs by RNA sequencing, the expressions of HIF1A, SYK, and PGK1 were analyzed for their expression in corresponding tumor tissues by RT-qPCR. RB tissues and retinas were thawed, and total RNA was recovered using TRIZIN reagent (GCC Biotech Pvt. Ltd., West Bengal, India). The samples were homogenized in trizin and incubated at room temperature for 5 min. Chloroform was added and the tubes were mixed vigorously for 15 S and incubated at room temperature for 2-3 min. Next, the samples were centrifuged at 12,000× g for 15 min at 4 • C. The aqueous phase was precipitated with isopropanol, followed by 75% ethanol washes. After the washes, the RNA pellet was air dried and dissolved in nuclease fee water. The isolated RNA was quantified using NanoDrop TM 2000 (ThermoFisher Scientific, Waltham, MA, USA). cDNA was prepared from 1 µg of RNA using RevertAid First Strand cDNA Synthesis kit (Thermo Scientific) using oligo (dT)18 primer and random hexamer primer according to the manufacturer instructions. qPCR was performed on a 7900HT Fast RT PCR system (Applied Biosystems, MA, USA) using DyNAmo™ Flash SYBR Green qPCR Kit (ThermoFisher Scientific MA, USA) as per manufacturer's instructions. The reaction was performed in 96-well transparent plates (Thermo Fisher Scientific) for real time in a final volume of 10 µL. For each gene, three technical replicates of each sample were analyzed along with negative controls and a 5-point relative standard curve and the non-template control. The following amplification conditions were used: 10 min at 95 • C, 40 cycles of 15 S at 95 • C (denaturation), and 1 min at 60 • C (annealing and extension). A dissociation protocol with incremental temperatures of 95 • C for 15 S plus 65 • C for 15 S was used to investigate the specificity of the qPCR reaction. The qPCR expression data for each reference gene were extracted in the form of crossing points. The data acquired was computed by SDS software v2.3 (Applied Biosystems, Waltham, MA, USA) and subjected to subsequent analysis. The specificity and integrity of the PCR product was confirmed by a single melt curve peak. Relative expression was normalized to βactin and calculated according to the 2 −∆∆Ct approach [47]. The primer sequences used for detecting the expression of HIF1A, SYK, PGK1, and β-actin are given in Table S1.

Statistical Analysis
Data are presented as the mean ± SEM and statistically significant differences were identified with Student's t test as indicated in the figure legends. The difference in size, concentration and ZP of RB and non-RB were analyzed by the Welch t-test. Student's t test was used to compare the normalized mRNA expression levels and two-sided p < 0.05 was considered to be statistically significant.

Results
The demographic and clinical features of RB and non-RB subjects are summarized in Table S2. All the nine RB patients had intraocular unilateral RB with two belonging to Group D and 7 to Group E based on ICRB classification. The mean (median; range) age of RB children and controls were 3.6 (4; 1 to 9 years) and 5.4 (5; 4 to 8 years) respectively. Of the nine RB children, six were male and three were female, and amongst the five controls, three were male and two were female.

Functional Enrichment Analysis of Differentially Expressed mRNAs and Protein-Protein Interaction-Network of Eye-Related Genes in RB Serum-Derived Small EVs
The top 20 up and down DE mRNAs in RB sEVs with their functions are listed in (Tables S5 and S6)

Functional Enrichment Analysis of Differentially Expressed mRNAs and Protein-Protein Interaction-Network of Eye-Related Genes in RB Serum-Derived Small EVs
The top 20 up and down DE mRNAs in RB sEVs with their functions are listed in (Tables S5 and S6)  We also identified 39 DE mRNAs associated with eye development in RB sEVs (Tables 1 and 2). A PPI network with these genes was examined for deciphering the functional role of these protein interactions in RB pathobiology. The PPI network of RB sEVs was divided into three significant clusters ( Figure S2). Cluster 1 proteins were shown to be involved in 9-cis retinoic acid biosynthesis (p = 0.003) and retinol metabolic process (FDR = 3.58 × 10 −5 ). Cluster 2 was enriched with retina development in camera-type eye (GO:0042127). Phosphatidylinositol 3′-kinase (PI3K)-Akt (hsa04151), transforming growth factor (TGF)-beta (hsa04350), Oxytocin (hsa04921), Insulin (hsa04910) and FoxO (hsa04068) signaling pathways, and phototransduction (hsa04744) are the enriched KEGG terms for RB sEVs DE mRNAs. GSEA results revealed that the top significant gene sets enriched for RB sEVs are associated with ether lipid metabolism (p = 0.01), alpha linolenic acid metabolism (p = 0.02), and regulation of insulin like growth factor receptor signaling pathway (p = 0.04) ( Figure 3J-L). Regulation of DNA damage response signal transduction by p53 class mediator gene set is enriched for non-RB sEVs (p = 0.005) ( Figure 3M).
The significance of DE miRNA targets in RB was determined by functional enrichment analysis. The top GO and KEGG terms were represented in bar graphs ( Figure 4A-D). Up regulated miRNA targets are associated with regulation of protein serine/threonine kinase activity (GO:0071900), covalent chromatin modification (GO:0016569), and down regulated miRNA targets are related to positive regulation of cellular catabolic process (GO:0031331). Both up and down miRNA targets were enriched with RNA catabolic process (GO:0006401), and transcription co-regulator activity (GO:0003712), MAPK signaling pathway (hsa04010), PI3K-Akt signaling pathway (hsa04151), and proteoglycans in cancer (hsa05205).

miRNA-mRNA Regulatory Network Results
Based on the miRNA enrichment results, the top DE miRNAs having the highest targets and mRNAs targeted by maximum number of DE miRNAs were used as input for constructing experimentally validated miRNA-target interaction (MTI) network. This co-expression network unveiled the expression status of 35 miRNAs and their targets 19 mRNAs in RB sEVs ( Figure 4E). A total of 562 functional MTI pairs were generated (Supplementary File S1). The key DE mRNAs VEGFA, CCND1, E2F3, WEE1, c-MYC, HIF1A, and XIAP in the network were found to be regulated by highly dysregulated miR-NAs such as 16-5p, 106a-5p, 24-3p, 17-5p, 15a-5p, 181a-5p, 181b-5p, and 15b-5p. miRNA-mRNA network enrichment results revealed that the target genes in the network are linked to pathways in cancer, microRNAs in cancer and cell cycle (KEGG, p < 0.05), Cyclin D associated events in G1, ubiquitin specific processing proteases (reactome pathways, p < 0.05), PI3K-Akt Signaling, G1 to S cell cycle control and retinoblastoma gene in cancer (wiki pathways, p < 0.05).
The significance of DE miRNA targets in RB was determined by functional enrichment analysis. The top GO and KEGG terms were represented in bar graphs ( Figure 4A

miRNA-mRNA Regulatory Network Results
Based on the miRNA enrichment results, the top DE miRNAs having the highest targets and mRNAs targeted by maximum number of DE miRNAs were used as input for constructing experimentally validated miRNA-target interaction (MTI) network. This coexpression network unveiled the expression status of 35 miRNAs and their targets 19 mRNAs in RB sEVs ( Figure 4E). A total of 562 functional MTI pairs were generated (Supplementary File S1). The key DE mRNAs VEGFA, CCND1, E2F3, WEE1, c-MYC, HIF1A, and XIAP in the network were found to be regulated by highly dysregulated miRNAs such as 16-5p, 106a-5p, 24-3p, 17-5p, 15a-5p, 181a-5p, 181b-5p, and 15b-5p. miRNA-mRNA network enrichment results revealed that the target genes in the network are linked to pathways in cancer, microRNAs in cancer and cell cycle (KEGG, p < 0.05), Cyclin D associated events in G1, ubiquitin specific processing proteases (reactome pathways, p < 0.05), PI3K-Akt Signaling, G1 to S cell cycle control and retinoblastoma gene in cancer (wiki pathways, p < 0.05).

DE lncRNA Analysis and Functional Enrichment Analysis of Target Genes
The majority of the DE lncRNAs identified in RB sEVs, including the top ones, belong to long intergenic non-coding RNAs (LINCRNAs): LINC02499, LINC02773, LINC01416 and LINC00994, whereas OXCT1-AS1, HDAC2-AS2, ACSL6-AS1, and SLC8A1-AS1 are some of the antisense DE lncRNAs that are able to control their own sense gene expression. The DE lncRNAs are represented in a heat map ( Figure 5A). LncRNA target gene analysis revealed that of 3634, only 242 lncRNAs in RB sEVs have targets. A total of 769 regulatory lncRNA-miRNAs pairs, and 541 lncRNA-mRNA pairs comprised of 242 lncRNAs, 323 miRNAs, and 332 mRNAs were predicted (Supplementary File S2). The significantly up regulated AFAP-AS1, and down regulated MALAT1, GAS5, ZFAS1, and SNHG16 have the highest number of target interactions. In addition, lncRNAs that directly target cell cycle specific genes were also detected in RB sEVs (Table 3).   The enriched GO and KEGG terms for lncRNA-target genes are depicted in bar graphs ( Figure 5B

DE lncRNA Analysis and Functional Enrichment Analysis of Target Genes
The majority of the DE lncRNAs identified in RB sEVs, including the top ones, belong to long intergenic non-coding RNAs (LINCRNAs): LINC02499, LINC02773, LINC01416 and LINC00994, whereas OXCT1-AS1, HDAC2-AS2, ACSL6-AS1, and SLC8A1-AS1 are some of the antisense DE lncRNAs that are able to control their own sense gene expression. The DE lncRNAs are represented in a heat map ( Figure 5A). LncRNA target gene analysis revealed that of 3634, only 242 lncRNAs in RB sEVs have targets. A total of 769 regulatory lncRNA-miRNAs pairs, and 541 lncRNA-mRNA pairs comprised of 242 lncRNAs, 323 miRNAs, and 332 mRNAs were predicted (Supplementary File S2). The significantly up regulated AFAP-AS1, and down regulated MALAT1, GAS5, ZFAS1, and SNHG16 have the highest number of target interactions. In addition, lncRNAs that directly target cell cycle specific genes were also detected in RB sEVs (Table 3).

LncRNA-miRNA-mRNA Network Results
To reveal the regulatory role of lncRNAs on miRNAs and protein-coding genes associated with RB tumorigenesis, a lncRNA-miRNA-mRNA network consisting of 23 lncRNA, 64 miRNA, and 46 mRNA nodes with 203 edges (101 lncRNA-miRNA pairs and 102 lncRNA-mRNA pairs) was constructed ( Figure 5F and Supplementary File S3). We found that the protein-coding genes targeted by lncRNAs and miRNAs in the network were associated with the cell cycle, such as RB1, CDK6, cyclins (CCND1 and CCNE1), CDK inhibitors (CDKN1A, CDKN1B, CDKN1C and CDKN2A). Based on the degree, closeness and betweenness centrality, the hub lncRNAs (MALAT1, AFAP1-AS1, HOTAIR, NEAT1 and MEG3), miR145 and miR 101, and mRNAs (CDKN1A, EZH2 and ZEB1) were identified from network analysis ( Table 4). Hub RNAs was found to be involved in more regulatory interactions having key roles in network organization. The CluGo enriched functional KEGG terms for protein coding genes in the network are cell cycle, apoptosis and tumor related signaling pathways such as HIF1A, ErbB, and P53 ( Figure 5G). In addition, the most relevant biological process related to RB "aberrant regulation of mitotic G1/S transition in cancer due to RB1 defects" was also enriched for target mRNAs in the network ( Figure 5H). These results suggest that lncRNAs play a crucial role in mediating RB progression.

Quantitative Reverse Transcriptase-Polymerase Chain Reaction Validation
The mRNA expression levels of HIF1A (hypoxia-inducible factor 1-alpha), PGK1 (phosphoglycerate kinase 1), and SYK (spleen associated tyrosine kinase) were tested in five primary RB tumors and two control retinas by RT-qPCR. There was a statistically significant (p < 0.05) overexpression of SYK in three RB tumors and HFIA in one RB tumor compared to control retina ( Figure S4). However, the expression levels of PGK1 remained same in RB tumors and control retina.

Discussion
Serum-derived small extracellular vesicles (sEVs) can be a possible potential source of liquid biopsy for tumors such as RB, where routine tissue biopsy for diagnosis is not done due to the risk of extraocular tumor spread [48]. In this comprehensive study, we characterized RB serum sEVs coding and non-coding RNA content, and explored the miRNA-mRNA, and lncRNA-miRNA-mRNA regulatory interactions. Until now, only one previous study had analyzed small non-coding miRNAs in serum sEVs of RB patients using NGS technology [21].
An average of 5-6 × 10 11 particles/mL, with size 120-135 nm, and ZP of 11.0-12.6 mV was recovered from RB and non-RB serum samples. There was no difference in sEV concentration detected between the two groups. However, smaller size and lower ZP values were observed for single freeze-thaw cycled RB sEVs compared to non-RB. All the isolated serum sEV samples were positive for well-known exosome markers CD9, CD81, and TSG 101 [49]. To our knowledge, the physical properties of sEVs of young children (<5 years) have not yet been reported. However, studies on adult sEVs revealed that they have a diameter of 109 ± 4 nm, 1.8 × 10 11 ± 3.1 particles/mL and ZP of −9.80 to −21.1 mV [50]. Another group documented the presence of 1-3 × 10 12 sEVs/1 mL in adult serum [51]. Physico-chemical properties of sEVs can provide insights about intricate pathological processes in the course of the disease and have been evolving as potential factors in cancer diagnosis and monitoring treatment response [5,7,9]. Elevated levels of sEVs have been found in cancer patients compared to healthy donors [52][53][54]. Plasma membranes of cells have a negative surface charge, which is known to influence various biological processes such as cellular uptake and cytotoxicity [25]. Zeta potential (ZP), an indicator of colloidal stability, is influenced by the surface charge. The difference in ZP of RB and non-RB sEVs could be due to aberrant biological processes in RB patients. The reason for smaller size of RB sEVs compared to non-RB patients is not known. The literature states that the difference in vesicle size is not a crucial criterion, as sEVs might aggregate into larger vesicles, or multi-vesicular bodies can split into smaller vesicles. Moreover, it remains unknown whether one cell produces EVs of different sizes, if the difference in size reflects EV production by different cells, or whether vesicles with a common size share the same or different composition [55].
Identification of the most commonly deregulated genes for a specific cancer is of potential diagnostic interest. Several studies have reported that tumor-derived exosomes from several cancers reflect original tumor molecular signatures, and they share distinct transcriptome and proteomic profiles with respect to healthy controls [7,79]. However, high throughput data of non-coding RNAs related to RB as well as RB EVs are not available in the literature. A large number of RB patients' sEVs, and corresponding RB tumors must be analyzed to obtain a sound conclusion. GO, KEGG and GSEA analysis of DE mRNAs, DE lncRNA-target genes, and DE miRNA target genes were significantly enriched in biological processes related to angiogenesis, EMT, mitotic cell cycle regulation, chromatin organization, metabolic pathways, and with cancer-related pathways, such as the PI3K-Ak signaling, TGF-beta signaling, p53 signaling, Insulin signaling, and HIF1 signaling pathways.
PPI network analysis of 39 DEMs related to eye and retina development were shown to be involved in biological processes associated with 9-cis-retinoic acid (RA) biosynthetic processes, retinal rod cell differentiation and retinal cone cell development. Surprisingly, cone specific phosphodiesterase 6C (PDE6C) expression, which is restricted to the retina, pineal gland, and retina-derived tumors [80], was also detected in sEVs. However, their levels were high in normal sEVs compared to RB. The presence of PDE6C in sEVs indicates that eye EVs can enter into the circulation system. In support of our findings, previous studies reported that EVs can cross the blood-retinal barrier (BRB) and blood-brain barrier (BBB) and serve as a source of communication between the central nervous system (CNS) and peripheral system [81]. In addition, consistent with our data, earlier studies reported that RA induces cone-specific and rod-specific gene inactivation, and cell cycle arrest during the differentiation of RB cells, which implicates the role of RA signaling in RB development [82]. Moreover, NGS data of normal retinas, RB tumors and RB cell lines revealed that abnormal retina development could be involved in RB origin and progression [83,84].
Exploration of regulatory interactions between ncRNA-target mRNAs is crucial for elucidating ncRNA-mediated gene regulation in RB, as emerging data has provided evidence that epigenetic molecular events driven by RB1 loss are necessary for malignant phenotype [85]. We identified several miRNAs (miR 101-3P, 16-5p, 155-5P, 181a-5p, miR17/192, LET-7 clusters, and lncRNAs (HOTAIR, KCNQ1OT1, MALAT1, AFAP-S1) that directly target common dysregulated cell cycle genes in RB including RB1, CCND1, E2F3, WEE1, XIAP, CDKN1A, c-MYC, MYCN, BCL2, VEGEF, HIF1A, and P53. Although mutations in RB1 gene are responsible for most cases of RB, identification of miRNAs and lncRNAs that directly target RB1 in unilateral sporadic cases with RB1 + / + genotype tumors aid in uncovering novel non-genetic mechanisms of RB1 inactivation. Similar to our findings, previous studies reported that miRNAs, and lncRNAs negatively regulate translation of target cell cycle-related gene cyclins, E2F members, cell cycle inhibitors, and TGFB2 [86][87][88]. However, in vitro and in vivo functional studies must be done to validate/confirm the individual coding-noncoding interactions identified in RB sEVs. Along with this, screening these genes in sEVs and corresponding tumor from different cohorts of RB patients would reveal potential EV diagnostic markers.
Based on in silico analysis and functional studies, several studies have investigated regulatory interactions between coding and ncRNAs in RB [89][90][91]. However, in most cases, only ncRNAs and their target genes that are aberrantly expressed in other cancers were selectively picked up and studied for their role in RB. However, the detailed analysis of direct interactions between mRNA-ncRNA regulatory axis has not been much explored as miRNA and lncRNA have very short lifespans and their expression pattern is very dynamic in a given cell. These properties make it difficult to understand ncRNA regulation by wet lab experiments in a single study. Computational analysis of lncRNA-miRNA-mRNA interactions would provide a theoretical basis for the molecular mechanism of disease, and for the identification of potential therapeutic targets [92]. However, to our knowledge only a single study has demonstrated autophagy-related lncRNA-miRNA-mRNA regulatory networks in RB [93]. Thus, we investigated the direct interactions between lncRNA-miRNA, miRNA-mRNA and lncRNA-mRNA, and constructed lncRNA-miRNA-mRNA network. Based on the degree, closeness and betweenness centrality, MALAT1, HOTAIR, NEAT1, AFAP1-AS1, miR 145, and miR101 were identified as hub ncRNAs that play central roles in RB pathogenesis. Target protein coding genes for these hub ncRNAs are associated with cell cycle (RB1, c-MYC, cyclins, CDKs), cellular senescence, remodeling of extracellular matrix (MMP1, MMP2, and MMP9), and epigenetics (EZH2 and ZEB1). PI3k-AKT, HIF1, ErbB, and P53 signaling pathways are enriched terms. These are the known signaling pathways that drive RB tumorigenesis [57,94,95]. The hub genes EZH2 and ZEB are well-studied epigenetic regulators in RB. ZEB1 is an EMT transcription factor, usually repressed by RB1. This facilitates the epigenetic silencing of CDH1 (E-cadherin). The XIST-miR101-ZEB1 axis has been shown to be responsible for malignant properties of RB cells [70]. Our RT-qPCR data showed a high expression of epigenetic regulator, SYK and hypoxia-inducible transcription factor, HIF1A in RB tissues compared to control retinas. However, the expression status of SYK and PGK1 in RB sEVs and corresponding RB tissues were not matching. It is worth noting that the study included low number of sEV RNA samples sequenced per group (n = 3), and the presence of minor contaminants of serum in the analyzed sEV data due to technical challenges of EV extraction with the available commercial kits.

Conclusions
The serum EV RNA profiling in non-invasive RB tumors suggests that sEVs have the signatures of RB tumors. The cargo in the small EVs is involved in the epigenetic regulation of cell cycle, metabolism, and tumor-associated signaling pathways. Further validation is warranted to use them as prognostic biomarkers.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14174179/s1. Table S1: Primers for RT-qPCR. Table S2: Clinical characteristics of subjects included in the study. Table S3: The concentration, size distribution and zeta potential of serum small extracellular vesicles from retinoblastoma (RB) and non-RB subjects. Table S4: RNA sequencing details of serum small EVs from RB and non-RB subjects.  Table S9: Dysregulated microRNAs identified in RB serum small EVs and their expression in RB tumor tissues from review of literature. Table S10: Dysregulated lncRNAs identified in retinoblastoma (RB) serum small EVs and their expression in RB tumor tissues from a review of the literature. Figure S1: Research flow chart. Figure S2: Protein-protein interaction network of differentially expressed mRNAs of RB serum small EVs related to eye development. Figure S3: miRNA-target interactions in RB serum small EVs. Figure S4: Real time expression data of HIF1A, PGK1, and SYK File S1: Exosomal miRNA-mRNA network in RB. File S2: RB small EV lncRNAs and lncRNA targets, and File S3: RB exosomal lncRNA-miRNA-mRNA network.  Informed Consent Statement: Informed consent was obtained from all subjects/legal guardians involved in the study.

Data Availability Statement:
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest. The sponsors had no role in the design, execution, interpretation, or writing of the study.