A novel intronic mutation of PDE6B is a major cause of autosomal recessive retinitis pigmentosa among Caucasus Jews

Patients

A total of 17 Israeli RP patients were recruited for this study, including 8 females and 9 males, with an age range of 20-78 years old. Nine of them were recruited through Hadassah-Hebrew University Medical Center in Jerusalem, six through Rambam Health Care Center in Haifa, one through Hillel Yaffe Medical Center in Hadera, and one through Tel Aviv Sourasky Medical Center. The tenets of the Declaration of Helsinki were followed, the study adhered to the ARVO statement on human subjects, the study was approved by institutional review boards in participating Medical Centers, and written informed consent was obtained from all participants. All patients reported Caucasus Jewish ancestry on both parental sides.

Ophthalmic examination included measurement of best corrected visual acuity (BCVA), fundus examination after pupillary dilatation, fundus photography, spectral domain–optical coherence tomography (SD-OCT; Heidelberg Engineering, Heidelberg, Germany), widefield autofluorescence and pseudocolor fundus imaging (Optos, Dunfermline, UK), full-field electroretinography (ff-ERG), and measurement of flash visual evoked potentials (VEP). ffERG was conducted according to International Society for Clinical Electrophysiology of Vision (ISCEV) guidelines [10], using UTAS Bigshot electrophysiological system (LKC Technologies, Gaithersburg, MD), with bipolar Burian-Allen corneal electrodes (Hansen Ophthalmic Development, Coralville, IA).

DNA analyses

Venous blood samples were obtained using K3EDTA vacuette tubes (Greiner Bio-One, Kremsmunster, Austria) and kept at 4 ºC for up to 48 h. Genomic DNA was extracted from venous blood samples using a high-salt solution according to a standard protocol [11]. Whole exome sequencing (WES) of patient A-1 was performed at Theragen Etex Bio Institute (Suwon, Korea) using SureSelectXT Human All Exon V5 (Agilent Technologies, Santa Clara, CA) and the HiSeq 2500 Sequencing System (Illumina, San Diego, CA). Sequence reads were aligned to the reference human genome (GRCh37/hg19), and variants were called via a customized pipeline as previously described [12]. For confirmation of the c.1921–9C>G mutation and screening of additional patients, specific primers were used to PCR-amplify PDE6B exon 16, including intron–exon boundaries (forward, 5′-GAG AGG CAC AGG CAG CCG AG-′3 and reverse, 5′-CCG TGG CGA TGA TGG CGA TG-′3). PCR was performed on 50 ng of genomic DNA in a 25 µl reaction volume in the presence of 5X Readymix (LAROVA GmbH, Teltow, Germany) and 10 pmol of each forward and reverse primers. Mutation screening was performed by direct sequencing with the Big Dye terminator cycle sequencing kit on an ABI 3130xl Genetic Analyzer (PE Applied Biosystems, Foster City, CA).

Splicing analysis

