Introduction

Vertebrate photoreceptors develop from optic vesicle retinal progenitor cells (RPCs) through progressive RPC lineage restriction, fate determination, and post-mitotic maturation (1, 2). While several transcription factors that govern these events have been identified, important aspects of photoreceptor fate determination and maturation remain unclear. For example, it is unclear if fate is determined in RPCs, where OTX2 and ONECUT1 are thought to control the post-mitotic expression of long- or medium-wavelength (L/M) cone determinant TRβ2 and rod determinant NRL (3), or is determined in post-mitotic photoreceptor precursors with concurrent TRβ2 and NRL expression (4). Following fate commitment, post-mitotic developmental stages have been defined based on morphologic features and phototransduction-related gene or protein expression (5, 6), but it is unclear if progression through such stages is subdivided into distinct cell states governed by unique transcription factor combinations or represents a developmental continuum.

An improved understanding of photoreceptor development may provide insight into the pathogenesis of retinal dystrophies, retinal degenerations, and the retinal cone precursor cancer, retinoblastoma (7, 8). In the latter case, L/M cone precursors lacking the retinoblastoma protein (pRB) were shown to proliferate in a manner dependent on the L/M-cone lineage factors RXRγ and TRβ2 and the intrinsically highly expressed MDM2 and MYCN oncoproteins, likely representing the first step of retinoblastoma tumorigenesis (8–11). Similarly, retinoblastoma cell proliferation depends on RXRγ, TRβ2, MDM2, and MYCN (9), implying that intrinsic L/M-cone factors contribute to the oncogenic state. However, retinoblastoma cells also express rod lineage factor NRL RNAs, which – along with other evidence – suggested a heretofore unexplained connection between rod gene expression and retinoblastoma development (12, 13). Improved discrimination of early photoreceptor states is needed to determine if co-expression of rod- and cone-related genes is adopted during tumorigenesis or reflects the co-expression of such genes in the retinoblastoma cell of origin.

The cone precursors’ propensity to form retinoblastoma is a human-specific feature whose study requires analysis of developing human retina (8). Single-cell RNA-sequencing (scRNA-seq) is well suited to such analyses as it enables discrimination of cell type-specific fate determination and maturation features. scRNA-seq studies employing 3’ end-counting have defined age-related post-mitotic transition populations, fate-determining features of post-mitotic photoreceptor precursors (14–16), and gene expression changes associated with the cone fate decision and early development (14–19). However, 3’ end-counting cannot be used to interrogate transcript isoforms, and the relatively low number of genes detected per cell limits the ability to distinguish closely related states in individual cells.

In this study, we sought to further define the transcriptomic underpinnings of human photoreceptor development and their relationship to retinoblastoma development. We generated deep, full-length single-cell transcriptomes of human retinal progenitor cells (RPCs) and developing photoreceptors from fetal week 13-19 retinae, with enrichment of rare cone precursor populations, and applied long-read cDNA sequencing, RNA velocity, pseudotemporal trajectory reconstruction, and single-cell regulatory network inference and clustering (SCENIC) to interrogate individual cell states. These analyses discriminated previously unrecognized photoreceptor developmental states, identified photoreceptor precursor states with cone and rod-related RNA co-expression, uncovered cell type-specific expression of RNA isoforms of photoreceptor fate-determining genes, elucidated post-mitotic photoreceptor developmental trajectories, and defined retinoblastoma cell-of-origin features that contribute to retinoblastoma genesis.

Results

Regulon-defined RPC and photoreceptor precursor states

To interrogate transcriptomic changes during human photoreceptor development, dissociated RPCs and photoreceptor precursors were FACS-enriched from 18 fetal week 13–19 retinae and isolated using microfluidics or direct sorting into microliter droplets, followed by full-length cDNA synthesis, paired-end sequencing, and alignment to Ensembl transcript isoforms (Figure 1A).

Figure 1:

The FACS enrichment was based on a prior cone isolation method (10) but with wider gating to include rods and RPCs. Exclusion of all cells with <100,000 reads and 18 cells expressing markers of retinal ganglion, amacrine, and/or horizontal cells (POU4F1, POU4F2, POU4F3, TFAP2A, TFAP2B, ISL1) and lacking photoreceptor lineage marker OTX2 yielded 794 single cells with averages of 3,750,417 uniquely aligned reads, 8,278 genes detected, and 20,343 Ensembl transcripts inferred (Figure S1A-C). Sequencing batches were normalized and transcriptomes clustered and visualized in uniform manifold approximation and projection (UMAP) plots that integrated cells across different retinae, isolation methods, and sequencing runs (Figure 1B, S1D-F).

Low resolution Louvain clustering generated six clusters that segregated into mostly distinct UMAP domains (Figure 1B). One cluster was comprised of RPCs and Mϋller glia (MG), with specific expression of LHX2, VSX2, SOX2, and SLC1A3, while five clusters were comprised of cells with photoreceptor features, with wide expression of OTX2 and CRX and cluster-specific rod- and cone gene expression (Figure 1C, S2). In UMAP space, the cluster designated immature photoreceptor precursors (iPRPs) intermixed with the RPC/MG population, extended towards early LM cone precursors, and was predominantly comprised of cells expressing the L/M cone determinant THRB (Figures 1B,C).

Two clusters expressing the rod determinant NR2E3 were designated early maturing rod (ER) and late maturing rod (LR) (Figure 1B), based on the latter’s increased expression of rod phototransduction genes GNAT1, CNGB1, PDE6G, GNGT1 and RHO (Figure 1D, S2C).

Differential expression analysis revealed other genes upregulated in late rods (GNB1, SAMD7, NT5E) (20, 21) as well as the downregulated CRABP2, DCT, and FABP7 (Figure S3A, Table S1). Genes upregulated in the LR cluster were enriched for photoreceptor and light sensing gene ontologies (Figure S3B) including spectrin binding, likely relating to proteins that control photoreceptor polarity and synapse formation (22, 23), and purine biosynthesis and ribonucleotide metabolic processes, potentially related to the developing photoreceptors’ high NAD+ requirements (24).

Cones segregated into distinct S- and L/M-cone clusters, with differential expression of cone subtype markers (OPN1SW, THRB), previously identified S-cone enriched genes (CCDC136, UPB1) (25–27), novel S cone enriched genes (MEGF10, NRXN3, ACKR3), and the L/M cone transcription factor ISL2 (15), among others (Figure 1B,C, S3C, Table S2). Gene ontology analysis did not reveal relevant terms enriched in S cones, whereas L/M cones were enriched for protein translation related ontologies (Figure S3D) due to increased expression of ribosomal protein genes RPL23A, RPLP0, RPS19, RPS27, RPS27A, RPS29, RPS3A (Table S2).

In UMAP space, the THRB+ L/M-cone cluster segregated into a large 363-cell proximal population and a 5-cell distal population inferred to represent early maturing and late maturing stages, respectively, based on the latter’s increased expression of cone phototransduction genes OPN1LW (encoding L-opsin), PDE6H, and GUCA1C (Figure 1D, S3E,F), analogous to RHO, PDE6G, and GNGT1 upregulation in late rods. Compared to the early L/M population, late L/M cones had upregulation of three M-opsin genes, TTR, encoding the retinol-binding protein transthyretin, PCP4, a small protein that binds calmodulin previously noted in foveal cones (28), and MYL4, a myosin light chain gene upregulated in retinal organoid L/M cones (27), among others (Figure S3E,F, Table S3). The low proportions of OPN1MW+ and OPN1LW+ late maturing L/M cones is consistent with a prior analysis of similar-age retinae and with the further upregulation of these proteins in later maturation (15).

To further interrogate cell identities, we used SCENIC to identify cluster-specific transcription factor regulons, which represent the overall expression of single transcription factors and their likely coregulated target genes (29). The highest specificity regulons defined major cell populations including RPC/MG-specific E2F2, E2F3, VSX2, and PAX6; pan-photoreceptor NEUROD1, OTX2, and CRX; rod-specific NRL; and L/M cone-specific THRB and ISL2 (Figure 2, Table S4). SCENIC also distinguished ER and LR states via the latter’s increased NRL, CRX, ATF4, and LHX3 and decreased HMX1 and RAX activities (p <0.0005 for each, Dunn test). RAX activity also decreased in the 5-cell late maturing L/M cone group (Figure 2B), supporting its distinct transcriptomic identity and suggesting a similar mode of late cone and late rod maturation. Additionally, iPRPs expressed their most specific regulons, LHX9 and OLIG2, at levels similar to RPC/MGs along with photoreceptor-related regulons, consistent with the transitional nature of this population (Figure 2A, Table S4). However, SCENIC did not identify S cone-specific regulons, in keeping with the notion that S cones represent a default photoreceptor state induced by pan-photoreceptor factors such as OTX2, CRX, and NEUROD1 in the absence of NRL and THRB (30). Thus, deep, full-length scRNA-seq enabled identification of regulons underlying RPC and developing photoreceptor states at the single cell level.

