Transcriptome-based molecular staging of human stem cell-derived retinal organoids uncovers accelerated photoreceptor differentiation by 9-cis retinal

Maintenance and differentiation of human pluripotent stem cells

CRX-GFP H9 is a subclone of the H9 human ESC line, carrying a green fluorescent protein (GFP) gene under the control of the cone‐rod homeobox (CRX; Gene ID: 1406, associated OMIM numbers: 120970, 613829, 268000) promoter as previously reported [26]. The human iPSC lines PEN8E and NEI377 were reprogrammed from skin biopsies using integration-free Sendai virus carrying the four Yamanaka factors, as described [33], and their genome integrity and pluripotency have been evaluated [33]. The ESP1 and ESP2 lines were reprogrammed by Oct4, Klf4, Sox2, c-Myc, and Lin-28 mRNAs [34] and the Sendai virus (Ctrl1 FiPS4F1, Spanish National Stem Cell Bank, Valencia, Spain), respectively. The H9, PEN8E, and NEI377 cell lines were maintained in Essential 8 medium (E8; ThermoFisher Scientific, Waltham, MA) under hypoxia (5% O₂), and the ESP1 and ESP2 lines in mTeSR1^® (Stem Cell Technologies, Vancouver, Canada) under normoxia (20% O₂). All lines were sustained on BD Matrigel^TM human embryonic stem cell-qualified Basement Membrane Matrix (Corning, New York, NY)-coated plates. The PSCs were passaged every 3–4 days at 70–80% confluency using the EDTA-based protocol [35].

