Improved Retinal Organoid Differentiation by Modulating Signaling Pathways Revealed by Comparative Transcriptome Analyses with Development In Vivo.

Stem cell-derived retinal organoids recapitulate many landmarks of in vivo differentiation but lack functional maturation of distinct cell types, especially photoreceptors. Using comprehensive temporal transcriptome analyses, we show that transcriptome shift from postnatal day 6 (P6) to P10, associated with morphogenesis and synapse formation during mouse retina development, was not evident in organoids, and co-expression clusters with similar patterns included different sets of genes. Furthermore, network analysis identified divergent regulatory dynamics between developing retina in vivo and in organoids, with temporal dysregulation of specific signaling pathways and delayed or reduced expression of genes involved in photoreceptor function(s) and survival. Accordingly, addition of docosahexaenoic acid and fibroblast growth factor 1 to organoid cultures specifically promoted the maturation of photoreceptors, including cones. Our study thus identifies regulatory signals deficient in developing retinal organoids and provides experimental validation by producing a more mature retina in vitro, thereby facilitating investigations in disease modeling and therapies.

Introduction

Generation of induced pluripotent stem cells (iPSCs) from somatic cells and their differentiation into diverse tissues in 3D organoid cultures have ushered a new era of renaissance in biology and medicine (Sasai, 2013). Self-organizing organoids contain multiple cell types patterned into a tissue-specific architecture and resemble in vivo counterparts (Lancaster and Knoblich, 2014). Stem cell-based technologies have enabled investigations into organogenesis, disease mechanisms, and design of novel therapies by recapitulating physiological events in vitro. However, morphological and/or functional maturity of tissues in vitro remains a major challenge (Soldner and Jaenisch, 2018), thereby limiting their potential for exploring biological complexities. Because organogenesis entails interaction(s) of intrinsic mechanisms with specific microenvironments (Yin et al., 2016), elucidation of molecular and regulatory networks that distinguish differentiation in vivo from that of organoids in vitro is critical for improvement of long-term 3D culture systems.

The mammalian retina is comprised of five major types of neurons: retinal ganglion cells (RGCs), and amacrine, horizontal, bipolar, and photoreceptor cells (rods and cones). These cell types are organized in a laminar structure consisting of three cell layers connected by two plexiform layers. All neurons and Müller glia originate from pools of retinal progenitor cells that pass through distinct states of competence in a spatiotemporal manner (Agathocleous and Harris, 2009, Bassett and Wallace, 2012, Hoon et al., 2014). This precise control is largely guided by intrinsic genetic programs in the retina (Cayouette et al., 2003); however, extrinsic signaling factors also contribute to functional differentiation (Cepko, 1999, Levine et al., 2000). Gene expression patterns that underlie retinal development and define cellular identity are stringently regulated by combinatorial actions of specific transcription factors (TFs) in concert with signaling molecules. For example, Pou4f1-3 (Brn3a-c) directs RGC differentiation (Sajgo et al., 2017), homeodomain protein CRX signifies the birth of photoreceptors (Gregory-Evans et al., 2013), and basic motif leucine zipper protein NRL specifies rod photoreceptor cell fate (Swaroop et al., 2010). Similarly, retinoic acid, fibroblast growth factor (FGF), and thyroid hormone modulate the expression of specific genes during photoreceptor differentiation (da Silva and Cepko, 2017, Khanna et al., 2006, McFarlane et al., 1998, Roberts et al., 2006).

Seminal studies have established fundamental principles for producing mouse and human retinal organoids from pluripotent stem cells (Eiraku et al., 2011, Meyer et al., 2011, Nakano et al., 2012). Retinal organoids in culture closely mimic fundamental stages in retinogenesis, demonstrating appropriate apical-basal and neuronal polarities, differentiation, and patterning of major cell types into a laminated structure (Volkner et al., 2016, Zhong et al., 2014). Furthermore, the use of bioreactors can facilitate the growth and differentiation of retinal organoids (DiStefano et al., 2018, Ovando-Roche et al., 2018). Transcriptome analyses of photoreceptors derived from organoids reveal cell-type-specific molecular signatures (Chen et al., 2016, Kaewkhaw et al., 2015, Welby et al., 2017). However, most studies have used RNA analysis and/or immunohistochemistry (IHC) profiles of markers to examine the presence of distinct cell types, and the precise temporal sequence of specification in retinal organoids is poorly understood. Despite technical advances, retinal organoids fail to reach functional maturation and eventually degenerate. Potential regulatory blocks of retinal maturation that may exist in vitro are currently unclear.

Advent of next-generation sequencing (NGS) technologies has expedited global profiling of retinal gene expression by RNA sequencing (RNA-seq) and permitted construction of regulatory networks (Aldiri et al., 2017, Chaitankar et al., 2016). In this study, we performed comprehensive and comparative transcriptome analyses to elucidate molecular distinctions between retinal development in vivo and in iPSC-derived organoids. We show that lineage specification in organoids is largely concordant with in vivo development, yet neural retina in organoids exhibits temporal distinctions in gene expression patterns and signaling pathways. Addition of docosahexaenoic acid (DHA) or FGF1, the two components missing in organoid cultures, facilitated photoreceptor differentiation and maturation. Temporal transcriptome dynamics can therefore provide valuable insights for improvements in long-term organoid cultures in vitro.

