RNA-seq data from CECs reveal CTG18.1 expansion-mediated FECD to be transcriptomically distinct from non-expanded FECD

To further increase our understanding of CTG18.1 expansion-associated splicing dysregulation events, and better understand the relationship between individual events to transcript-specific isoforms, we generated long-read RNA-seq data from biologically independent unpassaged primary CEC cultures derived from healthy control (n = 4, mean age 68.8 years) and Exp+ FECD (mono-allelic CTG18.1 allele repeat length ≥50) (n = 4, mean age 66.8 years) tissues (Tables 1 and S1). After alignment and stringent quality control, reads that displayed evidence of 5’ degradation, intrapriming and reverse transcription template switching were removed. Exp+ and control samples were subsequently analysed with SQANTI3 to assess isoform diversity (Fig 1A and 1B and Tables 2 and S2). We also generated short-read RNA-seq data from unpassaged primary CEC cultures derived from healthy control (n = 4, mean age 58.5 years), Exp+ FECD (mono-allelic CTG18.1 allele repeat length ≥50) (n = 3, mean age 63.3 years) and Exp- FECD (bi-allelic CTG18.1 allele repeat lengths ≤30) (n = 3, mean age 63.3 years) tissues (Table 1). Exome data generated from the Exp- FECD cases utilised in this transcriptomic work were analysed for the presence of variants in any previously FECD-associated genes, such as TCF4, COL8A2, SLC4A11, ZEB1, and AGBL1 [4], to identify any potential underlying genetic cause. No rare (MAF <0.005) or potentially deleterious (CADD >15) variants were identified. Hierarchical clustering using heatmap visualisation and PCA analysis including sex and ethnicity as covariates revealed that Exp+ samples grouped separately from not only the controls, but also the Exp- samples (Fig 1C and 1D). Driving genes and additional principal components are listed in S2 Table. Notably, Exp- samples also clustered together on the PCA plot, but the heatmap derived from hierarchical clustering highlighted that greater intra-group variability existed between Exp- samples compared to variability within control and Exp+ groups (Fig 1C). We hypothesise that the increased transcriptomic diversity within the Exp- group is due to potentially distinct, as yet unidentified, genetic causes driving disease for each FECD case-derived sample. We included all three groups in the downstream comparative pipelines using the following pairwise comparisons; (PWC1) control versus Exp+ FECD, (PWC2) Exp+ FECD versus Exp- FECD, and (PWC3) control versus Exp- FECD.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Exp+ FECD samples are transcriptomically distinct in both short- and long-read datasets when compared to controls.

A) Bar chart showing proportion of annotated genes in long read transcriptomic dataset with only one transcript detect vs 2 or more transcripts in Control and Exp+ samples. Exp+ long-read samples have a higher proportion of genes with more than one isoform identified. B) Proportion of annotated transcripts with full splice matches. Exp+ FECD and Control samples have similar percentages of full splice matches suggesting additional splicing diversity seen in Exp+ is not predominantly due to aberrant splicing. C) Hierarchical clustering heatmap of sample-to-sample distances within short-read RNA-Seq data. Darker colours denote higher similarities between samples, demonstrating that Exp+ and control are more similar within each subgroup than with other samples. Exp- shows a similar result but is less pronounced, likely due to differences in unidentified underlying genetic causes. D) Principal component plot (PCA) of all short-read samples with sex, disease state and ethnicity as covariates. All three subsets of data cluster separately in the PCA plot. Control samples appear as circles, FECD Exp+ appear as orange triangles, FECD Exp- appear as green triangles. Male individuals are represented as filled shapes while female individuals are outlined.

https://doi.org/10.1371/journal.pgen.1011230.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of FECD case and control-derived corneal endothelial cell (CEC) cultures analysed by short and long-read RNA-seq, and spatial transcriptomics.

https://doi.org/10.1371/journal.pgen.1011230.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of SQANTI3 analysis of IsoSeq data.

In-depth characterisation of isoforms in control and Exp+ long-read RNA-seq data utilising SQANTI3 pipeline following stringent quality control. SQANTI3 subdivides splice junctions into known (present in reference transcriptome) and novel (not in reference transcriptome). Canonical splice junctions refer to the most common type of splice junction (GTAG, GCAG, ATAC)), while non-canonical refers to other rarer and tissue specific splice junctions with different donor and acceptor sites. SQANTI3 defines novel as genes or splice junctions that are absent from the reference transcriptome.

https://doi.org/10.1371/journal.pgen.1011230.t002

Alternative splicing analysis of long-read and short-read transcriptomic data demonstrate the complexity of increased splicing events underlying CTG18.1 expansion-mediated disease in CECs

