A comprehensive genetic landscape of inherited retinal diseases in a large Pakistani cohort

Clinical and demographic data

Our cohort comprised a total of 722 IRD patients, belonging to 213 different and unrelated families. Of the index patients, 115 (54%) were female and 98 (46%) were male. Most of the patients (156, 73.3%) were from the Khyber Pakhtunkhwa province, 29 (13.6%) were from the Sindh province, and 28 (13.1%) from Punjab (Fig. S2 and Table 1), many of them residing in very remote areas.

All patients were broadly classified as having IRD based on descriptive clinical evaluation. For most of the patients, a clinical sub-diagnosis (RP, CRD, STGD, etc.) could not be obtained, mostly because of a lack of specialized ophthalmic examination facilities in the sampled areas. The only available ophthalmic evaluation data is provided in Fig. S1. A generic diagnosis was therefore inferred based on answers given to a medical questionnaire, filled out either by the patients themselves or by their legal guardians (Supplementary Material S1).

Molecular findings and mutational spectrum

To ensure an accurate representation of the molecular landscape of IRDs within the cohort and to avoid data inflation, we first performed a relatedness analysis using global WES data, which allowed us to exclude 34 genetically related pedigrees or duplicate index patients and analyze only the 213 unrelated pedigrees mentioned above (Fig. S3). Molecular diagnosis was achieved in 171 of them (80.3%), through the identification of pathogenic or likely pathogenic (PLP) variants in genes associated with IRDs. Out of these, 162 had homozygous mutations, two were compound heterozygotes, three had heterozygous dominant mutations, and four were hemizygotes for mutations in a gene on the X chromosome, for a total of 335 PLP alleles. In five of the 213 families (2.3%), molecular diagnosis was considered uncertain due to the presence of variants of uncertain significance (VUS), whereas diagnosis remained elusive in 37 (17.4%) pedigrees (Table 1, Fig. 1A), who were categorized as unsolved.

A Classification of patients based on genetic findings. Patients with pathogenic/likely pathogenic (PLP) variants were categorized as ‘solved’’ or ‘likely solved’’ according to the criteria detailed by Peter et al.¹⁶, those with variants of uncertain significance (VUS) as ‘uncertain’’, and all the others as ‘unsolved.’’ B Further stratification of solved and likely solved cases, according to disease genes. Percentages refer to the number of index patients found to be positive to mutations in a given gene, over the total number of patients from the solved or likely solved classes.

Overall, we identified 129 different PLP variants (65 P and 64 LP) in 60 IRD-associated genes. Of these, 54 (41.8%) were novel disease-causing DNA changes (Supplementary Data 1 and 2), which were assessed to co-segregate with disease in their respective pedigrees. Missense variants (n = 44 out of 129, 34.1%), were the most represented DNA changes and were detected in 67 individuals. Nonsense variants (n = 30, 23.2%) were present in 37 individuals, and small insertions/deletions (indels) leading to frameshifts were responsible for disease in 37 individuals as well (n = 30). Canonical splicing variants accounted for seven events (5.4%, in seven individuals), and structural variants/copy number variations (SVs/CNVs) were detected in five pedigrees (3.9%, in five individuals). Moreover, 13 variants (10.1%, in 20 individuals) were classified as “others.” These included non-canonical splicing mutations, non-frameshift indels, synonymous variants, and stop-loss variants (Supplementary Data 1, Fig. S4).

The seven most frequently mutated genes in our cohort were ABCA4 (9.4%, n = 16 index patients) and CRB1 (9.4%, n = 16), MYO7A (5.3%, n = 9), PDE6B (4.7%, n = 8), and RP1 (4.7%, n = 8), as well as RPE65 (4.1%, n = 7) and TULP1 (4.1%, n = 7) (Supplementary Data 1, Fig. 1B). RPGR and RHO, along with other genes, were among the least frequently mutated ones in our cohort, with each gene having been observed in only one family at a time. Unlike in patients from other cohorts, mutations in USH2A and EYS were also only marginally present in our families, accounting for disease in only three and four families, respectively. Collectively, 30.4% of the mutated genes were detected individually only once or twice (Supplementary Data 1, Fig. 1B).

In terms of disease inheritance, the largest majority of the pedigrees with a clear-cut molecular diagnosis (164 out of 171, or 95.9%) presented with disease segregating as an autosomal recessive trait, four (2.3%) of them displayed an X-linked inheritance, while three (1.8%) displayed an autosomal dominant inheritance pattern (Table 1).

Recurrent homozygous mutations and overall genomic autozygosity

In our cohort, three mutations were identified as highly recurrent disease-causing alleles, mostly in the homozygous state. These comprised two missense mutations (CRB1: NM_201253.3:c.1459T>C, p.Ser487Pro, and ABCA4: NM_000350.3:c.214G>A, p.Gly72Arg) and one frameshift deletion (MYO7A: NM_000260.4:c.4838del, p.Asp1613ValfsTer32) (Fig. 2). They were observed 16, 15, and eight times in eight, eight, and four genetically unrelated pedigrees, respectively (Fig. 2). The detection of these rare pathogenic alleles across ethnically and geographically matched families suggested the possibility that such variants represented founder mutations, inherited from common distant ancestors. Haplotype analysis confirmed this hypothesis, since all patients shared one distinct haplotype encompassing each of the three alleles. More specifically, the CRB1 missense was comprised within a 7.2 Mb common haplotype, whereas the ABCA4 and MYO7A variants lay within a 1.9 Mb and a 2.46 Mb interval, respectively (Supplementary Data 3). Additional recurrent alleles were observed six times (in three probands) and four times (in two probands), in as many as 17 genes within our cohort (Fig. 2). However, most of the pathogenic variants detected were not observed in more than one family (210 alleles out of 335, 62.7%) (Fig. 2). Of note, affected individuals homozygous for the MYO7A variants had a phenotype compatible with Usher syndrome type 1, including congenital deafness.

