Gene regulatory behaviors of mutant CRX proteins demonstrated with late-stage human retinal organoids

Multimodal single-cell sequencing of late-stage LCA7 retinal organoids

We used previously described iPSC lines [24] derived from two LCA7 patients with severe phenotypes (Table 1) to investigate the gene regulatory effects of two CRX variants when expressed in developing photoreceptors, as well as the consequence of knockout of the mutant allele for one variant as an approach to phenotypic rescue (position of variants shown in Fig. 1A). One variant studied (K88Q) causes an amino acid substitution in the DNA-binding homeodomain while the other (T155fs) causes a frameshift and early truncation within the transactivation domain. A CRISPR-corrected control with two intact CRX alleles was generated from the K88Q line (referred to as ‘control’). Organoids from each line were differentiated to 240 days at which time nuclei were isolated and barcoded for joint single-cell RNA and ATAC sequencing (Fig. 1B). After computational filtering for high-quality nuclei, 10 485 total nuclei from four samples were included in the final dataset (Table 1). A median of 8160 unique transcripts and a median of 4322 genes were detected across all cells included in the analysis.

Cell type annotation of single cells from retinal organoids

Gene expression and chromatin accessibility data for each cell in the dataset were used in combination to drive clustering of similar cells across samples with Weighted Nearest Neighbor (WNN) analysis as previously described [25]. Cell types were annotated based on expression of canonical marker genes within these clusters (Fig. 1C, Fig. S2). Clusters that represented likely doublets expressing markers of multiple cell types were excluded from further analysis. As expected, our late-stage retinal organoids yielded primarily photoreceptor cells and glia, consistent with previous studies at comparable timepoints [21, 26–30] (Fig. 1C). The T155fs line produced a distinct population of glia not observed in other lines (Fig. 1D). To confirm that cell type clustering was not an artifact resulting from low-quality nuclei, post-filtering quality control metrics across annotated clusters and samples were visualized (Fig. S1). These metrics show no systemic deviation across samples, including within the glia cluster from the T155fs line.

CRX patient lines and edited controls produce expected major inner and outer neural retina cell types following 240 days of culture

Using the cell type annotations shown in Fig. 1C, we next calculated the fraction of each cell type out of all cells captured from each sample, yielding a cell type proportion of analyzed cells (Fig. 1E). While all lines produced all annotated cell types (apart from K88Q which did not produce S cone photoreceptors), T155fs line yielded a higher proportion of glia. K88Q line yielded an increased proportion of cone photoreceptors at the apparent expense of rods, as compared to the isogenic control and K88Qko lines. No ganglion cells were observed in any lines, which is typical for 240-day-old organoids. All proportions are relative only to cells analyzed as absolute proportions and cell type counts are not attainable with scRNAseq.

Expression of CRX is not significantly altered by CRX allelic status

We examined the expression level of CRX mRNA within photoreceptor cells from each line (Fig. 1F) to determine whether pathogenic CRX alleles diminished expression of mutant or wildtype CRX at the transcript level. While a range of expression in individual cells was found across samples (likely due, in part, to sub-optimal sampling inherent to the single-cell sequencing approach), the median expression value was not significantly altered between samples carrying wildtype, mutant, or knock-out alleles (Fig. 1E). For the purposes of this analysis, all CRX transcripts are mapped and counted towards the common CRX feature, meaning that transcripts originating from the K88Qko allele were included in this analysis.

CRX T155fs and K88Q produce distinct changes in photoreceptor gene expression

We next focused on the gene expression profile of rod and cone photoreceptors due to CRX’s known function as an essential photoreceptor transcription factor. First, we assessed how CRX allelic status alters photoreceptor gene expression. Differential gene expression analysis was performed between the isogenic control line versus the T155fs, K88Q, and K88Qko lines across both the rod and cone clusters (Fig. 2A–C, and  G–I). For this analysis, genes with an absolute log₂(fold change) greater than 1, adjusted p-value < 0.05, and detected expression in at least 10% of cells in either condition were considered as differentially expressed genes (DEGs). Within the rod cluster, T155fs, K88Q, and K88Qko exhibited 974, 1100, and 278 DEGs respectively (Fig. 2A–C). Within the cone cluster, T155fs, K88Q, and K88Qko exhibited 1082, 1215, and 604 DEGs respectively (Fig. 2E–G). The number of differentially expressed genes between each sample and control line, as well as the intersection of these sets was calculated (Fig. 2D, E, J, and  K).