Iso-seq and subsequent SQANTI3 analysis (Tables 2 and S3) revealed a higher proportion of isoforms per gene in Exp+ samples (1.7 unique isoforms per unique gene in Exp+ vs 1.3 in control). Additionally, we see a higher proportion of known non-canonical splice junctions in Exp+ samples. When characterising transcripts based on splice junctions, we see similar proportions of isoform subtypes (S3 Table). This is consistent with previous work [6,10–12] and validated by our short read data, discussed below, (S1 and S2 Figs), whereby Exp+ samples showed greater isoform diversity compared to control and Exp- samples, but not an overall skew in the types of alternative splicing observed.

To further characterise and quantify alternative splicing events within our short-read dataset, we utilised rMATS, a junction-based method that employs a hierarchical framework, to model unpaired datasets and evaluate and quantify all major alternative splice types [23–25]. Comparing Exp+ to control samples (PWC1) we identified 2,118 significant alternative splicing events in 1,437 genes (S1 and S2 Figs and S4 Table). Alternative splicing comparisons including Exp+ (PWC1 and PWC2) yielded the most alternative splice types in every category (S1 and S2 Figs) further reinforcing the spliceopathy aspects of this genetic subtype of disease.

We then compared the rMATS and Iso-seq results to a subset of published datasets exploring differential splicing in FECD samples (S5 Table) and selected 24 published differential splicing events (from 23 genes) that had previously been identified by bioinformatic or in vitro methods to be strongly associated with CTG18.1 expansion-mediated FECD [4,6,7]. Alternative splicing analysis of our transcriptomic dataset identified that 75.0% (18/24) of these events were significant in PWC1 (control vs Exp+, S5 Table). In addition to these 24 selected events, rMATS identified another 10 novel and significant differential splicing events within Exp+ in these genes, highlighting the increasing complexity of splicing dysregulation within Exp+ FECD. Long-read sequencing confirmed the presence of these novel splicing events in Exp+ FECD (S6 Table). Also, it was interesting to note that FN1 was identified to have robust alternative splicing in Exp+, and further investigation showed that alternatively spliced exons in Exp+ corresponded with skipping of the Extra Domain A (EDA) and Extra Domain B (EDB) protein domains that play a critical role in signal transduction and fibrillar formation [26] (S7 Table).

Skipped exon alternative splicing events, where an additional exon is either included or excluded, were the most common in all three pairwise comparisons (S1 Fig). Prior work regarding alternative splicing in FECD has focused on MBNL1/2 sequestration [6,7,27] as MBNL proteins are key regulators of exon splicing [28,29]. Given their complex role, we attempted to determine if the flanking intronic sequences of alternatively spliced exons were enriched for common RNA-binding protein (RBP) binding motifs. We focussed our RBP motif enrichment analysis (rMAPS2 [30]) on the 250bp region upstream of alternatively spliced exons and isolated those areas where there was enrichment in PWC1 and PWC2, but not PWC3 (i.e. events enriched within Exp+ CEC lines). By filtering rMAPS results in this manner, a trend whereby skipped exons are significantly (p < 0.05) upregulated by RNA-binding proteins with TT-containing binding motifs and downregulated by RNA-binding proteins with CG-containing binding motifs was observed (S8 Table). As MBNL proteins favour YCGY binding motifs, this suggests that other RNA-binding proteins with TT-containing motifs could also play a role in driving the spliceopathy underpinning RNA toxicity, a key mechanism of disease in Exp+ FECD.

Evidence of TCF4 dysregulation appears enriched in CTG18.1-mediated FECD CECs, highlighting specific candidate biomarkers for this genetic subtype of disease

Differential gene expression was assessed with DESeq2 [31] in all three groups including sex and ethnicity as covariates. For all comparisons an FDR corrected P-value (padj) of <0.05 was used as a threshold for significance (S9 Table). A total of 4,288 genes were identified to be differentially expressed between the Exp+ and controls (PWC1; Fig 2A). Of these, 2,089 genes had a shrunkLFC value (moderated log₂fold values of expression accounting for high dispersion and low expression counts) of greater than 1 or less than -1. Overall, PWC2 and PWC3 groups resulted in lowered differential gene expression, with 1,695 and 1,911 significantly differentially expressed genes, respectively (Fig 2A and S9 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Differential gene expression analysis on Exp+ FECD, Exp- FECD, and control primary corneal endothelial cell (CEC) transcriptomes.

A) Differential gene expression was assessed via DESeq2. Three pairwise comparisons were conducted: PWC1 (Exp+ vs control), PWC2 (Exp+ vs Exp-), and PWC3 (Exp- vs control). Venn diagrams show overlap of significantly differentially expressed genes (adjusted p-value < 0.05) within these three comparisons. B) Total TCF4 expression showed no significant dysregulation between the different groups. C) Volcano plot shows the distribution of gene expression in Exp+ FECD compared to the control group. Horizontal dotted lines denote the boundary between significance (adjusted p-value < 0.05), whereas vertical dotted lines show a shrunkLFC cutoffs for -1 and 1. Dark grey solid markers show significantly expressed genes with a fold change magnitude greater than 1. Solid orange markers delineate uniquely Exp+ differentially expressed genes. D) Pathways enriched in Exp+ vs Control (g:profiler). Shaded grey boxes show direction of pathway enrichment in TCF4 knockdown models (Doostparast Torshizi et al. 2019) [33]. Green dots show pathways that are significantly enriched in PWC3 as well. Solid orange dots and open orange circle show that pathway enrichment persists in genes unique to Exp+ (1,028 genes and 3,292 genes, respectively).