Each slice of the chart corresponds to a specific disease-causing allele, except for variants occurring only once or twice. n represents the total allele count within a specific group.

Interestingly, only a small fraction of the PLP variants detected, specifically seven variants from five probands, were found in either a heterozygous or compound heterozygous state. The remaining ones were either in a homozygous or hemizygous (on the X chromosome) state, in noticeable contrast with observations in patients from other populations^12,16,17,18. We speculated that this low ratio of heterozygous variants may be attributed to the high cumulative autozygosity characteristic of the Pakistani population. Therefore, we calculated the cumulative genome-wide autozygosity of index patients from our families, comparing it to that of cohorts from various other ethnic groups. This analysis showed that our patients displayed a median cumulative autozygome of 285 Mb (for autosomes), which was significantly higher than aggregated autozygosity intervals detected in IRD individuals from other countries, such as Japan (120 Mb), Sweden (76 Mb), Hungary (70 Mb), Switzerland (67 Mb), Portugal (67 Mb), or Italy (66 Mb) (Fig. 3).

Boxplots refer to the sum of all runs of homozygosity (ROHs) within individual genomes of IRD patients from different cohorts. Data from individual patient are indicated by black circles, while median values are indicated by horizontal thick bars, with numbers (Mb: megabases). Standard notation for boxplots applies to other components of the graph.

Finally, in one individual, MAPK0313, a homozygous variant (p.His620Arg in PDE6B) was found to co-segregate with the condition in the proband’s branch of the family, but, likely because of high endogamy, was also identified at the heterozygous state in an affected member from another branch, whose condition was due to a yet to be identified genotype (Supplementary Data 1).

Copy number variations

Another class of prevalent pathogenic variants in IRDs is that of SVs/CNVs¹⁹. In our cohort, we could identify causative large deletions in five unrelated families, affecting EYS (two different events), CRB1, IMPG2, and RP1 (Figs. S5, S6, S7, S8, and 4). The most frequent of them, a large deletion in RP1, was also validated by PCR amplification and Sanger sequencing, which revealed the presence of a 11 kb gap (hg19, chr8:55532084_55543199del) (Fig. 4).

A Genotypes and variant segregation analysis in the two pedigrees carrying the deletion. These families are known to be related. B Visualization of the large structural deletion involving RP1, by OFF-PEAK⁵⁴. The shaded pink region indicates the deletion. C Molecular analysis of the breakpoint. Primers CR-7002 and CR-7003 were designed in the proximity of the borders of the deletion. In control subjects, such primers do not yield any PCR product in standard conditions, since the DNA stretch that separates them is too long (~ 11 kb). However, in subjects bearing the deletion, they resulted in a 343 bp amplification product, which was subsequently sequenced to reveal the precise breakpoint. Genomic coordinates are given with respect to the hg19 reference sequence (hg38 coordinates of the deletion are: chr8:54619524-54630639).

Variants affecting RNA splicing

Two of the variants identified were predicted by in silico tools to alter the splicing of the genes PDE6C and NYX, without affecting the canonical splicing motifs. The first was a synonymous change, NM_006204.4:c.633G>A, p.Glu211=, affecting the last base of exon 2. It was detected homozygously in three pedigrees that, according to both family history and molecular analysis, were unrelated, but originated from the same small town. The second was an intronic change (NM_001378477.3:c.22+5G>T), occurring in the second intron of NYX, on the X chromosome of multiple affected male individuals (Fig. 5). It was shared by two families, both genetically unrelated and residing in different and distant regions of the country.

The co-segregation analysis of the pathogenic variants detected in PDE6C (c.633G>A) and NYX (c.22+5G>T), across three and two unrelated families, respectively, is shown.

We assessed the predicted effect on RNA splicing of both variants by constructing minigenes derived from pCI-NEO-RHO, a plasmid routinely used for ex-vivo splicing assays²⁰. Upon transfection of HEK-293FT cells, we detected aberrant splicing for both the PDE6C and NYX constructs bearing the variants detected in our patients, therefore validating their pathogenic role (Fig. 6). Specifically, the minigene involving the synonymous change in PDE6C revealed the skipping of exon 2, in turn resulting in an in-frame deletion of 51 codons (NM_006204.4:c.481_633del, p.Asn161_Glu211del) (Fig. 6A, B). Similarly, the intronic variant detected in NYX (c.22+5G>T) also resulted in the skipping of exon 2 of this gene, leading to the loss of the start codon (Fig. 6A, B).

A Agarose gel electrophoresis of the RT-PCR-based amplification of transcripts from HEK-293FT cells, transfected with plasmids carrying wild-type and patient-derived sequences. n.t. = non-transfected cells; wt = cDNA from cells transfected with plasmids carrying wild-type sequences; mut = cDNA from cells transfected with plasmids carrying either c.633G>A in PDE6C or c.22+5G>T in NYX. B Electropherograms of the PCR products depicted in (A). Exons 3 and 5 of the rhodopsin gene derived from the backbone of plasmid pCI-NEO-RHO (see Methods).