Results

Transcriptomes of Developing Retina In Vivo and in Organoid Cultures

Fate specification and functional differentiation proceed sequentially over a long period during mouse retinal development, with neurogenesis beginning as early as embryonic day 11 (E11) and formation of a morphologically mature retina by or after postnatal day 21 (P21) (Figure 1A). The birth of rod photoreceptors (indicated by expression of Nrl) overlaps with all other cell types (Akimoto et al., 2006). Though no cell proliferation is detected after P10, synaptic integration and retinal function is complete around/after P21. Development of retinal organoids from mouse iPSCs reveals somewhat altered temporal pattern of cellular differentiation (Figure 1B). These organoids demonstrate appropriate apical-basal and neuronal polarities, laminated organization, and synaptic protein expression comparable with the mouse retina, yet RGCs are not observed after differentiation day 25 (D25) and the organoids begin to degenerate around D32 (Figures 1C and S1).

Figure 1

Study Design and Quality Control

(A) Timeline of major developmental events used for mouse retina transcriptomes. Colored bars correspond to cell genesis and peak of cell birth. Cartoon of the retina (adapted from Yang et al., 2015) depicts major cell types.

(B) Representative brightfield images and time points for analysis of major events during mouse retinal organoid differentiation.

(C) Marker genes for distinct cell types, shown for in vivo retina at P8 and retinal organoid at D32, with the exception of the retinal ganglion cell (RGC) marker shown at D25. Müller glia and rod photoreceptors were stained with SOX9 (red) and RHO (green), respectively. Bipolar cells are immunostained with CHX10 (green) and ON bipolar cells (including all rod bipolars) with PKCα (red). Cone and amacrine cells are labeled with OPN1SW (red) and PAX6 (green), respectively. Horizontal cells and RGCs are respectively indicated by staining with CALB (red) and BRN3A (green). Cell nuclei are visualized with DAPI (blue). Arrowheads indicate relevant immunostaining.

(D) Principal-component analysis (PCA) of the expressed CPM values for each sequencing sample. Percentages indicate experimental variance assigned to each PC.

(E) Pearson correlation plots calculated using the expressed gene log2 CPM values of DE genes from the in vivo dataset. See also Figures S1 and S2, and Table S1.

We generated transcriptome profiles of mouse retina at 12 developmental stages (from E11 to P28) and of iPSC-derived retinal organoids at 10 stages (D0 to D32) (see Figures 1 and S2 and Table S1). Principal-component analysis using normalized, log2 CPM (gene level counts per million mapped reads) values of the protein-coding genes revealed developmental time as the primary variance (PC1) in the two datasets (Figure 1D). Pearson correlation plots of the in vivo data showed a major shift in transcriptome dynamics between P6 and P10 (Figure 1E), consistent with rod photoreceptor transcriptome profiles (Kim et al., 2016). Curiously, transcriptome profiles of retinal organoids did not reveal any such shift in gene expression patterns during differentiation, suggesting impediments to a normal developmental trajectory observed in vivo.

Differentially Expressed Gene Clusters during Retinal Development In Vivo and in Organoids

We then identified co-expression patterns within the two datasets to streamline comparative analysis. We performed hierarchical clustering of 12,718 genes that exhibited significant change in expression during mouse retinal development and determined 36 optimal clusters (Figures 2A and S3A), which were further sequestered into seven super clusters (SCs). Functional profiling of each SC and individual clusters (Figure S3B) uncovered key biological pathways and genes associated with specific stages of retinal development (Figure 1A). SC1 and SC2 were enriched for genes that exhibited monotonic decreases in expression during development. SC1 was enriched for genes associated with signaling pathways (e.g., Wnt, Hippo) implicated in proliferation and lineage commitment (Aldiri et al., 2013, Lee et al., 2018, Liu et al., 2003), whereas SC2 contained cell-cycle genes that are suppressed when cell proliferation nears completion. SC3 was enriched for synapse genes showing increased expression starting at E12/E14 (inner retina differentiation), whereas higher expression of SC4 genes was observed postnatally (phototransduction, ligand receptors). The Reactome pathway “Transmission across chemical synapses” was enriched in both SC3 and SC4, specifically in clusters 15, 16, and 22 (Figure S3C). Cholinergic, GABAB, and GABAA (Gabrb3 and Gabrg2/3) receptors exhibited expression embryonically, but a subset of GABAA (Gabra1/4) and GABAC (Gabrr1/2) receptors was also detected postnatally. SC5 primarily relates to metabolism. SC6-SC7 included TFs involved in differentiation of RGCs (e.g., Pou4f1/2) and other early-born neurons (e.g., Ascl1, Neurog2, and Prdm1) indicative of their transient requirement. This may also reflect “dampening” due to rod dominance in late development.