https://doi.org/10.1371/journal.pgen.1011230.g002

Pairwise analysis revealed considerable overlap between dysregulated genes identified by all three pairwise comparisons (Fig 2A), however, highest convergence was observed by PWC1 (Exp+ versus control) and PWC2 (Exp+ v Exp-), suggesting these events are driven by CTG18.1 expansions. In total, both pairwise comparisons reveal 1,028 dysregulated genes uniquely and significantly dysregulated in Exp+ FECD case-derived samples. Conversely, by comparing the overlap between PWC1 and PWC3 (Exp+ vs control and Exp- vs control) we identified a further 873 genes that are dysregulated in FECD, irrespective of the underlying genetic cause, and thus representing more ubiquitous signatures of the disease. Dysregulated genes in PWC1, including Exp+ unique and common to FECD, did not skew directionally and were evenly up or down-regulated within our analysis (Fig 2C). We isolated 3,292 genes that are uniquely dysregulated in Exp+ FECD, 2,264 genes dysregulated solely in PWC1 and a further 1,028 dysregulated genes shared with PWC2; (highlighted in orange Fig 2A) and proceeded with downstream analysis (S9 Table).

We conducted enrichment analysis using g:profiler [32] (querying gene ontology (GO), KEGG pathway, and Reactome pathways). GO term enrichment and KEGG pathway analysis of PWC1 dysregulated genes was consistent with previously published transcriptomes [6,11,12], with several enriched terms and pathways relevant to RNA toxicity (i.e. GO:0000166 etc) (S10 Table). While total TCF4 gene expression was not significantly downregulated in Exp+ (Fig 2B), intriguingly, we observed strong parallels to pathway enrichment analyses highlighted from a study of TCF4 neural progenitor knockdown cell model [33] within the 3,292 genes uniquely dysregulated in Exp+, such as protein processing in endoplasmic reticulum, MAPK signalling, actin cytoskeletal regulation, endocytosis and apoptosis (Fig 2D). Critically, these pathways are not enriched in PWC3, suggesting that these are unique to Exp+ CECs. Furthermore, several enriched pathways and GO terms involving Wnt signalling were seen in the 3,292 dysregulated genes unique to Exp+ that are not present in PWC3, again suggesting TCF4 function may be modified within Exp+ FECD CECs (S10 Table). Furthermore, enrichment of dysregulated genes in pathways involving other repeat expansion-mediated diseases and misfolded protein aggregate diseases suggest mechanistic parallels which reinforce Exp+ CECs as a model system with broad relevance to repeatome biology (S10 Table).

The consistent presence of PI3K-Akt pathway enrichment within KEGG results in Exp+ dysregulated genes (and its consistent absence in PWC3 pathway enrichment) suggests that Exp+ FECD CECs may interact with growth factors and the extracellular matrix in a unique manner. Hence, components of the PI3K-Akt pathway may serve as informative biomarkers for Exp+ FECD in future clinical trial settings, and thus warrant further investigation. Additional pathways are enriched in this set of genes (and not in PWC3) that involve RNA processing and binding, other repeat-mediated disease related mechanisms, epithelial to mesenchymal transition, and multiple signal transduction cascades with extracellular matrix interactions (S10 Table).

To determine if other RNA-binding proteins may be implicated in the alternative splicing seen in Exp+ FECD, we cross-referenced our short-read data with a list of ~1500 human RNA-binding proteins [34], and discovered that there were 141 RBPs uniquely dysregulated in Exp+ (45 downregulated and 106 upregulated), including MBNL1, CELF1, ESRP2, and a variety of hnRNP and DDX encoding genes (S11 Table). Both the hnRNP and DDX family of RNA-binding proteins have previously been implicated in other repeat expansion diseases [32,35–37]. It is likely that a combination of these dysregulated RBPs plays a role in driving the alternative splicing and RNA toxicity seen in Exp+ FECD.

Differential TCF4 exon usage suggests isoform ratio disruption of TCF4 in CECs of CTG18.1-mediated FECD individuals

Given signatures of TCF4 dysregulation were apparent in the pathway enrichment analysis of Exp+ CECs (Fig 2D), we then aimed to determine if any tangible patterns of TCF4 transcript-specific dysregulation was observed in the Exp+ CECs, despite total TCF4 levels not being significantly different between Exp+, Exp- and control CECs (Fig 2B). Due to transcriptional diversity and complexity of TCF4 transcripts, DEXSeq [38] was used to assess differential exon usage utilising our short-read transcriptomic data.