Evaluation of the identified intronic variant was performed using the web tool Human Splicing Finder [13]. In vitro evaluation was performed using a minigene-based assay. To create wild-type (wt) and mutant minigene constructs, a 750 bp DNA fragment harboring PDE6B exon 15, intron 15, and exon 16 was PCR-amplified from genomic DNA of patients and controls. The fragment was inserted into a mammalian expression vector harboring the Cytomegalovirus (CMV) promoter (pCMV-Script). Constructs were transfected into COS-7 cells, using the jetPEI transfection reagent (Polyplus-transfection, Illkirch, France). Cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (both from Biological Industries, Beit Ha’emek, Israel) and maintained at 37 °C and 5% CO₂. Twenty-four hours following transfection, total RNA was extracted from cells with TRI Reagent (Sigma-Aldrich, St. Louis, MO) and treated with RQ1 RNase-free DNase (Promega, Madison, WI). Reverse transcription was performed with 1 μg of DNase-treated total RNA in a 20 μl reaction volume using 200 U of M-MLV Reverse Transcriptase and 100 ng of random primers (Stratagene, La Jolla, CA). Two microliters of cDNA were subjected to PCR amplification with primers located in exons 15 and 16. PCR was performed in a 25 µl reaction volume in the presence of 5X Readymix (LAROVA GmbH, Teltow, Germany) and 10 pmol of each forward and reverse primers. Primer sequences were detailed in the previous paragraph. PCR products were cloned into the pGEM-Teasy vector (Promega). COS-7 cells used for this experiment were authenticated by genotyping 4 STR markers (D17S1304, D5S1467, D4S2408, D19S245). Genotyping was performed by PCR amplification of each STR and direct sequencing. PCR was performed on 50 ng of genomic DNA in a 25 µl reaction volume in the presence of 5X Readymix (LAROVA GmbH, Teltow, Germany) and 10 pmol of each forward and reverse primers. PCR primer sequences are provided in Almeida et al [14]. Sequencing was performed with the Big Dye terminator cycle sequencing kit on an ABI 3130xl Genetic Analyzer (PE Applied Biosystems, Foster City, CA). Results (in terms of repeat number for each STR) were compatible with those described for COS-7 cells [14] (Appendix 1).

Linkage disequilibrium (LD) analysis

DNA samples were from five unrelated patients homozygous for the c.1921–9C>G mutation and from 45 unrelated healthy Jewish subjects of Eastern European origin (Caucasus and Bukhara), who served as controls. All control subjects reported Caucasus Jewish or Bukhara Jewish ancestry on both parental sides. Individuals of mixed ancestry were not included. Five highly polymorphic (expected heterozygosity: 0.63‒0.79) dinucleotide short tandem repeats (STRP: D4S3360/4PTEL04, D4S3038/AFMa060ze5, D4S43; D4S1614/AFM331xh5, and D4S127) selected through the GeneLoc integrated map of chromosome 4 (v.4.8; Weizmann Institute of Science) were genotyped by PCR amplification with fluorescent dye-labeled primers and capillary electrophoresis (ABI Prism 310; PE Applied Biosystems, Foster City, CA). PCR was performed on 50 ng of genomic DNA in a 15 µl reaction volume in the presence of True Allele PCR Premix (PE Applied Biosystems), and 5 pmol of each forward and reverse primers. Marker haplotypes encompassing the PDE6B locus were inferred either by homozygosity or with the use of the expectation-maximization algorithm as implemented in Haploview (v.4.0; Broad Institute, Cambridge, MA) [15]. LD was calculated according to Bengtsson and Thomson's method [16]: δ = (p₀ ‒ p_N)/(1 ‒ p_N), where p₀ represents the frequency of the allele among individuals with the PDE6B mutation, and p_N represents the frequency of the allele in the general population. Physical map distances were obtained from GeneLoc v.4.8 (based on NCBI GRCh38.p12 built), and genetic distances were calculated through a conversion factor derived from the Marshfield genetic map of human chromosome 4 [17]. Mutation age was first estimated with two moment methods [18,19], both of which are based on the “genetic clock.” Its equation, lnQ = –θg, relates the time (in generations [g]) since the most recent common ancestor (MRCA) of mutant chromosomes, the frequency of recombination between the disease locus and the marker (θ), and the probability that a marker’s allele on a disease chromosome is the ancestral one (Q). An unbiased estimate of Q is the proportion of observed haplotypes that are ancestral. The genetic clock was set following Luria and Delbrück’s approach [20] which takes into account the population growth rate (d). Accordingly, the estimated age was corrected by adding g₀ = –(1/d)ln(θ/d) to the g value [21]. With the moment method introduced by Reich and Goldstein [18], it was possible to correct the decay of LD over generations for the mutation rate (μ = 0.00056 for dinucleotide repeats [22]) at marker loci. LD data were also used to estimate the MRCA age with Rannala and Reeve’s Bayesian Markov chain–Monte Carlo method [23] implemented in the DMLE+ program, v.2.3.