To start differentiation, all lines were detached and dissociated into small clumps with the EDTA dissociation protocol [35], and cultured in polyHEMA-coated or Ultra Low Attachment Culture Dishes (Corning) in E8 with 10 μM Y-27632 (Tocris) to form embryoid bodies (EBs). Neural-induction medium (NIM) consisting of Dulbecco's Modified Eagle's medium (DMEM/Ham's F-12 Nutrient Mix (F-12) (1:1), 1% N2 supplement (ThermoFisher Scientific), 1X Minimum essential medium (MEM) non-essential amino acids (NEAA; MilliporeSigma), and 2 µg/ml heparin (MilliporeSigma) was added at differentiation day (D)1 and D2 to reach a final ratio of 1:3 and 1:1 NIM:E8, respectively. Starting at D3, EBs were cultured in 100% NIM for an additional 4 days. EBs were collected at D7 and plated onto BD Matrigel^TM Growth Factor Reduced (GFR) Basement Membrane Matrix (Corning)-coated plates. On D16, the medium was changed to photoreceptor induction medium (PIM) consisting of DMEM/F-12 (3:1) with 2% B27 without vitamin A, 1% antibiotic-antimycotic solution (ThermoFisher Scientific), and 1X NEAA. Medium was changed daily.

Upon the appearance of optic vesicles (OVs; typically, D21–28) with neuroepithelium morphology, the regions were excised with tungsten needles under magnification and transferred to suspension culture in polyHEMA-coated or Ultra Low Attachment 10 cm² Culture Dishes (Corning) in PIM. After D42, PIM was supplemented with 10% fetal bovine serum (FBS; ThermoFisher Scientific), 100 mM taurine (MilliporeSigma), and 2 mM GlutaMAX (ThermoFisher Scientific). Starting at D63, the PIM was supplemented with 1 µM retinoid (MilliporeSigma) three times a week and switched to 0.5 µM at D92. The 1% N2 supplement was included starting at D92. The medium for the H9 and PEN8E-derived organoids was supplemented with 20 ng/ml insulin-like growth factor 1 (IGF1; ThermoFisher Scientific) and 55 nM beta-mercaptoethanol (2-ME; ThermoFisher Scientific) from dissection until the end of differentiation. 9-cis retinal (9CRAL; MilliporeSigma) was added at the described concentrations and corresponding time periods during the media change. Organoids for the data set PEN8E_2 were cultured with 10 ng/ml IGF1 from D42 to D200, and ATRA (MilliporeSigma) was supplemented to the PIM at 1 µM 5 times per week and at 0.5 µM twice a week for D63–92 and D92–180, respectively. The ESP-derived organoid cultures were maintained in the same manner as the PEN8E_2 organoids, except no IGF1 was used. For all experiments directly comparing 9CRAL and ATRA, media changes were performed in the dark under dim red light to reduce isomerization of the retinoids. Retinal organoids were collected in the dark and processed in the light as described for all other experiments.

Immunofluorescence

hPSC-derived retinal organoids were collected throughout the differentiation process, fixed, and immunostained as previously described [36]. Briefly, the organoids were fixed in 2% or 4% paraformaldehyde in 1X PBS (D-PBS, 1X, 138 mM NaCl, 2.6 mM KCl, 8 mM Na₂HPO₄-7H₂O, 1.5 mM KH₂PO₄, pH 7.4) for 1 h at room temperature. Fixation was followed by cryoprotection in a sucrose gradient, 10–30%. Retinal tissues were then embedded in Shandon M1 embedding matrix and frozen on dry ice. Ten-micron sections were obtained using a Leica cryostat (Leica Biosystems, Buffalo Grove, IL) at −14 °C. The sections were placed on Superfrost Plus slides (Fisher Scientific, Hampton, NH) and stored at −20 °C until use. Specific antibodies and concentrations are summarized in Appendix 1.

Immunoblots

The protocol for immunoblot analysis was as previously described [37]. Individual organoids were lysed in 50 μl of 1% Triton X-100 (MilliporeSigma) in 1X PBS (ThermoFisher Scientific) supplemented with 1X protease inhibitors (Roche, Basel, Switzerland). Gentle trituration was then used to dissociate the organoid until a homogenous solution was obtained. The samples were then incubated on a nutator at 4 °C for 1 h and centrifuged at 4 °C at 4,000 ×g for 5 min. The protein concentration was measured via the Pierce bicinchoninic acid (BCA) protein assay (ThermoFisher). Approximately 15 μg supernatant protein was diluted 4:1 in reducing 4X Laemmeli buffer for 1 h at room temperature and separated at 100 V for 1 h on 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) gels, which were then transferred to polyvinylidene fluoride (PVDF) membranes at 100 V for 1 h. After blocking in 5% milk in 1X Tris-Buffered Saline, 0.1% Tween® 20 Detergent (TBST; Santa Cruz Biotech, Dallas, TX) for 1 h at room temperature, the blots were incubated in an antibody cocktail (1:1,000 1D4, 3A6, and 4D2, a generous gift from Dr. Robert Molday, Vancouver, Canada) overnight in 1% milk in 1X TBST. The following morning, the membranes were washed in 1X TBST and shaken at room temperature four times for 10 min each. The membranes were then incubated in appropriate secondary species immunoglobulin G (IgG) conjugated to horseradish peroxidase (1:10,000) in 1X TBST shaking at room temperature for 1.5 h. The membranes were then washed again four times in 1X TBST for 10 min each. Before imaging, the membranes were exposed to SuperSignal® West Pico enhanced chemiluminescence (ECL) solution (ThermoFisher) for 5 min, and chemiluminescence was captured using a Bio-Rad ChemiDoc™ touch (Bio-Rad, Hercules, CA). Membranes were next dried overnight at room temperature to remove antibody binding to the protein. Dry membranes were reactivated in 100% methanol and incubated with gamma-tubulin (1:1,000, ab11317, Abcam, Cambridge, United Kingdom) overnight as described above. Images were analyzed in Image Lab (Bio-Rad) and exported to Adobe Photoshop (San Jose, CA) for figure generation.

Transmission electron microscopy

Retinal organoids were processed for transmission electron microscopy (TEM) analysis as previously described [33]. Briefly, the organoids were initially fixed in 4% formaldehyde and 2% glutaraldehyde in 0.1 M cacodylate buffer (pH 7.4; Tousimis, Rockville, MD) for 2 h, washed in cacodylate buffer three times, and then fixed in osmium tetroxide (1% v/v in 0.1 M cacodylate buffer; Electron Microscope Science, Hatfield, PA) for 1 h in a room temperature. The organoids were washed in the same buffer three times, followed by acetate buffer (0.1 M, pH 4.2), and en-bloc staining in uranyl acetate (0.5% w/v; Electron Microscope Science) in acetate buffer for 1 h. The samples were dehydrated in ethanol solution (35%, 50%, 75%, 95%, and 100%), followed by propylene oxide. Subsequently, the samples were infiltrated in a mixture of propylene oxide and epoxy resin (1:1) overnight, embedded in a pure epoxy resin in a flat mold, and cured in a 55 °C oven for 48 h. Thin sections (70 to 80 nm) were made with an ultramicrotome (UC 7) and diamond knife (Diatome, Hatfield, PA), attached on 200-mesh copper grid, and counter-stained in aqueous solution of uranyl acetate (0.5% w/v) followed by lead citrate solutions. The thin sections were stabilized with carbon evaporation under a vacuum evaporator before the EM examination. The digital images were taken in the electron microscope (H7650)-equipped digital camera (AMT, Rosemont, IL).

