Protein modeling and in silico analysis to assess pathogenicity of ABCA4 variants in patients with inherited retinal disease

Patients

This retrospective study was conducted in accordance with Institutional Review Board approval (IRB#14–403E) of the Wills Eye Hospital and Thomas Jefferson University. A retrospective chart review was conducted for those patients with retinal disease–harboring mutations in the ABCA4 gene. All 13 patients were seen in the Ocular Genetics Clinic at Wills Eye Hospital from August 2010 to August 2015 by an ocular genetics specialist (AVL). All patients who were found to have only a heterozygous ABCA4 mutation or who had another suspected gene that might have modified the phenotype were excluded from this study. The medical records were reviewed, and the following data were collected: age of onset, age first seen, best-corrected visual acuity, and results of slit lamp examination and dilated fundus examination. Ophthalmic diagnostic investigations were reviewed, including fundus autofluorescence (AF), spectral-domain optical coherence tomography (SD-OCT), Goldmann visual field (GVF), and, if done, intravenous fluorescein angiography (IVFA). Full-field electroretinogram (ffERG) and multifocal electroretinogram (mfERG) were performed in accordance with the International Society of Clinical Electrophysiology and Vision standards [36]. Testing was customized for each patient as clinically indicated, based on age, visual acuity, and other factors. No clinical testing was done specifically for this study.

Genetic testing

The ABCA4 gene was screened for variants using a range of molecular diagnostic methods, in isolation or in parallel with other genes, as determined from the phenotype and considering the patients’ insurance coverage for testing (Appendix 1). When clinically indicated testing could not be performed—for example, due to a lack of insurance coverage—research testing was offered with the confirmation of detected variants performed in a Clinical Laboratory Improvement Amendments (CLIA)–certified laboratory.

Protein modeling and pathogenicity prediction

To predict the structural and functional consequences of missense ABCA4 variants, protein structure analysis was performed on AF2-predicted protein models and the existing experimental structures of the human ABCA4 protein [4,37,38]. Full-length (FL) and domain-specific protein modeling was performed using AF2 [35] for wild-type (WT) proteins and each variant protein separately, as described previously [34]. We used AF2 through the ColabFold notebook in the Google Colaboratory [39], which is reported to have a level of accuracy that is almost identical to that of AlphaFold v2.3.1 locally for monomeric structures [36]. The default parameters on the ColabFold notebooks (AlphaFold2.ipynb and AlphaFold2_advanced.ipynb) were used, and the highest confidence models (rank 1) were selected among the generated models. The predicted WT and variant structures were refined for energy minimization using the Amber-Relax option in ColabFold, and the “repair object” feature in the FoldX plugin in YASARA View (free version: 22.5.22) [40,41]. The refinement process was repeated several times until no movement of the side chains was observed in each model. PyMOL2 (The PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC) was used for model visualization and structure analysis.

The structure analysis of variants was performed using two approaches. First, we compared the WT and variant AF2 models to assess 3D conformational changes, to calculate the root-mean-square deviation (RMSD) alignment scores, and to identify potential impairments in secondary structure elements. Second, we evaluated the substitutions on the ABCA4 cryo-EM structures, considering various aspects, including stereochemical properties, binding, and molecular interactions (i.e., polar content, weak interactions, surface charge, and accessibility, cavities, and pockets). We measured the all-atom RMSD between the AF2-generated FL WT model and each variant, employing the “align” command without outlier rejection (cycles = 0) in PyMol2. We analyzed clashing interactions in available experimental structures and the refined AF models, as steric clashes resulting from the unnatural overlapping of newly positioned side chains with other residues can lead to conformational changes and structural destabilization due to a mutation.

To predict stability changes in protein variants, we used the FoldX plugin for the YASARA software. The ΔΔG free-energy change, calculated as the difference between the folded and unfolded free-energy states (ΔGmutant – ΔGwild-type), was determined using the “mutagenesis” feature of the software [40,41]. We reported the ΔΔG values obtained from all available ABCA4 experimental structures.

Furthermore, we evaluated the suitability of the AF2-predicted FL ABCA4 model for this stability analysis, by comparing the ΔΔG values per variant calculated using the AF2 model with those based on the experimentally determined structures. To facilitate this comparison, we conducted a Spearman’s rank correlation using the JMP statistical software package (JMP®, Version 17, SAS Institute, Inc., Cary, NC, 1989–2021). It is important to note that in silico ΔΔG values are approximate indicators of potential stability changes rather than absolute, experimentally determined values.

For the informatics-based pathogenicity prediction of missense variants, we used Combined Annotation Dependent Depletion (CADD) [42], Polyphen-2 [43], REVEL [44], and ConSurf [45] for the evolutionary conservation of the amino acids. We also used Human Splicing Finder (HSF-Pro) Version 4 to predict the possible splicing effects of all variants in the study [46].

Patient characteristics

We enrolled 13 patients diagnosed with retinal dystrophy caused by variants in the ABCA4 gene. Segregation analysis was performed for all patients to ascertain that the patients were biallelic for variants in the ABCA4 gene. To our knowledge, the patients did not harbor any other pathogenic variants in other confounding retinopathy genes. Genotype–phenotype comparisons were conducted for the ABCA4 variants occurring in these patients.

The clinical characteristics of the patient cohort are summarized in Table 1. The age first seen ranged from 6 to 65 years old, with a median age of 14 years old. The age of onset ranged from 6 to 15 years old, with a median age of 12 years. The best-corrected visual acuity varied from 20/30 to light perception. Two patients (A03M, A08M) had a visual acuity difference between the fellow eyes of more than two Snellen lines. In patients A02M and A05F, eccentric fixation was observed. The extraocular movements and applanation tonometry were normal in all patients. Nuclear sclerosis cataract was found in two patients (A02M and A04F).

OCT showed some degree of photoreceptor loss in 12 patients (92.3%), while 1 patient had a normal OCT. Seven patients were shown to have subretinal deposits with geographic atrophy, and five patients had only subretinal deposits. The foveal thickness in the right eye varied from 69 to 226 μm. Six patients had a normal GVF, three had a constricted visual field, and one had a central scotoma. For three patients, visual field testing was not performed. An ffERG was done in 12 patients, and the results were abnormal for 11 (91.66%). Seven patients had an abnormal response only in the photopic and four patients in both photopic and scotopic ERG. Following an isoelectric ERG, one patient had both abnormal scotopic and photopic responses. For all 12 patients who had an mfERG, abnormal results were produced. IVFA was done in 8 patients, and a silent choroid was found in 5 (62.5%).

There were 3 homozygous ABCA4 mutation patients and 10 compound heterozygous ABCA4 mutation patients. We identified 19 different ABCA4 variants among the 13 patients: 3 frameshift variants, 1 nonsense variant, 5 splicing variants, and 10 missense variants (Table 2). We also identified one VUS, four CIs, and one novel ABCA4 variant in this study, as reported in the ClinVar database (Table 3 and Table 4).

Protein modeling and structure analysis

To compare the deleterious impacts of the missense variants with the patients’ disease severities, we used AF2 protein modeling and all existing human ABCA4 cryo-EM structures (Protein Data Bank identifiers [PDB IDs]: 7lkp, 7lkz, 7e7i, 7e7q, 7e7o, 7e7q, 7m1p, 7m1q). ABCA4 WT and variant models were obtained in high confidence (pLDDT > 80%) and with high structural similarity with the experimental structures (Appendix 2, Appendix 3 and Appendix 4).

Variants inside the open reading frame (ORF) can have various impacts at the protein level. Some variants have clearly devastating consequences, such as those that lead to the misfolding and mislocalization of the ABCA4 or truncating mutations, which result in the complete loss of function (LoF). However, LoF can also result from local effects on the protein structure and stability that are more subtle and less obvious. Computational analysis aids in the evaluation of these types of effects. These types of local damage may lead to a spectrum of function loss, including affected intramolecular and ligand interactions. To interpret these types of functional damage, we used in silico protein analysis. For all missense variants, we examined the conformational changes, secondary structure elements, relative solvent accessibility, RMSD alignment scores, polar content alterations, clashing interactions, surface properties, and Gibbs ∆∆G free-energy changes. In silico analysis found a range of structural alterations in 10 missense ABCA4 variants carried by affected individuals (Table 3, Figure 2).

Of the 10 missense variants analyzed, 9 resulted in substantial structural damage. These variants caused clashing interactions with neighboring residues (p.C54Y, p.L541P, p.P1380L, p.A1598D, and p.G1961E), disrupted secondary structure elements (p.L541P and p.T1253M), replaced a hydrophobic buried residue with a hydrophilic one (p.A1598D and p.G1961E), or affected molecular interactions (p.R212C, p.R1098C, p.P1380L, and p.E2233V). The remaining variant, p.A1038V, led to a local conformational change in the backbone and reduced stability (Figure 2, Table 3). The stability analysis suggested that 6 of the missense variants (p.C54Y, p.L541P, p.A1038V, p.P1380L, p.A1598D, and p.G1961E) were particularly destabilizing with ∆∆G > 5 (Table 3 , Appendix 4).

We examined how the AF2 model determined ΔΔG values compared to those calculated using the ABCA4 cryo-EM structures. The ∆∆G profiles over 10 missense variants (and 2 complex alleles) were formed to compare these measures. Seven values were missing because the Arg-212 residue was unavailable in the 7lkp, 7m1p, and 7lkz structures, and the Gly-1961 and Glu-2233 residues were absent in the 7mlq structure. Appendix 5 shows agreement in the relative ∆∆G values for each missense variant among the experimental structures and the AF2 model. We formalized the similarity seen with a pairwise correlation analysis using Spearman’s rank correlation. Unadjusted p values showed 26 of 28 possible paired correlations to be significant (maximum p value = 0.0114); the ΔΔGs calculated based on AF2 specifically correlated with those done on the experimentally determined structures at approximately 0.89 or above for all except 7lkz (⍴ = 0.75; Appendix 6). Based on these findings, ΔΔG predictions from the AF2 model can be considered an alternative to those obtained from experimentally determined structures.

In silico protein structure analysis and phenotypic association

The missense variants were distributed across the ABCA4 protein sequence: there were 3 in ECD1, 3 in NBD1, 1 in ECD2, 1 in transmembrane domain-2 (TMD2), and 2 in NBD2 (Figure 1, Table 3). The variants identified in our cohort are described in terms of the clinical phenotype and assessment of the impact of the variant on protein structure.

The p.C54Y variant was identified in severely affected siblings, A04F and A05F, with ages of onset of 14 and 13 years, respectively. The two patients had the same biallelic combination, with p.E2233V on the opposite allele (in trans). The visual acuity of these patients was poor, ranging from 20/200 to 20/400. The predicted protein model and experimental structures showed disulfide bond breakage and clashing interactions (Figure 2A), along with a poor alignment score (RMSD = 1.26) and destabilization as reflected by the ∆∆G (+ 41.68 kcal/mol), supporting the deleterious effect of this variant (Table 3). The informatics tools also predicted the variant to be pathogenic (Table 3, Supplementary Table S2). The observed severity of this variant, based on the patient’s phenotype (Table 1), aligned well with the protein structural change.

Variant p.R212C was identified in a 6-year-old patient (A09M) with a rapidly progressing disease. The patient’s visual acuity at the time of the study ranged from 20/80 to 20/200, with constricted visual fields. OCT showed subretinal deposits and geographic atrophy (Table 1). The p.R212C variant was found in trans with the c.768G>T:p(V256=) synonymous variant that was predicted to be pathogenic due to its likely effect on splicing (HSF-Pro; Table 4) [46]. Given that one allele was predicted to be completely compromised (due to coding for a severely truncated polypeptide), the pathogenicity level of the p.R212C allele was expected to be directly related to the patient’s disease severity. The p.R212C variant was predicted to be pathogenic by the informatics tools (Supplementary Table S2), and structure analysis revealed a missing H-bond with Gln-190 due to the substitution (Figure 2B, Table 3). However, p.R212C is in the part of the ECD1 domain where the experimental structures have not been well resolved. This residue might be important for the structural flexibility of the domain.

Complex allele p.L541P/p.A1038V was detected in 3 patients (A06F, A07M, and A08M) in our study. The siblings A07M and A08M had the same biallelic combination with the c.5714+5G>A splicing variant on the trans allele, while the third patient had a p.Q1412* nonsense mutation in trans (Table 2). All three patients exhibited fundus flecks, subretinal deposits, and pigmentary clumps. The age of onset, visual acuity, and other findings were comparable in the siblings (Table 1). Protein structure analysis revealed that p.L541P had a more deleterious impact than did p.A1038V by causing significant impairment to the protein structure. This substitution is located within an α-helical region and is predicted to disrupt the α-helix structure and cause clashing interactions, leading to conformational changes (RMSD = 1.27; Figure 2C, Table 3). Although we did not observe a noticeable structural difference between the WT and p.A1038V variant, it showed an overall poor alignment score (RMSD = 1.25) and ∆∆G (+ 6.56 kcal/mol), indicating that it can be a destabilizing mutation (Figure 2D). The informatics-based tools suggested that the p.L541P variant is pathogenic while p.A1038V is benign (Table S2). For the complex allele L541P/A1038V, the computational analysis of the variant on the ABCA4 protein structure aligned well with patient phenotype.

Variant p.R1098C was found in patient A13M, who had an age of onset of 11. It was in trans with a frameshift mutation (Ala222Glnfs*19), which leads to a complete loss of function. The Arg1098 residue, localized inside the NBD1, was found to have an essential role in the NBD1–RD2 interdomain interaction by forming a salt bridge, which is disrupted by the Cys substitution (Figure 2E). Prediction tools classified p.R1098C as pathogenic (Table 3, Supplementary Table S2). For this variant, the computational analysis of the effects on the ABCA4 protein structure coincided with the disease phenotype.

Complex allele p.T1253M/p.G1961E was identified in patient A03M, with the c.5714+5G>A splicing variant in trans, which is predicted to affect the splicing pattern according to HSF-Pro [46]. The patient presented with a “beaten bronze” appearance in fundus imaging, with lipofuscin deposits centrally and a mild surrounding ring of hyperpigmentation (Table 1). Remarkably, 10 years after the diagnosis, the patient’s vision remained relatively preserved compared to the other patients in our cohort. In silico analysis predicted that both missense variants individually cause protein structural changes. The p.T1253M variant, located within RD1, led to a major structural alteration with a degenerated secondary structure element (Figure 2F). While, the p.G1961E variant replaced a hydrophobic buried residue with a hydrophilic amino acid, which was predicted to destabilize the structure, as indicated by the ∆∆G value (Table 3) and as shown in Figure 2I. The prediction tools unanimously supported the pathogenicity of both variants (Supplementary Table S2).

The p.P1380L variant was identified in a 10-year-old homozygous patient (A010M) with a highly progressive ABCA4 retinopathy. The patient’s severe phenotypic findings (Table 1) are consistent with protein modeling that indicates the Pro → Leu substitution within the transmembrane domain has significant structural consequences. The substitution causes steric clashes with neighboring residues in all possible positions, leading to altered protein conformation (RMSD = 1.26). Additionally, the variant protein model superimposed with the substrate-bound cryo-EM structure suggests that this proline residue likely plays a role in a kinked transmembrane helix, with implications of a functional role in substrate translocation (Figure 2G).

The p.A1598D variant, located in the ECD2 domain, was identified in patient A01F, who presented with macular flecks and alterations in the retinal pigment epithelium (RPE) subretinal deposits, with a visual acuity of 20/200 in both eyes. The patient was homozygous for the p.A1598D variant, with no other ABCA4 variants detected. In silico analysis revealed that by introducing a hydrophilic residue to a buried location, the p.A1598D variant resulted in destabilization (∆∆G = + 13.56 kcal/mol) and caused steric clashes (Figure 2H). However, the conventional pathogenicity assessment of this variant resulted in inconclusive predictions. PolyPhen-2 and REVEL both classified it as possibly damaging, with scores of 0.57 and 0.685, respectively. In contrast, the CADD score of 16.3 was below the threshold of 20, indicating a lower likelihood of pathogenicity. Ala-1598 is not evolutionarily conserved [45], but the variant has a low allele frequency (2.63 × 10⁻⁵; Supplementary Table S2) [47]. Overall, the patient’s severe phenotype was consistent with the substantial changes observed in the protein structure predicted by structure analysis.

Variant p.E2233V was identified in siblings A04F and A05F, who had comparable age of onset and visual acuity (Table 1). A04F presented with large retinal pigmentary clumps involving the macula and mild periphery, choroidal attenuation, but no fleck lipofuscin deposits. A05F showed geographic atrophy in the macular area with mid-periphery pigmentary clumping, with several subretinal flecks in the macula. A grayish elevation in the superior macula could represent either a choroidal nevus or subretinal neovascular membrane. Both patients had relatively healthy optic nerves, with minimal pallor as well as mild retinal vascular attenuation (Table 1). This novel variant is located in RD2. Computational analysis revealed a loss of an H-bond with a residue in the NBD1 domain, significantly impairing the protein structure (RMSD = 1.35; Figure 2J). This variant was also predicted to affect splicing by altering the auxiliary sequences’ ESE/ESS motifs ratio (−10; HSF-Pro [46]). Although the true consequence of this variant is not yet fully understood, we saw an agreement between the computational analysis of the missense effects on the ABCA4 protein structure and the patient phenotype.

Our cohort had three patients who did not have a missense ABCA4 variation. Patient A02M carried the c.5018+2T>C variant in a homozygous state. This intronic variant is predicted to affect the splicing pattern by causing a broken WT splice site (HSF-Pro [46],), resulting in translated product lacking the last quarter of the protein. Patient A11F was compound heterozygous with the c.859–13T>C and c.4539+2028C>T intronic variants, present in the ClinVar database as a CI and VUS, respectively. This variant was also predicted to affect splicing by altering the ratio of Exonic Splicing Enhancer (ESE) to Exonic Splicing Silencer (ESS) motifs, with a score of −10 as predicted by the HSF-Pro [46]. The c.4539+2028C>T deep intronic variant has previously been shown to cause 345-nt pseudoexon insertion with a likely effect on the protein level (p.R1514fs) [48]. Unless there are other undetermined variations in the patient, the phenotype suggests that these intronic variants are pathogenic. Patient A12F carried two different frameshift variants—p.I284Vfs*33 and p.V2062Tfs*51—both with high levels of pathogenicity.

Conclusion

In silico methods are essential elements of biomedical research. Computational protein structure analysis is a practical tool used to understand protein structure and function, as well as to infer the pathogenic impacts of sequence variations. This study illustrated the likely structural consequences of 10 missense variants. Our study demonstrated that a disease severity that is closely defined by the age of onset and by visual acuity is well associated with in silico–predicted protein structural changes caused by missense ABCA4 variants. Determining disease-causing variants is a crucial part of clinical management, and we believe that our research will contribute to elucidating the molecular pathology of variants associated with ABCA4 retinal diseases.