Clinical findings in the studied cohort of patients with RP

A male patient from a non-consanguineous Jewish family who originated from Caucasia was diagnosed with RP at the age of 40 years (patient A-1). He had had night blindness since childhood. At the age of 39 years, his scotopic ERG was non-recordable, and his photopic ERG was significantly reduced. At this age, his BCVA was compromised, and fundus examination revealed pigmentation at the macular area in both eyes (Appendix 2). To further study the genetic basis for RP in Caucasus Jews, we ascertained additional patients with RP of the same ethnic background (Figure 1A).

The diagnosis of RP was primarily based on ERG testing in 6/8 patients (A-1, C-1, C-2, D-1, E-1, and F-1), and was further supported by funduscopic findings and/or visual field testing. In two additional patients (B-1 and E-2), ERG was not performed, and the RP diagnosis was based on fundus findings, visual field constriction, and self-report of night blindness as the initial symptom. All patients had classic symptoms of RP, including night blindness as the initial symptom, concentric restriction of the visual field, and a rod>cone disease process, as indicated by the ERG recordings. Typical funduscopic findings (including narrowed retinal blood vessels and bone spicule-like pigmentation in the periphery) were observed in six patients (A-1, B-1, D-1, E-1, E-2, and F-1). For two patients (C-1 and C-2), funduscopic exam records were not available. Interestingly, significant macular involvement (including macular atrophy, pigmentation at the macular region, and cystic macular edema) was noted in most patients (Appendix 2, Appendix 3, and Figure 2). Age of onset (as reported subjectively by the patients) varied between childhood and 49 years. However, in most patients (6/8), onset was during the fourth decade or later, indicating a relatively late onset of the disease.

Identification of an intronic mutation in PDE6B in a patient with RP with whole exome sequencing