RNA-seq alignment and quantitation

High-quality total RNA (100 ng, RIN >7) from at least two independent replicates for hPSCs at various differentiation time points were subjected to mRNA directional library construction as described previously [26,38]. Paired-end sequencing was performed to a length of 125 bases using HiSeq2500 (Illumina, San Diego, CA). Genome reference sequence GRCm38.p7 and Ensembl v82 annotation was used for alignment and quantitation. Quality control, sequence alignment, transcript, and gene-level quantification of primary RNA sequence (RNA-seq) data were accomplished using an established bioinformatics pipeline [39] Appendix 2. The transcriptome of adult human retina samples was obtained from NEI Commons (Brooks, MJ and Swaroop, A).

Differential expression and clustering analysis

Data generated for the five different cell lines, transcriptome comparison and 9CRAL/ATRA analysis were analyzed separately. For each data set, genes were kept for further analysis only if there were ≥5 count per million (CPM) in all replicates of at least one group of the data sets. The data were subjected to trimmed mean of log expression ratios values (TMM) normalization with edgeR v3.20.9 [40,41], and principle component analysis (PCA) and Pearson correlation were performed with normalized CPM (log₂) values. Differential expression analysis was performed using limma v3.34.9 [42], and genes that had ≥2-fold change at least between two time points in either transcriptome of comparison and a false discovery rate (FDR) ≤0.01 were considered statistically significantly differentially expressed.

Gene expression clustering was performed on gene-wise Z-scores using AP (apcluster v1.4.7 package) [43,44] with the k-means option, where k=12 by choosing the corSimMat function as the similarity measure. Determination of the cluster number k was accomplished by empirical observation of the heatmap produced from the PCA projections of the genes-wise z-score.

Gene Ontology analysis

GO enrichment analysis was performed using clusterProfiler v3.6.0 [45]. To reduce redundancy associated with GO analysis, we performed Wang semantic similarity comparisons of the enriched GO terms using GOSemSim v2.4.1 [46] and semantic similarities clustered with AP. The term with the highest level in each cluster was taken as the enriched GO term. Heatmaps of GO term gene lists were generated using the log2 CPM values. The GO enrichment ratio was calculated with reference to the number of pathway genes found overlapping with the analysis.

Time-course stage comparison of retinal organoids with human fetal retina

Open-ended dynamic time warping (OE-DTW) analysis was used to compare the maturation state of the different retina organoid time-course stage expression data to the human fetal retina (GSE104827) and the human adult retina, as previously described [47]. Retina organoid time-course stages were defined as Stage 1: 25–36 days, Stage 2: 50–75 days, Stage 3: 80–90 days, Stage 4: 105–125 days, Stage 5: 145–172 days, and Stage 6: 186–205. Mean gene expression values at each stage were used for OE-DTW analysis. The gene set used for the OE-DTW analysis consisted of 118 well-known retina-centric transcription factors and cell-type markers (see Appendix 3 and Appendix 4). Co-inertia analysis was performed using MADE4 v1.52.0 package in R.

Comparable morphology between hESC- and hiPSC-derived retinal organoids

Two human pluripotent stem cell lines, H9 (hESC) and PEN8E (hiPSC), were differentiated into retinal organoids using a widely used robust protocol [29], with minor modifications based on our mouse organoid differentiation method that included 9CRAL [36]. Both lines could generate phase-bright neural retina with cone and rod photoreceptors of comparable morphology (Figure 1A–C). S opsin (OPN1SW) was polarized to the apical side of the organoids (outer surface of organoid, exposed to media) as early as D90, with a further increase in expression from D120 to D200. L/M cone photoreceptors (labeled with L/M opsin, OPN1L/MW) were barely evident at D120 but increased at D200 (Figure 1C, left). Rhodopsin (RHO) immunostaining was observed at the apical side as early as D120 in the H9- and PEN8E-derived organoids, and elongated as the photoreceptors matured (Figure 1C, right). Retinal ganglion cells (RGCs) were evident at early stages, but could not be maintained through the end of differentiation, while all other neural retina cell types showed similar morphology and development in the H9- and PEN8E-organoids (Appendix 5).