DEXSeq utilises generalised linear models to detect differential exon usage while controlling for biological variation. This analysis allowed us to query gene expression and splicing at the exon level comparing PWC1 (with and without (PWC1_nocov) sex and ethnicity as covariates), PWC2, PWC3, FECD compared to control (FvC) and Exp+ compared to all expansion negative samples (Exp+ v Ctr+Exp-). In concordance with our rMATs and DESeq2 analysis, higher levels of transcriptome-wide differential exon usage events were observed in PWC1 compared to PWC3 (S3 Fig), suggesting these patterns of differential exon usage are likely attributed to CTG18.1-related pathology.

Our short-read transcriptomic data reveals that multiple TCF4 exons (n = 5) were significantly dysregulated in Exp+ compared to control and Exp- FECD samples (Fig 3A and S12 and S13 Tables), whereas no significant differences were observed between Exp- and control (PWC3, Fig 3B and S12 and S13 Tables), demonstrating for the first time that differential TCF4 exon usage is driven by CTG18.1 expansions and not a downstream feature of FECD unrelated to the primary genetic driver. A subset of exons contained within a small group of TCF4 transcripts (4 out of 93) are upregulated in Exp+, while an increased number of more broadly expressed exons are downregulated (appearing in 56 out of 93 transcripts) (Fig 3A and S12 and S13 Tables). Transcripts containing downregulated exons are primarily included in longer isoforms of the gene which contain 3 AD domains, while the upregulated exons (excluding those surrounding the repeat) appear in shorter TCF4 isoforms containing only two AD domains. This suggests that Exp+ CECs may be employing certain compensatory measures in response to the transcript-specific dysregulation that is induced by the presence of CTG18.1 repeat expansions. Interestingly, our rMATS analysis independently detects this same phenomenon (S14 Table), which showed a significant mutually exclusive exon splicing event whereby the last exon in common with most 3 AD domain-containing isoforms (Exon 6 in ENST000000354452.8) was favoured in the control group, while the first exon in common with all isoforms (Exon 7 in ENST000000354452.8) was favoured in the Exp+ group.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. CTG18.1 expansion-mediated FECD results in alternative TCF4 exon usage within corneal endothelial cells (CECs).

A) Differences in exon usage between control (grey dotted line) and Exp+ (orange dotted line) CECs. Significantly upregulated and downregulated exons are highlighted with blue and red vertical lines, respectively. The location of CTG18.1 is depicted with an orange line. Inset boxes magnify dysregulated regions. A genomic schematic of TCF4 with transcripts aligned is presented below the axis. Key protein domains highlighted (green—basic helix loop helix (bHLH) domain, blue—activation domains (ADs), red—main nuclear localization signal (NLS). B) Differences in exon usage between control (grey dotted line) and Exp- (green dotted line) CECs. No exons were statistically dysregulated.

https://doi.org/10.1371/journal.pgen.1011230.g003

We also detect retention of the intron flanking the CTG18.1 repeat, previously described by Sznajder et al. 2018 [39] in Exp+ FECD (S12 and S13 Tables). The intronic region is significantly upregulated in Exp+ CECs (E173 and E174) (S12 and S13 Tables). While this intron represents the 5’ untranslated region (UTR) of several transcripts of TCF4, it is also included within an intronic region in 32 out of 93 currently annotated TCF4 transcripts.

RNAScope shows robust detection of altered tissue-specific TCF4 isoform expression in CTG18.1-mediated FECD case-derived corneal endothelial cells