Figure 2:

Differential expression of NRL and THRB isoforms in rod and cone precursors

Although cone and rod precursors segregated into distinct clusters, mRNAs encoding rod-determining factor NRL, L/M cone-determining factor TRβ, and cone marker RXRγ were co-expressed in both rod and cone precursor populations, with mean NRL expression only 4.3-fold higher in the ER vs LM cluster, mean THRB expression 5.0-fold higher in LM vs LR, and mean RXRG expression 4.3-fold higher in LM vs ER (Figure S2C,D and below). Cone NRL expression was unexpected given NRL’s role in rod fate determination (27) and rod-specific NRL regulon activity (Figure 2A,D). Similarly, rod RXRG and THRB expression were unexpected given their roles in cone gene expression and fate determination (31, 32) and L/M cone-specific THRB regulon activity (Figure 2A,D). Accordingly, we used full-length scRNA-seq data to determine if cone and rod precursors differentially express NRL, THRB, and RXRG transcript isoforms.

For NRL, the most highly assigned transcript isoform (ENST00000397002) encoded the canonical full-length NRL (FL-NRL) (RefSeq NP_001341697.1), while the second-most assigned isoform (ENST00000560550) encoded a previously uncharacterized transcript that used an alternative “P2” promoter and first exon, here termed exon 1T (Figure 3A-C). The novel transcript was predicted to encode an N-terminally truncated NRL (Tr-NRL) retaining the leucine zipper DNA binding domain but lacking the minimal transactivation domain (33) (Figure 3C).

Figure 3:

While both isoforms were inferred to be more highly expressed in the ER vs LM cluster, the average FL-NRL (ENST00000397002): Tr-NRL (ENST00000560550) ratio was 2.8:1 in early rods and 0.72:1 in L/M cones (Figure 3B). Consistent with the assigned isoform ratios, mean Tr-NRL exon 1T coverage was higher in S and L/M cones while FL-NRL exon 1 coverage was higher in rods (Figure 3C, red vs. black arrowheads). Comparing the reads mapped to each first exon relative to total reads further confirmed that the Tr-NRL exon 1T predominated in individual cones whereas the FL-NRL exon 1 predominated in rods (Figure 3D). The cone cells’ higher proportional expression of Tr-NRL first exon sequences was validated by RNA fluorescence in situ hybridization (FISH) of FW16 fetal retina (Figure 3E,F).

While the Tr-NRL isoform was not to our knowledge previously described, another NRL isoform that initiated within exon 2 and lacked the NRL transactivation domain due to alternative exon 2 splicing, termed DD10 (Figure 3C), was previously identified in adult retina (34). Concordantly, we detected reads spanning the unique DD10 splice junction, yet at lower levels than the unique Tr-NRL junction (5,942 versus 57,048). As both Tr-NRL and DD10 lack the NRL transactivation domain and as DD10 interferes with FL-NRL transactivation (35), we examined if Tr-NRL similarly opposes FL-NRL transcriptional activity. In luciferase assays, Tr-NRL suppressed FL-NRL activation of a PDE6B promoter (Figure 3G), demonstrating that the cone-enriched Tr-NRL antagonized FL-NRL-mediated transactivation.

To further define NRL isoforms, we performed NRL 5’ rapid amplification of cDNA ends (RACE) on the already generated single cell cDNA libraries from 23 early rod (ER) and 21 early L/M cone (LM) cells, followed by nanopore long-read sequencing of the pooled RACE products. The long-read sequencing revealed isoforms consistent with FL-NRL, Tr-NRL, DD10, and several other NRL isoforms with alternative transcription initiation and alternative splicing within exon 2 as well as within the Tr-NRL exon 1T (Figure 3H, Figure S4). In keeping with the computationally assigned isoforms and differential exon usage, ER libraries had a higher proportion of FL-NRL exon 1 and exon 2 reads and LM libraries had a higher proportion of Tr-NRL exon 1T reads (Figure 3H, Figure S4). Moreover, alternative splicing within NRL exon 2 was more prevalent in LM libraries, affecting 6,788 of 29,177 (23 %) of exon 2 reads compared to 5,378 of 85,860 (6 %) in ER cells (p<0.0001; Chi-square with Yates correction), resulting in rarer full-length (exon 1-2-3) transcripts than inferred from short-read sequencing. Thus, long-read sequencing revealed cell type-biased NRL isoform expression with a preponderance of FL-NRL in early rods and disrupted FL-NRL and Tr-NRL isoforms in LM cones.

Despite our detection of L/M cone RNAs encoding Tr-NRL and FL-NRL, cone expression of NRL protein has not been reported. To assess endogenous Tr-NRL expression, we performed immunoblot analysis of CHLA-VC-RB-31 retinoblastoma cells (36), which express FL-NRL- and Tr-NRL transcripts in a cone-like 0.57:1 ratio, with an antibody that recognizes both FL- and Tr-NRL proteins (Figure S5A,B). As in other retinoblastoma cells (13), FL-NRL increased in response to retinoic acid and proteasome inhibition, whereas Tr-NRL was not detected (Figure 5C). Similarly, Tr-NRL was not detected in EGFP+, RXRγ+ cones following lentiviral transduction of an explanted fetal retina with an EGFP-P2A-Tr-NRL cassette (Figure S5D-G). These findings suggest that Tr-NRL is poorly translated or unstable in most cone precursors.

For THRB, the most highly assigned transcript isoforms included those encoding the L/M cone-specific TRβ2 (ENST00000280696) and the more widely expressed TRβ1 (ENST00000396671 and others) (Figure 4A-C, S6). While both TRβ1 and TRβ2 promote L/M cone fate determination (37), isoforms encoding TRβ2 predominated in L/M cones while those encoding TRβ1 predominated in early rods (Figure 4B). Moreover, late rods preferentially expressed THRB exons 1-6 and the THRB 3’ untranslated region (UTR), implying that RNAs encoding full-length TRβ proteins were rare (Figure 4B,C). Notably, a higher percentage of reads extended from the exon 4 and exon 6 splice donor sequences into the subsequent introns in LR versus LM cells (Figure 4D,E) (p ≤ 0.001 for both, two-tailed Chi square test), suggesting that premature transcription termination (PTT) in introns 4 and 6 preferentially limits TRβ expression in the LR population. The inferred PTT events are consistent with structures of the assigned Ensembl isoforms (Figure 4B, S6).

Figure 4:

Figure 5:

To further evaluate PTT events, we used THRB 3’ RACE and long-read sequencing on single cell cDNA libraries from 23 L/M cone cells and five LR cells selected for high THRB expression. RACE reactions were performed separately with primers complementary to the TRβ1-specific exon 4 and to the TRβ2-specific first exon (Figure 4F). Sequencing of LM and LR RACE products corroborated pronounced PTT in introns 4 and 6, with greater intron 6 PTT in LR versus LM cells (Figure 4F, Figure S7). However, we did not corroborate the late rods’ proportionately higher intron 4 PTT, likely due to the small number of LR cells examined and heterogeneity in PTT frequency in individual cells. Long read sequencing also revealed PTT following the TRβ2-specific exon and a novel transcription-terminating exon following the canonical exon 5 solely in TRβ2 transcripts (Figure 4F, S7). 3’ RACE transcripts rarely extended into the 3’ UTR in LM or LR cells, suggesting that reads mapping to this region do not reflect protein-coding RNAs and confound the assessment of THRB expression. These analyses demonstrate that THRB is regulated by multiple PTT events in rod and cone precursors as well as by cell type-specific promoter utilization and independent 3’ UTR expression.

For RXRG, short read sequencing reads were assigned to several isoforms that differed in their 5’ promoter position and exon utilization (Figure S8A-C). However, 5’ RACE and long-read sequencing did not support differential isoform expression (Figure S8D), and quantitative imaging revealed an average ∼3.5-fold higher RXRγ protein in cones compared to rods (Figure S8E). Thus, long-read sequencing clarified that RXRG expression is moderately higher in human cone versus rod precursors without evidence of cell-type specific isoforms.

