ABCA4 c.859-25A>G, a Frequent Palestinian Founder Mutation Affecting the Intron 7 Branchpoint, Is Associated With Early-Onset Stargardt Disease

Abstract

Purpose: The effect of noncoding variants is often unknown in the absence of functional assays. Here, we characterized an ABCA4 intron 7 variant, c.859-25A>G, identified in Palestinian probands with Stargardt disease (STGD) or cone-rod dystrophy (CRD). We investigated the effect of this variant on the ABCA4 mRNA and retinal phenotype, and its prevalence in Palestine.

Methods: The ABCA4 gene was sequenced completely or partially in 1998 cases with STGD or CRD. The effect of c.859-25A>G on splicing was investigated in silico using SpliceAI and in vitro using splice assays. Homozygosity mapping was performed for 16 affected individuals homozygous for c.859-25A>G. The clinical phenotype was assessed using functional and structural analyses including visual acuity, full-field electroretinography, and multimodal imaging.

Results: The smMIPs-based ABCA4 sequencing revealed c.859-25A>G in 10 Palestinian probands from Hebron and Jerusalem. SpliceAI predicted a significant effect of this putative branchpoint-inactivating variant on the nearby intron 7 splice acceptor site. Splice assays revealed exon 8 skipping and two partial inclusions of intron 7, each having a deleterious effect. Additional genotyping revealed another 46 affected homozygous or compound heterozygous individuals carrying variant c.859-25A>G. Homozygotes shared a genomic segment of 59.6 to 87.9 kb and showed severe retinal defects on ophthalmoscopic evaluation.

Conclusions: The ABCA4 variant c.859-25A>G disrupts a predicted branchpoint, resulting in protein truncation because of different splice defects, and is associated with early-onset STGD1 when present in homozygosity. This variant was found in 25/525 Palestinian inherited retinal dystrophy probands, representing one of the most frequent inherited retinal disease-causing variants in West-Bank Palestine.