Knockout of CRX K88Q partially restores normal photoreceptor gene expression

In both photoreceptor clusters, the K88Qko line was more transcriptionally similar to the isogenic control line than the unedited K88Q line (Fig. 2D, E, J, and  K), implying that K88Q DEGs were mostly due to a gain of abnormal function due to the CRX K88Q variant. However, a minority of genes remain differentially expressed following knockout, suggesting a set of genes sensitive to CRX haploinsufficiency. The majority of DEGs found in the K88Qko line were shared with one or both mutant lines (Fig. 2D and J).

CRX variants impair expression of photoreceptor maturity marker genes in an allele-specific manner

CRX is known to regulate expression of essential photoreceptor genes during retinal development. Using the rod and cone populations, we next examined the expression of genes involved in photoreceptor maturation and function across the four lines. Within the rod cluster (Fig. 2D), two genes that were well expressed in the K88Q and control lines were poorly expressed in the T155fs line (TULP1 and IMPG1). Two genes were poorly expressed in both mutant lines compared to controls (SAG and GNAT1). Within the cone cluster (Fig. 2H), PODXL, GNGT2, and RDH10 were all expressed at a lower level in K88Q than the other mutant or control lines. Both T155fs and K88Q exhibited lower GNAT2 and SAMD7 expression compared to control. Most curated photoreceptor maturity markers of both rod and cone maturation were well expressed in mutant and control lines.

CRX dysfunction alters photoreceptor chromatin accessibility

Regions of chromatin accessibility (‘peaks’) were defined by overlapped ATAC reads above a background level as previously described [31]. All cells of each cell type were pooled across lines for this analysis to reduce calling of genotype-specific peaks which would be difficult to distinguish from technical artifact. Accessibility at these peaks was then compared between samples across each cell type cluster (Fig. 3A and B). Aside from the large number of differentially accessible regions in the glia from the T155fs line (Fig. 3B), most of the variation in chromatin accessibility was restricted to rod and cone clusters (Fig. 3A). The K88Qko line displayed reduced deviation of accessibility compared to the isogenic control line in both clusters, highlighting the dominant negative functionality of this variant and corroborating what was observed at the transcript-level.

We next compared differentially accessible regions (DARs) of chromatin from each sample to determine whether these were largely overlapping or allele specific. Within the photoreceptor clusters (i.e. rod and cone clusters), most variation in chromatin accessibility was restricted to a single sample. The T155fs and K88Q lines each respectively had 353 and 544 differentially accessible regions compared to the same control (Fig. 3C). A total of 148 regions of differential accessibility were shared between these lines (Fig. 3C, in yellow). All but 2 of these peaks exhibited the same directionality of differential expression in both T155fs and K88Q lines (Fig. 3D), implying a core set of responsive regions downstream of normal CRX function. DARs identified in each comparison are highlighted in Fig. 3E and F, with common DARs shown in yellow. Both CRX T155fs and K88Q variants lead to both increased and decreased accessibility at several sites across the genome compared to the control line.

To evaluate whether abnormal CRX function leads to a stalling of normal photoreceptor development, we compared chromatin accessibility between rod and cone clusters and the photoreceptor progenitor cluster within each line (Fig. 3G). We found far more DARs in the control and K88Qko comparisons than within the K88Q and T155fs (Fig. 3H), indicating that the chromatin accessibility landscape of photoreceptors with mutant CRX may fail to mature in a typical fashion. Knockout of the K88Q allele resulted in DARs largely overlapping with those observed in the isogenic control.

The K88Q substitution in CRX causes increased accessibility of VSX2-family motifs

We next asked how the K88Q variant alters normal CRX target gene regulation at the transcript and chromatin level. This variant leads to an amino acid substitution from lysine to glutamine at relative position 50 within the DNA-binding homeodomain of CRX [32, 33]. Substitution at this position is known to discriminate CRX and OTX2 factors from RAX/VSX2/ARX family factors [34] and is important to the differentiation of photoreceptor and bipolar cell gene expression wherein Q50 homeodomain factors normally suppress expression of photoreceptor-specific genes [34]. The consensus motifs for CRX (Fig. 4B), VSX2 (Fig. 4B) differ based on this homeodomain difference.