Two post-mitotic photoreceptor precursor populations

To further define photoreceptor developmental states, we subdivided the initial six cell populations with higher resolution clustering, identified high resolution cluster-associated genes and regulons, and inferred each cell’s rate and direction of transcriptomic change using RNA velocity (38) (Figure 5A,B). Increased clustering resolution divided the RPC/MG cluster into separate RPC and Müller glia (MG) groups, divided L/M cones into four subgroups (LM1 – LM4), which partially overlapped in UMAP space, and yielded a new cluster, here designated transitional rods (TR), derived from the low resolution iPRP and rod clusters (Figure 5C). Similar clusters were observed at reduced k.param values that defines nearest neighbors.

The distinction between RPCs and MGs was corroborated by the expression of known marker genes, with the RPC cluster having increased expression of cell cycle markers (CCNE2, CCNA2, CCNB2) and the cell cycle-associated PBK and E2F7 (Figure 5D, S9A,B, Table S5). The MG cluster was spatially segregated, lacked expression of most cell cycle genes, and expressed genes known to be expressed in both populations (CCND1, SLC1A3, PAX6, LHX2, VSX2) as well as the MG marker RLBP1 (39–41) (Figure 5C,D, S9A). Differential gene expression analyses revealed upregulation of cell cycle-related genes and G2/M checkpoint and E2F target ontologies in the RPC cluster but no significantly upregulated ontologies in the MG population (Figure S9B,C), consistent with early MGs resembling quiescent RPCs (42).

High resolution clustering also resolved iPRP and TR cells positioned near the RPCs in UMAP space (Figure 5C). RNA velocity suggested that these ‘RPC localized’ cell groups flowed towards the larger iPRP population and towards early maturing rods, respectively (Figure 5B,C). Compared to RPCs and MG, RPC-localized iPRP and TR cells had minimal expression of cyclin RNAs or RPC markers PAX6, LHX2, and VSX2, and had increased expression of photoreceptor determinants OTX2 and CRX (Figure 5D), consistent with their being immediately post-mitotic photoreceptor precursors. The iPRP cells in this region also upregulated the early cell fate determinant ONECUT1 (3, 43), the L/M cone determinant THRB, and the cone differentiation marker GNAT2, but lacked the rod determinant NR2E3, indicative of a cone-bias (Figure 5D). In contrast, TR cells in the same region trended towards higher expression of the rod-related NRL and the rod-specific NR2E3 (p= 0.19, p=0.054, respectively). Prior droplet scRNA-seq analyses did not resolve early post-mitotic cone versus rod precursors (15, 16), nor could we discriminate early post-mitotic THRB+ iPRP and NR2E3+ TR populations through a similar analysis of single cell transcriptomes generated by 3’ end-counting. Thus, deep sequencing enabled detection of THRB, NR2E3, and perhaps other low-expressed RNAs that are important to discrimination of early iPRP and TR populations.

SCENIC regulons further clarified the identities of the four RPC-localized clusters (Figure 5E, S10). The MG cluster was best specified (i.e., had highest regulon specificity scores) by PAX6 and SOX9, both previously described in RPCs and MGs (39, 44, 45), whereas RPCs were best specified by E2F2 and E2F3. These RPC and MG regulons were low or absent in the early post-mitotic iPRP and TR cells, suggesting they undergo an abrupt cell state change. Moreover, immediately post-mitotic iPRPs and TRs that colocalized with RPCs showed differential THRB and NRL regulon activities. iPRP cells had greater THRB regulon signal than TR cells and TR cells had higher NRL regulon signal than iPRPs (p<0.05 for both), consistent with the bifurcated direction of gene expression change implied by RNA velocity (Figure 5C,E). TR cells also had lower activity of the pan-photoreceptor regulons LHX3, OTX2, CRX, and NEUROD1 (Figure S10), consistent with their immature nature.

Overall, high resolution clustering resolved two post-mitotic PR precursor populations with iPRPs having more cone-related and TRs having more rod-related gene expression and regulon activities, consistent with RPC-localized TR cells proceeding towards the ER cluster and iPRP cells continuing towards the larger iPRP population. To our knowledge, this is the first indication of two immediately post-mitotic photoreceptor precursor populations with cone versus rod-biased gene expression (15, 16, 19).

Early cone precursors with cone- and rod-related RNA expression

In addition to the immediately post-mitotic iPRP cells, the iPRP cluster included cells bridging the UMAP region between early maturing cones and early maturing rods (Figure 6A). Cells in this region expressed cone markers (GNAT2, THRB), rod markers (GNAT1, NR2E3), or, in many cases, both (Figure 6B), resembling a proposed hybrid cone/rod precursor (4). To determine whether the intermixed cone and rod gene expression reflect a hybrid cone-rod precursor state, we examined bridge region NRL and THRB regulon activities, which embody the overall transcriptomic effects of these factors. Indeed, in the iPRP bridge and adjacent ER regions, cells with NRL and THRB regulon signals intermixed, and some cells showed strong signals for both (Figure 6C, arrows). However, bridge region iPRP cells had low NRL regulon z-scores similar to L/M cones and had THRB z-scores above that of adjacent ER cells (Figure 6D). The difference between regulon signals increased in rod and cone clusters outside the bridge as NRL or THRB regulons became more dominant (Figure 6D).

Figure 6:

To determine if cells with early cone and rod RNA co-expression can be detected in the developing retina, we performed multiplex RNA FISH for early cone marker GNAT2 and rod marker NR2E3. To identify cell types with such co-expression, FISH analyses were combined with immunostaining for RXRγ, which has high expression in outer neuroblastic layer (NBL) cone precursors and low expression in inner NBL rod precursors (Figure S8E), and for NR2E3, which is detected solely in rods. To infer the spatiotemporal pattern of such expression across human retinal development, we examined GNAT2 and NR2E3 RNA co-expression in outer NBL RXRγ+ cones and in inner NBL RXRγ^lo rods across 13 peripheral fetal retina regions (Figure 6E,F, S11). The analyses revealed that most of the far peripheral (hence, nascent) outer NBL RXRγ+ cones (Figure 6E, regions 2 and 13) lacked detectable GNAT2 and NR2E3 RNA, whereas those in the more mature central retina (regions 3-12) were uniformly GNAT2+ (Figure 6G). However, starting with regions 3 and 11, some outer NBL GNAT2+ cones were also NR2E3+ (Figure 6G). The proportion of GNAT2+, NR2E3+ outer NBL cells increased in the more central retina, yet all of these cells had strong RXRγ staining and lacked NR2E3 protein, consistent with their having a cone identity. In contrast, most of the inner NBL RXRγ^lo cells had either low NR2E3 RNA without NR2E3 protein, in nascent rods, or prominent NR2E3 RNA with NR2E3 protein, in maturing rods (Figure 6H). However, we did not detect GNAT2 puncta in inner NBL RXRγ^lo, NR2E3+ cells, including the most peripheral RXRγ^lo cells with low-level NR2E3 RNA. In summary, most, if not all, photoreceptor precursors with GNAT2 and NR2E3 RNA co-expression had high RXRγ, no detectable NR2E3 protein, and outer NBL positions expected of cone precursors. This supports the notion that early bridge region iPRP cells with combined cone- and rod-related RNA expression are likely cone-directed.

To further assess whether iPRP bridge region cells comprise a distinct subpopulation, we sought to identify RNA markers of this region. Among four iPRP marker genes (Figure S9A), CHRNA1 was mainly expressed in the bridge region and adjacent rods and cones, whereas ONECUT1 was biased to cone-directed iPRPs, lncRNA CTC-378H22.2 was expressed in a narrow zone of THRB+, GNAT2+ cone precursors, and S100A6 was more widely expressed (Figure 7A, S12A). Bridge region cells also expressed ATOH7 and DLL3, which were previously proposed to define transitional photoreceptor precursors and promote cone fate determination (14–16). ATOH7 was largely restricted to cone-directed cells whereas DLL3 was broadly expressed similar to S100A6 (Figure 7A, S12A). As CHRNA1 was the most specific marker of bridge region cone and rod precursors, we evaluated whether CHRNA1 RNA marked specific cone and rod populations using FISH combined with RXRγ and NRL immunofluorescent staining. In a FW12 retina, we detected highest CHRNA1 in the earliest, most peripheral NRL+ rods and RXRγ+ cones and fewer CHRNA1+ cells in more mature, central regions (p <0.005 in cones) (Figure 7B-D), consistent with its expression in bridge region photoreceptor precursors.