Differentially expressed gene clusters during organoid development

We then performed RNA-seq analysis of the developing retinal organoids to decipher the major differentiation stages in vitro. The H9- and PEN8E-derived retinal organoids showed a substantial overlap of a total of 3,851 differentially expressed (DE) genes that exhibited statistically significant changes in expression during differentiation (Figure 2A). The 1,041 DE genes unique to H9 and the 1,749 genes unique to PEN8E showed the same general trend in both data sets (Appendix 6). We noted that many of these genes were statistically significantly differentially expressed in both data sets under a less stringent cutoff of 1.5-fold and 5% FDR (data not shown). AP clustering analysis of the common DE genes between the two datasets yielded eight clusters (Cs) with either monotonically increasing or decreasing expression (Figure 2B). GO Biologic Process enrichment analysis uncovered key biologic pathways and genes that could be associated with specific stages of retinal development, and identified statistically significantly enriched terms for all clusters except C1 and C7 (Figure 2C). The C2 through C4 genes displayed a gradual decrease throughout differentiation, and were associated with negative regulation of neuronal differentiation, mitotic cell division, and neurodevelopmental processes and signaling (Figure 2D, Appendix 7). C4 also included axon guidance genes, which likely reflected the loss of ganglion cells in the organoid cultures. In contrast, the C5, C6, and C8 genes exhibited progressively increasing expression during organoid development, and included genes associated with neuronal or retinal differentiation and functions, including visual perception, synaptogenesis, and phototransduction.

Transcriptomes of organoids from different iPSCs and protocols

To evaluate the impact of iPSC lines or modifications in the differentiation protocol, we compared the expression of the 3,851 DE genes with organoids produced from different hiPSC lines (ESP1, ESP2, and NEI377) and protocols (PEN8E_2). PCA revealed a similar temporal progression in all organoid transcriptomes with increments on the differentiation days, as evident from PC1 (Appendix 4). As the harvesting time points varied among individuals and laboratories, we organized the organoid samples into six groups, for ease of comparison, based on the collection day (chronological age) and state of differentiation (Figure 3A) and performed PCA (Figure 3B) using a selection of established, retina-centric genes (Appendix 3). Although all organoids clustered together at an early stage of differentiation (group 1), the H9- and PEN8E-derived organoids (which included 9CRAL supplementation) showed accelerated differentiation from group 3 and onward, compared to the ESP, NEI377, and PEN8E_2 organoids (which used ATRA as described by [29]). In concordance, the ESP, NEI377, and PEN8E_2 organoids displayed lower expression of genes associated with phototransduction and outer segments compared to the H9- and PEN8E-derived organoids (Appendix 4). The NEI377 and PEN8E_2 organoids were generated by an identical protocol, and they displayed a similar developmental pattern, except group 6 of PEN8E_2, which is likely due to iPSC line differences.

Molecular staging of retinal organoids based on developing human retina in vivo

To determine the maturity of organoid cultures using objective parameters, we performed OE-DTW analysis to align the time-series expression data between the organoids and the human fetal and adult retina gene profiles (Figure 3C) using a set of retina-centric genes (Appendix 3). Across all stem cell lines, group 1 organoids (D25–37) matched with the earliest human fetal time points (D52–67). Subsequently, the distinct PSC organoids showed differences in their highest correlation to corresponding human retinal stages. While organoids derived from ESP1, ESP2, NEI377, and PEN8E_2 lines revealed higher correlation with fetal retina samples; groups 4–6 H9 and PEN8E organoids were more concordant with the late fetal and adult retina (Figure 3C). For example, the group 4 H9 and PEN8E-derived organoids (D145–172) were highly correlated with adult samples (17, 24, and 28 years of age), whereas the ESP1, ESP2, NEI377, and PEN8E_2 organoids corresponded with late fetal samples (D105–136) at a similar stage (Figure 3C).

To confirm the result of the DTW analysis, we compared the expression of mature retinal genes associated with phototransduction and synaptic function between organoids and in vivo retinas. The late stage (group 5 and 6) H9 and PEN8E organoids showed similar expression of phototransduction genes as that of the adult retina (Figure 3D). However, most genes from “Transmission across chemical synapses” pathways were downregulated in the organoid cultures, suggesting the retinal organoids were not yet fully mature (Appendix 4).

Accelerated rod photoreceptor differentiation with 9CRAL