Figure 2

Molecular Characterization by Intra-dataset DE Analysis

(A) A total of 12,718 DE genes from the in vivo dataset were clustered by their Z score (36 clusters, showing an average representative from each) and grouped into seven super clusters (SCs) of broad expression similarity (left panel). Pathway analysis of the top enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of each SC is shown in the center panel. Expression heatmaps (right panel) show selected genes from representative pathways (bolded in center panel) for each SC enrichment.

(B) A total of 12,979 DE genes from the retinal organoids dataset were clustered by their Z score as in (A) and grouped into seven SCs of broad expression similarity. The heatmap uses the same Z score scale as in (A) (left panel).

(C) Top enriched KEGG pathways of each retinal organoid SC are shown in the left panel. Expression (log2 CPM) heatmaps (right panel) show selected genes from a representative pathway (bolded in the left panel) for SC2, SC4, and SC7 enrichment.

(D) Expression heatmaps of selected marker genes for photoreceptors, bipolar cells, and the early eye-field TFs in the retina. Expression values are represented as log2 CPM values and use the same scale as in (A, right panel).

Functional gene enrichment uses cumulative hypergeometric p values. See also Figures S3 and S4, and Tables S2 and S3.

Transcriptome analysis of developing organoids identified a similar number of differentially expressed (DE) genes (12,979) as in the retina in vivo. We performed analogous hierarchical clustering (Figure 2A) and generated 36 clusters sequestered into seven SCs showing comparable expression patterns (Figures 2B and S4A). Many pathways identified in the functional enrichment analysis from the in vivo retina were recapitulated in the organoid data (e.g., Wnt and Hippo signaling, cell cycle) (Figures 2C, S4B, and S4C). We then looked at cell-type-enriched genes to assess how well the retinal organoids recapitulated in vivo development (Figures 2D and S4D). As predicted based on IHC data, expression of RGC-specific genes including Pou4f1/2/3 was barely detectable after D22 in organoids. The horizontal and amacrine cell genes demonstrated a significant shift in organoids relative to the in vivo retina, whereas cone genes were detected at lower levels and transiently in organoids. Most other early-born cell marker genes (e.g., Isl1, Onecut1/2, and Prox1) showed high expression in retinal organoids at or after D10, equivalent to in vivo retina at E12/E14. Expression of marker genes for the late-born cells (e.g., Rho, Prkca, Gabrr1, and Grm6) in D32 organoids appeared similar to P6/P10 in vivo. The homeobox eye-field TF genes are highly expressed at the earliest times in the in vivo data and are among the first genes upregulated in the D4 and later organoids, with the exception of Tbx3. Curiously, Tbx3 is expressed in the undifferentiated iPSCs (D0) at 96 CPM, going down at D4 (2.9 CPM) and then upregulated at D7 (19.8 CPM), suggesting a “reset” in expression during early organoid differentiation.

Parallel Developmental Trajectory of Retina In Vivo and in Organoids

A comparative analysis of retinal transcriptomes in vivo and in organoids revealed a significant overlap of 9,599 genes (Figure 3A, top), showing similar patterns of expression (Figure 3A, bottom). However, we also identified a significant proportion of DE genes in one or the other dataset. Further examination of DE genes in the in vivo retina but not the organoids (3,119 genes) indicated a similar pattern of expression (Figure 3A, left) but with high variability among replicates, resulting in higher p value and elimination from the organoid analysis. These genes are enriched in the catabolism and mitochondria-related pathways (Figure 3B). DE genes unique to retinal organoids (3,380 genes) also showed high variability across replicates (Figure 3A, right). Many of these are enriched for early embryonic development and transcription regulation pathways that were altered at the beginning of organoid differentiation (SC1-SC3) and remained relatively unchanged at later stages (SC7).

Figure 3

Comparison of the Intra-dataset DE Analysis

(A) Euler diagram of the significantly DE genes within each dataset (top) and expression heatmaps of the genes in each section of the diagram in both datasets. Heatmaps for the 3,119 genes in the in vivo section (left) and the 9,599 genes in the overlap section (bottom) are arranged relative to the in vivo retina (Figure S3A), whereas the 3,380 genes in the retinal organoid section are arranged relative to the retinal organoid (Figure S4A). SC membership is indicated relative to the dataset from which the clustering order was employed.

(B) Top enriched KEGG pathways for genes belonging to the in vivo retina and retinal organoid sections of the Euler diagram from (A).

(C) Gene membership of the retinal organoid significantly DE genes relative to the in vivo retina SCs. Genes from the Euler diagram overlap in (A) are indicated in red, whereas genes found only in the in vivo retina analysis are in blue.

(D) Gene expression heatmaps of MAPK cascade and lipid localization pathways, which are enriched in the SC4 cluster of in vivo retina but not in the organoid analysis.