Figure 7.

To identify factors that regulate transcriptomic states during iPRP cell fate determination, we examined the distribution of the most iPRP-specific regulons identified after high-resolution clustering, OLIG2 and LHX9 (Figure S12B, Table S6). As Olig2 was previously detected in mouse RPCs in which Onecut1 enabled Thrb expression and cone fate determination (3, 46), the OLIG2 regulon was expected to be most active in RPCs and in early, immediately post-mitotic iPRPs. While the OLIG2 and LHX9 regulons were indeed active in RPCs, both were also active in a narrow zone of ONECUT1+ cone-directed iPRP cells distal to the bridge region and immediately preceding the upregulation of the THRB regulon (Figure S12B). These data indicate that the OLIG2 and LHX9 regulons are active in and potentially relevant to human cone fate determination in post-mitotic photoreceptor precursors.

An early L/M cone trajectory marked by successive lncRNA expression

After fate commitment, L/M cones initiate a maturation process with upregulation of RNAs and proteins related to phototransduction, axonogenesis, synaptogenesis, and outer segment morphogenesis (5, 6). To evaluate whether early L/M cone maturation is comprised of distinct transcriptomic cell states, we assessed marker gene and regulon differences between high resolution clusters LM1, LM2, LM3, and LM4. While clusters LM1-4 were distinguished by sequential increasing expression of ACOT7, RTN1, PDE6H, OLAH, and NPFF, they showed only subtle differences among the three regulons with highest LM1-4 specificity scores, THRB, ISL2, and LHX3 (Figure 8A,B). The lack of cluster-specific marker genes and regulons suggest that LM1-LM4 represent different stages of a graded maturation process.

Figure 8:

We further evaluated gene expression changes in maturing cones by pseudotemporally ordering iPRP and L/M cones (Figure 8C). This identified 967 pseudotime-correlated genes (q-value < 0.05, expression > 0.05 in >5% of cells; Table S7) in seven gene modules (Figure S13). Among the top 20 pseudotime-correlated genes in each module, we identified four lncRNAs that were sequentially expressed (Figure 8D). To determine if these lncRNAs distinguish developmentally distinct maturing cones in vivo, we probed their expression via multiplex RNA FISH in co-stained RXRγ+ cone precursors across a FW16 retina (Figure 8E,F). Quantitation of FISH puncta defined four cone maturation zones based on expression peaks and significant count differences for each lncRNA: the most peripheral cones with high CTC-378H22.2 and HOTAIRM1, peripheral cones with high HOTAIRM1 only, cones with low expression of all four lncRNAs, and parafoveal/foveal cones with upregulated RP13143G15.4, CTD-2034I21.2, and CTC-378H22.2 (Figure 8G). The detection of foveal CTC-378H22 by ISH but not in late LM transcriptomes may relate to a lack of the most mature cones in our scRNA-seq analyses. These data support the concept that maturing L/M cones sequentially express specific lncRNAs as they develop.

Cone intrinsic SYK contributions to the proliferative response to pRB loss

We next assessed whether early maturing cone transcriptomic features are conducive to the proliferative response to pRB loss (10, 11). Given the high transcriptomic similarity across the L/M cone population, we compared gene expression between all early maturing L/M cones and early maturing rods. This identified 422 genes upregulated and 119 downregulated in cones (p<0.05, log2FC>|0.4|) (Figure S14A, Table S8). Among cone-enriched genes, the top three enriched ontologies related to translation initiation, protein localization to membrane, and MYC targets (Figure S14B). Although MYC RNA was not detected, the upregulation of MYC target genes was of interest given the cone precursors’ high MYCN protein expression (Figure 9A) and the role of MYCN in the cone precursor response to pRB loss (8–10). Concordantly, the LM cluster had increased MYCN RNA (log₂FC = 0.54) and MYCN regulon activity, representing the seventh highest LM cluster regulon specificity score (Figure 9B,C, Table S4).

Figure 9:

Among other differentially expressed genes, we noted cone-specific upregulation of SYK (Figure 9D, Table S8), which encodes a non-receptor tyrosine kinase. Whereas SYK was previously implicated in retinoblastoma genesis and proposed to be induced in response to pRB loss (47), its expression was not previously reported in developing fetal retina. Indeed, SYK RNA expression increased at the iPRP stage and into clusters LM3 and 4 and was abolished during the transition from early-maturing to late-maturing L/M cones, in contrast to its minimal expression in the corresponding rod TR, ER, and LR stages (Figure 9E). Similarly, immunohistochemical staining revealed high SYK protein expression in immature (ARR3-) and early maturing (ARR3+) cones from the retinal periphery to the maturing foveal cones at FW16, while SYK was not detected in the most mature foveal cones at FW18 (Figure 9F). The loss of SYK protein and RNA expression with cone maturation is consistent with the lack of SYK in cones of normal retina adjacent to a retinoblastoma tumors (47) and implies that SYK is a defining feature of the human early cone precursor state. While MYCN was also preferentially expressed in early cones, it did not increase as much relative to RPCs nor decline as much in late maturing cone precursors when compared to SYK RNA dynamics (Figure 9E).

To determine if SYK contributes to retinoblastoma initiation, dissociated fetal retinal cells were transduced with an RB1-directed shRNA (shRB1-733) known to induce cone precursor proliferation (10), treated with the highly specific SYK inhibitor GS-9876 (48) for 12 days, and examined for cone precursor cell cycle entry by co-staining for RXRγ and Ki67. GS-9876 treatment suppressed the proliferative response to pRB knockdown at all concentrations from 1.0 – 5.0 μM (Figure 9G), indicating that cone precursor intrinsic SYK activity is critical to the proliferative response to pRB loss and is retained in retinoblastoma cells (Figure 9H).

Discussion

This study evaluated cell state changes associated with human cone and rod photoreceptor development using deep full-length scRNA-seq. Whereas prior scRNA-seq studies provided insights into RPC fate determination and trajectories using 3’ end counting (14–18, 25), deep full-length sequencing enabled more precise resolution of cell states and identification of cell type-specific regulon activities and transcript isoforms. Additionally, our cell enrichment strategy provided insight into the transcriptomic profiles and cancer-predisposing features of developing cones, a rare population with unique cancer cell-of-origin properties (8, 10, 11).

One advantage of full-length scRNA-seq is that it enables the detection of differential transcript isoform expression. Here, we show that the rod determinant NRL, the L/M-cone determinant THRB, and the cone maker RXRG are all co-expressed in cone and rod precursors at the RNA level yet use distinct cell type-specific mechanisms to enable appropriate differential protein expression. Long-read cDNA sequencing revealed that L/M cones preferentially express novel NRL transcripts predicted to encode a truncated NRL (Tr-NRL) lacking a transactivation domain as well as NRL isoforms with intra-exon-2 transcription initiation and alternative splicing. As with DD10 (35), Tr-NRL opposed transactivation by full-length NRL and thus may suppress rod-related transcription. However, our inability to detect intrinsic Tr-NRL or to overexpress Tr-NRL in cone precursors suggests there are additional layers of regulation, possible effects of untranslated Tr-NRL RNA, and alternative contexts in which these isoforms act. A similar assortment of NRL isoforms were expressed in rods, yet a far higher proportion of rod NRL transcripts encoded full-length NRL protein. Thus, L/M cone precursors’ greater use of the Tr-NRL first exon and NRL exon 2 disruptions reveal a multipronged approach to suppress FL-NRL function while downregulating overall NRL RNA by only ∼ 4-fold.

Similarly, we uncovered cell type-specific differences in THRB transcript isoforms, albeit generated through premature transcription termination (PTT) and 3’ UTR expression in late rod precursors. PTT is widely used to govern expression of transcription regulators (49), and while PTT-shortened THRB transcripts had been identified (50), their retinal cell specificity was not previously recognized. In contrast, the expression of THRB 3’ UTR sequences independent of THRB protein-coding sequences was not previously reported and may enable additional regulation (51). As a general matter, independent 3’ UTR transcripts may be inferred from full-length scRNA-seq but may confound the interpretation of scRNA-seq performed with 3’ end-counting. To address this issue, our publicly available Shiny app displays exon usage for all genes expressed in the RPC and photoreceptor clusters in this study.

