Human retinal organoids release extracellular vesicles that regulate gene expression in target human retinal progenitor cells

The mechanisms underlying retinal development have not been completely elucidated. Extracellular vesicles (EVs) are novel essential mediators of cell-to-cell communication with emerging roles in developmental processes. Nevertheless, the identification of EVs in human retinal tissue, characterization of their cargo, and analysis of their potential role in retina development has not been accomplished. Three-dimensional retinal tissue derived from human induced pluripotent stem cells (hiPSC) provide an ideal developmental system to achieve this goal. Here we report that hiPSC-derived retinal organoids release exosomes and microvesicles with small noncoding RNA cargo. EV miRNA cargo-predicted targetome correlates with Gene Ontology (GO) pathways involved in mechanisms of retinogenesis relevant to specific developmental stages corresponding to hallmarks of native human retina development. Furthermore, uptake of EVs by human retinal progenitor cells leads to changes in gene expression correlated with EV miRNA cargo predicted gene targets, and mechanisms involved in retinal development, ganglion cell and photoreceptor differentiation and function.


hiPSC-derived 3D retinal organoids as a model of human retinal development. Retinal orga-
noids provide an in vitro model of human retina development that recapitulates well defined stages of retinogenesis, including early formation of a pseudostratified retinal neuroepithelium with its characteristic RPC proliferation and interkinetic nuclear migration, followed by lamination and temporally ordered fate specification 10 . For this study, we selected three developmental time points (day (D) 42, 63, and 90) based on the dynamics of cell differentiation and migration of retinal organoids, and that represent distinctive stages during retinal cell fate specification and lamination (Fig. 1). These stages were thoroughly characterized in our previous report 10  hiPSC-derived 3D retinal organoids release extracellular vesicles including exosomes and microvesicles. In order to determine whether hiPSC-derived 3D retinal organoids are capable of releasing EVs throughout their developmental process, we analyzed retinal organoid conditioned media at D42, D63 and D90 of differentiation. EVs were successfully isolated from the three developmental stages and analyzed using nanoparticle tracking analysis technology (NTA) to determine their size and release rate (Fig. 2). EV diameters were consistently within a range of 30-570 nm ( Fig. 2A-D) throughout the developmental stages, with average peak diameters shown in Fig. 2E. Diameter analysis demonstrated that, while EVs from retinal organoids contained both exosomes and microvesicles 1 , the majority of EVs (69.4 % +/− 6.7) across timepoints corresponded to exosomes with diameters of 30-150 nm 13 . Next, the concentration of released EVs for each time point was calculated, and the corresponding concentrations after normalization shown in Fig. 2F. A trend of increased EV concentration above control levels was visible for each time point suggesting ubiquitous release 14 .
EVs isolated from hiPSC-derived 3D retinal organoids were further characterized by immunogold detection of EV proteins on standard transmission electron microscopy (TEM) and scanning electron microscopy (SEM). Immunogold TEM images confirmed proper localization of canonical exosome (CD63) and microvesicle (Tsg101) transmembrane proteins at the surface of isolated EVs (Fig, 3A,C) 15,16 . Standard TEM images revealed a heterogenous population of spheroid and cup-shaped EV morphologies, with an average diameter near 110 nm, consistent with the known EV size and morphology and within the range size of exosomes ( Fig. 3D-F) 1 . EV confirmation was also performed via immunoblotting against the EV specific markers TSG-101 and CD63 in Supplemental Figure SF1. determine their profile in retinal organoids and corresponding released EVs. Hierarchical cluster analysis of miRNAs present in EVs versus their respective donor retinal organoids revealed clear differences in the miRNA profile between cells and EVs at all developmental time points analyzed (Fig. 4). Once again, this outcome is in agreement with previous studies, including the original report identifying the presence of miRNA cargo in extracellular vesicles 18 . To illustrate unique and common miRNAs in 3D retinal organoids and released EVs at the three developmental time points, Venn diagrams were generated (Fig. 5A). Released EVs at D42, D63 and D90 contained 12, 3 and 2 miRNAs with >5 reads per million (RPM) 19 , respectively. The majority of these miRNAs were also present in their respective retinal organoids, while a small number was detected exclusively in the EVs (Fig. 5A,B); of note, miR-4516 was found exclusively within EVs at all developmental times analyzed with no detectable levels in 3D retinal organoids (Fig. 5B). 3D retinal organoids at D63 contained the highest number of miRNAs identified, while the highest number of miRNA species in EV samples was found at D42 (Fig. 5A). miRNAs showed differential enrichment in EVs compared to their respective retinal organoids, with the most abundant miRNAs contained in EV cargo being generally different to the most abundant miRNAs detected in the 3D retinal organoids at the corresponding developmental time points (Fig. 5C). In addition, the diversity and expression level of miRNAs in both, retinal organoids and released EVs, showed a decreasing trend concomitant with developmental progress (Fig. 5B,C). A comprehensive list of all miRNAs identified in hiPSC-derived 3D retinal organoids at each developmental time point is presented in Supplemental Table ST1A.
In addition to miRNA, piRNA and tRNA species were also identified in 3D retinal organoids and derived EVs. Venn diagrams and comprehensive lists of all piRNA and tRNAs present in our samples are described in s 1 and 2, and Supplemental Tables 1B and 1C. The full dataset of RNAseq performed on iPSC-derived 3D retinal organoids and EVs at D42, D63 and D90 is available through GEO under the accession number GSE18252.
Confirmation of EV miRNA cargo and EV-mediated RNA transfer to human retinal progenitor cells. Considering the potential biological relevance of EV miRNA cargo, we set to further validate their expression in both hiPSC-derived retinal organoids and corresponding released EVs. For this purpose, we selected 6 out of the 12 miRNAs identified by NGS for validation by qPCR. Both, 3D retinal organoids and EVs were confirmed to contain the 6 selected miRNAs (Fig. 6A). To demonstrate cell-to-cell cargo transfer by EVs released from retinal organoids, we labeled EV RNA with the green (FITC) fluorescent dye Syto RNASelect and, EV lipid envelopes with red (TRITC) fluorescent dye PKH26, and performed live imaging on their uptake into non-labeled human retinal progenitor cells (hRPCs) 20,21 . Following 15 hours of co-incubation of unlabeled hRPCs and fluorescently labeled EVs, media was replaced with fresh media without labeled EVs to quantify  www.nature.com/scientificreports/ uptake, retention and reduction of EV-transferred RNA in hRPCs. At 30hrs (15 hours of co-incubation with labeled EVs followed by rinse and 15 hours incubation in fresh EV-free media) taken up and retained labeled EVs were visualized internalized in the cytoplasm of hRPCs, often with highest concentrations in regions near the putative endoplasmic reticulum and nuclear membrane ( Fig. 6B-F). Furthermore, live-cell fluorescence microscopy demonstrated active transport in central and peripheral cytoplasmic regions of internalized EVs in hRPCs (Supplemental Video SV2). EV-labeled RNA fluorescence increased significantly in hRPCs from 2 to15 hours, consistently with active EV uptake (Fig. 6G,H). After 15 additional hours in the presence of fresh EV-free media (30 hours total) transferred EV RNA signal began to decrease but remained significantly higher than 2 hours. Finally, EV RNA fluorescence in hRPCs significantly decreased from 15 to 72 hours (Fig. 6H).