We used gene expression and chromatin accessibility data from photoreceptors carrying the CRX K88Q allele to ask how this disease-causing variant leads to malformation of photoreceptors. We first identified the set of DARs between each line and control within rod and cone clusters as described above (i.e. those in Fig. 3A). We next computed enriched and de-enriched transcription factor binding motifs within each set of peaks compared to a set of background peaks matched for GC-content (Fig. 4C and D). In both mutant lines (T155fs and K88Q) we found that CRX and OTX2 motifs were de-enriched, while the K88Qko line did not display de-enrichment of these motifs within its differentially accessible peaks (Fig. 4C). We examined motifs that were enriched in each differentially accessible peak set. RAX/ARX/VSX2 motifs were significantly enriched within the K88Q line but not enriched in the other tested lines (Fig. 4D).

CRX variant presence does not alter CRX or OTX2 expression in developing photoreceptors

To confirm that our organoids exhibited the expected developmental pattern of CRX, OTX2, and VSX2 gene expression, we examined expression of these genes, as well as ARX and RAX within cone, rod, and bipolar clusters for all 4 lines (Fig. 4E–I, Fig. S4). OTX2 was expressed in photoreceptors and bipolars in all lines, while CRX was primarily restricted to photoreceptors (Fig. 4E and F). Notably, VSX2 was expressed only in bipolars in all lines including K88Q (Fig. 4G). As VSX2 transcript was not detected in photoreceptor cells, the aberrant enrichment of VSX2 motifs observed in this cluster (Fig. 4D) is likely not due to the presence of VSX2 itself but instead may be due to inappropriate binding of mutant CRX to this motif.

VSX2 motif activity and gene expression in single cells

We next applied a per-cell approached to investigate the co-occurrence of transcription factor gene expression and motif accessibility leveraging our multiomic dataset. Motif activity scores were generated for each cell in the dataset using ChromVAR [35]. The per-cell motif activity score for VSX2 as well as the magnitude of VSX2 RNA expression are shown (Fig. 4J–M). Normally, neither VSX2 expression nor motif activity are detectable in photoreceptor cells (Fig. 4J). In the K88Q line, VSX2 motif activity is seen in both rod and cone clusters (Fig. 4K). Notably, these clusters do not demonstrate expression of VSX2 RNA. No VSX2 motif activity is seen in T155fs photoreceptors (Fig. 4L), and knockout of K88Q also abolishes presumed aberrant VSX2 motif activity (Fig. 4M). CRX motif chromVAR activity was decreased in cones in mutant lines and decreased in rods in the T155fs line (Fig. S3). K88Q rods and K88Qko rods and cones exhibited comparable activity levels to the control line. Bipolar cells from all samples exhibited lower CRX activity compared to photoreceptors.

Retinal organoid differentiation

iPSC lines were cultured as previously described [24] and cryopreserved in mTeSR Plus medium (Stemcell Technologies, 5825) with 10% DMSO. For each line (control, T155fs, K88Q, and K88Qko), retinal organoid differentiation was initiated using cells at passage 56, 35, 30, and 65 respectively. iPSC quality control was confirmed as previously described [24]. Organoids were generated using a previously published protocol [24]. At day 240 (D240) of differentiation, organoids were snap frozen (without culture medium) and stored at −80°C until processing for single nucleus sequencing, as described below.

Nucleus isolation from organoids

Nuclei were isolated from snap-frozen organoids using a protocol based on 10× Genomics’ demonstrated protocol CG000366. In brief, chilled 0.1× Lysis Buffer (Tris–HCl (10 mM), NaCl (10 mM), MgCl2 (3 mM), Tween-20 (0.01%), Nonidet P40 Substitute (0.01%), Digitonin (0.001%), BSA (1%), DTT (1 mM), RNase inhibitor (1 U/μL)) was added to each sample. Samples were homogenized using a sterile pestle. The resulting mixture was incubated on ice for 5 minutes, then mixed 10 times with a P1000 pipette tip. The mixture was incubated on ice for 8.5 minutes. 500 μL chilled Wash Buffer (Tris–HCl (10 mM), NaCl (10 mM), MgCl2 (3 mM), BSA (1%), Tween-20 (0.1%), DTT (1 mM), RNase inhibitor (1 U/μL)) was then added and mixed with a P1000 pipette five times. The nuclei were pelleted at 500× g for 5 minutes at 4°C. Washing was repeated in this manner for a total of three times. Following the final centrifugation, wash buffer was removed, and nuclei were resuspended in Diluted Nuclei Buffer (10× Genomics Nuclei Buffer (1× from 20× stock)), DTT (1 mM), RNase inhibitor (1 U/μL). Nuclei were counted manually using a hemocytometer and processed immediately using the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression v1 system (10× Genomics).