(E) Comparison of enriched Reactome pathways from SC functional gene analysis. Nodes represent pathways shared in both analyses (gold), exclusive to in vivo retina (red), or exclusive to retinal organoids (blue). The circular hierarchical pathway structure represents each root node as a primary functional group in the Reactome. For increased clarity, two plots were generated.

Functional gene enrichment uses cumulative hypergeometric p values. See also Figure S5, Tables S2 and S3.

To further understand distinctions between retinal development in vivo and in organoids, we focused on the representation of in vivo SCs of DE genes in organoids (Figure 3C). We identified greater under-representation in organoids of SC2 and SC4 genes, which are associated with the transcriptome transition at ages P6–P10 in vivo. Further examination of gene ontology (GO) pathways enriched in SC4 of the in vivo retina, but not in the organoid data, revealed mitogen-activated protein kinase (MAPK) cascade and lipid localization (Figure 3D). Reactome pathway analysis of in vivo SCs revealed enrichment of most biological pathways in the retinal organoid data as well (Figure 3E). The pathways uniquely detected in the in vivo retina transcriptome primarily centered on metabolism and extracellular matrix organization, whereas those only in retinal organoids were associated with RNA metabolism and transport of small molecules. Interestingly, signal transduction-associated pathways exhibited most variability between the two transcriptome datasets.

Altered Dynamics of Transcriptional Regulatory Network(s) in Organoids

We then performed open-ended dynamic time warping (DTW) to align time series expression data between the organoid and in vivo datasets using DE genes from the in vivo transcriptomes. As indicated by cell-type-specific genes (Figures 2D and S4D), temporal alignment of the two datasets revealed a match of D10–D15 organoid transcriptome to the E14 data in vivo and thereafter a linear equivalence of D15–D32 organoids with E14 through P6/P10 in vivo (Figure 4A). Negligible similarity was evident between iPSCs (labeled as D0) and the examined in vivo transcriptome data.

Figure 4

Temporal Alignment of Datasets and Direct, Paired DE Analysis

(A) DTW analysis comparing the temporal alignment of each dataset. CPM values from the DE genes of the in vivo dataset were used to generate the local cost matrix (LCM) depicted in the heatmap. The bold line indicates the minimum global dissimilarity value at each retinal organoid time point.

(B) Significantly DE TFs (red dots) identified by the paired DE analysis. Significance thresholds being fold change ≥1.5 and false discovery rate <1% (dashed lines). Gene labels indicate those having fold change ≥2.

(C) Enriched KEGG pathways identified in the paired DE analysis of retinal organoids with respect to the in vivo dataset.

(D) Functional gene set enrichment analysis (FGSEA) enrichment for four enriched pathways.

(E) Gene expression values (see Table S2) in the “Glycolysis” (top panel) and “Oxidative phosphorylation” (bottom panel) pathways.

(F) Genes from “Transmission across chemical synapse” KEGG pathway.

(G) Gene expression heatmap of MAPK pathway KEGG pathway downregulated in retinal organoids. (H) Gene expression of signaling molecule Fgf1 and the MAPK pathway small GTPase, Kras.

(I) Schematic of MAPK signaling pathway, showing significantly downregulated genes identified in the retinal organoids dataset (blue).

Functional gene enrichment uses cumulative hypergeometric p values. See also Figure S6 and Tables S4 and S5.

To gain insights into arrested development in retinal organoids at equivalent days P6–P10 in vivo, we performed a paired (blocking for the time factor) DE analysis of the linearly aligned time points revealed by DTW. This analysis allowed a direct comparison of datasets without confounding genes associated with early iPSC differentiation or maturation of the retina in vivo. We took five time points from each dataset, D15 through D32 in the organoid dataset and the corresponding DTW-matched time points in the in vivo retina (E14–P10). We focused first on TFs (Figure 4B) and identified significantly lower expression of Onecut2/3, Esrrb, Epas1, Foxg1, Irx2/4/5/6, and Rorb, among others, in retinal organoids. RGC TFs (Pou4f1/2/3 and Isl1) showed low expression at later stages of organoid differentiation, whereas several nuclear receptor TFs (Rarb/g, Rxrg, Nr2f2, and Nr4a1) are dramatically reduced in vivo and demonstrated consistently higher expression values in organoids.

We then performed functional genomic mRNA (FGM) profiling using all expressed TFs in each dataset and generated independent co-regulation networks. We focused on the top 50 most connected hubs from each network and identified retina-specific TF subnetworks (Figure S5A). We observed retina-specific connectivity of Plagl1 and Hes1 factors in vivo (Figure S5B). Plagl1 showed higher expression in later stages of organoid development. Hes1 expression was sharply decreased between P2 and P6 in vivo, but the reduction was more gradual in organoids. Zinc-finger proteins of the Zic family showed discrete connectivity patterns in the two networks (Figure S5C), displaying high expression in late stages of retinal organoids and low levels in the maturing retina in vivo. The POK (POZ and Krüppel)/ZBTB (zing finger and BTB) family of TFs also demonstrated distinct expression dynamics (Figure S5D). In retinal organoids, Zbtb7c generated a specific hub and was upregulated at D25–D32, whereas little to no expression was detectable in vivo.