Predicted retina targetome of EV-miRNA cargo correlates with GO pathways involved in mechanisms of retinogenesis.
Having confirmed the presence of miRNAs identified by NGS within hiPSC-derived retinal organoid released EV cargo, as well as the ability of these EVs to transfer their RNA cargo to neighboring cells, we pursued bioinformatics prediction of EV miRNA targetome as a first approximation to the potential role of EVs from retinal organoids. Putative miRNA target mRNAs were extracted from the MicroRNA Target Prediction And Functional Study Database (miRDB) 22 and filtered by the expression level of the target genes in retina according to published data 23 . A total of 3,523 microRNA-target mRNA pairs were predicted corresponding to the 12 miRNAs identified in EVs released from 3D retinal organoids. The number of genes targeted by each of the 12 miRNAs ranged from 9 to 1,069 (Fig. 7A). Interestingly, the majority of these genes (2,322) are putatively targeted by only one of the 12 identified miRNAs, while 467 of them are targeted by two miRNAs and 86 by 3-5 miRNAs (Fig. 7B). We then asked whether the mRNAs targeted by 4 and 5 miRNAs could represent key gene regulators of retinal cell differentiation at a particular developmental time point. As shown in Fig. 7C, network analysis revealed that most of these genes are putatively targeted by sets of different miRNAs in EVs at different developmental time points. This suggests that EV miRNA may act in a timedependent combinatorial manner, facilitating post-translational modification during retinal development 4 . We performed a second analysis including all the genes targeted by 2 or more miRNAs to evaluate the extent of this predicted combinatorial targeting pattern. The data in Fig. 7D suggested that targeting most commonly follows a combinatorial pattern, with different combinations of 2 to 5 miRNAs targeting different genes, rather than the same set of 2, 3, or more miRNAs targeting different genes. www.nature.com/scientificreports/ As a first approximation to the potential role that retinal EV miRNA cargo may exert during retina development, we performed Gene Ontology (GO) and pathway enrichment analyses. At a first level, unbiased GO analysis showed correlation with mRNA targets involved in mechanisms of general cell function, homeostasis, and differentiation (Supplemental Figure SF5). When analyzed by developmental time point, GO categories exhibited significantly higher enrichment at D42, for transcription, RNA binding, protein binding, positive regulation of cell proliferation and nuclear transcription factor binding. When analyzed by individual miRNAs, individual miRNAs showed significant enrichment for GO terms associated to specific cellular mechanisms (Supplemental Figure SF5). A second level of analysis revealed that predicted EV-miRNA targets correlated with GO-terms related to biological processes involved in retina development, including regulation of extracellular vesicles, nervous system development, eye and retina development, retina layer formation, neurogenesis, cell proliferation, and neuron differentiation, morphogenesis and migration (Supplemental Figure SF4). Notably, GO-term "neuronal stem cell population maintenance" showed a pronounced decreasing trend from D42 to D63 and D90 (values: 1.5, 0.72, and 0.78, respectively). Interestingly, another GO-term showing pronounced decreasing trend was "retinoic acid receptor signaling pathway" (D42, D63, and D90 values: 1.7, 0.67, and 0.63), a pathway well known to play an important role in retina development and differentiation 24 . Additional GO-terms that showed decreasing trends, although less pronounced included: regulation of cell cycle (D42, D63, and D90 Conversely, several GO-terms showed increasing trends of enrichment from D42 to D90. The most marked increasing trends were observed for retina development (D42, D63, and D90 values: 1.2, 2, and 2.1), negative regulation of neuron differentiation (D42, D63, and D90 values: 1.1, 1.9, and 2.1), cell morphogenesis (D42, D63, and D90 values: 1.4, 1.8, and 2), retinal layer formation (D42, D63, and D90 values: 1.7, 1.8, and 1.9), retinal ganglion cell axon guidance (D42, D63, and D90 values: 1.8, 2.2, and 2.3), cellular response to retinoic acid (D42, D63, and D90 values: 1.4, 1.8, and 1.9), and photoreceptor cell development (D42, D63, and D90 values: 1.5, 2.9, and 3.1). Remarkably, the GO-terms for retinal ganglion cell axon guidance and retinal photoreceptor development exhibited the highest enrichment, consistent with the dynamics of ganglion and photoreceptor cell generation and differentiation in human retinal development. Additional Venn Diagrams describing retinal organoid EV miRNA target GO function overlap with associated pathways is shown in Supplemental Figure SF6.
Several EV miRNAs showed high enrichment for target genes associated with developmental processes (Supplemental Figure SF4). hsa-miR-1298-5p appears involved in regulating development of neuronal cell projections together with negative regulation of cell proliferation, two mechanisms known to act conjointly in directing cell differentiation and morphogenesis. Notably, hsa-miR-1298-5p showed a unique enrichment for the GO category cytoskeleton (Supplemental Figure SF5). On the other hand, hsa-miR-99b-3p showed specific association to cell proliferation, cell division, cell differentiation, and retinoic acid receptor signaling pathway. Four of the EV miRNAs showed enrichment for retinal ganglion cell axon guidance and eye photoreceptor cell development. Of these, hsa-miR-4488 was highly enriched for eye photoreceptor cell development (score value of 11) compared to retinal ganglion cell axon guidance (score value of 0); hsa-miR-92b-3p and hsa-miR-204-3p were enriched for retinal ganglion cell axon compared to eye photoreceptor cell development (score values of 2.3/0 and 5.7/1.9 respectively); while hsa-miR-4516 showed similar enrichment for both GO terms (score value of 2.5) (Supplemental Figure SF4). Of note, these were not the only identified EV miRNAs with predicted target genes involved in photoreceptor and ganglion cell mechanisms. Ten out of the 12 miRNAs identified in EVs released by retinal organoids showed predicted targeted genes associated with specific aspects of photoreceptor and ganglion cell differentiation and function. Supplemental Table ST2 shows the temporal expression of EV miRNAs and targeted genes associated with ganglion and photoreceptor cell development. As shown in Supplemental Table ST2, hsa-miR-4488 was identified at the three differentiation time points analyzed, having a single predicted target gene, TULP1, which was associated to several GO categories relevant to photoreceptor differentiation and development. Interestingly, among the 12 retinal organoid miRNAs identified in this study, hsa-miR-4488 was the only one predicted to target TULP1. hsa-miR-92b-3p was identified in EVs isolated from retinal organoids only at D42, with one predicted target gene, associated to GO category retinal ganglion cell axon guidance (ROBO2) and five targeted genes associated to GO categories relevant to photoreceptor cellular mechanisms (MYO5A, PTPRK, GNAQ, MAP1B, PHLPP2). ROBO2 was also predicted to be targeted by miR204-5p from all 12 identified miRNAs. Coincidentally, hsa-miR-204-5p was also identified only at D42 but showed seven predicted target genes. Three of the predicted target genes were associated to the retinal ganglion cell axon guidance (ALCAM, ROBO2, and EPHB2), while PRKCI, GUCA1B, HCN1, and DNM2 were associated to categories relevant to photoreceptor cell development. Finally, hsa-miR-4516 was detected at all time points (D42, D63, and D90), and targeted genes involved in eye photoreceptor cell development (STAT3, GNAT1, CRB1) and retinal ganglion cell axon guidance (SEMA4F, EFNA5, EPHB3). The majority of EV miRNA predicted target genes were predicted to be targeted by hsa-miR-4516.
EVs from 3D Retinal Organoids mediate changes in hRPC differential gene expression that are associated to mechanisms of retinal development, ganglion and photoreceptor cell differentiation. To confirm the ability of EVs released from 3D retinal organoids to regulate gene expression in target retinal cells, we co-cultured EVs isolated from 3D retinal organoids at D42 with multipotent hRPCs in vitro for 96 hours, followed by RNAseq analysis of differential gene expression. D42 3D retinal organoids are characterized by early retinogenesis including the emergence of the first post-mitotic ganglion and photoreceptor cell precursors ( Fig. 1 and 10 ); their released EVs contain all the miRNA species identified in this study (described in Fig. 5) and likely additional early retinogenic factors. On the other hand, the hRPC population is mitotic, multipotent and undifferentiated at the time of co-culture with D42 EVs, thus with the potential to reveal EV mediated changes in gene expression associated with early retinogenesis. In EV-treated hRPCs, nine of the ten genes predicted to be targeted by at least one of the 3D retinal organoid derived EV-miRNAs (Fig. 7C) showed a trend of down-regulation (Fig. 8A). This downregulation trend was confirmed by qPCR for at least four out of five selected targeted genes, CCSER2, PVRL1, FAM117B, ILDR2 and CTDSPL (Supplemental Table ST3). Across all EV-treated hRPC genes downregulated at least a 2-fold change or greater (0.05 p-value) 33 genes were predicted to be targeted by one or more miRNA present in 3D retinal organoid EV-cargo (Fig. 8B). To map the predicted EV-miRNA downregulated gene targets in Fig. 8B, four databases were used -mirDB, targetscan, mirWalk and mirGate-with 11 out of the 12 3D retinal organoid EV-miRNAs having at least 1 predicted match among the EV-treated hRPC downregulated genes. Including genes predicted to be targeted by EV-miRNA, there were a total of 216 genes downregulated in EV treated hRPCs at a 2-fold change, 0.05 P-value (Fig. 8C, Supplemental Table ST3).
The 33 EV-miRNA targeted and downregulated genes correlated to GO functions including nuclear transport, translation, transcription, metabolism, development and homeostasis (Fig. 8D), and showed significant overlap to GO pathways identified through our predicted targetome analysis (Supplemental Figure SF5) 25,26 . Noteworthy, this is in agreement with our predicted targetome analysis which identified the NSL-bearing protein import into nucleus as a GO pathway enriched for D42 EV miRNAs (Supplemental Figure SF5). In addition to RGPD, the downregulation of TBC1 Domain Family Member (TBC1D3) significantly correlated to regulation of GTPase activity (also a GO pathway enriched for D42 EV miRNA predicted targetome); importantly, in neurons regulation of Rho GTPases directs differentiation, neurite growth, axogenesis and migration 27,28 . Additional relevant functions correlated to the EV-miRNA targeted downregulated gene POTE Ankyrin Domain Family Member E (POTEE, an ATP binding protein with roles in cell proliferation and regulation of metabolic activity 29 ) include neural nuclear development and retinal homeostasis. Furthermore, two of the downregulated genes, mitochondrially encoded dehydrogenase (MT-ND) and mitochondrially encoded ATP synthase (MT-ATP), have been shown to be differentially expressed in Notch1 positive progenitor cells in the developing retina 30 , where Notch1 regulates photoreceptor differentiation and contributes to ganglion cell fate specification.
In addition, a total of 38 genes showed increased expression at least a twofold change or greater, 0.05 P-value (Fig. 8E, Supplemental Table ST4). Interestingly, the cell function pathway most significantly correlated to the upregulated expression patterns of EV treated hRPCs is the Visual Cycle. The visual cycle in photoreceptors transduces photon detection into modification of electrochemical signaling at the first retinal synapse, thereby initiating the process of light perception 31 . An additional activated pathway correlated to upregulated hRPC genes is the nuclear receptor vitamin D receptor (VDR) and retinoid X receptor (RXR) (VDR/RXR) complex which remodels chromatin and modulates gene transcription 32 . Additionally, the Retinoate Biosynthesis II pathway is activated and correlated to upregulated genes in EV treated hRPCs (Fig. 8F). The full dataset of RNAseq performed on iPSC-derived 3D retinal organoid EV treated hRPCs is available through GEO under the accession number GSE182520.

Discussion
EV-mediated transfer of small RNAs is increasingly seen as a novel mechanism of genetic exchange between cells 33 modulating a range of cell behaviors. While relatively small (30 nm-1 μm), steady release rates and enrichment of small RNA species suggests EVs are capable of facilitating robust signaling 18 . Our group previously described functional EV mRNA transfer between mouse retina progenitor cells in vitro 34 . In this study, the characterization of hiPSC-derived retinal organoid EVs as carriers of miRNAs, piRNAs and tRNAs suggests a novel molecular mechanism for intercellular communication during differentiation of 3D retinal organoids specifically, and retinal development in general.
The use of hiPSCs to generate 3D retinal tissue is a novel approach for studying developmental processes. In this work, a well characterized and highly reproducible in vitro strategy is used for inducing hiPSCs into 3D retina, which mimics in vivo developmental processes, where normal retinal cell specification occurs 10 . In this system, hiPSCs recapitulate each of the main steps during human retinal development and form a laminated 3D retina tissue containing all major retinal cell types. Photoreceptors can reach an advanced stage of maturation, begin to form outer segments and develop photosensitivity. To better understand mechanisms regulating development of retinal organoids in vitro, we selected three developmental time points that represent distinctive stages of cell fate specification and lamination 10 and characterized EV release, morphology, small RNA cargo, cell-to-cell transfer, modulation of differential gene expression and predicted miRNA target genes correlated to retinal development.
Post-translational miRNA targeting regulates the expression of approximately 60% of mammalian genes, contributing to a range of cellular processes including proliferation, neuronal development and fate specification 35 . We identified twelve miRNA species in 3D retinal organoid EVs, some of which were present across different developmental stages. Only a subset of miRNAs identified in the retinal organoids was enriched in isolated EVs, and, though rarely, some miRNAs were detected exclusively in the EVs. These results are consistent with previous studies showing that miRNA enriched in EVs differ from those of the parent cells 36 . A number of studies www.nature.com/scientificreports/ suggest that cellular loading of miRNA species into EVs is not random, rather a mechanism of selective miRNA enrichment and export 37,38 .
Although a few studies have reported the presence and function of piRNA in the nervous system, the exact roles of piRNAs in developing nervous system tissues remain to be fully elucidated 39,40 . In the mammalian brain, piRNA is involved in silencing of retro transposons and neuronal piRNAs in Aplysia have been shown to modulate synaptic plasticity 41 . A number of studies suggest that piRNAs facilitate genome-wide DNA methylation, contributing to epigenetic modulation 42 . Some piRNA species have been detected in the mouse brain throughout postnatal brain development, whereas others are present primarily in adult stages of development 40 . To the best of our knowledge, this study is the first to detect piRNAs in EVs from 3D retinal organoids or native retinal tissue and suggests the potential for their involvement in retinal development.
An additional finding in this work was that hiPSC-derived 3D retinal organoid EVs were highly enriched for several tRNA species. This result is consistent with those of Bellingham et al. who indicated that neuronal cellderived EVs are enriched in tRNA species 43 . tRNAs have been detected in the mammalian brain and hypothesized to function similarly to miRNAs in post-transcriptional regulation 44 . Interestingly, effects in mitochondria tRNA genes are associated with pigmentary retinopathies and other retinal degenerations such as Usher syndrome 45,46 . In addition, mutations in a nucleotidyltransferase involved in tRNA processing, cause retinitis pigmentosa 47 .
Communication between retinal cells, via secreted factors, is essential for proper development and functioning of the eye. Neurotransmitters, chemokines, cytokines, hormones and other factors could propagate signals in a paracrine/endocrine manner. EVs released from retinal organoids provide an alternative route of cell-to-cell communication, introducing a new paradigm of retinal development. As novel mediators of intercellular communication, EVs act as paracrine effectors. We demonstrated that released EVs from human retinal organoids transferred RNA cargo to human retinal progenitor cells and, in turn, regulated the expression of genes involved in different mechanisms relevant to retinal cells.
The visualization of transfer and uptake of 3D retinal organoids EV lipid envelopes and enclosed RNA, support a model of EV mediated transfer of genetic information during retinal development. This aligns with our recent findings showing that mRNA and miRNA are packed into retinal progenitor cell EVs, and the content is functional following transfer to target retinal cells 34 . Additional recent work suggests that EV transfer of genetic material occurs between differentiated adult neural retina and multipotent retinal progenitor cells 20 . In this study, early differentiating 3D retinal organoid EVs were taken up into multipotent hRPCs. Neural EV uptake involves properties of the EVs and target cells, occurring via clathrin-mediated endocytosis 6,48 . Furthermore, the localization of internalized 3D retinal organoid EVs in hRPCs appears to align with other EV RNA transfer studies showing initial endosomal targeting followed by RNA content release into the cytoplasm 49 .
Differential expression analysis of hRPCs treated with retinal organoid EVs revealed changes in genes and signaling pathways associated with cell cycle, differentiation, retina homeostasis, and neuroprotection. Among downregulated genes in hRPCs that are targeted by one or more of the miRNAs identified in 3D retinal organoid EVs, those of special interest include CDK11, RGPD, MT-ND, MT-ATP, and EIF4EBP3. CDK11 is active in embryonic development involved in control of proper cell proliferation, cell-cycle rates and RNA-processing events 50,51 . The downregulation of CDK is also associated with cell-cycle exit and may indicate a transition of EV treated hRPCs toward a post-mitotic state 52 . Consistently, changes observed in RGPD protein expression levels may be associated to regulation of expression of downstream genes and fate determination 25,26 . Furthermore, changes in expression of MT-ND1 and MT-ATP6 may be related to photoreceptor and ganglion cell differentiation, as differential expression of these genes has been observed in retinal progenitor cells expressing Notch1, a transcription factor involved in photoreceptor differentiation and ganglion cell fate specification 30 . Worth noting, photoreceptors and ganglion cells are the two first cell types generated during the process of retina development and, what is more, are actively being generated in 3D retinal organoids at D42 of differentiation 10 . In addition, downregulation of MT-ND3 has been shown to provide neuroprotection for retinal cells by elevating the level of stress-related factors and reducing factors contributing to apoptosis 53 .
Interestingly, a number of genes upregulated 2-fold or greater in hRPCs treated with retinal organoid EVs are also associated to retinal differentiation and function. Upregulated genes of interest include RLBP1, FABP5, and APOC1. RLBP1 encodes cellular retinaldehyde binding protein (CRALBP) essential for visual cycle enzymatic pathway in rods and cones 54 . FABP5 is present in ganglion cells with predicted roles in neurite and axon growth 55 . APOC1 is expressed in retina and facilitates lipoprotein-mediated cholesterol transport and homeostasis 56 .
Finally, in our bioinformatic modeling of potential in vivo EV miRNA activity, several target gene pathways and functions were related to the eye and retina. Selected GO categories correlated to retinal organoid EV miRNA target genes included eye photoreceptor cell differentiation (GO:0001754), eye photoreceptor cell development (GO:0042462), photoreceptor outer segment (GO:0001750), regulation of retinal ganglion cell axon guidance (GO:0090259), retinal cone cell development (GO:0046549), retinal metabolic process (GO:0042574), photoreceptor cell maintenance (GO:0045494), negative regulation of photoreceptor cell differentiation (GO:0046533), photoreceptor outer segment membrane (GO:0042622) and photoreceptor connecting cilium (GO:0032391). Interestingly, our analysis showed that a number of the identified miRNAs have several predicted targets, as well as the redundancy of multiple miRNAs targeting the same gene at one or multiple time points. This suggests that the miRNAs might act through a combinatorial targeting mechanism to exert a synergistic modulation of differentiation and development.

Conclusion
This study demonstrates for the first time that human iPSC-derived 3D retinal organoids consistently release extracellular vesicles containing genetic cargo. It also provides the first sequencing analysis of sncRNAs contained in hiPSC-derived 3D retinal organoids and their EVs. Furthermore, this analysis includes three developmental Scientific Reports | (2021) 11:21128 | https://doi.org/10.1038/s41598-021-00542-w www.nature.com/scientificreports/ time points, providing snapshots of changing expression profiles. Small RNAs in EVs including miRNAs, tRNAs and piRNAs may serve as novel modulators of post-translational modification and biomarkers of retinal development and disease 57 . The uptake of EV RNA and consequent changes to differential gene expression in hRPCs suggests that cell-to-cell communication may occur in vivo as a novel mechanism of differentiation and development. This is further indicated by in silico modeling showing that 3D retinal organoid EV miRNA has predicted targets associated with retina cell differentiation and function.

Methods
hiPSC-derived 3D retinal organoid culture and conditioned media sample collection. A human episomal iPSC line derived from CD34+ cord blood was used in this study (A18945, Thermo Fisher Scientific; 58 ). This cell line has already been validated in its ability to generate light-sensitive 3D retinal organoids, and the process of differentiation thoroughly characterized 10 . The cell line was verified for normal karyotype and routinely tested for mycoplasma contamination by PCR. hiPSCs were maintained on Matrigel (growth-factorreduced; BD Biosciences)-coated plates with mTeSR1 medium (Stemcell Technologies) according to WiCell protocols. hiPSC-derived 3D retinal organoids were generated and maintained in culture as previously described 10 . RPE tissue growing attached to the 3D retinal organoids was dissected out to avoid the presence of RPE-derived EVs. Ten organoids per dish (a minimum of 5 dishes per conditions) were cultured in suspension. Cell culture media was normally changed twice a week. For the experimental EVs collections, cell culture media was fully replaced with 15ml of fresh media at D42, D63 and D90, and both conditioned media and 3D retinal organoids collected after 24h and 96h. Conditioned media was centrifuged (10 min at 2000 rpm) to remove cell debris and frozen at -80 °C until further analysis.
Nanosight analysis of EV size and concentration. EV size and number were assessed with a NanoSight NS500 system at D42, D63 and D90 in hiPSC-derived 3D retina cultures. After retinal organoids were cultured, conditioned media was collected and transferred to centrifuge tubes. The supernatant collection was identical to the previously described isolation of EVs except that the supernatant was used for NanoSight analysis before ultracentrifugation at 100,000x g. Non-conditioned control media from D42 and D63/90 were processed at the same time. Non-conditioned control media was used because it contains serum exosome concentrations and cargo that have to be removed from experimental conditions. Next, the supernatant was diluted 1:20 in phosphate-buffered saline (PBS), and 1 ml was used for NanoSight analysis. The NanoSight system uses a laser light source to illuminate nano-scale particles, which are individually detected as light-scattered points moving via Brownian motion. Polydispersity was quantified, and we used Nanoparticle Tracking Analysis (NTA) software 2.3 to track and determine the size of nanoparticles on an individual basis. The results are displayed as a frequency-size distribution graph describing the number of particles per ml.

Extracellular vesicle isolation.
EVs were isolated using a modified protocol based on previous work by Yuanet et al. [59]. hiPSC-derived 3D retinal organoids were cultured as previously described to collect EVs 10 . Briefly, conditioned media were collected, transferred to centrifuge tubes (polypropylene conical bottom) and centrifuged at 500x g for 10 min to pellet cells at room temperature; the supernatant was collected and centrifuged at 10,000x g for 20 min at 4°C (Beckman Ultracentrifuge, Rotor 60ti) to remove cell fragments and debris. The final supernatant was spun at 100,000 x g for 70 min at 4°C to pellet the EVs. EVs were re-suspended in PBS and stored at -80 °C for further analysis.
EV Immunogold TEM. EVs released from hiPSC-derived 3D retinal organoids were isolated from 25 ml of culture supernatant using differential centrifugation. Immunogold TEM was performed as previously described 34 . Briefly, 5 μl of resuspended 2% paraformaldehyde fixed EVs were put on glow discharged formvarcarbon coated nickel grids. After washing with PBS, the grids were incubated with 50 mM glycine/PBS for 3 min. The grids were blocked for 10 min with either 1% coldwater fish skin gelatin (Sigma) for the surface immunolabeling (CD63, Tsg101, Abcam), or 5% BSA, 5% goat serum, 0.1% cold-water fish skin gelatin and 0.1% saponin in PBS. Primary antibodies in blocking solution) were applied for 2 hours at room temperature. Control was prepared in the absence of primary antibodies. After washing with PBS, Nanogold-labeled Fab' anti-rabbit or antimouse (Nanoprobes, NY), or 5 nm, 10 nm gold conjugated goat anti-mouse antibodies (Ted Pella Inc. Redding, CA) were applied in the correlated antibody incubation buffer for 1 hour. The grids were then washed with PBS, fixed in 1% glutaraldehyde for 5 min. After thoroughly washed with distilled water, the grids were either directly go to methylcellulose embedding for 5 nm or 10 nm gold, or continue with silver enhancement for nanogold. For the silver enhancement, the grids were washed with 0.02 M sodium citrate (pH 7.0), and performed silver enhancement in the dark using HQ Silver enhancement kit (Nanoprobes, NY) at room temperature for 8 min.
After washing with distilled water, the grids were contrasted and embedded in a mixture of 3% uranyl acetate and 2% methylcellulose in a ratio of 1-9. Stained grids were examined under Philips CM-12 electron microscope and photographed with a Gatan (1 k × 1 k) digital camera. All antibodies were purchased from Abcam, USA, Anti-Tsg101 (1:200) and Anti-CD63 (1:200) were diluted according to manufacturer's instruction.
EV SEM. EVs were isolated from 25 ml of culture supernatant using differential centrifugation. EV SEM was performed as previously described 34