One limitation of these studies is that individual cells may express a small and varying spectrum of transcript isoforms. While this variability was mitigated by combining cDNAs of cells deemed to be in similar states based on short-read sequencing, analysis of more cells might better reveal the spectrum of isoforms expressed by specific cell populations.

Understanding human photoreceptor development requires the identification of photoreceptor lineage developmental states and elucidation of mechanisms underlying their transitions. While cell states may be defined in various ways, at the transcriptome level they are perhaps best defined by unique combinations of transcription factor regulon activities (29). Our deep sequencing approach uncovered discrete cell states distinguished by gene expression as well as by regulon activities. For example, early and late maturing rod populations were distinguished by the latter’s increased expression of the NRL, ATF4, CRX, and LHX3 regulons and decreased HMX1 and RAX regulons. Similarly, L/M cones formed an early maturing cone cluster characterized by high THRB and ISL2 regulons, consistent with THRB and ISL2 driving L/M cone identity (4, 15, 52), and a late maturing cone group with decreased RAX activity. Downregulation of the RAX regulon in late maturing rod and cone precursors is consistent with decreasing RAX protein during photoreceptor maturation and with RAX-mediated suppression of the cone opsin and rhodopsin promoters (53, 54).

Increased clustering resolution further divided the L/M cone cluster into four subgroups with subtle maturation-associated changes in marker gene expression and regulon activity, consistent with LM1-4 comprising different stages of a graded maturation process. However, trajectory analyses revealed successive expression of lncRNAs that was validated in developing retinal tissue. Three of these lncRNAs (HOTAIRM1, CTD-2034I21.1 (neighbor gene to CTD-2034I21.2), and the mouse ortholog of CTC-378H22.2 (D930028M14Rik)) were previously observed in cone scRNA-seq (15, 17, 55). Their sequential expression in the absence of discrete regulon changes suggests that they demarcate early L/M cone precursor substates.

High resolution clustering also segregated the immature photoreceptor precursor (iPRP) and transitional rod (TR) populations. iPRP and TR cells positioned adjacent to RPC/MG cells in UMAP plots showed higher cone-like or rod-like gene expression and regulon activities, respectively, and had RNA velocities directed toward distinct cone and rod populations. Prior droplet-based scRNA-seq analyses have not to our knowledge distinguished immediately post-mitotic cone- and rod-directed photoreceptor precursors. However, scRNA-seq and scATAC-seq analyses showed distinct early neurogenic precursors directed in part towards cones and late neurogenic precursors directed in part towards rods (14, 15, 19). Thus, as one possibility, early iPRPs may comprise a cone-directed subset of early neurogenic cells while early TRs comprise a rod-directed subset of late neurogenic cells.

Although not identified as a distinct cluster, our analyses also revealed iPRP ‘bridge region’ cells that co-expressed rod and cone genes and regulons. This population resembled the T3 transition and early neurogenic precursor states identified via 3’ end-counting and marked by expression of ATOH7 and DLL3 (15, 16), although the bridge region was more precisely delineated by the new marker CHRNA1. Bridge region cells with velocity directed towards early cones expressed THRB and ONECUT1 RNAs and had OLIG2 regulon activity; as these elements were proposed to promote cone fate in lineage-restricted RPCs (3), their expression in the iPRP population suggests a similar role in post-mitotic photoreceptor precursors. Conversely, the status of bridge region cells with velocity directed towards early rods is less clear; while the velocity vectors and concurrent cone- and rod-related gene expression are consistent with iPRP cells that follow an alternative route to rods, our finding that only outer NBL, RXRγ^hi cells co-express cone marker GNAT2 and rod marker NR2E3 RNAs suggests that most cells with combined cone- and rod-related gene expression are cone-committed. Moreover, the early iPRP and LM cone precursors’ expression of NR2E3 and NRL RNAs suggest that their presence in retinoblastomas (12, 13) reflects their expression in the L/M cone precursor cells of origin.

We also mined gene expression differences in early cone versus rod photoreceptors to identify factors that impact human cone sensitivity to RB1 loss (10, 11). We detected upregulation of genes targeted by MYC along with upregulated MYCN RNA, MYCN protein, and MYCN regulon activity. As MYCN is required for the proliferative response to pRB loss (10) and triggers retinoblastoma when amplified and overexpressed (8), these findings suggest that increased MYCN RNA expression and regulon activity contribute to pRB-deficient cone precursor proliferation and retinoblastoma genesis. We further found that cone precursors highly express SYK, an oncoprotein previously detected in human retinoblastomas and RB1-null retinal organoids but not in healthy tumor-associated retina (47, 56). The high SYK expression preceding early cone precursor maturation and its loss during late maturation implies that it is a retinoblastoma cell-of-origin-specific feature. SYK inhibition impaired pRB-depleted cone precursor cell cycle entry, implying that native SYK expression rather than de novo induction contributes to the cone precursors’ initial proliferation (Figure 9H).

In summary, through deep, full-length RNA sequencing, we identified photoreceptor cell type-specific differences in transcript isoform expression and post-mitotic photoreceptor precursor states with distinctive gene expression and regulon activities. The discrimination of these states enables the identification of developmental stage-specific cone precursor features that underlie retinoblastoma predisposition.

Materials and methods

Primary human tissue

Human fetal samples were provided by Advanced Bioscience Resources, Inc., Alameda, CA, or collected under Institutional Review Board approval at the University of Southern California and Children’s Hospital Los Angeles. Following the patient decision for pregnancy termination, patients were offered the option of donation of the products of conception for research purposes, and those that agreed signed an informed consent. This did not alter the choice of termination procedure, and the products of conception from those who declined participation were disposed of in a standard fashion. Fetal age was determined according to the American College of Obstetrics and Gynecology guidelines (57). Each retina subjected to sequencing was from a different donor.

Retinal dissociation and RNA-sequencing

Single cell isolation

Retinae were dissected and dissociated in ice-cold phosphate-buffered saline (PBS) as described (10). Briefly, retina was removed and placed in 200 µl Earle’s balanced salt solution (EBSS) in a 6-well plate on ice. 2 ml of 10 u/ml papain (Worthington LK003176) solution was added and then incubated at 37°C for 10 min, followed by pipetting with a 1000 µl pipet tip to break up large pieces and additional 5 min incubations until cells were dissociated to single cell level. An additional 1-2 ml papain was added after 20 min if dissociation was insufficient (∼40 min total). Retina culture media [Iscove’s Modified Dulbecco’s Medium (IMDM; Corning), 10% FBS, 0.28 U/mL insulin, 55 µM β-mercaptoethanol (Sigma-Aldrich), glutamine, penicillin, and streptomycin] was used to stop enzyme activity and cells were centrifuged in 14 ml round bottom tubes at 300 x g, 4°C for 10 min. Supernatants were centrifuged at 1100 x g for an additional 3 min. After resuspension in HBSS, tubes were recombined and cells counted.

FACS isolation of single cells and cDNA synthesis

Single cells were FACS-isolated as described (10) with differences as follows. Dissociated cells were centrifuged and resuspended in a 5% fetal bovine serum in PBS (FBS/PBS) with CD133-PE (Miltenyi Biotec 130-080-801), CD44-FITC (BD Pharmingen 555478) and CD49b-FITC antibodies (BD Pharmingen 555498) to a final concentration of 10,000 cells/µl. After incubation at room temperature (RT) for 1 h, samples were diluted in 5% FBS/PBS, centrifuged as above, and resuspended in 300-400 µl 1% FBS/PBS containing 10 µg/ml 4′,6-diamidino-2-phenylindole (DAPI). Single cells were sorted on a BD FACSAria I at 4°C using 100 µm nozzle in single-cell mode into 1.2 µl lysis buffer droplets on parafilm-covered glass slides (S. Lee et al., in preparation). Gating removed low forward- and side-scatter cells before collecting a CD133-high, CD44/49b-low population. Upon collection of eight cells, droplets were transferred to individual low-retention PCR tubes (Bioplastics K69901, B57801) on ice. cDNA was prepared and amplified using the SMART-Seq V4 Ultra Low Input RNA Kit (Takara Bio 634891) in 10x reduced volume reactions using a Mantis liquid handling system (S. Lee et al., in preparation). Samples were stored at −20°C until quantitation and library preparation. All RNA/cDNA volumes during FACS and processing were handled using low retention siliconized pipette tips (VWR 53503-800, 535093-794).

C1 isolation of single cells and cDNA synthesis