We performed spatial transcriptomics, via RNAScope, to further validate our hypothesis that a subset of TCF4 isoforms containing the third AD3 domain are being selectively downregulated due to the presence of an expanded CTG18.1 allele (≥50 repeats). For the first time, we visualised TCF4 transcripts within FECD and control CECs and fibroblasts using two distinct sets of probes, termed A and B (Fig 4A). Probe set A comprises a set of probes that allow the visualisation of all TCF4 transcripts in the cell (probe A covers all exons in the transcript Isoform B+ which contains 3 AD domains, TCF4-201, ENST00000354452.8, NM_001083962.2). Probe set B targets exons in the canonical longer TCF4-201 TCF4 transcript which do not overlap with the canonical shorter TCF4 transcript (Isoform A+/J which contains 2 AD domains, AD1 and AD2, TCF4-204, ENST00000457482.7, NM_001243234.2). Hence, these probe sets enable the detection, visualisation, and quantification of the ratio of longer AD3-containing transcripts, previously identified to be downregulated in Exp+ CECs via our differential exon usage analysis, compared to total TCF4 transcripts (Fig 3A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. RNAScope probes targeting TCF4 transcripts demonstrate that longer transcripts containing activation domain 3 (AD3) account for a smaller proportion of total transcripts specifically in Exp+ FECD corneal endothelial cells (CECs).

A) Probe targeting schematic for TCF4. Probe A targets all exons in NM_001083962.2 and visualises total TCF4 transcripts, while Probe B targets only the 5’ exons, excluding those that are contained within shorter transcripts (NM_001243234.2). Probe B visualises the longer subset of AD3-containing transcripts within the cell. CTG18.1 is highlighted in orange. Green bar denotes TCF4 bHLH region, blue regions highlight TCF4 activation domains, red box shows the bipartite TCF4 NLS signal and orange shows the genomic region containing the CTG18.1 repeat. Key domains highlighted (green—basic helix loop helix (bHLH) domain, blue—activation domains (ADs), red—main nuclear localization signal (NLS)). B) Representative confocal images of Exp+, Exp- and control CECs and Exp+ and control fibroblasts with nuclei (blue, DAPI) and TCF4 puncta (visualised with Cy3 but represented here with white). Scale bars represent 20 μm. C) Plotting the proportion of longer AD3-containing TCF4 transcripts in Exp+ (orange), control (black), and Exp- (green) shows that Exp+ have significantly lower proportion of AD3-containing transcripts compared to both control (p = 0.0145) and Exp- (p = 0.0213) CECs, while there is no significant difference between control and Exp- (p = 0.362). D) Plotting the proportion of longer AD3-containing TCF4 transcripts in Exp+ (orange), and control (black) shows that there is no statistically significant (p = 0.61) difference in TCF4 isoform ratios between the two groups. Open circle denotes mean and error bars show standard deviation. Diamond datapoint shows bi-allelic expanded FECD cells. Significance determined with student T-test.

https://doi.org/10.1371/journal.pgen.1011230.g004

The TCF4 transcripts were visualised in control, Exp+ and Exp- FECD case-derived CECs (n = 3 for each group) (Fig 4B and Tables 1, S1, and S15). Each independent primary CEC culture was also probed with a negative control (targeting bacterial genes) and a (CAG)₇-Cy3 probe, using a standard FISH protocol, to validate the presence or absence of CUG-specific RNA foci derived from expanded copies of CTG18.1 (S15 Table). RNA foci were consistently detected in all CTG18.1 Exp+ CECs, and were absent from control and Exp- CECs. Using CellProfiler (version 4), puncta representing TCF4 transcripts, detected with both probe set A and B, were quantified and puncta/cell values were established for each image and averaged across images with at least 100 nuclei imaged per condition per line. Finally, a ratio was calculated establishing the proportion of longer AD3 containing transcripts compared to total TCF4 transcripts for each biological replicate. This showed that, on average, Exp+ CECs had a significantly smaller proportion of the longer AD3-containing transcripts (16.33% ± 4.53) compared to unaffected control CECs (36% ± 6.9, p = 0.0145) and non-expanded FECD CECs (30.9% ± 5.17, p = 0.0213) (Fig 4C). To explore potential cell-type specificity of the observed skewed TCF4 transcript ratios we also analysed a series of dermal fibroblast lines generated from unrelated Exp+ patients with FECD (n = 6) [7] and control adult dermal fibroblast lines (n = 3) (S16 Table). No significant difference (p = 0.61) was observed between the Exp+ (44.2% ± 5.68) and control fibroblast (41.17% ± 12.1) lines, suggesting the dysregulation of TCF4 is cell-type specific. This result suggests that the presence of the expanded repeat alone consistently skews TCF4 transcript ratios to favour shorter isoforms, specifically in CECs. As this subgroup of isoforms typically contains only two out of three activation domains, a higher prevalence of shorter TCF4 isoforms may affect potential dimerization and TCF4 transcriptional function in Exp+ CECs.

Rare TCF4 variants identified within a genetically refined CTG18.1 expansion-negative FECD cohort suggest isoform-specific TCF4 dysregulation may be a risk factors for FECD in the absence of CTG18.1 expansions

We next wanted to explore if rare and potentially deleterious TCF4 variants could be contributing to disease in FECD cases without a CTG18.1 expansion, given that Exp+ CECs display alternative TCF4 exon usage and skewed ratios of TCF4 isoforms (Figs 3 and 4), with detectable downstream features of TCF4-driven dysregulation (Fig 2D). To test this hypothesis, we interrogated exome data generated from 134 Exp-FECD cases, acquired from our ongoing inherited corneal disease research program, representing approximately 22% of total FECD cases recruited to date. Overall, 8 rare and potentially deleterious (MAF <0.005; CADD score >15) TCF4 variants were identified in 7/134 probands (Table 3). Sanger sequencing confirmed the presence of all 8 variants and exome data from these cases were also interrogated using the same thresholds (MAF <0.005, CADD >15) to determine if any other rare variants in genes previously associated with FECD could also be contributing to disease [4] (S17 Table). The large number of rare and potentially deleterious TCF4 variants detected in this cohort seemed intriguing given the high levels of evolutionary constraint on this essential and ubiquitously expressed transcription factor encoding gene (pLI = 1, missense Z = 4.1 gnomAD). We therefore decided to apply CoCoRV, a rare variant analysis framework to determine if the Exp-FECD cohort was indeed enriched for rare and potentially deleterious TCF4 variants compared to the gnomAD dataset. CoCoRV utilises the publicly available genotype summary counts to prioritise disease-predisposition genes in case cohorts [40]. Under a dominant model, applying the ethnically-stratified (S18 Table) analysis pipeline, and filtering for rare (MAF <0.005) and potentially deleterious (CADD score threshold of >15) variants, an inflation factor estimate (λemp) 1.21 was achieved. At the exome wide level TCF4 was ranked as 34/17,025 genes with a p-value of 0.001 suggesting the Exp-FECD cohort could be enriched for rare and potentially deleterious TCF4 variants. However, after applying CoCoRV’s stringent false discovery rate correction methodology TCF4 no longer remained significant (corrected p-value = 0.249).