Library preparation, sequencing, and data processing

Final libraries were generated according to Chromium Next GEM Single Cell Multiome ATAC + Gene Expression User Guide CG000338 (10× Genomics, Rev E) without deviation. The libraries were quantified using the Qubit dsDNA HS Assay Kit and diluted to 3 ng/μL. Library quality was checked using the Bioanalyzer or Tapestation (Agilent). Raw sequencing reads were mapped to the GRCh38 reference genome provided by 10× Genomics (GENCODE Release 32 (GRCh38.p13)) without modification using CellRanger ARC v2.0.0 (10× Genomics) using default settings.

Sample integration and clustering

Gene expression and chromatin accessibility datasets from the four samples were filtered for high quality nuclei and integrated count matrices were generated. Gene expression data were integrated using canonical correlation analysis using Seurat [41] and chromatin accessibility data were integrated using latent semantic indexing using Signac [42]. Dimensionality reduction was performed on both integrated datasets separately. Only nuclei passing quality control thresholds in both gene expression and ATAC modalities were included in the final, multimodal dataset.

RNA and ATAC data for each nucleus was used to drive clustering using weighted nearest neighbor analysis (Signac). Dimensions 1:30 of the gene expression PCA and dimensions 2:50 of the chromatin accessibility LSI analysis were taken as WNN input. The resulting ‘weighted’ shared nearest neighbor graph was used to identify clusters of nuclei using the smart local moving (SLM) algorithm [41] with a resolution of 0.5.

Cluster cell type annotation

Annotation of the above clusters was determined manually based on marker gene expression and comparison with previously published single-cell RNA sequencing data from developing human fetal retina [21]. The identity of three clusters was not determined. DoubletFinder [43] was used to identify one of these clusters as likely doublets originating from photoreceptors and glia cells. This and the other two unidentified clusters were dropped from further analysis.

ATAC peak calling, differential chromatin accessibility analysis, and motif enrichment

ATAC peaks to be used in motif enrichment analysis were first called using MACS2 [31, 44] as implemented by the CallPeaks function in Signac. Cells were grouped by cell type for peak calling. Cells from all lines were pooled for peak calling. To identify transcription factor binding motif enrichment, the FindMarkers function was used with the following parameters: only.pos = FALSE, test.use = LR, latent.vars = nCount_macs2_peaks. Differentially accessible regions were filtered for an adjusted p-value less than 0.05 and absolute Log2(Fold change) > 1. Enriched motifs were identified in differentially accessible regions using the FindMotifs function, with a set of 50 000 GC-matched peaks accessible in the same cell type used as background. Enrichment over background and significance of motifs were plotted, highlighting the motifs of interest. ChromVAR was used to generate per-cell motif activity scores as implemented by Signac [42].

Differential gene expression analysis

The FindMarkers function (Seurat) identified differentially expressed genes between samples within each cluster tested. Genes with a log2(fold change) of at least 1 and an adjusted p-value below 0.05 were considered significantly differentially expressed.

Trajectory construction with slingshot

Slingshot [45] was used to quantify differing maturity within the photoreceptor and photoreceptor progenitor clusters. The Slingshot function was run taking WNN UMAP-derived embeddings as ‘input.’ The starting (progenitor) and end (cone, rod) clusters were given, and the following parameters were used: stretch = 0.1, thresh = 0.01. Two trajectories were identified. Curves were drawn based on these trajectories from which pseudotime values were derived and given to each cell.

Statistics

All statistical testing for differential gene expression and chromatin accessibility was implemented using the Seurat/Signac package in R [41, 42]. For gene expression analysis, significantly differentially expressed genes were defined as those with absolute log2(fold change) > 1 and adjusted P-value < 0.05 based on the Wilcoxon Rank Sum test with Bonferroni correction using all genes. Only genes with detectable expression in at least 10% of cells in either comparison group were considered. For ATAC-seq analysis, significantly differentially accessible regions were defined as those with absolute log2(fold change) > 1 and adjusted P-value < 0.05 based on logistic regression model with likelihood ratio testing adjusted with Bonferroni correction using all peaks. Only peaks identified in at least 1% of cells in either comparison group were considered. For all experiments, 4 organoids were used per sample. Organoids were pooled prior to library preparation resulting in one sequencing library per sample. For statistical tests, individual cells were used as replicates.