Much like its physiologically relevant isomer 11-cis retinal, 9CRAL can also bind opsin to form functioning rhodopsin, and retinoic acid (RA) can be generated as an oxidation product. Therefore, we hypothesized that replacing RA with 9CRAL in the modified protocol would result in improved rod differentiation. To evaluate this hypothesis, we performed a direct comparison of the D90 and D120 H9 organoids cultured with ATRA or 9CRAL from D64 onward. We then applied coinertia analysis to the transcriptomes of the ATRA- and 9CRAL-supplemented organoids and the developing human fetal retina (D53 to D136). The D90 and D120 of the 9CRAL organoid transcriptomes corresponded with the D80–94 and D125–136 human fetal retina gene profiles, respectively. However, the D90 and D120 ATRA transcriptomes matched with earlier fetal retina time points (D57–67 and D105–115, respectively; Figure 4A). Thus, treatment with 9CRAL expedited retinal development in organoid cultures and showed closer temporal transcriptome dynamics to the human fetal retina. The DE genes between D90 and D120 in the ATRA (1284 genes) and 9CRAL (852 genes) data presented an overlap of 594 genes (Figure 4B), which could, in turn, be grouped into five clusters with hierarchical clustering analysis (Figure 4C). Redundancy-reduced GO analysis of the common DE genes revealed a monotonic decrease in the C1 and C2 genes involved in the cell cycle, suggesting neural progenitor cells exit the cell cycle and become mature in retinal organoids (Figure 4D). The C3 through C5 genes were involved in energy metabolism and visual perception, which further characterized their progressive increase from D90 to D120. Although the ATRA and 9CRAL organoids showed comparable developmental patterns, changes in the DE genes during development varied between the two groups (Figure 4C). Day-matched comparison of DE genes revealed lower expression of genes involved in regulation of the RA signaling pathway, RA metabolic processing, and energy metabolism, whereas higher expression of visual perception and function genes was evident in the D90 and D120 9CRAL-treated organoids (Figure 4E). The heatmap of early eye field transcription factors and photoreceptor genes consistently showed expedited photoreceptor differentiation in the 9CRAL organoids, as revealed by the higher expression of these genes, compared to the ATRA organoids (Figure 4F).

To validate the findings that 9CRAL accelerated photoreceptor differentiation, we performed immunohistochemistry on the D90 and D120 organoids. The D90 9CRAL organoids showed lower expression of retinal progenitor cell marker VSX2 (also called CHX10) and higher expression of a pan-photoreceptor marker (recoverin, RCVRN) compared to the ATRA organoids (Figure 5A, upper panel). At D120, the 9CRAL-treated organoids showed distinctively more RHO+ cells, which were barely detected in the ATRA organoids (Figure 5A, bottom panel), consistent with the higher rhodopsin expression in 9CRAL (14.0±9.00 CPM) compared to the ATRA (3.8±1.7 CPM) RNA-seq data. No marked difference was observed in the expression or morphology of the L/M-opsin and S-opsin+ cells (see Figure 5A, bottom panel), as well as other neural retina cell types (Appendix 8). In concordance, immunoblot analysis demonstrated increased rhodopsin expression in the 9CRAL organoids compared to the ATRA organoids at D120 (Figure 5B). Ultrastructural analysis of D130 retinal organoids using transmission electron microscopy (TEM) showed an advanced stage of maturation of the 9CRAL photoreceptors compared to the ATRA photoreceptors (Figure 5C). The ellipsoid (apical) side of the photoreceptor inner segments in the 9CRAL organoids contained a higher number of mitochondria with typical morphology, which could explain the divergent expression of genes involved in energy metabolism between the 9CRAL and ATRA organoids. No differences were evident in the morphology of the photoreceptor connecting cilia between the two groups. Taken together, these data demonstrate that 9CRAL supplementation expedited differentiation and maturation of rod photoreceptors in retinal organoids.

Summary statement

Three-dimensional organoids derived from human pluripotent stem cells have been extensively applied for investigating organogenesis, modeling diseases, and developing therapies. However, substantial variations within organoids pose challenges for comparison among different cultures and studies. We generated transcriptomes of multiple distinct retinal organoids and compared them to human fetal and adult retina gene profiles for molecular staging of differentiation state of the cultures. The analysis revealed the advantage of using 9-cis retinal, instead of the widely used all-trans retinoic acid, in facilitating rod photoreceptor differentiation. Thus, a transcriptome-based comparison can provide an objective method for uncovering the maturity of organoid cultures across different lines and in various study platforms.