Interestingly rare variants identified in Proband A (Variant 1a/1b), B (Variant 2) and D (Variant 4) and E (Variant 5) are all predicted to only affect a small fraction of total TCF4 transcripts (Tables 3 and S19 and Fig 5). For Proband A, two consecutive in cis heterozygous missense and nonsense variants were identified. These occur within the coding region of only one (out of 93) TCF4 transcripts (ENST00000566286.5; c.[57G>T; 58A>T]; p.[(Arg19Ser; p.Lys20*)]) and in the 5’ UTR of a further 5/93 transcripts. The variant identified in Proband B is only encompassed by the same 6 TCF4 transcripts (ENST00000566286.5; c.66G>A, p.(Glu22 =)). Despite being synonymous, Splice AI predicts the variant introduces loss of the splice donor site (SpliceAI Δ score 0.78). Furthermore, SpliceRover also predicts that this could result in the activation of a cryptic splice donor site downstream (from 0.098 to 0.233) of the wildtype donor site, which would introduce a short frameshift insertion followed by a premature termination codon (PTC) c.66_67insGTGCTCGATGAATTTTC, p.(Arg23Valfs*12). The variant identified in Proband D (Variant 4) is within the protein coding region of three (out of 93) TCF4 transcripts and in the 5’ UTR of one additional transcript. Variant 5 is only encompassed by the 5’ UTR of a single TCF4 transcript. Collectively, all of these identified variants have the potential to exert a functional impact on TCF4, but importantly only on a small subset of overall transcripts.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Schematic map of rare and potentially deleterious variants identified in individuals with molecularly unsolved FECD within TCF4, without CTG18.1 expansions.

All 8 rare (minor allele frequency (MAF)<0.005) and potentially pathogenic (CADD>15) variants in 7 Exp- Probands are mapped along TCF4 genomic region and transcriptomic locations contextualised. Proband A variants are found in cis within the proband. Conserved intronic transcripts refer to variants located in intronic regions identified as conserved via Genomic Evolutionary Rate Profiling (GERP) conservation scores on Ensembl. Splice modifying refers to transcripts where the variant identified would induce alternative splicing. Transcript numbers where the variant occurred within introns were not included. UTR: untranslated region. NLS: nuclear localization signal. bHLH: basic helix loop helix.

https://doi.org/10.1371/journal.pgen.1011230.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Summary of rare and potentially deleterious TCF4 variants identified in a FECD CTG18.1 Exp- cohort (n = 134).

https://doi.org/10.1371/journal.pgen.1011230.t003

Ethics statement

The study adhered to the tenets of the Declaration of Helsinki and was approved by the Research Ethics Committees of University College London (UCL) (22/EE/0090), Moorfields Eye Hospital (MEH) London (13/LO/1084), or the General University Hospital (GUH) Prague (2/19 GACR). All participants provided written informed consent.

Subject recruitment, phenotyping and bio-sample processing

All subjects had either clinical signs of FECD or had a previous corneal transplant for FECD. For participants who had planned endothelial keratoplasty, we collected the diseased corneal endothelium and Descemet membrane (DM) removed as part of the surgical procedure (MEH). All individuals recruited with FECD undergoing DMEK surgery had confluent guttae affecting the central cornea with visual symptoms and were of similar age. Importantly, all individuals at the time of bio-sample recruitment had similarly advanced disease (Krachmer Grade 5).

We applied strict criteria in the selection for control tissues. We obtained control tissue as donor (age >45 years) corneoscleral discs classed as unsuitable for clinical transplantation due to non-CEC related reasons (Miracles in Sight Eye Bank, Texas, USA). We selected the age of the donor tissues to ensure there was no statistically significant differences in the ages between control and case at the time of bio-sample recruitment. All control corneas were examined and were shown to have a high CEC counts (≥2,600 cells/mm²) with no guttae or other endothelial abnormality. To standardise handling, the control DM and corneal endothelium were excised by a consultant-grade surgical clinician in accordance with standardised clinical practice. All excised case and control tissues were stored in Lebovitz L-15 media (Life Technologies) before transferring to the laboratory for processing within 24 hours of excision from either donor corneas or patients undergoing elective posterior corneal transplantation surgery.