Dysregulation of Maturation Pathways during Organoid Differentiation

Functional profiling of 1,790 genes from the paired DE analysis revealed enrichment of axon guidance (Figure 4C), synaptic transmission, phototransduction, and metabolism. The cAMP, calcium signaling, and circadian entrainment pathways were similarly enriched (Figure S5E). An unbiased approach using functional gene set enrichment analysis further validated dysregulation of genes associated with phototransduction, metabolism, and synaptic transmission (Figure 4D) in developing organoids. Notably, key genes within glycolysis and oxidative phosphorylation pathways were enriched at all time points of retinal differentiation in organoid culture (Figure 4E), suggesting a higher metabolic requirement. In concordance, the “Ribosome” and “Biosynthesis of Amino Acids” pathways were also upregulated in retinal organoids.

Major developmental events during P6–P10 transition in vivo include formation of the synapses, especially in the outer plexiform layer, and photoreceptor outer segment morphogenesis. By D32, genes associated with visual perception, specifically phototransduction, demonstrated delayed expression in retinal organoids, along with lower expression of cone-specific genes (Opn1sw, Gnat2, and Pde6h) (Figure S4F). Many genes involved in synaptic transmission showed higher expression in the in vivo dataset compared with the organoids (Figure 4F). Taken together, our analyses suggest arrested maturation of photoreceptors and synaptogenesis in organoids before their degeneration after D32.

Functional enrichment analysis exhibited significantly reduced expression of genes involved in calcium signaling, cAMP signaling, circadian entrainment, and MAPK pathways in organoids. Lower expression of many DE genes from these linked pathways included channels, receptors, and signaling molecules (Cacna1a/f, Cnga1, Grin2a, Grin3a, Gnai1, Gnb4, Gng4, Adcy7, and Camk2a/d) (Figure S5E). Furthermore, we observed downregulation of two key TFs (Gli1 and Hhip) involved in the Sonic Hedgehog pathway (Figure 4B). In concordance, downregulation of the MAPK/ERK pathway, an important target of cAMP/Ca2+/circadian signaling, was evident in retinal organoids (Figure 4G). We also noticed delayed expression of voltage-gated calcium channel genes (Cacna2d3, Cacna1f, and Cang4/5) and dysregulation of Ras genes Kras and Fgf1 (Figure 4H). Interestingly, Fgf1 is a known stimulator of the MAPK pathway and its expression is specifically upregulated at the critical transition period in developing retina at P6–P10 but is undetectable in retinal organoids. Thus, our data suggest aberrant MAPK and FGF1 signaling as a potential cause of arrested maturation in retinal organoids (Figure 4I).

Recapitulation of Altered Pathways in Human Organoids

To examine whether the altered pathways observed in the mouse organoids are unique to mouse or applicable to human organoid maturation, we analyzed recently published developmental human RNA-seq datasets consisting of fetal (Hoshino et al., 2017) and adult retina (Ratnapriya et al., 2019), and iPSC-derived organoids (Eldred et al., 2018). We examined whether the pathways perturbed in the mouse paired DE analysis were also altered in the human data (Figure S6A). Some of the MAPK pathway genes showed higher expression and others lower in the human retinal organoids compared with the in vivo data. However, unlike mouse organoids, KRAS and FGF1 were expressed in the human retinal organoids. The cAMP, circadian rhythm, and calcium signaling pathways were downregulated in the human retinal organoids as compared with the in vivo data. Visual perception genes demonstrated higher, and axon guidance genes lower, expression in the human organoids. The genes associated with glycolysis and oxidative phosphorylation showed higher expression in organoids, as in the mouse organoids (Figure S6B).

A paired DE analysis was performed between the human fetal retina and the organoid data using DTW of 118 retina-centric genes (Table S3) to align the time points between the datasets (Figure S6C). We then clustered 3,376 significantly DE protein-coding genes (Figure S6D) and determined dysregulated pathways in the retinal organoids (Figure S6E). Glutamatergic synapse and circadian entrainment pathways were significantly lower in expression in the retinal organoids (cluster 5). Cluster 12, which contained upregulated genes in organoids, was enriched in glycolysis and hypoxia-inducible factor-1 signaling genes, as in mouse retinal organoids. Although the in vivo samples suffer from time coverage between the oldest fetal (D136 post conception) and the adult (55 years old), many of the pathways dysregulated in the mouse organoids were also identified in the human organoid data.

Supplementation of DHA to Retinal Organoids

