Targeted next generation sequencing reveals genetic defects underlying inherited retinal disease in Iranian families.

Purpose: Inherited retinal diseases (IRDs) are clinically and genetically heterogeneous showing progressive retinal cell death which results in vision loss. IRDs include a wide spectrum of disorders, such as retinitis pigmentosa (RP), Leber congenital amaurosis (LCA), cone-rod dystrophy (CRD), and Stargardt disease (STGD1).
Methods: In this study, we performed targeted next-generation sequencing based on molecular inversion probes (MIPs) that allowed the sequence analysis of 108 IRD-associated genes in 50 Iranian IRD probands.
Results: The sequencing and variant filtering led to the identification of putative pathogenic variants in 36 out of 50 (72%) probands. Among 36 unique variants, we identified 20 novel variants in 15 genes. Four out of 36 probands carry compound heterozygous variants, and 32 probands carry homozygous variants. Conclusions: Employing a cost-effective targeted next-generation sequencing procedure, we identified the genetic causes of different retinal disorders in the majority of Iranian families in this study.

Introduction

Inherited retinal dystrophies (IRDs) are a group of clinically and genetically heterogeneous disorders leading to vision loss because of degeneration of retinal cells [1-5]. Mutations in more than 250 genes have been reported to cause IRDs. Most identified mutations are inherited in an autosomal recessive fashion (50–60%) although some can be autosomal dominant (30–40%), X-linked (2–20%), or mitochondrial [6]. IRDs are divided into non-syndromic forms, such as retinitis pigmentosa (RP), cone-rod dystrophy (CRD), Leber congenital amaurosis (LCA), and Stargardt disease (STGD1), as well as syndromic forms [7]. These inherited retinal dystrophies have an estimated incidence of 1:2,000 individuals [4].

The earliest clinical symptoms of RP (OMIM 268000) are night blindness and peripheral vision loss (tunnel vision) due to rod cell degeneration. In later stages of the disease, central vision can also be lost once cone cells degenerate, and eventually, patients with RP can become (legally) blind [6,8]. In CRD (OMIM 120970), the main manifestations are loss of visual acuity, color vision, and photophobia followed by night blindness and tunnel vision [9]. The rarest and most severe retinal dystrophy form is LCA (OMIM 204000), in which patients have severe visual impairment in the first year of life [10]. LCA appears with severe and early visual loss, nystagmus, and sluggish pupils [11]. STGD1 (OMIM 248200) is one of the most common macular degenerative disorders with the features of macular atrophy, impaired central vision, and decreased visual acuity. Color vision, especially a red-green defect, is apparent in most affected individuals [12].

Because of the complexity of IRDs, it is often necessary to reconcile the clinical diagnosis, family diagnosis, and molecular diagnosis. This can be done in a multidisciplinary team which ideally includes trained ophthalmologists, genetic counselors, and molecular geneticists. Thus, identification of the pathogenic DNA variants can help clinicians and health professionals to improve and expedite the diagnosis of the disease and provide patients and their families with more accurate genetic counseling. Moreover, increasing knowledge of the molecular mechanisms involved in IRD pathogenesis may also lead to the development of new treatment options for affected individuals in the future [13].

Since the introduction of next-generation sequencing (NGS), several different targeted gene panel tests that sequence a large set of disease-associated genes have been developed for specific inherited eye disorders [14]. In this study, we applied a targeted gene panel–based NGS approach based on molecular inversion probes (MIPs) to obtain a molecular diagnosis for 50 patients who had been diagnosed with IRDs and for whom no genetic tests had been performed previously. Using this method, variants in 108 genes associated with non-syndromic IRDs were screened.

Methods

Subjects