3D retinal organoid and EV NGS expression analysis.
Because of the sequence homology between some bovine and human miRNA, piRNA and tRNA species, the quantification of each in the cell culture or EV samples could be affected. Therefore, we removed reads from the cell culture and EV libraries if the same reads were observed in the control samples. Then, for each sample, we used the same number of reads used for alignment as the library size to normalize the expression of RNA. We calculated RPM for each RNA by dividing the read count mapped to each RNA by the library size and scaled the library size to one million reads. To determine which RNAs were expressed above a certain threshold (e.g., > 5 RPM) in cell culture or EVs at a time point, we used the following criteria: (1) the RNA is expressed >5 RPM in both replicates, and (2) the RNA is expressed ≤5 RPM in both control replicates. These criteria were applied in all analyses and described as "expressed >5 RPM" in the text. To visualize the expression patterns between samples, we log2 transformed (log2 (RPM +1), with a pseudo count of one added to avoid log-zero) the RPM values. We combined the replicates by calculating the average log2-RPM of each RNA between the two replicate samples. These values were then standardized across samples (mean 0 and s.d. 1) for visualization in heatmaps. In each heatmap, we ordered the RNAs by completelinkage hierarchical clustering (Euclidean distance). For each Venn diagram (Figs. 5, 6, 7), the following thresholds were used to define "expressed": >5. Those RNAs expressed in the corresponding control samples were removed. Only the RNAs that were expressed in both batches and not expressed in controls were used to generate the Venn diagram.
For tables of Figs. 5, 6, 7,: for each RNA type (miRNA, piRNA, or tRNA), RNAs that were expressed in EVs but not in controls were extracted to generate tables. These RNAs corresponded to those expressed in both cells and EVs (the overlap) and expressed only in EVs in the Venn diagrams. NGS expression analysis was performed by the New York Genome Center.
For heatmaps of Figs. 5, 6, and 7: The RPM values of the two batches were log2 transformed (with a pseudo count of 1 added) and averaged. Then, the data were standardized across time and library types for each gene. We selected RNAs expressed in EVs (shown in the Venn diagrams) at any time point. The standardized data of these RNAs were then clustered on the gene dimension (hierarchical clustering based on Euclidean distance) and rendered as heatmaps. qPCR analysis of miRNA in 3D retinal organoid EVs. EVs 22 and filtered by the expression level of the target genes in retina according to published data 23 . The functional enrichment analysis for miRNA target genes was done using DAVID tool. The fold change was used to measure the enrichment of the specific gene ontology term in the target gene list compared to the retina transcriptome. The heatmap was drawn using matlab to show the enrichment fold change of the interested terms in different conditions. Statistical analysis. Statistical analyses were performed using GraphPad Prism 7.0 software and shown as mean ± SD. Three biological replicates were used for analysis of cell EV uptake. The statistical significance of the difference was determined using unpaired one-sided Student's t-tests. Statistical significance was defined as P less than 0.05. *P<0.05. Bioinformatics and NGS expression analysis are described in their corresponding sections.