To identify the genetic cause of disease in patient A-1, WES was performed. Because patient A-1 is an isolated case, AR and XL patterns of inheritance were considered. WES resulted in 474,641 variants. Data were initially filtered for rare variants (minor allele frequency ≤1% in the (MAF ≤1% in the 1000 genomes database, the Exome Variant Server and dbSNP) leading to 1,580 missense, nonsense, deletions/insertions, or candidate splicing mutations (based on being located up to 10 bp from an intron–exon junction). We first focused on variants located in known IRD-related genes (RetNet - Retinal Information Network). No putative pathogenic variants (missense, nonsense, deletions/insertions, or canonical (GU-AG) splice-site alterations) in known IRD genes were identified. However, we noted a homozygous C>G transversion, located at position −9 of PDE6B intron 15 (NM_000283: c.1921–9C>G; IVS15–9C>G; Figure 1B). This novel variant is not present in the 1000 Genomes database, the Exome Variant Server, dbSNP, or the Genome Aggregation Database (gnomAD).

c.1921–9C>G mutation affects splicing

In silico analysis predicted that c.1921–9C>G generates a novel acceptor splice site, 9 bp upstream from the original splice site of intron 15. PDE6B is endogenously expressed mainly in the human retina (and specifically in rod photoreceptors), and to a much lesser extent in other tissues, such as brain, heart, spleen, intestine, prostate, ovary, kidney, and muscle. As these human tissues are not easily accessible, we could not evaluate the effect of c.1921–9C>G on splicing in patient-derived RNA. Alternatively, we used an in vitro splicing assay approach. For this purpose, we created two minigene constructs (wt and mutant) harboring PDE6B exons 15 and 16, and the intron between them, downstream of a CMV promoter (Figure 3A). For transfection, we needed a human or primate cell line (which provides the mammalian splicing machinery), which does not express PDE6B endogenously (so we can easily detect transcripts obtained from our minigene construct, without an endogenic background). We chose the COS-7 cell line, which is derived from the monkey kidney. Constructs were transfected into COS-7 cells, followed by RNA extraction and reverse transcription (RT)–PCR with primers located in exons 15 and 16. Both constructs yielded similar products, although the mutant product was slightly bigger (Figure 3B). Both products were subcloned, and eight independent clones of each product were sequenced. Sequencing results revealed that cDNA derived from the wt construct yielded the expected product of 172 bp, in which intron 15 was removed and exons 15 and 16 were correctly spliced. In contrast, cDNA derived from the mutant construct yielded an aberrantly spliced product of 180 bp, in which the newly created acceptor site was used instead of the wt site, leading to an 8-bp insertion into exon 16 (Figure 3A). This insertion is predicted to generate a frameshift and premature termination of translation of the PDE6B protein product (p.T641Lfs*5).

Identification of c.1921–9C>G in additional patients

We used PCR amplification and direct sequencing to test 16 additional patients with RP from 14 unrelated Caucasus Jewish families. c.1921–9C>G was identified homozygously in seven additional patients from five families. Only one of the families was consanguineous (Figure 1A). In total, c.1921–9C>G was identified as the cause of arRP in 40% (6/15) of the Caucasus Jewish families tested.

To test the possibility that the earlier age of onset reported in patients A-1 (childhood) and C-2 (18 years) is due to additional genetic factors, we reviewed the next-generation sequencing results of both patients for additional variants in IRD-related genes. In patient A-1, we observed a heterozygous variant in the ROM1 gene (OMIM 180721; NM_000327, c. 671C>T, p.Pro224Leu). This rare variant (gnomAD AF=0.00004) is predicted to be pathogenic by several prediction tools. Interestingly, ROM1 is involved in digenic inheritance of RP, together with PRPH2 (OMIM 179605; NM_000322), and acts as a modifier of the PRPH2-induced retinal phenotype [24,25]. Whether ROM1 acts as a modifier of PDE6B as well is still to be determined.

Age of the founder mutation

The occurrence of the same biallelic PDE6B c.1921‒9C>G mutation in eight patients with arRP belonging to six independent Caucasus Jewish families is suggestive of a common ancestor who introduced the mutation in the population by a founder effect. We addressed this intriguing question by performing haplotype analysis with the use of five STRPs and the c.1921‒9C>G mutation. The presence of an ancestral six-marker haplotype (5-G-8-4-3-5) and several related haplotypes derived from recombination events in the mutation-bearing chromosomes is compatible with the hypothesis that the pathogenic variant c.1921‒9C>G can be traced to an ancient common founder belonging to the Caucasus Jewish population ancestry. Thereafter, we estimated the age of the MRCA of the PDE6B c.1921‒9C>G chromosomes. The presence of significant LD (δ>0.1) at five highly polymorphic STRP markers (Table 1) allowed the estimation of the coalescence time with two moment methods [18,19], both of which are based on the “genetic clock” [21]. Due to the uncertainty in inferring the historical growth rate of the Caucasus Jewish population, the Luria-Delbrück adjustment of the age [20] was applied twice, under the assumption of the extremes of a reasonable interval of the d value (0.2‒0.4).

The mean (± standard deviation [SD]) overall age estimate for the PDE6B c.1921‒9C>G mutation is 58.4±27.1 g (95% confidence interval [CI]: 39.1–77.8 g) with use of the first algorithm [19] and is 53.9±24.1 g (95% CI: 36.7–71.2 g) with use of the second method [18]. Because simple parametric age estimators based on the genetic-clock equation suffer from the uncertainty about the intra-allelic genealogy, the LD data were reanalyzed from a Bayesian perspective. Applying a mean population growth rate of 0.2‒0.4, a Markov chain–Monte Carlo algorithm [23] provided an estimation of 57.3 g (95% CI: 45.4–76.2 g), an age similar to that obtained by the moment methods. However, the lower limit of the CI is higher than the corresponding confidence limits provided by these methods.