In this study, we included 50 IRD probands with Iranian origin and the clinical diagnosis of RP (n=31), CRD (n=5), LCA (n=11), or STGD1 (n=3) who had not been genetically investigated previously. The clinical diagnosis was established by a trained ophthalmologist including the patient pedigree and history of the disease, electroretinography (MonPack3; Metrovision, Perenchies, France), evaluation of visual acuity using slit-lamp biomicroscopy, and fundoscopy (TRC-50EX retinal camera; Topcon, Tokyo, Japan). Written informed consent was obtained from all study subjects. This study was approved by the institutional review board of the ethics committee of Isfahan Medical University, Iran. The study adhered to the tenets of the Declaration of Helsinki and the ARVO statement on human subjects was approved by the institutional review board of the ethics committee of Isfahan Medical University, Iran.

Targeted sequencing and genetic analysis

Five cc of peripheral blood samples of affected and unaffected members were collected and mixed with EDTA anticoagulant (Merck KGaA, Darmstadt, Germany). Genomic DNA of all samples was extracted using the phenol-chloroform method [15].

Targeted NGS was performed to identify the genetic etiology of IRDs in 50 patients with different types of IRDs. The 108 gene panel MIP analysis designed at the Department of Human Genetics at Radboud University Medical Center (Nijmegen, The Netherlands) consisted of genes implicated in non-syndromic achromatopsia, cone dystrophy (CD) and CRD (autosomal dominant, autosomal recessive, and X-linked), congenital stationary night blindness (autosomal recessive and X-linked), macular degeneration (autosomal recessive and autosomal dominant), LCA (autosomal dominant and autosomal recessive), RP (autosomal recessive, autosomal dominant, and X-linked), and choroideremia. MIPs with 5-bp molecular tags were used to perform targeted NGS of 108 genes associated with IRD (Appendix 1). The 1,524 coding exons and at least 10 bp flanking each exon were targeted using 6,129 MIPs for an overall target size of 647,574 bp. On average, four to six MIPs cover one exon. When known single nucleotide polymorphisms (SNPs) resided in the genomic regions complementary to the annealing or ligation arm, SNP-specific MIPs were designed. The panel also included the frequent LCA-associated pathogenic intronic variant c.2992+1655A>G in CEP290 (Gene ID 80184, OMIM 610142) [16], as well as five deep-intronic variants in ABCA4 (Gene ID 496442, OMIM 248200) [17]. Pooled and phosphorylated probes were added to the capture reactions with 100 ng of genomic DNA from each individual to produce a sequence library for each individual. The libraries were amplified with 21 cycles of PCR (first denaturation step at 95 ºC for 4 min, followed by 21 cycles at 94 ºC for 30 s, 57 ºC for 30 s, 72 ºC for 1 min and then a final extension step at 72 ºC for 5 min), during which an 8-bp sample barcode was introduced. The barcoded libraries were then pooled and purified with AMPPureXP beads (Beckman-Coulter, Indianapolis, IN). Sequencing was performed on an Illumina NextSeq 500 system (San Diego, CA). Demultiplexed BAM files were aligned to a human reference sequence (UCSC Genome Browser hg19) via the Burrows-Wheeler Aligner (BWA) v.0.6.2 [18]. In-house automated data analysis pipeline and variant interpretation tools were used for variant calling. Only non-synonymous single nucleotide variants (nsSNVs), nonsense variants, putative splice site (±10 bps) variants, insertions, duplications, and deletions represented by more than 20 sequence reads were considered for further analysis. In addition, variants with a minor allele frequency (MAF) >0.5% in the Exome Aggregation Consortium Database (ExAC) Version 0.3.1 [19] were excluded from further investigation, with the exception of a few frequent variants in ABCA4 (i.e., c.2588G>C, c.5882G>A, and c.5603A>T [20-25]) that were selected separately.

In silico predictions

Several criteria were taken into account to evaluate the potential pathogenicity of variants: 1) The variants had previously been reported to be pathogenic, 2) the variants were absent in the ExAC database or showed a low MAF (<0.5%), 3) the variants caused a loss-of-function (LOF) of the encoded protein, such as nonsense, canonical splice site, and frameshift mutations, and 4) missense variants were predicted to be damaging by in silico prediction algorithms. The pathogenicity prediction tools used for missense variants were PolyPhen-2, Sorting Intolerant From Tolerant (SIFT), PhyloP score, the Grantham score, Mutation Taster, and Combined Annotation Dependent Depletion (CADD) score assessment. For the selection of non-canonical splice site variants, in silico predictions were performed by using five algorithms (SpliceSiteFinder-like, MaxEntScan, NNSPLICE, GeneSplicer, and Human Splicing Finder) via Alamut Visual software version 2.7 [26-30], and Biosoftware (Interactive Biosoftware, Rouen, France), by comparing the splicing scores for wild-type and variant nucleotides.