Based on comparative transcriptome analysis, we supplied some of the missing nutrients/growth factors to evaluate their impact on retinal, and specifically photoreceptor, development. Fatty acids (FAs), provided through the retinal pigment epithelium (RPE) to photoreceptors during development in vivo, constitute a major component of photoreceptor outer segments (Shindou et al., 2017). DHA reportedly activates the MAPK pathway and prevents photoreceptor apoptosis (German et al., 2006). We therefore included DHA in organoid cultures at D10 in an attempt to ameliorate its deficiency in vitro. BSA and linoleic acid (LA) served as a vehicle control and FA specificity control, respectively, because FAs can also be utilized by mitochondria and their metabolic intermediates can act as co-factors for gene regulation (Zhang et al., 2018). Immunostaining of rhodopsin and S-opsin revealed a gradually increased expression of both marker proteins in all three groups (BSA, LA, and DHA) from D18 to D32, indicating the progressive maturation of photoreceptors under all conditions (Figure 5A). No significant differences were observed among the three groups at D18 and D26. At D32, DHA-treated neural retina showed a moderately higher expression. Immunoblot analysis of rhodopsin confirmed the IHC analysis, and DHA-treated organoids had a 30% increase at D32 (Figure 5B). Interestingly, LA-treated organoids demonstrated a decrease in rhodopsin expression, suggesting that increased rhodopsin in DHA-treated organoids was not due to the addition of FAs. No significant differences were evident in cell types other than rods (Figure S7A). Gene profiles of BSA- and DHA-treated neural retina consistently revealed a moderate increase in rod-specific genes including Rho, Pde6a, and Gnat1 (Figure 5C). Cone-specific genes appear to be somewhat decreased in the DHA-treated group, but immunostaining of S-opsin was undetectable (Figure 5A). We then assessed the morphology of the photoreceptor cilia, which is an important feature in photoreceptor maturation. Immunostaining of ciliary axoneme markers ARL13B and RDS demonstrated improved ciliogenesis in DHA-treated photoreceptors (Figure 5D, upper). Transmission electron microscopy (TEM) analysis further revealed more developed inner segments, with mature mitochondria polarized to the apical side, in DHA-supplied organoid cultures (Figure 5D, lower). Although photoreceptor cilia showed comparable morphology in both groups, several cilia (7 out of 15 cilia captured in TEM) of >2 μm were evident in DHA-treated organoids compared with none in the control organoids (Figures 5D, lower; S7E). Notably, we frequently observed structures resembling nascent outer segment spirals in the inner segments of DHA-treated photoreceptors, and these structures were barely detected in the control organoids.

Figure 5

Photoreceptor Biogenesis in Fatty Acid-Treated Retinal Organoids

(A) RHO (green) and OPN1SW (red) showing rod and cone photoreceptors, respectively. Cell nuclei are visualized with DAPI (blue). Arrowheads indicate relevant immunostaining. D, differentiation day; LA, linoleic acid; DHA, docosahexaenoic acid.

(B) Rhodopsin expression in organoid cultures under different culture conditions revealed by immunoblotting. γ-Tubulin was used as housekeeping gene (upper). The histogram (lower) summarizes the data from three independent biological experiments (n = 3, at least three organoids were quantified in each experiment) and represented as mean ± SD. ∗p < 0.05.

(C) Gene expression heatmap of rod- and cone-specific genes for DHA-treated organoids at D32.

(D) Ciliogenesis of photoreceptors in organoid cultures shown by immunostaining of RHO (red), ciliary axoneme (ARL13B, magenta), and outer segment marker (RDS, green) (upper) and transmission electron microscopy (lower). Hollow, solid, and v-shaped arrowheads indicate relevant structures of photoreceptor cilium, mitochondria, and nascent outer segment spirals, respectively.

See also Figure S7.

Treatment of Retinal Organoids with FGF1

FGF1 is another growth factor that was missing in transcriptome data from in vitro retinal organoids (Figures 4H and 4I). In the modified protocol, we included FGF1 at D26, a stage corresponding to the upregulation of Fgf1 in vivo (Figure 4H). IHC analyses revealed significantly higher expression of S-opsin, which was polarized to the apical side of the neural retina in FGF1-treated organoids at D32 as compared with the control (Figure 6A). Rods as indicated by rhodopsin expression and other cell types were comparable between the two groups (Figure S7A). Transcriptome analysis of rod- and cone-specific genes confirmed the IHC results, with significant increase in several cone maturation genes including Opn1sw, Gngt2, Gnat2, Arr3, and Pde6h (Figure 6B). Except for a few genes, such as Rho, Gucy2f, and Pde6g, most rod-specific genes appear to be downregulated. As in case of DHA-treated organoids, longer photoreceptor cilia of >2 μm (five out of eight cilia captured in TEM) and nascent outer segment spirals were observed in FGF1-treated organoids (Figures 6C and S7E). Redundancy-removed GO analysis of the 5,375 DE genes identified between the two groups showed downregulation of cell cycle and cell proliferation genes and upregulation of genes related to vesicle transport, axon formation, and neuronal development in the FGF1-treated group (Figures S7B and S7C). Consistently, the presynaptic vesicle marker synaptophysin and the synaptic ribbon marker bassoon showed a more robust immunostaining in FGF1-treated organoids compared with the control groups (Figure 6D). Transcriptome profiling of organoid cultures supplemented with both DHA and FGF1 demonstrated no synergistic effect of the two factors, with rod-expressed genes showing similar expression as in FGF1-treated organoids and cone genes exhibiting similarity with DHA-treatment (Figure S7D).