Participant genomic DNA for CTG18.1 genotyping was extracted from whole blood (QIAgen Gentra Puregene Blood Kit), whereas control DNA was extracted from the residual corneoscleral disc (QIAgen Blood and Tissue kits).

CTG18.1 genotyping

We determined the CTG18.1 repeat length using a short tandem repeat (STR) genotyping assay as previously described in Zarouchlioti et al. 2018 [7], originally adapted from Wieben et al. 2012 [55]. Briefly, STR-PCR was performed using a 6’-FAM-conjugated primer (5’ CAGATGAGTTTGGTGTAAGAT 3’) upstream of the repeat and an unlabeled primer (5’ ACAAGCAGAAAGGGGGCTGCAA 3’) downstream. Post-PCR, products were combined with GeneMarker ROX500 ladder and separated on an ABI 3730 Electrophoresis capillary DNA analyzer (Applied Biosystems). Sizing and data analysis was performed on GeneMarker software (SoftGenetics). To distinguish between samples with presumed bi-allelic non-expanded alleles of equal length and samples with a larger repeat length beyond the detection limit of the STR-PCR, we performed a further triplet-primed PCR assay to confirm the presence or absence of a CTG18.1 allele above the detection range of the STR assay (approximately 125 repeats), following previously published methods [56].

Primary corneal endothelial cell (CEC) culture

We isolated and cultured FECD and control-derived CECs as described by Zarouchlioti et al. 2019 [7], adapted from Peh et al. 2015 [57]. Briefly, to dissolve the DM and detach the CECs, the tissue sample was incubated for 2 to 4 hours in 0.2% collagenase type I. After incubation, the CECs were resuspended in stabilisation media (containing Human Endothelial-SFM (Life technologies) supplemented with 5% foetal bovine serum (FBS), 1% antibiotic/antimycotic, and 0.1% selective ROCK inhibitor Y-27632 (AdooQ BioSciences)), and then seeded into 6-well plates pre-coated with FNC coating mixture (Stratech). After 24 hours, we replaced the culture media with expansion media (Ham’s F-12 Nutrient Mix with GlutaMAX Supplement (Life Technologies)/Medium 199 GlutaMAX Supplement (Life Technologies), 20 μg/mL ascorbic acid, 1% insulin-transferrin-selenium (Life Technologies), 5% FBS, 1% antibiotic/antimycotic, 10 ng/mL bFGF (Life Technologies) and 0.1% selective ROCK inhibitor Y-27632 (AdooQ BioSciences)). All cells were cultured at 37°C and 5% CO₂ with twice-weekly media changes.

Fibroblast cell culture

Fibroblast cell lines from dermal skin biopsies were cultured DMEM/F-12, GlutaMAX (Life Technologies) supplemented with 10% FBS and 1% penicillin/streptomycin. Throughout culture, cells were kept in an incubator at 37°C, 5% CO₂ and medium was refreshed every 48 hr until the cells showed appropriate confluence for experimentation or passage.

RNA extraction

All primary CEC monolayer cultures selected for downstream analysis displayed a consistent morphology characteristic of CECs in vivo. Total RNA was extracted from the primary CEC lines using a NucleoSpin RNA XS kit (Macherey-Nagel), in accordance with the manufacturer’s protocol. The RNA was eluted and stored at -80°C prior to sequencing. To assess RNA integrity, RNA integrity numbers (RIN) were determined for all RNA samples using a Bioanalyzer 2100 (Agilent). All samples analysed as part of this study had a RIN value ≥9.6.

Short-read RNA-seq

Stranded RNA-seq libraries were prepared by BGI Genomics using the TruSeq RNA Library Prep Kit (Illumina) following the manufacturer’s protocol. We sequenced these with a HiSeq 4000 platform (Illumina) with short-read strand-specific library prep with poly A enrichment (Illumina). Paired-end reads were filtered to remove low-quality reads (Q ≤5) or those with >10% of N base calls. Adaptor sequences were also removed during data quality analysis. We checked the quality of the resulting fastq files for the short-read data using FastQC Quantification [58]. Reads were aligned to the human genome (hg38) patch 13 using STAR [59] (v2.7.10a) and Salmon (V1.4.0) [60]. After alignment, hierarchical clustering was performed using PoiClaClu [61]. We also performed principal component analysis (PCA) using PCAtools [62] (v2.2.0) to determine global sample clustering. Sex and ethnicity have been factored in as covariates in all applicable downstream pipelines.

We analysed differential gene expression using DESeq2 [31] (version 1.8.2), with independent hypothesis weighting (IHW) testing of the results to generate FDR values [63,64], and log₂fold change was calculated using the apeglm (v 1.2) package [65]. We then determined significance as FDR corrected-p value (padj) <0.05. In addition, we analysed different exon usage amongst the subgroups with DEXSeq [38] (v 1.44.0). Significant results were filtered using an FDR-corrected P-value of <0.05. Pathway analysis was performed by g:Profiler on significantly dysregulated genes using g:SCS statistics for thresholding [66].