Segregation analysis

Candidate variants were validated and tested for cosegregation within the families using Sanger sequencing. For this purpose, specific primers were designed for PCR amplification of fragments containing the identified candidate causal mutations in each patient. PCR amplicons were sequenced using Sanger sequencing on an ABI 3730 DNA analyzer (Cologne, Germany) with BigDye chemistry v3.1.

Results

Study population

Fifty unrelated Iranian patients with non-syndromic IRD were included in this study, and their parents were first (34 families) or second cousins (16 families). All clinical information for the 50 families is summarized in Table 1. Among the 50 patients, 31 have RP, five have CRD, 11 have LCA, and three have STGD1.

Table 1

Clinical information overview for 50 Iranian families with inherited retinal dystrophies.

Diagnosis	Symptoms	Number of patients	Sex		Mean age of onset (±SD)	Mean age (±SD)	Family history of IRDs
			Male	Female			Yes	No
RP	NB, Tunnel vision, near blindness	31	15	16	11.6±4.2 years	31.5±6.3 years	21	10
CRD	Photophobia, defect in color vision, near blindness	5	1	4	11.7± 4.2 years	31.7±6.3 years	4	1
LCA	Nystagmus and total blindness	11	7	4	9.6±2.08 months	31.3±6.7 years	7	4
STGD1	Photophobia, defect in color vision, difficulty in adaption to the dark after sunlight exposure, central vision loss	3	0	3	11.2±5.1 years	30.6±2.5 years	2	1

NB: Night blindness, RP: Retinitis pigmentosa, CRD: Cone-rod dystrophy, LCA: Leber congenital amaurosis, STGD1: Stargardt disease, IRD: Inherited retinal dystrophies, SD: Standard deviation

Variant prioritization to identify causal variants

One affected individual from each of the 50 families was subjected to targeted NGS of 108 known IRD-associated genes (Appendix 1). The averages of 1.2 million reads on target per sample were obtained using our capture panel technology, with an average coverage of 213 reads per MIP. Moreover, an average of 88% of targeted regions had 10X coverage or more, which was sufficient for accurate variant calling. The pipeline initially called an average of 532 single nucleotide variants and 64 insertions or deletions for each sample. In this study, putative pathogenic variants were identified in 36 index cases (Table 2), and we could not find the pathogenic variants for the remaining 14 families. All putative disease-associated variants were validated with conventional Sanger sequencing, and segregation analysis was performed for all 36 families with identified pathogenic variants. Homozygosity was observed for 32/36 patients (88%), and four patients were compound heterozygous.

Table 2

Pathogenic variants identified by targeted next-generation sequencing in 36 Iranian probands affected with IRDs.