Figure 6

Photoreceptor Biogenesis in FGF1-Treated Retinal Organoids

(A) Rod and cone photoreceptors revealed by RHO (green) and OPN1SW (red), respectively. Cell nuclei are visualized with DAPI (blue). Arrowheads indicate relevant immunostaining. D, differentiation day.

(B) Gene expression heatmap of rod- and cone-specific genes of D32 organoids under different treatment conditions.

(C) Immunostaining of rod photoreceptors (RHO, red), ciliary axoneme (ARL13B, magenta) and outer segment marker (RDS, green) (left) and transmission electron microscopy (right). Hollow, solid, and v-shaped arrowheads indicate relevant structures of photoreceptor cilium, mitochondria, and nascent outer segment spirals, respectively.

(D) Synapses of neural retina shown by presynaptic vesicles (synaptophysin, green) and ribbon synapses (bassoon, red). Cell nuclei are visualized with DAPI (blue).

See also Figure S7 and Table S7.

Discussion

Generation of neuronal tissues, such as the retina, involves differentiation of discrete cell types and their functional assemblies that are guided by complex interactions between the intrinsic genetic program and the microenvironment. Formation of the optic vesicle, and subsequently the optic cup, from neural ectoderm is accompanied by a defined set of eye-field TFs that stringently control patterns of gene expression to produce distinct retinal cell types. Cell-cell interactions and environmental factors provide further spatiotemporal precision to construct a functional retinal architecture capable of capturing and transmitting visual information. Neural retina with multiple cell types in organoid cultures has now permitted in vitro manipulation of specific physiological cues to better elucidate differentiation mechanisms and evaluate therapeutic paradigms. Here, we directly compare temporal transcriptome dynamics of developing mouse retina in vivo and in organoid cultures to better recognize commonalities as well as distinctions between the two systems. Our studies present detailed analyses of genes that specify different retinal cell types, regulatory networks that guide their expression, and cellular pathways that can be modulated further to create improved long-term modeling of retinal organoids in vitro.

Transcriptomes of developing retina in vivo and in vitro demonstrated a broad concordance in temporal gene profile, but delayed expression of cell-type-specific genes was observed in retinal organoids. Furthermore, organoid transcriptomes did not show the marked shift in gene expression pattern that normally occurs between P6 and P10 in vivo and is correlated to functional maturation of the retina, including photoreceptor outer segment and outer plexiform layer formation in vivo. Thus, cellular differentiation in retinal organoid cultures proceeds at a divergent pace despite the expression of early eye-field and morphogenic markers. Comparative and targeted analysis of gene regulatory factors in organoid transcriptome identified temporal dysregulation of multiple transcriptional repressors and reduced or aberrant expression of positive regulators of cell differentiation. We can specifically point to altered expression patterns of Zbtb, Zic, Plegl1, and Hes genes that control progenitor proliferation and timing of cell fate decisions in coordination with MAP kinase, Hedgehog, and other signaling pathways (Ma et al., 2007, Wall et al., 2009, Watabe et al., 2011). In organoid transcriptomes, we also identified lower expression of many TF genes involved in governing progenitor competence (e.g., Vsx2, Hes1, Hmga2, Onecut1, and Bclaf1), RGC (such as Isl1, Atf2, Pou4f1, and Pou4f3), or interneuron (e.g., Sox2, Nr4a2, Irx6, Sal11, Bclaf1, Nhlh1, and Zeb2) differentiation. We also noted reduced expression of genes that direct axon growth (such as Foxg1 and Irx4) and higher expression of those repressing axon growth (e.g., Nfil3). We must mention that transcriptome analysis alone does not establish cause versus effect relationship between gene expression and cell type specification. The lack of physiological cues and/or inappropriate organoid culture conditions can lead to deficiency of certain cells (e.g., specific lineage-restricted progenitors, RGCs), which in turn could be reflected in the transcriptome data. It is also conceivable that aberrant gene expression causes death of specific cell types.

Morphogenesis of distinct compartments, such as outer segment membrane discs and ribbon synapses, are essential for photoreceptor function. Lack of well-developed outer segment discs in organoid photoreceptors is likely because of the absence of RPE (Nasonkin et al., 2013). In addition to relevant molecular signals, RPE is also the source of nutrients during photoreceptor development. One such molecule is the polyunsaturated FA DHA, which is reported to maintain photoreceptor homeostasis (Bazan, 2007) and activates the MAPK pathway (German et al., 2006). MAPK signaling is implicated in maintaining long-term synaptic plasticity (Kelleher et al., 2004) and cellular responses to oxidative stress (Allebrandt et al., 2005). Many components of MAPK signaling, including FGF1, are upregulated from P6 onward in vivo but show reduced or no expression in retinal organoids. FGF signaling is associated with regulation of NRL and rod maintenance (Siffroi-Fernandez et al., 2008). Thus, inadequate MAPK-mediated stress response, probably due to an avascular system, can offer another plausible explanation for degeneration of mouse retinal organoids after short-term cultures.