Stargardt disease (STGD1) is the most frequent inherited maculopathy with an estimated prevalence of one in 10000.¹ It is caused by biallelic variants in the gene encoding the transmembrane ATP-binding cassette transporter type A4 (ABCA4).² More than 2200 unique variants have been reported for ABCA4 (http://www.lovd.nl/ABCA4),^3,4 which have been linked to a spectrum of autosomal recessive maculopathies, such as classical STGD1, fundus flavimaculatus, cone-rod dystrophy (CRD), and atypical retinitis pigmentosa (RP).^2,5–8 

The observed clinical heterogeneity has led to a genotype-phenotype correlation model in which the severity of retinal dystrophy correlates with the residual activity of the ABCA4 protein because of the effect of the two causal alleles.⁹ In this model, variants can be categorized as deleterious (or severe), moderately severe, mild or hypomorphic.^8,10,11 The identification of two pathogenic alleles is important to provide an accurate clinical prognosis and for the application of possible future gene- or variant-specific therapies.^8,12 

A conspicuously high proportion of reported ABCA4 variants (∼25%) affect RNA splicing,¹³ either resulting in exon skipping or intron inclusion when they are located at exon-intron junctions,^11,14 or resulting in pseudo-exon inclusion in the case of deep-intronic variants.^{12,13,15–18} Notably, previous studies using whole ABCA4 gene (consisting of 50 exons spanning 128 kb) sequencing has uncovered pathogenic non-coding variants in ∼21% of STGD1 probands and in vitro splice assays have allowed a detailed characterization of putative splice variants.^12,17,19,20 

The prediction of the presence of pseudo-exons in ABCA4 was first shown by Braun et al.,²¹ who performed deep mRNA sequencing. In this way, low-abundant transcripts that contained putative pseudo-exons were revealed. Sequence analysis of these pseudo-exons in genetically unexplained STGD1 cases revealed the first causal deep-intronic variants in ABCA4. These low-abundant transcripts, in the absence of splice-site strengthening or creating variants in normal individuals, are often not present. In those cases, the discovery and selection of putative splice variants relies on the use of in silico prediction tools. Splicing defects occur when a variant alters one of the key splicing signals. These include the consensus sequences at the splice acceptor site (SAS) or splice donor site (SDS), the branchpoint site, as well as exonic or intronic cis regulators that act as enhancers or silencers.^22,23 Several in silico tools are available to predict canonical and cryptic splice sites based on these key elements. The approaches utilized by these tools vary from position weight matrix, Markov model, and maximum entropy principle^24–29 to more complex methods based on artificial neural networks and machine learning.^30,31 Nevertheless, the majority of prediction tools are lacking accuracy in predicting the effect of novel deep-intronic variants outside of the known constitutive splice sites. There are several aspects of splicing not considered by prediction tools, such as the concerted role of different splicing signals and splice-site interdependency,³² chromatin organization³³ and tissue-specific splicing events.^16,34,35 Some of these shortcomings can be addressed by using recently developed tools based on deep residual neural networks such as SpliceAI,³⁶ which is self-trained on pre-mRNA sequences to predict the likelihood that a certain position will be used as a splice site, and does not rely on a priori knowledge, allowing the incorporation of splicing determinants that are still not well understood. 

In the context of a large-scale single molecule molecular inversion probes (smMIPs)-based ABCA4 sequencing study of STGD and STGD-like probands, we found, through the deep learning prediction tool SpliceAI,³⁶ a putative pathogenic variant just upstream of exon 8, which potentially could result in a splice defect. Its frequent presence in other Palestinian STGD1 cases prompted us to investigate its effect on ABCA4 splicing, its associated ophthalmological characteristics, its founder characteristics, and its prevalence in the Palestinian population. 

Materials and Methods

Patients and Genotyping Methods

The “discovery” cohort (cohort 1) consisted of 876 probands with a diagnosis of STGD, CRD, or macular dystrophy (MD) based on the clinical phenotype. They were ascertained in 21 centers from all over the world, and include 175 cases that originate from Palestine. For these probands, the entire 128-kb ABCA4 gene was sequenced using smMIPs, as previously described.¹⁸ In brief, series of ∼210 probands were sequenced using the Illumina NextSeq 500 platform, followed by variant calling and annotation through an in-house bioinformatics pipeline. The second cohort in this study (cohort 2) was described previously¹⁷ and originally consisted of 1,054 probands with genetically unsolved STGD, CRD and MD cases. The smMIPs-based ABCA4 sequencing genetically explained 448 probands,¹⁷ whereas 606 probands remained genetically unsolved, of which 38 probands were submitted by the St. John of Jerusalem Eye Hospital Group and Hadassah Medical Center, Jerusalem, Palestine. Whole exome sequencing (WES) data (M. Salameh, A. AlTalbishi, D. Sharon, unpublished data) and MIPs-based sequencing data (M. Salameh, A. AlTalbishi, D. Sharon, F.P.M. Cremers, unpublished data) from another 312 Palestinian inherited retinal dystrophy (IRD) probands (cohort 3) were re-analyzed for the presence of the c.859-25A>G variant. Similarly, 201 Palestinian IRD probands and affected siblings were analyzed through Sanger sequencing of ABCA4 exon 8 and flanking sequences (cohort 4). Three additional individuals from the National Eye Institute, National Institutes of Health in the United States, were identified to carry the variant by assessment of custom capture next-generation sequencing data from commercial clinical genetic testing panels for STGD, MD or retinal dystrophy. For allele frequency (AF) calculations, an overall cohort of 525 Palestinian probands was considered. This final cohort includes the Palestinian probands present in cohorts 1 and 2 (175 and 38, respectively) and all 312 probands in cohort 3. Because cohort 4 is largely composed of affected siblings and relatives, it was not included in the calculation of the AF. A group of 1400 non-IRD Palestinian probands, for which WES data was available, were used as controls. Statistical analysis of the results was performed using a one-tailed Fisher's exact test. 

Variant Selection by in Silico Analysis

Variant selection from smMIPs sequencing data started with quality filtering on the overall coverage of the variants (>20 reads, given an average of 340 reads per nucleotide) and the percentage of variation at the position (variant present in ≥25% of all reads). Further variant prioritization was performed considering the AF, being <0.005 in the dbSNP database and <0.01 in general population databases, such as the Genome Aggregation Database (gnomAD; http://gnomad.broadinstitute.org/). In silico analysis was performed using the deep-learning tool SpliceAI to identify putative variants causing splice defects.³⁶ SpliceAI predictions were considered for each nucleotide position in ABCA4, and the script was set to return the five highest delta scores for acceptor gain, acceptor loss, donor gain, and donor loss, in a window of 5000 nucleotides (nt) up- and downstream from the selected position. Likely pathogenic deep-intronic or noncanonical splice site (NCSS) variants were selected if the delta score was >0.1 (range 0–1). Variants were considered to be part of the NCSS when located between −14 to −3 upstream of an exon and +3 to +6 downstream of an exon, or positioned within the first and last two nucleotides of an exon. Alamut Visual 2.13 software and its splice defect prediction tools (SpliceSiteFinder-like,³⁷ MaxEntScan,²⁶ NNSPLICE³⁰ and GeneSplicer²⁷) were used as a visual aid to identify the position in which SpliceAI delta scores were predicted and to characterize the genomic context around the variant. 

Midigene-Based Splice Assay

The splicing effect of c.859-25A>G was assessed using in vitro splice assays, based on a previously established wild-type midigene (BA7) containing ABCA4 exons 7 to 11.¹⁴ A mutant construct was generated for the variant through site-directed mutagenesis followed by Gateway Cloning. Subsequently, wild-type and mutant constructs were transfected separately in Human Embryonic Kidney (HEK293T, ATCC# CRL-3216) cells. HEK293T cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum, 1% penicillin-streptomycin, and 1% sodium pyruvate at 37°C and 5% CO₂. For transfection, cells were seeded in a six-well plate and transfected at 70% confluency with 600 nanograms (ng) of plasmid using FuGENE HD reagent (Promega, Madison, WI, USA), as specified in the manufacturers protocol. Transfection of the mutant construct was performed in duplicate. After 48 hours, RNA was obtained from the cells using the Nucleospin RNA kit (Machery-Nagel, Düren, Germany), and cDNA was synthesized from 1000 ng of RNA through the iScript cDNA Synthesis kit (Bio-Rad, Hercules, CA, USA). Reverse-transcription–polymerase chain reaction (RT-PCR), agarose gel analysis and Sanger sequencing were performed to assess the nature of splicing defects. RT-PCR conditions were as follows: 94°C for two minutes, followed by 35 cycles of 30 seconds at 94°C, 30 seconds at 58°C, and five minutes at 72°C, with a final extension step of two minutes at 72°C. Details on the primers used for mutagenesis, PCR, and Sanger sequencing can be found in Supplementary Table S1. Ratios between different RNA products were assessed by semi-quantification using Fiji software³⁸ after agarose gel electrophoresis. 

WES and Targeted Sanger Sequencing

WES was performed using 3billion Inc. (Seoul, South Korea) and Variantyx Inc. (Framingham, MA, USA) sequencing services for 250 and 63 additional samples, respectively. Sanger sequencing (Macrogen Europe B.V., Amsterdam, the Netherlands; analysis using Chromas software) was performed to confirm the presence of the c.859-25A>G variant in this cohort and for 200 additional Palestinian IRD probands. PCR was performed using the following conditions: 95°C for two minutes; 40 cycles of 15 seconds at 95°C, 15 seconds at 60°C, and 45 seconds at 72°C; final extension at 72°C for five minutes. Details on the primers used for PCR and Sanger sequencing can be found in Supplementary Table S1. 

Haplotype Analysis for the ABCA4 Locus

Haplotype analysis was performed on WES data for 16 probands using Automap tool (https://automap.iob.ch/) to show homozygous regions in chromosome 1. For haplotype identification within these regions, BAM files (generated using https://usegalaxy.org/) were uploaded into Integrative Genomics Viewer software (IGV v.2.9.4)^39,40 and were manually screened for each genomic location within the homozygous regions defined by Automap. The ABCA4 locus was screened manually. Variant calling format files (vcf) provided by 3billion were used to generate annotated genetic variants spreadsheets using ANNOVAR online tool (https://wannovar.wglab.org/). Within the ABCA4 region, haplotype information for 23 single nucleotide polymorphisms spanning over 125 kb was obtained. 

Clinical Phenotyping

All individuals were examined by an experienced ophthalmologist for clinical analysis. The phenotype data collected included the following: anamnestic information on disease onset, symptoms and progression, age at last examination, best-corrected Early Treatment of Diabetic Retinopathy Study (ETDRS) chart visual acuity converted to logMAR (31/33 probands), clinical ocular exam by slit lamp biomicroscopy, kinetic and/or static perimetry, fundus photos, fundus autofluorescence (FAF) and optical coherence tomography (OCT). When available, age at electroretinography (ERG) and ERG respecting International Society for Clinical Electrophysiology of Vision standards were collected for both probands and affected relatives (20/56).⁴¹ For most patients, information on a specific age of onset was not available; thus age at last examination was reported instead. In most cases first examination took place in early childhood when the first symptoms appeared. 

ERG for patients from the St John of Jerusalem Eye Hospital Group was performed using the Metrovision machine. The Topcon Triton 3D was used for OCT, fundus photos and FAF. The recently recruited patients had FAF and fundus photos taken using the Zeiss Clarus camera (Zeiss, Germany). In the National Eye Institute, ERG was performed using an LKC console and Burian-Allen contact lens electrodes. In addition, Cirrus-HD-OCT scans were performed as well as color fundus photography and fundus autofluorescence imaging. Full-field ERG (ffERG) information from Hadassah Medical Center was obtained using corneal electrodes and a computerized system (UTAS 3,000, LKC, MD). Briefly, in the dark-adapted state, a rod response to a dim blue flash and a mixed cone-rod response to a white flash were acquired. Cone responses to 30-Hz flashes of white light were acquired under a background light of 21 cd/m². All responses were filtered at 0.3–500 Hz and signal averaging was used. OCT was performed using the Heidelberg Spectralis system, while color, infrared and FAF imaging were assessed using a Zeiss and/or Optos fundus camera and the Heidelberg Spectralis system. At Sickkids in Toronto phenotyping included ERG testing using the Diagnosys system whereas OCT imaging used the Heidelberg Spectralis system. 

Results

Identification of ABCA4 Intron 7 Variant c.859-25A>G and In Silico Splice Defect Predictions

Incorporation of SpliceAI predictions in the annotation of smMIPs sequencing data led to the identification of the c.859-25A>G (chr1(GRCh37):g.94,546,299T>C) variant in 10 patients of Palestinian descent in the 876 STGD and STGD-like probands of cohort 1 (Table 1). This novel variant is absent in control populations in gnomAD. The position of the variant outside the SAS consensus sequence, but still close to the intron-exon junction, suggested a putative effect on other splicing elements, such as the branchpoint sequence (BPS). Therefore, an in-depth in silico analysis of the variant was performed using modified settings for the deep learning tool SpliceAI. Using a “scanning” window of 5000 nt upstream and 5000 nt downstream of c.859-25A>G, the strongest prediction entailed the skipping of exon 8 (SAS loss with a delta score of 0.39), together with a new SAS gain 138 nt upstream of exon 8 (delta score 0.31) (Supplementary Table S2). Interestingly, a donor loss with a delta score of 0.09 was predicted at the SDS of intron 7. Although this delta score did not reach the threshold of significance, it might be a contributing factor to exon 8 skipping (Supplementary Fig. S1A). 

Table 1.

 

Geographic Origin and Genotypic Details of Probands Carrying Variant c.859-25A>G

 

Considering that previous smMIPs sequencing datasets were not analyzed with the aid of SpliceAI, we re-examined the sequencing results of 606 genetically unsolved STGD or STGD-like probands (cohort 2)¹⁷ and identified c.859-25A>G in nine additional cases (Table 1). Similarly, analysis of cohorts 3 and 4 lead to the identification of seven and five additional probands carrying this variant, respectively. Additionally, 22 affected relatives were found to carry this variant (Supplementary Table S3). Finally, two additional probands and one affected relative carrying the variant were independently identified in the National Eye Institute in the United States (Cohort 5). In total, we found 56 cases carrying c.859-25G>A, comprising 41 homozygotes and 15 compound heterozygotes, belonging to 33 families. Most of the cases carrying c.859-25A>G reside in or in the vicinity of Hebron with the following distribution: Hebron Governorate (n = 39 cases), Bethlehem (n = 7), Jerusalem (n = 6) (Supplementary Fig. S2). The four cases who carry this variant and reside in the US or Canada originated from Jerusalem and Bethlehem (Table 1 and Supplementary Figs. S3 and S4). 

Allele Frequency Comparison Between Patient and Control Cohort

An overall cohort of 525 Palestinian probands was considered to calculate the AF of c.859-25A>G in the patient population. As 20 homozygous and 5 compound heterozygous cases were identified, we calculated an AF of 0.04. Screening of the WES control cohort (1400 non-IRD probands) returned an AF of 0.003, since the variant was identified in only 9 alleles. The difference between the two AF was assessed for statistical significance with a one-sided Fisher's exact test, that returned a p-value of 0.0001. Additionally, an odds ratio (OR) of 13.89 (confidence interval: 6.76–28.51) was obtained, which is in line with the OR>5 value ACMG guidelines report as a strong indicator of pathogenicity.⁴² These results show that c.859-25A>G is significantly enriched in the patient population. 

In Vitro Splicing Assay

The effect of c.859-25A>G was assessed using in vitro splice assays in HEK293T cells (Fig. 1). The cells were transfected either with a wild-type midigene construct (spanning ABCA4 exon 7 to exon 11) or a mutant construct carrying c.859-25A>G. Sanger sequencing analysis of RT-PCR products revealed three splicing defects and no remaining wild-type product in HEK293T cells. Skipping of exon 8 was observed, which resulted in a shift of the reading frame and the introduction of a premature stop codon (p.Phe287Hisfs*7). In addition, two distinct partial intron inclusions were observed among the RT-PCR products obtained from HEK293T cells transfected with the previously described midigene construct. These are likely due to the abolition of the BPS and lead to the activation of two cryptic SASs in intron 7. A partial intron inclusion of 138 nt at the 5ʹ end of exon 8 corresponds with the activation of a cryptic SAS at position c.859-138, as predicted by SpliceAI (Supplementary Table S2 and Fig. S1A), and results in a frameshift leading to the introduction of a premature stop codon (p.Phe287Leufs*3). A second partial intron inclusion of 685 nt is the result of the use of a second cryptic SAS at position c.859-685, also leading to a frameshift and a premature stop codon (p.Phe287Tyrfs*33). Interestingly, this partial inclusion was predicted by the second strongest delta score, which yielded a SAS gain with a delta score of 0.06 at position c.859-685 (Supplementary Table S2 and Supplementary Fig. S1B). Considering that all three mutated products lead to a nonfunctional ABCA4 protein, c.858-25A>G can be considered to be deleterious and will act as a severe allele. The protein notation of this defect, based on semi-quantification of the cDNA products (Supplementary Table S4), is expected to be p.[Phe287Hisfs*7,Phe287Tyrfs*33,Phe287Leufs*3]. 

Figure 1.

 

Overview of splice defects caused by variant c.859-25A>G in HEK293T cells. (A) Wild-type and mutant midigenes assay results. Rhodopsin exon 5 (RHO ex5) RT-PCR was used as a control for transfection efficiency. The housekeeping gene Beta-Actin (ACTB) was selected as a control for transcription. To the right, schematic representation of WT midigene (BA7_WT), in which the position of the variant is indicated with an arrow. Beneath, schematic representation of the four RT-PCR products identified in panel (A). c.859-25A>G leads to exon 8 skipping (Fragment 4), a partial intron 7 inclusion of 138-nt 5ʹ (Fragment 2) and a partial intron 7 inclusion of 685-nt 5ʹ elongation of exon 8 (Fragment 1). WT product (Fragment 3) was not detected in the mutant construct. (B) The chromatograms show the exact exonic and intronic breakpoints in the four fragments as confirmed by Sanger sequencing.

c.859-25A>G Resides on a Founder Haplotype

Homozygosity mapping was performed using WES data from 16 c.859-25G>A homozygous probands. Homozygous regions in probands carrying c.859-25G>A varied in size from 85 kb to 108 Mb. The identity-by-descent (IBD) region starts at g.94,546,299 (c.859-25G>A) and ends at position g.94,486,667 (rs933073), resulting in a shared homozygous region of 59.6 kb (Fig. 2 and Supplementary Table S5). Because WES data does not cover the ABCA4 locus completely, the exact boundaries of the IBD region cannot be determined. Intervals of 18.2 kb and 10.1 kb lacking SNP data to assess the haplotype are present at the 5ʹ (telomeric) and 3ʹ (centromeric) side, respectively. Thus the IBD region spans between 59.6 and 87.9 kb. In the 5ʹ region of the ABCA4 gene, four haplotypes can be distinguished (Supplementary Table S5), suggesting that at least four different recombination events occurred in the ancestors of the STGD1 cases between variant c.859-25G>A at position g.94,546,299 and the most 5’ (centromeric) ABCA4 SNP analyzed (rs3789451; g.94,586,328). Of these, haplotype 1 was the most prevalent in the analyzed cases (8/16), and proband 1 carrying this haplotype was thus considered as a reference to determine shared homozygosity regions between probands (Fig. 2 and Supplementary Table S5). At the 3ʹ end of the ABCA4 gene, all probands except proband 3 are homozygous and carry the same haplotype. 

Figure 2.

 

Founder haplotype analysis for c.859-25A>G in 16 homozygous probands. On the left, the genomic positions on chromosome 1 and at the top patient ID numbers are reported. The haplotype of patient 1 was taken as a reference for comparison. Black boxes represent identical homozygous segments between individuals, while grey boxes represent regions of homozygosity in single individuals. IBD: identity-by-descent, encompassing between 59.6 kb and 87.9 kb. The homozygous region in 37 ends at 12.8 Mb.

Ophthalmic Characteristics of Biallelic c.859-25G>A Cases

Phenotypic details for a representative selection of homozygous and compound heterozygous patients harboring c.859-25A>G are provided in Table 2 and Figures 3 and 4. The mean visual acuity of homozygous c.859-25G>A patients was 1.1 logMAR, indicating legal blindness (at mean age of 12 years). Examples of FAF, OCT, and fundus images are depicted for cases aged 10 years, 14 years and 42 years (Figs. 3A–C), which show different stages of STGD1. All cases showed foveal photoreceptor atrophy. None exhibit the yellow flecks that are typically seen in intermediate STGD1 cases. The 10-year-old presents with bull's eye maculopathy and the 14-year-old shows severe macular retinal pigment epithelium (RPE) and photoreceptor (PR) atrophy with localized hyperpigmentation deposits. The 42-year-old displays diffuse retinal atrophy and hypo-autofluorescence involving the whole retinal posterior pole and mid periphery, as observed in fundus images and FAF. OCT shows advanced atrophy of RPE and outer retinal layers of the macula. ERG information was available for 15 homozygous individuals. All of these showed abnormal cone responses and 13/15 showed abnormal rod responses suggesting a diagnosis of CRD (Supplementary Table S6). 

Table 2.

 

Clinical Characteristics of Selected ABCA4-Retinopathy Patients With c.859-25A>G

 

Figure 3.

 

Ophthalmic features of homozygous retinopathy cases carrying c.859-25A>G. Fundus autofluorescence (upper panels), optical coherence tomography (OCT) (middle panels) and color fundus (lower panels) for left (OS) and right (OD) eyes of three homozygous cases. The pictures are ordered by age to show different stages of STGD1 phenotypic characteristics. (A) Proband 18, 10 years of age. FAF images show enlarged foveal hypo-autofluorescence with hypo-autofluorescent spots in the posterior pole and mid periphery. OCT images show central macular atrophy of the RPE, photoreceptors and the outer retinal layers with preserved parafoveal retina, fundus images show bull's eye maculopathy (beaten-bronze pattern) with no clear flecks and few pigmentary clumps in the temporal retina. (B) Proband 27, 14 years of age. FAF images show patchy and pinpoint hypo-autofluorescence in the posterior pole and mid periphery, OCT images show diffused RPE and photoreceptor atrophy with outer retinal layers disruption. Fundus images show beaten-bronze foveal appearance with orange hue and pigment clumps. (C) Proband 24, 42 years of age. FAF images show a geographic pattern of diffuse hypo-autofluorescence involving the whole retinal posterior pole and mid-periphery. OCT show advanced atrophy of the RPE and outer retinal layers of the whole macula. Fundus images show diffuse retinal atrophy in the posterior pole with pigment clumps and orange macular hue.

Figure 4.

 

Ophthalmic features of compound heterozygous retinopathy cases carrying c.859-25A>G. Fundus autofluorescence (upper panels), OCT (middle panels), and color fundus (lower panels) for left (OS) and right (OD) eyes of five compound heterozygous cases. (A) Proband 6, carrying c.6816+2T>A p.(?) as the second allele. FAF images show enlarged foveal hypo-autofluorescence, OCTs show central macular atrophy of RPE, photoreceptors and the outer retinal layers with preserved parafoveal retina. Fundus images show beaten-bronze foveal atrophy pattern with no flecks. (B) Proband 20, carrying c.5882G>C p.(Gly1961Ala) as the second allele. FAF images show enlarged foveal hypo-autofluorescence, OCTs show central macular atrophy of RPE and photoreceptors with preserved parafoveal retina. Fundus images show beaten-bronze foveal atrophy pattern with no flecks. (C) Proband 32, carrying c.5882G>A p.(Gly1961Glu) as the second allele. FAF images abnormal foveal hypo-autofluorescence signal with hypo-hyper autofluorescence flecks in the posterior pole of the retina, OCTs show images with a narrow area of foveal RPE and photoreceptor atrophy. Fundus images show blunt foveal reflex with posterior pole yellowish flecks. (D–E) Probands 55 and 25, both carrying c.4979C>T p.(Pro1660Leu) as the second allele. (D) FAF images show patchy hypo-autofluorescence islands in the posterior pole with diffused pinpoint hypo-autofluorescence involving the posterior pole and mid-periphery, OCT images show diffused RPE and outer retinal layers atrophy involving the whole macula and fundus images show diffused atrophy of the posterior pole with pigment clumps and yellow hue of the macula. (E) FAF images show large islands of complete hypo-autofluorescence with remaining small island of hyper-autofluorescence, OCT images show diffused atrophy of the outer retina, RPE and choroid with posterior staphyloma appearance. Fundus images show diffused posterior pole chorioretinal atrophy with pigment clumps.

The compound heterozygous case carrying c.859-25A>G and the splice site variant c.6816+2T>A (no. 6; Fig. 4A) shows increased foveal hypo-autofluorescence in FAF images. Additionally, OCT shows central macular atrophy of RPE, photoreceptors and the outer retinal layers with preserved parafoveal retina. Fundus images show a beaten-bronze foveal atrophy pattern with no flecks. The variant c.6816+2T>A disrupts the exon 49 SDS and, in the absence of the activation of a novel or wild-type cryptic SDS in exon 49 or intron 49, is predicted to result in exon 49 skipping. The predicted mutant protein (p.(Val2244_Gln2272del)) would lack the C-terminus of the protein. 

Compound heterozygous cases carrying c.5882G>C;p.(Gly1961Ala) as a second allele have a mean visual acuity of 0.9 logMAR. Fundus photos and FAF in most cases show diffuse patchy retinal atrophy of the posterior pole with pigmentary changes and diffuse macular atrophy and retinal thinning in OCT (no. 20; Fig. 4B). Based on our data, we classified variant c.5882G>C as moderately severe. Patient 32 (Fig. 4C) carries c.5882G>A;p.(Gly1961Glu) as a second variant. Fundus and FAF images show an abnormal foveal hypo-autofluorescence signal with some hypo-autofluorescent and hyper-autofluorescent flecks in the posterior pole of the retina and blunt foveal reflex with posterior pole yellowish flecks. OCT images show a narrow area of foveal RPE and photoreceptor atrophy. The retinal images show defects that are confined to the vascular arcade, which is typical for STGD1 cases carrying p.(Gly1961Glu) in trans with any deleterious allele.⁴³ 

Six patients from three separate families (Table 1 and Supplementary Fig. S4) carry c.4979C>T, p.(Pro1660Leu) as the second allele. This is a frequent pathogenic ABCA4 variant in the Israeli population.⁴⁴ The two unrelated probands carrying p.(Pro1660Leu) in Figure 4D and 4E show a severe CRD phenotype with visual acuities ranging from 1.3 logMAR to no light perception, with an average of 1.5 logMAR. In patient 55 (Fig. 4D), living in the United States, color fundus pictures at 31 years have a beaten-bronze macular atrophic appearance with a few flecks extending beyond the nasal aspect of the optic nerve head. FAF imaging shows a combination of macular hyper-autofluorescence (corresponding to flecks) and some hypo-autofluorescent lesions (RPE atrophy), extending beyond the nasal edge of the optic nerve head. OCT confirms the macular atrophy with the affected area showing loss of ellipsoid zone and retinal pigment epithelium. ERG for the patient showed amplitude reductions in both scotopic and photopic responses, measuring at around 50% of the lower range amplitude value compared to normal, with asymmetric responses (worse in the right eye) and notable delay in photopic response implicit times. The abnormal cone and rod responses are suggestive of a CRD diagnosis. In patient 25 (Fig. 4E) FAF and OCT at 51 years show large islands of complete hypo-autofluorescence with remaining small island of hyper-autofluorescence, an diffuse atrophy of the outer retina, RPE and choroid with posterior staphyloma appearance. Fundus images show diffuse posterior pole chorioretinal atrophy with pigment clumps, however no ERG data was available for this patient. The suggested diagnosis for the patient is CRD. In the third family carrying c.4979C>T as second allele, ERG results were available for patient 14 who had a clinical diagnosis of STGD and presented with normal ERG values for both cones and rods. 

Discussion

In this study, we identified a novel Palestinian founder mutation, c.859-25A>G, in ABCA4. The variant affects the putative branchpoint of intron 7, leading to a complex in vitro splice defect in which exon 8 is skipped and two partial inclusions of intron 7 occur. All splice defects are likely to induce nonsense-mediated RNA decay and result in the absence of ABCA4 activity. Additionally, no other pathogenic variants were identified in the probands in the shared homozygous region, offering c.859-25A>G as a compelling pathogenic variant. Patients homozygous for the variant show early-onset STGD1 and a phenotype that strongly points to a deleterious effect of the variant. This is further supported by the genotype-phenotype correlation observed in compound heterozygous probands carrying p.(Gly1961Glu) as second allele, which is in line with previously reported cases of p.(Gly1961Glu) being present in trans with a deleterious allele.⁴³ 

Genotype-phenotype correlations in our study cohort also suggests the classification of variant c.5882G>C;p.(Gly1961Ala) as moderately severe. Regarding c.4979C>T, different phenotypes were observed in the three families carrying the variant in trans with c.859-25A>G. In particular, two families showed a more severe phenotype with a proposed diagnosis of CRD, suggesting a classification severe, whereas in the third family, all affected individuals presented with a milder phenotype associated to classical STGD, pointing to a moderately severe effect. No other pathogenic variants in ABCA4 that would explain the different phenotypes were identified. It is possible that currently unknown modifiers, either in ABCA4 itself or in other genes, play a role in the difference of severity observed. From the data collected in this study, we propose a classification of moderately severe or severe. Finally, phenotypic assessment of the patient carrying c.6816+2T>A suggests a mild or moderately severe effect of the variant. To better classify the effect, further analysis is required to assess the percentage of correct residual ABCA4 transcript. 

ABCA4 founder mutations have been extensively described in literature^45–48 and are known to widely differ between ethnic groups and geographic regions. The Palestinian population has a particularly high incidence of IRD-associated founder mutations, due to the presence of geographically and culturally isolated settlements that originate from small founder populations. Consequently, many of these villages have high consanguinity rates, which results in an enrichment of haplotype homozygosity and a high number of founder mutations.^44,49 Statistical analysis of the AFs calculated for the patient and control cohorts showed a significantly higher frequency of the variant in probands compared to non-IRD individuals, supporting the pathogenicity of the founder variant c.859-25A>G. Considering that we found 27 homozygous c.859-25A>G STGD1 cases in the Hebron Governorate that has approximately 1,000,000 inhabitants, and that not all IRD cases from this area were investigated for the presence of this variant, we can deduce a minimal AF of 0.0052 for c.859-25A>G, and a minimal carrier rate of 0.0103 for this variant in healthy individuals (i.e., 1/97 individuals). This is likely an overestimate because of the high consanguinity levels observed in the population. Nevertheless, knowledge of the carriership of c.859-25A>G is relevant for genetic counseling in Palestinian families living in Hebron. 

It is of interest that despite its high frequency the variant was not previously detected, which is likely due to the position of the variant in a BPS. Although the recognition of the BPS by the splicing machinery is critical in RNA splicing, because it leads to the formation of the lariat structure that precedes intron excision from pre-mRNA, BPS variants are scarce in literature compared to other variants affecting splicing. Nevertheless, several examples of variants abolishing the BP have been reported in hereditary diseases, such as LCAT-associated fish-eye disease,⁵⁰ Ehlers-Danlos syndrome due to COL5A1 variants,⁵¹ hypophosphatasia due to variants in ALPL⁵² and in hereditary cancer-associated genes, such as XPC mutated in patients with xeroderma pigmentosum,⁵³ and RB1 mutated in patients with retinoblastoma.⁵⁴ Earlier this year, the first BPS variant associated with non-syndromic retinitis pigmentosa was described for BBS1.⁵⁵ The rarity of published pathogenic BPS variants is likely due to the difficulties in identifying the BPS sites. Correct identification is a challenge for two main reasons: the conserved motif and the localization of BPS. Firstly, the consensus motif has been identified as the extremely short and degenerate sequence “yUnAy,”⁵⁶ where y = a pyrimidine and n = any nucleotide. In this sequence, only the A residue (position 0, corresponding to the point of formation of the lariat structure) and the U residue (at position −2) are highly conserved. Mutations of these two residues cause BPS abolishment, leading to splice defects, such as exon skipping or pseudo-exon inclusion.^50,57 Second, the majority of BPS have been identified in a window of 18 to 44 nt upstream of the SAS, although BPSs as far away as 400 nt are known.⁵⁸ 

The sparse number of experimentally proven BPSs, both wild-type and mutated, has resulted in additional challenges in the development of effective tools to predict the impact of variants upstream of SASs. The currently available prediction tools perform poorly, and are mostly focused on the recognition of wild-type BPSs.⁵⁹ In this study we show that SpliceAI, given the right parameters, was able to correctly predict the effect of the BPS variants on splicing. These results suggest that SpliceAI is one of the best tools for putative splice affecting variant prioritization, because it is able to flag variants independently of the biological cause of the expected splicing disruption. Our results are in line with the recently reported performance assessment of SpliceAI for intronic ABCA4 splicing variants.⁶⁰ Further confirmation and implementation of visualization tools to check the genetic context around the predicted effects is still advisable. Nevertheless, use of SpliceAI allows one to bypass more lengthy and spurious research of variants, which would otherwise entail the use of multiple tools focused on predicting single aspects of splicing, such as effects on NCSS, BPS, enhancer and silencer motifs separately. As for analysis parameters, the window of flanking sequences considered by the program has a significant impact on the predicted splicing defects of a variant, thus should be carefully taken into account. Both partial intron inclusions observed in the in vitro splice assay were predicted by SpliceAI, but a small window of analysis (e.g., ±100 nt) would not have identified either of these inclusions as they are located 113 and 660 nt upstream of the variant. Notably, the considerable distance of the cryptic SASs also prevented their detection in all other programs used in this study. In particular, analysis of the effect of c.859-25A>G using Alamut Visual prediction tools resulted in a reduction of 1.9 in the predicted intron 7 SAS strength (wild-type, 3.7; mutant 1.8 [range 0–21]) on GeneSplicer, which on its own would not be considered sufficient to select the variant for further testing. 

In this study, the effect of the variant was assessed at the RNA level using an in vitro splice assay using HEK293T cells. There is ample evidence that this in vitro model is able to correctly recapitulate splice defects affecting the consensus splice site sequences at the junctions of exons and introns.⁶¹ This is also true for almost all variants that create new splice sites in introns or strengthen cryptic splice sites in introns, resulting in the recognition by the splicing machinery and pseudo-exon inclusion, exon elongation or intron retention.^12,13,17,62 Nevertheless, considering the absence of retina-specific factors and the artificial nature of the midigene system, it remains interesting to further validate the results obtained in HEK293T cells in an alternative model. The use of induced pluripotent stem cells to obtain patient-derived photoreceptor precursor cells or retinal organoids, allows observation of the variant effect in a context that is more similar to that of the patient. Retina-specific pseudo-exon recognition was found for intronic variants in ABCA4, creating new or strengthening existing exonic splice enhancer motifs.¹⁶ The most frequent variant associated with Leber congenital amaurosis, an intron 26 variant in CEP290, c.2991+1655A>G, was shown to enhance a pseudo-exon RNA insertion in retinal organoids compared to fibroblasts or EBV-immortalized B lymphocytes.^34,35 

Intronic variants resulting in aberrant splicing are good targets for antisense oligonucleotide (AON)-based therapeutic strategies. The use of AONs to modulate splicing and obtain correct transcripts has been proven to be effective in ABCA4.^{13,16,63–65} Regarding the treatment of the effects of c.859-25A>G, the use of AONs might prove challenging, since binding to the region upstream of the canonical SAS comes with the risk of disrupting regulatory motifs and the binding of auxiliary splice proteins. Considering the involvement of two cryptic splice sites upstream of the BPS variant, a potential rescue approach would be to target the new SASs of the partial intron retention with AONs, which are unlikely to affect exon 8 recognition by the splicing machinery. On the other hand, it remains to be seen whether blocking these new SASs would lead to a correct transcript. We hypothesize, partially based on SpliceAI predictions, that this approach may lead to more exon 8 skipping. Experimental testing of potential AONs is needed to investigate this further. 

Outside of AON-based rescue, other therapeutic approaches that could be implemented are gene augmentation or gene editing. Gene augmentation therapy is successful for some retinal disorders (e.g. RPE65-associated Leber congenital amaurosis or early onset RP)⁶⁶ but remains challenging for STGD1. The relatively large size of the ABCA4 transcript prevents delivery with a single adeno-associated virus vector, which has a maximum packaging capacity of about 5 kb. Multiple approaches are being developed to bypass this issue (e.g. the use of a dual vector system)^67–70 and the implementation of nanoparticles, such as non-viral gene therapy delivery system.^67–70 An attractive alternative to gene replacement for ABCA4 is RNA-editing, which allows the correction of point mutations at the transcript level with low risk of permanent off-target effects. Unfortunately, currently only A>G and C>U corrections can be obtained using engineered deaminase enzymes (such as ADAR and APOBEC), neither of which are viable options for our A>G variant.^71–73 

In conclusion, we identified a severe and frequent Palestinian founder variant, c.859-25A>G, leading to a complex splicing defect very likely due to the disruption of the BPS in intron 7 of ABCA4. The high frequency of the variant in the Palestinian population makes it highly relevant in terms of genetic counselling, in particular in the Hebron Governorate and in Jerusalem. Our study highlights the relevance of the use of correct prediction tools to identify elusive splice-altering variants, such as BPS variants, that have the potential to strongly impact genetic diagnosis and genetic counselling in patients affected by STGD1. 

Acknowledgments

The authors thank Moien Kanaan for allele frequency data from 1,400 Palestinian persons who do not show inherited retinal dystrophies. We thank Bin Guan, Ehsan Ullah, and Chelsea Bender for technical assistance, and Stéphanie S. Cornelis for help with Hardy-Weinberg calculations. 

Supported by the RetinaUK, grant no. GR591 (to FPMC), a Horizon 2020, Marie Sklodowska-Curie Innovative Training Network entitled European Training Network to Diagnose, Understand and Treat Stargardt Disease, a Frequent Inherited Blinding Disorder-StarT (813490) (to FPMC), the Foundation Fighting Blindness USA, grant no. PPA-0517-0717-RAD (to FPMC and SR), the Rotterdamse Stichting Blindenbelangen, the Stichting Blindenhulp, and the Stichting tot Verbetering van het Lot der Blinden (to FPMC), and by the Landelijke Stichting voor Blinden en Slechtzienden, the Macula Degeneratie fonds and the Stichting Blinden-Penning that contributed through Uitzicht 2016-12 (to FPMC). This work was also supported by the ProRetina foundation Germany, the Stichting Blindenhulp, the Stichting ter Verbetering van het Lot der Blinden, the Gelderse Blindenstichting, the Stichting voor Ooglijders, and the Oogfonds that contributed through Uitzicht 2020-17 (FPMC and SR). This work was also supported by EJPRD19-234 (Solve-RET) (FPMC and SR). This work was also supported by “Groupement de Coopération Sanitaire Interrégional G4 qui réunit les Centres Hospitaliers Universitaires Amiens, Caen, Lille et Rouen (GCS G4)” and by the Fondation Stargardt France (to CMD), the Foundation Fighting Blindness, Grant/Award Number: BR‐GE‐0214–0639‐TECH and BRGE‐0518–0734‐TECH TECH (to DS and EB), the Israel Science Foundation (grant No. 1778/20) within the Israel Precision Medicine Partnership program (DS and EB), the STEP-GTP fellowship (to MS), and EuropeAid (Peace of Sight project; DS, EB, and AAT). BB, WMZ and RBH were supported from the intramural research fund of the National Eye Institute, National Institutes of Health (NIH). The Henry Brent Chair in Innovative Pediatric Ophthalmology Research (EH). 

Disclosure: Z. Corradi, None; M. Salameh, None; M. Khan, None; E. Héon, Novartis (C), Jansen (C), Atsena Therapeutics (E); K. Mishra, None; R.J. Hitti-Malin, None; Y. AlSwaiti, None; A. Aslanian, None; E. Banin, None; B.P. Brooks, None; W.M. Zein, None; R.B. Hufnagel, None; S. Roosing, None; C‐M. Dhaenens, None; D. Sharon, None; F.P.M. Cremers, None; A. AlTalbishi, None