Index ID	Gene	Refseq	Phenotype	Exon/Intron	Nucleotide change	Protein change	Zygosity	Reference
066,592	ABCA4	NM_000350	STGD1	E1	c.1A>G	p.(Met1?)	Homo	[31]
066,595	ABCA4	NM_000350	STGD1	I43	c.6005+1G>A	p.?	Compound het	[32]
066,595	ABCA4	NM_000350	STGD1	E6	c.634C>T	p.(Arg212Cys)	Compound het	[33]
066,599	ABCA4	NM_000350	STGD1	E30	c.4462T>C	p.(Cys1488Arg)	Compound het	[34]
066,599	ABCA4	NM_000350	STGD1	E8	c.913C>T	p.(Gln305*)	Compound het	[35]
066,854	ABCA4	NM_000350	RP	E46	c.6385A>G	p.(Ser2129Gly)	Homo	[36]
066,863	ABCA4	NM_000350	RP	I20	c.3051–2A>G	p.?	Homo	Novel
066,876	ABCA4	NM_000350	RP	I27	c.4128+1G>A	p.?	Homo	Novel
066,877	AIPL1	NM_014336	LCA	E6	c.834G>A	p.(Trp278*)	Homo	[37]
066,872	BBS1	NM_024649	RP	E5	c.479G>A	p.(Arg160Gln)	Homo	[38]
066,568	BBS2	NM_031885	RP	E14	c.1705C>T	p.(Gln569*)	Compound het	Novel
066,568	BBS2	NM_031885	RP	I1	c.117+1G>T	p.?	Compound het	Novel
066,574	BBS2	NM_031885	RP	E2	c.224T>G	p.(Val75Gly)	Homo	[39]
066,891	CEP290	NM_025114	LCA	E54	c.7304T>G	p.(Val2435Gly)	Homo	Novel
066,871	CERKL	NM_001030311	RP	E6	c.847C>T	p.(Arg283*)	Homo	[40]
066,889;066892	CERKL	NM_001030311	CRD, RP	E13	c.1453_1454insA	p.(Pro485Hisfs*4)	Homo	Novel
066,880	CRB1	NM_201253	RP	E7	c.2234C>T	p.(Thr745Met)	Homo	[41]
066,886	CRB1	NM_201253	RP	E2	c.548G>A	p.(Cys183Tyr)	Homo	Novel
066,607	GUCY2D	NM_000180	LCA	E12	c.2383C>T	p.(Arg795Trp)	Homo	Novel
066,858	GUCY2D	NM_000180	LCA	E4	c.1216G>C	p.(Asp406His)	Homo	Novel
066,855	LCA5	NM_001122769	RP	E4	c.720+1G>A	p.?	Homo	[42]
066,852	LRAT	NM_004744	RP	E2	c.487C>G	p.(His163Asp)	Homo	Novel
066,853	LRAT	NM_004744	RP	E2	c.157_159dup	p.(Val53dup)	Homo	Novel
066,882	MERTK	NM_006343	RP	E2	c.390G>A	p.(Trp130*)	Homo	[43]
066,867	PDE6B	NM_000283	RP	E5	c.782_784del	p.(Phe261del)	Homo	Novel
066,884	PDE6C	NM_006204	CRD	E12	c.1612T>A	p.(Phe538Ile)	Homo	Novel
066,887	PDZD7	NM_001195263	RP	E6	c.793A>G	p.(Arg265Gly)	Homo	Novel
066,865	PROM1	NM_006017	RP	E25	c.2461C>T	p.(Arg821*)	Homo	Novel
066,859	RP1	NM_006269	RP	E4	c.2749C>T	p.(Gln917*)	Homo	Novel
066,885	RP1	NM_006269	RP	E4	c.1498_1499del	p.(Met500Valfs*7)	Compound het	[44]
066,885	RP1	NM_006269	RP	E4	c.4105C>T	p.(Gln1369*)	Compound het	[40]
066,860	RPE65	NM_000329	RP	E12	c.1299T>A	p.(Tyr433*)	Homo	Novel
066,864;066869	RPE65	NM_000329	RP	E10	c.1102T>C	p.(Tyr368His)	Homo	[45]
066,879	RPE65	NM_000329	RP	I4	c.354–1G>A	p.?	Homo	Novel
066,584;066866;066868	RPGRIP1	NM_020366	CRD, LCA	E10	c.1306G>T	p.(Ala436Ser)	Homo	Novel
066,878	SPATA7	NM_018418	RP	E6	c.506del	p.(Thr169Lysfs*25)	Homo	Novel