Addition of DHA or FGF1 to organoid cultures had a significant impact on photoreceptor development. Both DHA- and FGF1-treated organoids demonstrated more mature photoreceptors with longer cilia (by TEM) and formation of nascent outer segment-like structures compared with controls. FGF1-treated organoids also showed robust immunostaining of synaptic markers. Interestingly, rod cells in DHA-treated organoids showed higher expression of rhodopsin. Cone photoreceptors, not prominent in mouse retinal organoids (Chen et al., 2016, Gonzalez-Cordero et al., 2013, Volkner et al., 2016), showed exuberant biogenesis and maintenance in FGF1-treated groups; however, the number of rods were reduced. Thus, DHA and FGF1 appear to have distinct effects on rod and cone differentiation, respectively. Transcriptome analysis of organoids treated with both DHA and FGF1 supported this hypothesis and indicated the involvement of independent mechanisms. Further investigations are warranted to delineate the role of DHA and FGF1 on photoreceptor development and survival.

Comparative and temporal transcriptomic studies, reported here, have permitted the development of methodologies for improved photoreceptor maturation and maintenance of S-opsin-expressing cones in retinal organoids. We note that DHA and FGF1 are not sufficient to produce functionally mature photoreceptors and that the organoids still do not survive long after D32 in culture. Manipulation of additional signaling pathways might therefore be required for organoid maturation, including the need for co-culture systems with RPE or use of small molecules to modulate specific pathways. Furthermore, the photoreceptors in mouse organoids are located at the inner side of neural retina, a configuration that may hinder sufficient delivery of nutrients and oxygen needed for their maturation. Application of efficient bioreactors and other bioengineering tools might help facilitate long-term culture and maturation of retinal organoids (DiStefano et al., 2018). Regardless, NGS-based comparative temporal dynamics offers a powerful tool for examining distinct approaches to produce functional retina in vitro. In that direction, our studies provide a valuable framework and transcriptomic resource for modeling development and disease in mouse retinal organoids. We also note that transcriptome analyses of human retinal organoids has provided an objective system to evaluate their differentiation status and shown accelerated rod maturation by using 9-cis retinal (Kaya et al., 2019).

Experimental Procedures

Animals, Tissue Collection, and RNA Isolation

C57Bl/6J mice without rd8 or rd1 mutation were used in the study, following the protocols approved by NEI Animal Care and Use Committee (ASP650) and adhered to the ARVO Statement for the Use of Animals. The retinas were collected from mice at the indicated stages (Figure 1A; Table S1).

RNA Sequencing

RNA-seq and data analysis were described previously (Kim et al., 2016). In brief, the in vivo dataset library was single-end sequenced to 76 bases and the retinal organoid dataset was paired-end sequenced to 125 bases. The bioinformatic pipeline (Figure S2A) used GRCm38.p4 and Ensembl v.84 annotation. Sequencing depth and alignment statistics are shown in Figure S2B and Table S1. For DE, we performed three separate analyses including intra-dataset ANOVA-like DE, paired-factor on five DTW-matched time points blocked for time, and D32 FGF1-treated versus untreated.

Human fetal retina (GEO: GSE104827), adult retina (GEO: GSE115828), and developing iPSC retinal organoids (GEO: GSE119320) fastq files were obtained from GEO (www.ncbi.nlm.nih.gov/geo), and the analysis pipeline (Figure S2A) was run using GRCh38.p12 with Ensembl v.94 for gene-level analysis.

Cluster Analysis and Functional Gene Analysis

Hierarchical clustering was performed on DE genes identified from the intra-dataset DE analysis and tree dendrogram was empirically cut to generate 36 clusters and grouped to yield 7 SCs. Functional gene analysis was performed using the gProfiler v.0.6.6 and enrichment analysis was performed using fgsea v.1.8.0 packages in R.

Gene Co-regulation Network Generation Using FGM Profiling

FGM profiles were generated from DTW-matched time points of the two datasets, separately, and pairwise Pearson's correlation (≥ abs(0.8)) of PCs was selected for downstream analysis and visualization.

Author Contributions

Overall Conceptualization, M.J.B., H.Y.C., and A.S.; Methodology and Investigation, H.Y.C., R.A.K., and K.N.; Resources, N.D.V. and T.L.; NGS Data Analysis, M.J.B., A.K.M., and V.C.; Data Curation, M.J.B.; Writing – Original Draft, M.J.B., H.Y.C., and A.S.; Writing – Review & Editing, all authors; Funding Acquisition, Supervision, and Project Administration, A.S.