We assessed differential splicing at a junction-based level with rMATS-turbo 4.1.2 [23] with the Ensembl hg38 release 105 GTF file as reference. Default rMATS settings were used. For this study, we focused on splicing events that are fully/partly annotated by Ensembl. Significant alternative splicing events were identified as an absolute delta percent spliced in (psi) magnitude ≥ 0.1 and FDR ≤ 0.05. rMATs output was subsequently analysed using rMAPS2 [30] and maser 1.12.1 [67]. All script utilised to generate data can be found at: https://github.com/InheritedCornealDisease/TCF4-transcriptomic-paper-2024

Long-read RNA-seq

PacBio Sequel long-read transcriptomic sequencing (Iso-Seq) [68] was applied to generate full-length cDNA sequences from control and FECD Exp+ primary CEC-derived RNA samples. Raw reads were processed and de-multiplexed with Iso-Seq (v3) and merged for transcript collapse using Cupcake [68] (parameters: “-c 0.85 -i 0.95—dun-merge-5-shorter”). High-quality, full-length transcripts from the merged dataset were then aligned to the human reference genome (hg38) using Minimap2 [69] (v2.17, parameters: “-ax splice -uf—secondary = no -C5 -O6,24 -B4”) and annotated using SQANTI3 [70] (v7.4). SQANTI filtering of artefacts (intra-priming and reverse transcriptase template switching) was then applied. Github repository for the analysis described here can be found at: https://github.com/SziKayLeung/UCL_FECD

RNAScope detection of TCF4 transcripts

We seeded primary FECD and control CECs and fibroblasts on chamber slides (Lab-tek) and confluent cell cultures were stained using standard RNAScope protocols with minor adjustments (ACD Bio). We used the Hs-TCF4-C2 (cat. 557411) probe set to detect all TCF4 isoforms, as well as a custom set of TCF4 probes to detect all 5’ exons in the TCF4 transcript NM_001083962.2 excluding all 3’ exons contained in the shorter transcript NM_001243234.2. Akoya TSA Cy3 fluorophores were used to visualise probe staining. Negative control probes (cat. 320871) were provided by ACD-Bio to verify specificity. All cell lines analysed by RNAScope were also tested for RNA foci in parallel via fluorescence in situ hybridization (FISH) using a CUG-specific probe (Cy3-(CAG)₇) [7]. A minimum of 100 cells and four images were acquired for each cell line and each probe.

Exome sequencing and rare variant analysis

Exome sequencing libraries were generated from 134 whole blood derived gDNA samples from individuals with Exp- FECD, using a SureSelect Human All Exome V6 capture kit (Agilent, USA) or SeqCap EZ MedExome Enrichment Kit (Roche) and sequenced on a HiSeq 4000 or HiSeq 2500 platform (Illumina). All raw FECD case sequencing data was aligned and annotated following our previously published methods [71]. In brief, sequencing reads were aligned using Novoalign (version 3.02.08) and variants and indels were called according to Genome Analysis Toolkit (GATK) (version 4.4.0.0) best practices (joint variant calling followed by variant quality score recalibration) were normalised and quality-controlled, as recommended by CoCoRV [40]. CoCoRV, a rare variant analysis framework, was applied to prioritise disease-predisposition genes in the Exp-FECD cohort, utilising publicly available genotype summary counts from gnomAD (genomes version 2.1.1) as a control dataset. For consistency, both case and control variants were annotated using the Variant Effect Predictor (VEP) (release 110) [72]. The CoCoRV pipeline was then run, applying a dominant model and ethnically stratified in accordance with ethnicity predictions generated from the case cohort using Fraposa [22]. Fraposa is based on principal component analysis, using 1000 Genomes project as the reference panel (with known ancestral background, categorised into 5 distinctive superpopulations: Africans, admixed Americans, East Asians, Europeans and South Asians). Common variants between the control and exome cohort were identified and filtered (based on minor allele frequency, linkage equilibrium and missing genotype call rate) by plink and the principal components were computed with the default OADP approach (Online Augmentation, Decomposition, and Procrustes Transformation). Variants of interest were defined as having a CADD score (CADD) >15 [73], and a minor allele frequency (MAF) <0.005 in the Genome Aggregation Database (version 2.1.1) [74]. FECD cases identified to harbour rare TCF4 variants were additionally interrogated for rare and potentially disease-associated variants (MAF<0.005 in gnomAD and Kaviar Genomic Variant Database; CADD>15) in previously reported FECD-associated genes, including COL8A2, SLC4A11, AGBL1 and ZEB1 [4]. Splice site predictions were subsequently made using SpliceAI [75] and SpliceRover [76] for all qualifying variants in TCF4 and any other FECD-associated genes.