Overall, we identified 36 unique pathogenic variants based on the low allele frequency or absence of variants in ExAC, their high evolutionary conservation, the type of the amino acid substitutions, and the high predicted impacts of an amino acid substitution on the structure and function of a human protein in 19 genes. Of these, 20 were novel with the remaining 16 previously reported as IRD-causing variants [31-45]. Furthermore, 14 out of the 50 families remained genetically unexplained (Figure 1). Thus, the sensitivity of the gene panel for detection rate was 72%. The unique variants comprised 15 missense variants, nine nonsense variants, three deletions or insertions leading to a frameshift, six canonical splice site variants, two in-frame deletion and duplications, and one start loss. Pathogenicity was interpreted in accordance with different predictions (Table 3). According to the prediction tools in this study, all variants were annotated as pathogenic except one missense variant in ABCA4 (c.6385A>G; p.(Ser2129Gly)) that is predicted to be benign by PolyPhen-2 and tolerated by SIFT. However, this variant was predicted to be disease causing by MutationTaster and a high CADD score. As shown in Appendix 1, this variant is situated at the penultimate position of exon 46 and is part of the splice donor site consensus sequence. The canonical splice site nucleotides are GC, and therefore, the corresponding splice site was recognized only by SpliceSiteFinder-like. It predicted a strength of 91.8 (on a scale from 0 to 100) for the wild-type sequence and 80.3 for the mutant c.6385A>G, which is considered a large decrease in splice strength. Accordingly, we included the p.? in the predicted protein change. In addition, this homozygous variant has not been reported in ExAC.

Discussion

Although consanguineous marriages are common in the Iranian population [46], which makes the Iranian gene pool a valuable asset for genetic studies, little is known about inherited retinal degeneration in the Iranian community. Overall, 37.4% of Iranian marriages are consanguineous. Among them, 19.3% are first cousins and 18.1% second cousins [47]. Thus far, three small studies on genetic causes of Iranian patients with IRD have been reported [41,48,49].

Compared to the multiple targeted gene panel tests for different subtypes of IRDs that have been developed in the past few years with a detection rate of approximately 50–70% [50], we found similar percentages of pathogenic variants using our panel testing. Haer-Wigman et al. showed that using an in silico gene panel (based on whole exome sequencing) that was performed for 266 Dutch patients with visual impairment, in 52% of the cases, the genetic cause was identified [51]. However, in some studies the detection rate using panel testing was low. For instance, Fu et al. recruited 31 families with autosomal recessive RP and identified the pathogenic variant through a gene panel comprising 163 retinal disease genes. The authors achieved a detection rate of approximately 40% [52].

Variant c.6385A>G, apart from a predicted amino acid exchange, p.(Ser2129Gly), lowers the strength of a splice donor site of exon 46. A splice assay can shed light on its effect. Moreover, most of the variants detected in the patients are homozygous except four families that we identified as carrying compound heterozygous variants in RP1 (Gene ID 19888, OMIM 180100), ABCA4, and BBS2 (Gene ID 67378, OMIM 615981). The most frequently mutated gene was RPE65 (Gene ID 6121, OMIM 180069), mutated in 4/50 (8%) of the families (Figure 1). Sundaresan et al. reported that although the mutations in some genes such as RPE65 are a common cause of LCA in North America, they are rare in southern India [53]. To determine whether mutations in RPE65 are a common cause of RP in the Iranian population, we need to screen a larger cohort of RP probands in the future.

Our genetic findings, in particular the recurrence for variants in CERKL (Gene ID 375298, OMIM 608381), RPE65, and RPGRIP1 (Gene ID 57096, OMIM 605446) will aid genetic testing for IRD in the Iranian population. Furthermore, different variants in one gene can cause the same phenotype as, for example, different homozygous missense variants in GUCY2D (Gene ID 3000, OMIM 600179) in two unrelated patients cause LCA. Similarly, various pathogenic variants in LRAT (Gene ID 9227, OMIM 604863) and CRB1 (Gene ID 23418, OMIM 604210 ) in unrelated patients cause RP. However, the same variant in one gene can cause different phenotypes, such as a homozygous frameshift variant in the CERKL gene that led to CRD and RP in two unrelated families. Another example of this category is a homozygous missense variant in RPGRIP1 that caused LCA in two unrelated patients and CRD in one patient. In addition, different variants in one gene can cause different phenotypes as exemplified by ABCA4 that causes STGD1 and autosomal recessive (arRP).