Dissociated cells were FACS isolated as above and collected into a single 1.5 ml tube. Cells were centrifuged and resuspended at 400-800 cells/µl, combined with C1 suspension reagent, and loaded onto the C1 chip and imaged at each site to define cell number in each well. cDNA was synthesized using SMARTer chemistry (SMARTer Ultra Low RNA Kit for the Fluidigm C1 System, Clontech 634835/634935) (three retinae, Seq 1, FW17 and FW13-1), or SMART-Seq V4 (all others). Final cDNA was harvested into low-retention tube strips and stored at −20°C.

Quality control, library preparation and sequencing

DNA was quantitated using Quant-iT PicoGreen dsDNA assay and quality checked by Bioanalyzer. Libraries > 0.05 ng/µl were prepared using the Nextera XT DNA Library Preparation Kit (Illumina, FC-131-1096) and sequenced on the Illumina NextSeq 500 (2×75) (Seq 1) or on Illumina HiSeq 4000 (2×75) (all others).]

5’ and 3’ RACE and nanopore sequencing

1 μl of amplified single cell cDNAs libraries produced as above (23 ER cells, 21 LM cells, and 5 LR cells) were re-amplified using the 5’ PCR Primer II A with CloneAmp HiFi PCR Premix (Clontech Laboratories, Inc. Takara Bio #639298, Table S9) by heating to 98°C x 2 min followed by 20 cycles of 98°C x 10 s, 65 °C x 15 s, and 72 °C x 35 s, followed by 72 °C x 5 min. Re-amplified cDNAs were quantified using Qubit^TMFlex Fluorometer (InVitrogen) and 50 pg of each re-amplified library used for separate RACE PCR reactions with gene-specific NRL, RXRγ, TRβ1 and TRβ2 primers and universal 5’ RACE and 3’RACE primers (Table S9) using CloneAmp HiFi PCR Premix by heating to 98°C x 2 min followed by 30 cycles of 98°C x 10 s, 60 °C x 15 s, and 72 °C x 35 s, followed by 72 °C x 5 min. Equal volumes of RACE PCR products were pooled according to cell type, and DNA concentrations of LM (32.8 ng/μl), ER (22.6 ng/μl), and LR (16.5ng/μl) samples determined using Qubit fluorometer. Long read sequencing libraries were prepared and samples barcoded using Oxford Nanopore Technology (ONT) Native Barcoding Kit 24 V14 as specified by the manufacturer for average DNA fragment lengths of ∼ 2kb. 200-300ng of each DNA pool was subjected to repair/dA-tailing (End-Prep) with NEB Ultra II End-Prep Enzyme Mix (New England Biolabs, Ipswich, MA) followed by Native Barcode Ligation as specified by ONT with NB01 assigned to the ER pool, NB02 assigned to the LM pool, and NB03 assigned to the LR pool. Ligation reactions were quenched by addition of 4 μl EDTA and each reaction mixture (24 μl) pooled for a 72 μl volume. The final pool was incubated with Native Adapter (NA) and NEB Quick T4 ligase. The adapter-ligated pool was recovered by AMPure XP magnetic beads and then loaded onto a PromethION flowcell (R10.4.1) as specified by the manufacturer. Sequencing was performed on a PromethION 24 system using MinKNOW software for 72 hours. Raw FAST5 data files were converted to FASTQ files using ONT MinKNOW base calling software.

Computational Methods

Software packages used in this study are described in Table S10.

Read processing and dimensionality reduction

Adapter sequences were removed using the trimgalore wrapper for Cutadapt (58). Trimmed FASTQ files were used as input for HISAT2 (59) with a non-canonical splice penalty of 20, maximum and minimum penalties for mismatch of 0 and 1, and maximum and minimum penalties for soft-clipping set to 3 and 1. Aligned bam files were quantified with StringTie (60) and Tximport (61), yielding cell-by-transcript and cell-by-gene count matrices annotated according to Ensembl build 87. All analyses can be reproduced using a custom snakemake (62) workflow accessible at https://github.com/whtns/ARMOR.

For short-read sequencing, dimensionality reduction and data visualization were performed using the Seurat (V3) toolset (63, 64). Expression counts from all sequencing datasets were integrated using Seurat’s standard integration workflow. Seven sample sequencing batches were integrated after default normalization and scaling using the top 2000 most variable features in each set to identify anchor features. Normalized read counts for gene expression are reported as raw feature read counts divided by total cell read counts, then multiplied by a scaling factor of 10,000 and natural log transformed. These features of the integrated dataset were used to calculate principal component analysis (PCA) and UMAP embeddings (65, 66). A nearest-neighbors graph was constructed from the PCA embedding and clusters were identified using a Louvain algorithm at low and high resolutions (0.4 and 1.6) (67).

Read coverage plots for genes of interest were generated using wiggleplotr (68) to visualize BigWig files with ENSEMBL isoforms. Individual exon counts were identified using DEXseq (69), which takes exons from available isoforms and bins them into intervals before assigning any read counts that overlap a bin to that same bin. These counts were normalized for length for calculation of fold change or relative difference of exon use. Reads crossing target splice donor sites were evaluated as spliced or unspliced from BAM files to determine intron readthrough.

Analysis of long-read sequencing data was performed with the nf-core nanoseq version 3.1.0 pipeline (70) (https://nf-co.re/nanoseq/3.1.0). Quality control on raw reads was conducted using FastQC. Reads were aligned using minimap2 (71). SAM files were converted to coordinate-sorted BAM files and mapping metrics were obtained using samtools. bigWig coverage tracks were created for visualization using BEDTools and bedGraphToBigWig. bigBed coverage tracks were created using BEDTools and bedToBigBed. Transcripts were reconstructed and quantified using bambu. Quality control results for raw reads and alignment results were presented using MultiQC. All reads whose ends overlapped gene-specific RACE primers were displayed and quantitated using the Integrated Genome Browser (IGV). For THRB 3’ RACE reactions, reads inferred to initiate through exonic internal oligo(dT) priming were removed based on the following features adjacent to the 3′ end of the alignments: six continuous adenines (As), more than seven As in a 10 nucleotide window, AG-runs of six or more nucleotides, eight or more As or Gs in a 10 nucleotide window, eight or more As or high A/T content (27 out of 30 bases), or 12 or more adenines present in an 18 nt window (72). The proportion of alternative NRL exon 2 splicing was calculated based on the maximum read coverage between the exon 2 splice acceptor and the first splice donor and the minimum coverage after the internal splice donor as reported in IGV.

Differential expression and overrepresentation analyses

Marker features for each cluster were identified using the Wilcoxon-rank sum tests and specificity scores were computed using the genesorteR (73) R package. Only genes with an adjusted p-value (pAdj) <0.05 were considered for further analyses. Differential expression was performed using the Seurat FindMarkers function with the Wilcoxon-rank sum test on log-normalized counts. Results were displayed as volcano plots made with EnhancedVolcano (74).

Overrepresentation analyses were performed with WebGestalt (75). Analyses were run with the default settings (5-2000 genes per category, Benjamini-Hochberg multiple-testing correction), gene list enrichment was compared to the genome reference set and ontologies were displayed with a false discovery rate (FDR) <0.05 and ‘weighted set cover’ redundancy reduction. All gene sets evaluated were provided within WebGestalt and the same three were used unless otherwise indicated: GO – Biological Process, KEGG, and Hallmark50.

Transcription factor regulons

Transcription factor regulons were identified using pySCENIC, the python-based implementation of Single-Cell Regulatory Network Inference and Clustering (SCENIC) (29, 76). The tool was run using the basic settings as shown (https://pyscenic.readthedocs.io/en/latest/tutorial.html). Initial filtering required a gene to have a minimum of 3 raw counts in 1% of cells to be considered for inclusion in a regulon. AUC scores were z-score normalized for comparison between regulons.

Trajectory analysis

RNA velocity analysis was performed with R package velocyto (38). Spliced and unspliced count matrices were assembled using the run_smartseq2 subcommand with repeat-masked ENSEMBL build 87. RNA velocity was calculated using cell kNN pooling with a gamma fit based on a 95% quantile and velocities were overlaid on integrated UMAP visualizations using velocyto. Pseudotemporal trajectories were calculated using Monocle 3 (77). The SeuratWrappers package provided integration between Seurat’s final integrated dimensional reduction and Monocle. The data was subset to include only cone-directed iPRP and L/M cones, then a principal graph was fit across the UMAP plot and a ‘root cell’ chosen to represent the latest maturation point in the OPN1LW+ late maturing cone group. After assigning pseudotime values to each cell, genes that significantly changed as a function of pseudotime were identified and those with a correlation q-value of <0.05 with expression > 0.5 in at least 5% of the cells present in the pseudotime were grouped into modules of co-regulation by performing UMAP and Louvain community analysis.

Histology

Fixation and cryosectioning of fetal retina

Retinae were procured and dissected in cold 1x PBS. The cornea and lens were removed with a cut around the limbus of the iris leaving the front of the eye open. The tissue was submerged in ∼25 ml of 4% paraformaldehyde (PFA) in 1x TBS and placed at 4°C on a rocker at low speed for ∼16-18 hrs. The tissue was washed 3 times with 1x PBS before incubation in 30% sucrose [sucrose + 1x PBS] overnight at 4°C. A mold was made by cutting off the top 5 ml of a 50 ml conical tube and placed on dry ice. A solution of 1:2 OCT Compound:30% sucrose was mixed and centrifuged to remove bubbles. The mold was partially filled before adding the dissected eye and then covering fully. 10 µm frozen sections were cut at ∼-24°C with the front of the eye facing side on to the blade. Explanted retinae were fixed for 15 min in 4% PFA, washed with 1x PBS, incubated in 30% sucrose for 15 min, after which the membrane was cut from its frame and transferred to a mold containing 1:2 OCT:30% sucrose and frozen.

Immunohistochemical staining

For SYK and ARR3 co-staining, tissue sections were warmed at RT until dry and then washed in Coplin jars twice with 1x PBS for 5 min at RT. Samples were permeabilized in 1x PBS-T for 5 min, washed again in 1x PBS twice for 5 min, blocked for 1 hr in blocking solution [1x PBS, 0.1% Triton X-100, 5% normal donkey serum, 5% normal horse serum]. Primary mouse anti-SYK and rabbit anti-ARR3 (LUMIF-hCAR) antibodies were mixed in blocking solution and applied to samples overnight (16-18 h) at 4°C. Slides were washed three times with 1x PBS for 10 min each, probed with secondary antibodies (Key Resources Table) in blocking solution for 1 hr at RT, washed twice 1x PBS for 15 min each, incubated in1x PBS containing 10 µg/ml DAPI for 10 min, and mounted with Mowiol solution. Additional immunohistological staining of tissue sections and cells on coverslips was performed as described (11) but without initial EDTA wash.

Combined in situ hybridization and immunohistological staining

Human retinal sections were prepared as above but using RNase-free reagents for all stages (PBS, TBS) and ultrapure sucrose (JT Baker, 4097-04) and cryostats cleaned with ELIMINase (Decon Labs, Inc) and 70% ethanol prior to use. RNA FISH hybridization chain reaction probes were designed by and obtained from Molecular Instruments, Inc, with 20 20-bp single strand DNA probes per target RNA sequence except Tr-NRL exon 1 (4 probes), FL-NRL exons 1 and (8 probes), and HOTAIRM1 (14 probes). The accession numbers for target sequences are shown in Table S11, except for RP13-143G15.4 which included regions of several isoforms for maximum gene coverage. The in situ protocol was performed with fluorescent hairpins and buffers from Molecular Instruments using the manufacturer’s instructions for mammalian cells on a slide (https://files.molecularinstruments.com/MI-Protocol-RNAFISH-MammalianCellSlide-Rev7.pdf), beginning with a 4% PFA treatment after thawing and drying slides at RT for 10 min, with the following changes: samples were incubated in 70% ethanol overnight in a Coplin jar, probes were added to the final hybridization volume at 4x the recommended concentration for tissue sections (2.4 pmol in 150 µl of hybridization buffer), tissues were pre-hybridized at 42°C and then lowered to 37°C once probes were added, and hairpins were used at a volume of 1 µl (6 pmol) in 50 µl amplification buffer per section. After removing amplification solution and washing, samples were directly used in the immunofluorescence protocol. Briefly, samples were blocked in a basic RNase-free solution (1% BSA, 0.1% Triton X-100, 0.05% Tween-20 in PBS) for 1 hr and primary antibodies mixed in blocking solution and applied to samples overnight at 4°C. Samples were washed three times for 10 min with 1x PBS, then incubated with secondary antibodies in blocking solution for 40 min at RT. Samples were washed as above with DAPI in the third and final wash, and mounted with Mowiol before imaging using a Leica STELLARIS 5 inverted confocal microscope. For GNAT2 and NRL FISH, samples were washed with DAPI in the third and final wash, mounted and then imaged as above. After this first imaging, coverslip was removed, slides were washed in PBS, and then were used in the immunofluorescence protocol with an alternate blocking solution) (2.5% Horse, donkey, human sera, 1% BSA, 0.1% Triton X-100, 0.05% Tween-20 in TBS).

Quantitation of FISH

For NRL, CHRNA1, and lncRNAs, QuPath (78) was used to identify nuclei, predict cell expansion, classify cells, and count FISH puncta. Apical retina regions containing RXRγ+ and/or NRL+ cells were visually selected for evaluation. The StarDist extension (79) was used to outline nuclei and cell expansion zones with deep learning model ‘dsb2018_heavy_augment.pb’ on the DAPI image channel (threshold = 0.6, pixelSize = 0.15). Model results were reviewed and cells with partial, multiple, or overlapping nuclei removed from consideration. FISH puncta were thresholded then output, using the subcellular detection function (expected spot size = 0.5 µm², min spot size = 0.2 µm², max spot size = 0.7 µm² (NRL isoforms) or 1 µm² (CHRNA1 and lncRNAs) to capture the most puncta above background. Image positions were determined by tracing the retina apical edge with line segments in FIJI, identifying the closest point on the line to the midpoint of each image, then measuring the distance from the line starting position. For GNAT2 and NR2E3 FISH, images were analyzed in FIJI. The outer layer of cells was delimited manually based on the RXRγ+ cells and the number of RXRγ+ and/or GNAT2+ and/or NR2E3+ cells were manually counted. For the inner layer, the number of NR2E3 RNA positive cells and NR2E3 protein positive cells were counted.

NRL isoform analyses

Proteasome inhibition of CHLA-VC-RB-31

10⁶ CHLA-VC-RB-31 cells (36) were cultured in 2 ml retina culture medium in a 24-well plate with or without 10 µM retinoic acid (Cayman Chemical 11017) and incubated at 37°C in 5% CO₂. MG132 (Sigma-Aldrich 474790) was added at 23 h to a final concentration of 10 µM and cells incubated for 4 h. Cells were collected and centrifuged at 300 x g for 4 min, resuspended in 60 µl cold 1x RIPA buffer (Millipore 20-188), 10% sodium dodecyl sulfate (SDS), protease and phosphatase inhibitors diluted per manufacturer instructions (Sigma-Aldrich 5892970001, 4906837001)], incubated on ice for 30 min, then centrifuged at 4°C for 20 min at 14,000 RPM. 15 µl of supernatant were mixed with 3 µl 6X sample buffer (300 mM Tris, 60% glycerol,12% SDS, 86 mM B-mercaptoethanol, 0.6 mg/ml bromophenol blue), heated to 95°C for 5 min, and separated on a 4-12% Bis-Tris gel (Invitrogen, NP0321) using 1xMOPS running buffer (Life Technologies NP0001). Protein was wet transferred to a PVDF membrane in Towbin transfer buffer (25 mM Tris, 192 mM glycine pH 8.3,10% methanol) at 20 V ON at RT. The membrane was blocked in 1xTBS-T containing 5% w/v nonfat dry milk for 1 hr, washed with TBS-T three times for 5 min, incubated overnight in 0.1% dry milk TBS-T solution containing primary antibody at 4°C with slow rocking, washed three times for 15 min in TBS-T, incubated with HRP-conjugated secondary antibody for 1 hr at RT in the 0.1% milk solution, washed six times for 15 min, incubated in chemiluminescent substrate (Thermo Fisher, 34094) and imaged.

NRL isoform and reporter constructs

The pcDNA4-His-Max-C-Nrl (80) (gift of A. Swaroop) CMV promoter was replaced with EF1α promoter to produce pcDNA4-His-Max-C-EF1a-FL-NRL. Using In-Fusion (Takara Bio) the His-Max tag was removed to produce full-length NRL expression plasmid pCDNA4-C-EF1a-FL-NRL (primers #1,2). The His-Max tag plus first 417 bp of NRL N-terminal coding sequence were removed to produce truncated NRL expression plasmid pCDNA4-C-EF1a-Tr-NRL (primers #1,2,3,4). pcDNA-C-EF1a-empty transfection carrier plasmid was made by removing NRL coding sequence from pCDNA4-C-EF1a-Fl-NRL (primers #5,6). pGL3-PDE6B-146-luc was produced by amplifying a 146 bp promoter region of PDE6B (81) from H9 hESC genomic DNA and inserting into a pGL3-SV40 (Promega, primers #7,8), then removing the SV40 promoter (primers #9,10) to generate pGL3-PDE6B-146. The promoter-less pGL3-empty was produced by removing SV40 sequences from pGL3-SV40 (primers #9,11). pUltra-EGFP-FL-NRL was produced by amplification of the FL-NRL coding sequence from pcDNA4-C-EF1α-FL-NRL and inserting into BamHI digested pUltra-EGFP (Addgene 24129) (primers #12, 13). All primers are listed in the Key Resources Table.

NRL expression and analysis

pcDNA4-C-EF1α-NRL constructs were transfected into HEK293T cells. 4 µg plasmid in 1 ml DMEM was combined with 1 ml DMEM + 12 µg PEI and combined with 10⁶ trypsinized HEK293T cells, plated into two wells of a 6-well plate with an additional 1 ml of DMEM + 10% FBS, then cultured overnight in a 37°C incubator. Media was replaced with fresh 293T media the next day and after two additional days cells were collected, lysed, and protein quantitated via BCA assay (Thermo Scientific PI23227). 30 µg of each well was used for western blot as above, with membranes probed with antibody raised against the NRL N-terminus (Santa Cruz SC-374277) or against full-length NRL (R&D AF2945).

Luciferase assays

575 ng of total DNA (50 ng reporter, 100ng FL-NRL and/or 100-300 ng Tr-NRL expression vector, and the remainder pcDNA4C-EF1α empty vector) in 5 µl DMEM was mixed with 2.3 µg PEI in another 5 µl DMEM (per condition, in triplicate) and then plated into one well of a white bottom 96-well cell culture plate with 50,000 NIH3T3 cells in 100 µl 293T media. The luciferase assay was carried out 48 h after transfection with the Nano-Glo Dual-Luciferase Reporter Assay System (Promega N1610) per manufacturer instructions. Luminescence was recorded using the Promega Glomax Multi+ Detection System.

Truncated NRL overexpression in fetal retina

Concentrated pUltra-EGFP-NRL-205 and control pUltra-EGFP lentivirus were produced as described (82). At 16-24 h after 293T transfection in 15 cm dishes, media was exchanged for 20 mL UltraCULTURE (Lonza Inc.) containing 1% HEPES, 1% GlutaMax and 1% Penicillin-Streptomycin. After 60-64 h, media was harvested, centrifuged at 3000 x g for 10 min at 4°C, filtered through 0.45 µm PVDF flask filter, concentrated via tangential flow filtration using a Midikros 20 cm 500 KD 0.5 mm column (D02-5500-05-S), with both input and final supernatant containers kept in ice. The concentrate was re-filtered using a 0.45 µm PVDF syringe filter (Millipore), stored in aliquots at −80°C, and titered using a p24 ELISA kit (ZeptoMetrix 801002).

As described (11), fetal eyes were washed in 70% ethanol, washed three times in 1x sterile PBS, and dissected in cold 1x PBS. The cornea was removed and tissue cut ∼one-third of the distance towards the posterior pole along 4 lines approximately equidistant and between attached tendons. The retina was removed in this flattened state after cutting the optic nerve, transferred to a polytetrafluoroethylene culture insert (Millipore, PICM0RG50) with photoreceptor side down, and placed in a 6-well plate with 1,200 µl retinal culture medium. Retinae were infected as described, cultured for 7 days with half media changes every other day, then fixed in 4% PFA for 15 min, washed three times in 1x PBS, equilibrated in sucrose (30% in PBS) for 15 min, and embedded and sectioned as above. Immunostaining and imaging were performed as above.

RB1 knockdown and SYK inhibitor treatment

FW18.5 retina was partially dissociated with papain, cultured overnight in retinal culture media at 37°C in a 6-well plate, then frozen in 10% DMSO solution and stored in liquid nitrogen. Samples were thawed and revived in retinal culture media overnight and infected with lentivirus carrying shRNAs targeting RB1 (shRB1-733) (10) or a scrambled control (shSCR) and diluted β-mercaptoethanol, insulin, glutamine, penicillin, and streptomycin to the final retina culture medium concentrations along with 5 µg/ml Polybrene, for 24 hrs. After replacing 2/3 of media in each well with fresh media, cells were treated with GS-9876 (MedChemExpress, CAT# HY-109091) in DMSO at or with the same volume of DMSO vehicle control. After 12 d, cells were attached to poly-L-Lysine-coated coverslips for 3 h, fixed in 4% PFA for 10 min, washed in 1x PBS three times, and stored at −20°C until immunocytological staining.

Statistical Analysis

Statistical methods and packages are cited in the Methods text and figure legends. Statistical comparison between more than two groups of cells for gene expression, regulon activity, or exon count proportions used the non-parametric Kruskal-Wallis test followed by pairwise post-hoc Dunn tests with Benjamini-Hochberg correction. Mean expression levels of NRL and RXRG isoforms across cell clusters were evaluated by Welch’s t-tests for groups with unequal variance; p-values for these tests were estimated using bootstrapping with 1,000,000 replications per comparison. Changes in FISH puncta counts were examined within each channel individually using zero-inflated negative binomial regression models with robust standard error terms. Distance across the sample tissue was used to model degenerate zero inflation. P values reported in the manuscript text and legends are Wald Chi-Squared tests that evaluate differences in count between inflection points. Analyses were performed using Microsoft Excel, the R statistical language in RStudio, or Stata SE v14.2.

Acknowledgements

We thank Melissa L. Wilson (USC Department of Preventive Medicine) and Family Planning Associates for assistance in obtaining fetal tissue, the Stem Cell Analytics Core Facility, FACS Core Facility, and Imaging Core Facility of The Saban Research Institute of Children’s Hospital Los Angeles, and Biraj Mahato and Sumitha Bharathan for critical reading an earlier version of the manuscript.

Additional information

Funding

National Institutes of Health grant R01EY026661 (DC)

National Institutes of Health grant R01CA137124 (DC)

National Institutes of Health grant 5T32HD060549 to USC Department of Development, Stem Cells, and Regenerative Medicine (DWHS)

USC Provost Fellowship (DWHS)

Saban Research Institute of Children’s Hospital Los Angeles fellowship (JB)

Research to Prevent Blindness (unrestricted grant to USC Dept. of Ophthalmology)

Larry and Celia Moh Foundation (DC)

Neonatal Blindness Research Fund (DC)

A.B. Reins Foundation (DC)

Knights Templar Eye Foundation (DC)

Author contributions

Conceptualization: DWHS, DC

Methodology: DWHS, LC, DJW, HPS, DC

Investigation: DWHS, LC, SL, JB, DJW, JGA, YK, HPS, ZF

Computational analyses: KS, BB, MS

Statistical analyses: DWHS, MWR

Bioinformatic approaches: KS, MB, MAB

Computational infrastructure: ED, SGE

Tissue procurement: MT, BHG

Key reagent: CMC

Writing --- original draft: DWHS, DC

Writing --- review and editing: DWHS, KS, JGA, MAB, BHG, CMC, DC

Competing interests

The authors declare that they have no competing interests.

Data and materials availability

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. A web portal providing access to a Shiny app for manipulation of the presented single-cell data (including UMAP or tSNE display of clusters, gene or isoform expression, and sample metadata; cluster marker gene and violin plots; exon coverage plots; Ensembl isoform assignments) is available at https://docker.saban.chla.usc.edu/cobrinik/app/seuratApp/. Single-cell RNA-sequencing files and processed final gene and transcript matrices are publicly available at the GEO database under accession number GSE207802.

Note

This reviewed preprint has been updated to use the correct text; previously, a version prior to the one reviewed was presented here.

Supplementary Materials

Supplementary Figures

Figure S1:

Figure S2:

Figure S3.

Figure S4:

Figure S5:

Figure S6:

Figure S7:

Figure S8:

Figure S9:

Figure S10:

Figure S11:

Figure S12:

Figure S13:

Figure S14:

Key Resources Table

Information related to strains, cell lines, biological samples, antibodies, recombinant DNA reagents, sequence-based reagents, commercial assays or kits, and software.