Mutations in genes associated with Bardet-Biedl syndrome (BBS), such as BBS1 (Gene ID 582, OMIM 209900) and BBS2, can also be found in non-syndromic arRP [54-56]. The variants within these two genes identified in this study are associated with non-syndromic RP as the probands did not show extraocular features.

We identified homozygous CRB1 variants in two RP probands. One of these missense mutations, p.(Thr745Met), was previously detected in three Iranian RP cases by Ghofrani et al. [41]. Mutations in CRB1 have been shown to be responsible for 10% of LCA cases in northwestern Europe and 2.5% of arRP cases in the French population [10,57].

Although finding known and new mutations in known IRD-associated genes through cost-effective gene panel testing is important for patients who cannot afford the price of whole exome sequencing (WES) and whole genome sequencing (WGS), WES and WGS provide the opportunity to discover novel genes, as well as copy number variations (CNVs) and noncoding variants in IRD-associated genes [58-60]. The sequencing costs for this gene panel testing (MIP analysis), excluding the MIP synthesis costs and Sanger sequencing validation, were about €40 per sample, which is 10–20 times lower than the commercial price for WES.

Despite the high yield in this study, for 28% of the patients with IRDs, we still do not know the causal gene defects. The MIP design was based on the non-syndromic IRD-associated genes known in October 2013. In the meantime, another 60 non-syndromic IRD genes have been discovered that together may well explain another 10–15% of IRD cases. In conclusion, we demonstrated that this panel testing, employing MIPs, is a cost-effective targeted NGS-based method with a high detection rate for pathogenic variants in patients with IRDs.

Relative involvement of 19 inherited retinal disease-associated genes among the 50 probands.

Table 3

In-silico analysis of the non-truncating pathogenic variants identified in this study.

Gene	Amino acid change	Number of Probands	Mutationtaster	Polyphen-2	Grantham score	PhyloP score	SIFT	CADD score	ExAC MAF	Known/ Novel variant
ABCA4	p.(Arg212Cys)	1	Disease causing	Probably damaging	180	5.2	Deleterious	26.4	0.000117	Known
ABCA4	p.(Cys1488Arg)	1	Disease causing	Possibly damaging	180	8.9	Deleterious	21	1.42E-05	Known
ABCA4	p.(Ser2129Gly)	1	Disease causing	benign	56	7.9	Tolerated	19.96	N/A	Known
BBS1	p.(Arg160Gln)	1	Disease causing	Probably damaging	43	7.9	Deleterious	33	4.12E-05	Known
BBS2	p.(Val75Gly)	1	Disease causing	Probably damaging	109	9.1	Deleterious	27.2	N/A	Known
CEP290	p.(Val2435Gly)	1	Disease causing	Probably damaging	109	8.4	Deleterious	24.8	N/A	Novel
CRB1	p.(Thr745Met)	1	Disease causing	Probably damaging	81	5.5	Deleterious	15	8.26E-05	Known
CRB1	p.(Cys183Tyr)	1	Disease causing	Probably damaging	194	7.4	Deleterious	15.9	N/A	Novel
GUCY2D	p.(Arg795Trp)	1	Disease causing	Probably damaging	101	1.37	Deleterious	21.7	8.26E-06	Novel
GUCY2D	p.(Asp406His)	1	Disease causing	Probably damaging	81	9.1	Deleterious	20.9	N/A	Novel
LRAT	p.(His163Asp)	1	Disease causing	Probably damaging	81	7.6	Deleterious	20	N/A	Novel
PDE6C	p.(Phe538Ile)	1	Disease causing	Probably damaging	21	7.4	Deleterious	25.1	N/A	Novel
PDZD7	p.(Arg265Gly)	1	Disease causing	Possibly damaging	125	2.4	Deleterious	15.1	8.24E-06	Novel
RPE65	p.(Tyr368His)	2	Disease causing	Probably damaging	83	8.8	Deleterious	24.4	6.6E-05	Known
RPGRIP1	p.(Ala436Ser)	3	Disease causing	Probably damaging	99	5.4	Deleterious	15	N/A	Novel

MAF: Minor allele frequency, N/A: Not available