Pathogenicity of new BEST1 variants identified in Italian patients with best vitelliform macular dystrophy assessed by computational structural biology

Abstract

Background

Best vitelliform macular dystrophy (BVMD) is an autosomal dominant macular degeneration. The typical central yellowish yolk-like lesion usually appears in childhood and gradually worsens. Most cases are caused by variants in the BEST1 gene which encodes bestrophin-1, an integral membrane protein found primarily in the retinal pigment epithelium.

Methods

Here we describe the spectrum of BEST1 variants identified in a cohort of 57 Italian patients analyzed by Sanger sequencing. In 13 cases, the study also included segregation analysis in affected and unaffected relatives. We used molecular mechanics to calculate two quantitative parameters related to calcium-activated chloride channel (CaCC composed of 5 BEST1 subunits) stability and calcium-dependent activation and related them to the potential pathogenicity of individual missense variants detected in the probands.

Results

Thirty-six out of 57 probands (63% positivity) and 16 out of 18 relatives proved positive to genetic testing. Family study confirmed the variable penetrance and expressivity of the disease. Six of the 27 genetic variants discovered were novel: p.(Val9Gly), p.(Ser108Arg), p.(Asn179Asp), p.(Trp182Arg), p.(Glu292Gln) and p.(Asn296Lys). All BEST1 variants were assessed in silico for potential pathogenicity. Our computational structural biology approach based on 3D model structure of the CaCC showed that individual amino acid replacements may affect channel shape, stability, activation, gating, selectivity and throughput, and possibly also other features, depending on where the individual mutated amino acid residues are located in the tertiary structure of BEST1. Statistically significant correlations between mean logMAR best-corrected visual acuity (BCVA), age and modulus of computed BEST1 dimerization energies, which reflect variations in the in CaCC stability due to amino acid changes, permitted us to assess the pathogenicity of individual BEST1 variants.

Conclusions

Using this computational approach, we designed a method for estimating BCVA progression in patients with BEST1 variants.

Background

Best vitelliform macular dystrophy (BVMD) (OMIM #153700), also known as Best’s disease, is an autosomal dominant slowly progressive form of retinal macular degeneration. The typical central yellowish yolk-like lesion due to accumulation of lipofuscin in the retinal pigment epithelium (RPE) usually appears in childhood. The lesion gradually worsens [1] causing progressive macular atrophy or fibrosis and subsequent loss of visual acuity [2].

Five stages of the disease have been described: in the first or pre-vitelliform stage, usually discovered incidentally, subtle RPE alterations of the macula cause no symptoms. The second or vitelliform stage is characterized by a well-defined 0.5 to 2 mm diameter “egg-yolk” lesion in the macula, and patients may experience symptoms such as metamorphopsia, blurred vision and a decrease in visual acuity. In stage 3, this deposit is partially reabsorbed and deposited in a layer in the macula, known as “pseudohypopyon”. In the fourth or vitelliruptive stage, there is partial reabsorption of the material (scrambled-egg lesion) and macular atrophy. Macular fibrosis develops in stage 5 [3].

The vast majority of cases are caused by a variant in the BEST1 gene (also known as VMD2) which encodes bestrophin-1 (BEST1) [4].

Bestrophin-1 protein is expressed in the basolateral membrane of the RPE [4]. It is essential for normal eye development during embryogenesis and for retinal homeostasis throughout life [5]. Bestrophin-1 belongs to the family of calcium-activated chloride channels (CaCCs) that regulate the flow of chloride and other monovalent anions across cell membranes in response to intracellular Ca²⁺ levels [6,7,8]. These channels occur in various types of cell. The mechanisms of CaCC anion selectivity, calcium-dependent gating and the bestrophin-1 domains involved in channel regulation are not fully understood. The X-ray structure of BEST1 from Gallus gallus, which shares 74% sequence identity with human BEST1, showed that five bestrophin-1 molecules compose a CaCC, which contains a single long pore along a five-fold symmetry axis (Fig. 1) [9, 10]. The pore extends from the extracellular side across the cell membrane and through the cytosol. Its variable diameter defines two narrow hydrophobic regions: the neck and the aperture, which are thought to be related to channel gating and anion selectivity [9, 10]. The X-ray structure of BEST1 CaCC also shows five calcium-binding sites in the cytosolic portions of the subunits (Ca²⁺ clasps) and three binding sites for Cl⁻ for each subunit in the channel pore [9, 10].

Fig. 1

X-ray structure of bestrophin-1 from Gallus gallus [9, 10] forming a calcium-activated chloride transmembrane channel (CaCC) composed of five identical subunits and containing a single long pore along the axis of symmetry. a Top view (extracellular side), b side view, c bottom view (cytosol side, channel aperture). The CaCC is shown in Corey–Pauling–Koltun (CPK) representation, dark green spheres—chloride ions (A), orange spheres—potassium ions (not visible), blue spheres—calcium ions (not visible)

The N-terminal region of bestrophin-1 (residues 1–390 in human BEST1) is highly conserved in eukaryotic CaCCs and has been shown sufficient for CaCC activity and Cl⁻ ion transport [11]. The C-terminal cytosolic portion of the protein (amino acids 391–585) is predicted to be unstructured with an as yet undefined role in channel regulation [10]. A water molecule, acidic residues Asp301, Asp304 and backbone carbonyl groups of two neighboring bestrophin-1 subunits cooperate in calcium ion coordination in the Ca²⁺ clasps [9, 10], Fig. 2.

Fig. 2

Crystal structure of the “Ca²⁺ clasp” of BEST1 from Gallus gallus [9, 10]. The calcium binding site is formed by two bestrophin-1 subunits represented by the red (subunit A) and blue ribbons (subunit B) and is involved in the CaCC activation. Ca²⁺ is shown as a dark blue sphere. Yellow lines indicate calcium ion coordination by residues side chains and backbone carbonyl groups of both subunits. A water molecule coordinating the calcium ion in the crystal structure is not shown

Most of the approximatively 200 known disease-causing variants in BEST1 associated with retinal degenerative disorders have been reported to induce CaCC aberrations [5, 12,13,14,15]. Alterations in BEST1 CaCC functions lead to a diminished electrooculogram (EOG) light peak to dark trough ratio typical of BVMD [16]. BVMD is a heterogeneous pleomorphic disease underling different phenotypes and with alternative modes of inheritance. Although most cases show autosomal dominant inheritance, recessive inheritance has also been described [7, 9, 10]. Patients with recessive inheritance may have the classic features of Best disease [17], including a central vitelliform lesion, or may have extramacular punctate flecks without any notable central lesion [18]. Variants in BEST1 are also responsible for other clinically distinct human diseases: autosomal dominant vitreoretinochoroidopathy (ADVIRC) (OMIM #193220), autosomal recessive bestrophinopathy (ARB) (OMIM #611809) and retinitis pigmentosa (RP) (OMIM #613194) [19].

Here we report the results of the characterization of BEST1 variants in an Italian population consisting of 57 probands and 13 families, thus contributing to the molecular epidemiology of Best disease in our country. Six of the 27 BEST1 variants detected are novel.

All variants were assessed by computational structural biology. The conclusions of computational modelling led to the formulation of simple approximate quantitative structure-pathogenicity relationships (QSPR) of the BEST1 variants and a tool for predicting visual acuity progression in individuals with inherited BEST1 variants.

Methods

Genetic testing

Fifty-seven probands with clinical features suggesting Best’s disease were examined in different eye clinics and hospitals. Most of the patients underwent complete ophthalmic examination, including logMAR best-corrected visual acuity (BCVA), anterior and posterior segment examination and retinal imaging. Imaging included fundus photography (Zeiss, Visucam, Oberkochen, Germany), spectral domain optical coherence tomography (SD-OCT, Spectralis HRA + OCT, Heidelberg Engineering, Heidelberg, Germany), infra-red (IR) imaging and blue fundus autofluorescence (B-FAF) (Spectralis HRA and HRA II, Heidelberg Engineering). Some patients also underwent electrophysiological testing including electro-oculogram and multifocal electroretinogram [EOG, mfERG, VERIS Clinic 4.9, Electro-Diagnostic Imaging, San Mateo, CA and Retimax instrument (CSO, Firenze, Italy)] recordings and ultra-structural morphology evaluation by adaptive optics (AO, rtx1, Imagine Eyes, Orsay, France) retinal imaging. A small number of patients were only addressed to our laboratories for genetic testing from different institutions, therefore with partial clinical data.

A genetic test was performed to confirm the diagnosis of Best disease; whenever possible, ophthalmic examination and genetic analysis were extended to members of the proband’s family.

A total of 57 blood samples from patients and 18 from family members were received and analyzed by MAGI Laboratory (MAGI’s Lab, Rovereto, Italy). DNA was extracted by commercial kit (Blood DNA Kit E.Z.N.A.; Omega Bio-Tek Inc., Norcross, GA, USA). Extracted DNA underwent polymerase chain reaction (PCR) to amplify all coding regions (exons 2–11) and the intron/exon junctions of BEST1 (NM_004183.3). The products were purified and sequenced with a Beckman Coulter CEQ 8000 sequencer (Beckmann Coulter, Milan, Italy). All laboratory protocols are available upon request.

The electropherograms of all amplified fragments were analyzed with ChromasPro 1.5 software (Technelysium Pty Ltd, Australia; http://technelysium.com.au) and then compared with the reference sequences using Basic Local Alignment Search Tool 2 (BLAST2) sequences (http://blast.ncbi.nlm.nih.gov).

For the identification of variants, the Retina International database (http://www.retina-international.org), the Human Gene Mutation Database Professional version 2017.2 (https://portal.biobase-international.com/hgmd/pro/), the Exome Variant Server (http://evs.gs.washington.edu/EVS/) and dbSNP (https://www.ncbi.nlm.nih.gov/snp/) were also consulted.

New variants were analyzed for their putative pathogenicity using three on-line softwares: Polyphen 2 (Polymorphism Phenotyping v2; http://genetics.bwh.harvard.edu/pph2), SIFT (Sorting Intolerant From Tolerant; http://sift.bii.a-star.edu.sg/) and Mutation Taster (http://mutationtaster.org/). Pedigrees were designed according to Bennett et al. [20] using HaploPainter software (http://haplopainter.sourceforge.net/). New variants were classified according to the American College of Medical Genetics and Genomics standards and guidelines [21].

Molecular modelling

A refined 3D template model of human CaCC was prepared from the X-ray structure of chicken BEST1 (PDB entries 4RDQ and 5T5 N, 409 residues, resolution 3.1 Å) [9, 10] by protein homology modelling and molecular geometry optimization. Amino acid sequences of the N-terminal domain of Homo sapiens bestrophin-1 (UniProtKB O76090, residues 1-366) and of bestrophin-1 from Gallus gallus were aligned by the EMBOSS Needle sequence alignment server (73.6% identity, 88.3% similarity, gaps 0.5%) (https://www.ebi.ac.uk/Tools/psa/emboss_needle/). The 3D homology model of human bestrophin-1 channel subunit A was built by comparative modelling from the chicken CaCC template structure with the help of the Prime protein structure prediction approach [22] using standard parameters (Schrödinger, LLC, New York, 2014). The loops of the human bestrophin-1 model were perfected using the refine loops tool of Prime. The homology model was inspected with help of Prime to ensure that there was no violation of protein stereochemistry. The human BEST1 model contained one Ca²⁺ cation, two K⁺, three Cl⁻ and one water molecule coordinating the calcium ion. After building the 3D model of BEST1 protein, five protein subunits were superimposed on the chicken CaCC template structure using the Prime structural alignment module (Schrödinger, LLC, New York, 2014) to obtain the quaternary structure of the human chloride channel. A dimer formed by two BEST1 subunits was separated from the CaCC model and further refined to convergence by energy-minimization using molecular mechanics (MM). During minimization, the OPLS-2005 force field [23,24,25] and the generalized Born implicit solvation model (GB/SA) [26] were employed with MacroModel software (Schrödinger LLC, New York, NY, 2014). At first, the protein backbone atoms and all ions were constrained in their initial positions, while all other atoms were unrestrained. In the final minimization step, all atoms were set free to find their relaxed positions in the model structure of human bestrophin-1 subunits.

To study the effect of individual variants detected in probands, residue replacement with optimal side chain rotamer selection were conducted in the BEST1 dimer using MacroModel, followed by extensive dimer and monomer structure minimization. The consequence of variants for the bestrophin-1 subunit binding to form the pentameric channel and of subunit binding capacity to Ca²⁺ ions needed for channel regulation were estimated with the help of relative dimerization and calcium binding energies (ΔΔE_dim and ΔΔE_Cabin) computed with respect to the native BEST1 protein. The ΔΔE_dim energy estimates the extent of possible damage or change in CaCC formation, expressed as the change in binding energy between neighboring bestrophin-1 subunits A and B in CaCC, caused by a variant (var), as compared to the reference native BEST1 [27,28,29]:

$$\Delta \Delta {\text{E}}_{ \dim } = \Delta {\text{E}}_{ \dim } \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{AB}}}} } \right\}_{\text{aq}} - \Delta {\text{E}}_{ \dim } \left\{ {{\text{BEST}}1_{{{\text{native}},{\text{AB}}}} } \right\}_{\text{aq}}$$

(1)

$$\Delta {\text{E}}_{ \dim } \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{AB}}}} } \right\}_{\text{aq}} = {\text{E}}_{\text{tot}} \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{AB}}}} } \right\}_{\text{aq}} - {\text{E}}_{\text{tot}} \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{A}}}} } \right\}_{\text{aq}} - {\text{E}}_{\text{tot}} \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{B}}}} } \right\}_{\text{aq}}$$

(2)

where E_tot{Z}_aq is the total molecular mechanics (MM) energy of hydrated {}_aq protein Z (dimer AB or monomers A and B) computed by MacroModel (Schrödinger LLC, New York, NY, 2014). The relative Ca²⁺ binding energy (ΔΔE_Cabin) reflects the extent of possible damage or change to channel regulation expressed as the altered ability of the BEST1 dimer AB to bind calcium ions due to a variant, with respect to the reference native BEST1:

$$\Delta \Delta {\text{E}}_{\text{Cabin}} = \Delta {\text{E}}_{\text{Cabin}} \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{AB}},{\text{Ca}}}} } \right\}_{\text{aq}} - \Delta {\text{E}}_{\text{Cabin}} \left\{ {{\text{BEST}}1_{{{\text{native}},{\text{AB}},{\text{Ca}}}} } \right\}_{\text{aq}}$$

(3)

$$\Delta {\text{E}}_{\text{Cabin}} \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{AB}},{\text{Ca}}}} } \right\}_{\text{aq}} = \frac{1}{2}\left[ {{\text{E}}_{\text{tot}} \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{AB}},2{\text{Ca}}}} } \right\}_{\text{aq}} - {\text{E}}_{\text{tot}} \left\{ {{\text{BEST}}1_{{{\text{var}},{\text{AB}}}} } \right\}_{\text{aq}} - 2{\text{E}}_{\text{sol}} \left\{ {{\text{Ca}}^{2 + } } \right\}_{\text{aq}} } \right]$$

(4)

where E_tot{Z}_aq and E_sol{Ca²⁺}_aq are the total MM energy of the hydrated bestrophin-1 dimer with or without bound calcium ions, or the solvation energy of the Ca²⁺ ion.

QSPR of BEST1 variants affecting the visual acuity of the probands were elaborated with help of the Cerius² software package (Cerius² Life Sciences, version 4.5, 2000. Accelrys, San Diego, CA, USA) [30, 31]. Clinically evaluated stages and visual acuities averaged over both eyes of a proband and averaged over probands carrying the same BEST1 variant were correlated by linear regression with the computed subunit dimerization and calcium binding energies (ΔΔE_dim, ΔΔE_Cabin) as well as the modulus (absolute value) of dimerization energy, |ΔΔE_dim|. Due to the lack of other functional parameters such as the Arden Index or Mf-ERG for all the patients other correlations were not possible. A set of 20 distinct BEST1 variants in 27 probands, for whom BCVA data was available (Table 1), was entered in the QSPR model. Statistical techniques for validation of the regression (leave-one-out cross-validation) were used to identify outliers (data points modelled poorly by the regression equation). In the Cerius², an outlier is defined as a data point with a residual more than twice the standard deviation of the residuals generated in the validation procedure.

Table 1 Demographic, genetic and clinical characteristics of probands

Results

Genotype

The subjects of the study were 30 males and 27 females (mean age ± SD 40.7 ± 18.6 years; range 9-85) with clinical features suggesting Best’s disease. Genetic testing revealed 26 different variants in 36 patients who tested positive (63% positivity); 20 variants are already known to be associated with Best disease, and six were novel and associated with BVMD for the first time (Table 1). One of the six new variants regards a nucleotide substitution that causes an amino acid change already established as pathogenic, and we therefore do not discuss it here. We describe the clinical features of probands carrying the other five new variants in Additional files 1, 2.

Some family members of 13 patients testing positive for a variant in BEST1 were also studied by target sequencing: 16 out of 18 carried a variant in BEST1, eight of whom were clinically healthy and eight affected (Fig. 3). Penetrance in our families was estimated at 72.4% with no sex differences (71.4% and 73.3% for females and males, respectively).

Fig. 3

Pedigrees of 13 probands with Best disease. P, proband; E−, negative to genetic test; square, male subject; circle, female subject; black symbol, affected subject; white symbol, healthy subject; * documented clinical evaluation; ● obligate carrier; yr, age in years; maternal and paternal alleles are identified by square brackets; =, indicates no change in the second protein allele

Extension of study to family members was only possible for 2 out of 6 new variants: p.(Val9Gly) found in the asymptomatic father of proband 3, and p.(Asn296Lys) found in the affected father and grandmother of proband 31 (Fig. 3).

Pathogenic BEST1 variants were scattered throughout the gene (Fig. 4).

Fig. 4

Distribution of BEST1 variants found in our patients. Graphical representation of the distribution of BEST1 variants throughout the gene. New variants are shown in rectangular boxes

Structural biology of BEST1 variants

A variant can affect the structure and function of BEST1 CaCC in a number of ways. Depending on its site in bestrophin-1 molecular structure, an amino acid replacement can lead to: (i) altered Ca²⁺ binding and defective channel activation; (ii) conformational changes affecting channel gating; (iii) anomalous binding between bestrophin-1 subunits and modified channel stability; (iv) variations in channel pore size and shape and altered ion throughput rate or specificity; and (v) other changes (Fig. 5). We studied structural changes linked to the 36 variants (Table 2) in a 3D model of human bestrophin-1 prepared by comparative modelling from the crystal structure of chicken BEST1 [9, 10]. The locations of the amino acid variants in the 3D structure of bestrophin-1 of the probands are shown in Fig. 6.

Fig. 5

(figure adapted from: Ref. [10])

Subunits C and E of the X-ray structure of chicken bestrophin-1 (in ribbon representation) forming part of the calcium-activated chloride channel [9, 10]. The axial cross-section of the ion pore (yellow tube) shows the approximate size and shape of the pore and indicates the positions of the pore neck and aperture as well as calcium binding sites involved in channel gating. Boundaries of the cell membrane are indicated by the horizontal red lines

Table 2 Computational biology characterization of detected variants

Fig. 6

3D model of subunit A of human bestrophin-1 (residues 1-366; CPK representation) prepared by homology modelling and relaxed by geometrical optimization in water by molecular mechanics from crystal structure of chicken BEST1 [9, 10]. Ca²⁺ and K⁺ ions of neighbouring subunit B in the chloride channel coordinated by the C-terminal domain residues of subunit A are marked with (B). The locations of individual variant amino acids identified in the probands (Table 1) are shown in purple. a Front view of subunit A with “tail” wrapped around subunit B (channel pore inner side). b Side view (binding site to neighbouring subunit E). c Rear view (membrane side). d Side view (binding site to neighbouring subunit B)

Besides locating variant residues, we used molecular mechanics to calculate the relative energies of dimerization (ΔΔE_dim, Eqs. 1 and 2) of two bestrophin-1 subunits bearing a variant. The dimerization energies approximate the probability that assembled bestrophin-1 variants display an alteration of the ideal quaternary structure of CaCC affecting its function. We also computed relative energies of calcium binding (ΔΔE_Cabin, Eqs. 3 and 4) by the variant bestrophin-1 subunits as an approximation of the possibility of altered calcium-dependent regulation of CaCC for all the variants detected in the probands (Table 2). The relative energies therefore characterize deviation from native bestrophin-1 structure and/or function, presumably reflecting potential pathogenicity.

To explore the possibility of a link between the computed energies (ΔΔE_dim, ΔΔE_Cabin and |ΔΔE_dim|) and clinical symptoms of probands with the BEST1 variants in question, we established a QSPR [29]. In the QSPR model we used the age-adjusted mean Best’s Disease Severity Index (BDSI), averaged over all N probands with the same variant, to represent the severity of observed clinical symptoms on a quantitative scale (0% to 100%):

$${\text{BDSI}} = \frac{1}{\text{N}}\mathop \sum \limits_{\text{i}}^{\text{N}} \left[ {\frac{{{\text{logMAR}}_{{{\text{LE}},{\text{i}}}} + {\text{logMAR}}_{{{\text{RE}},{\text{i}}}} }}{2.6} \cdot 100{\text{\% }} \cdot {\text{e}}^{{ - \frac{{{\text{A}}_{\text{i}} }}{100}}} } \right]$$

(5)

The summation in Eq. (5) includes probands with the same variant, where logMAR_LE,i and logMAR_RE,i are left (LE) and right eye (RE) best-corrected visual acuities (BCVA) of probands in the log Minimum Angle of Resolution scale (logMAR—Table 1) determined by clinical evaluation. In Eq. (5), the sum of BCVA was normalized by a factor of 2.6 to express loss of vision as a percentage (logMAR = 0/0 means perfect vision of LE/RE, while logMAR ≥ 1.3/1.3 means complete blindness of both eyes). In Eq. (5), A_i is the age of proband i at the time of clinical evaluation and the corresponding time factor of the BDSI: ${\text{f}} = {\text{e}}^{{ - \frac{{{\text{A}}_{\text{i}} }}{100}}}$ adjusts the severity index for age of disease onset (${\text{f}} \cong 1$ for a patient diagnosed at age ${\text{A}}_{\text{i}} = 1$ year while ${\text{f}} \cong 0.368$ for a patient diagnosed with BVMD at age ${\text{A}}_{\text{i}} = 100$ years). BDSI = 0% means perfect vision (no symptoms), BDSI = 100% means complete loss of vision. See legend of Table 2 for examples of BDSI index values in relation to BCVA and age of proband. The BDSI was analyzed for correlations with computed energies (ΔΔE_dim, ΔΔE_Cabin) and with |ΔΔE_dim|. The modulus of dimerization energy was used in the QSPR models because we assumed that elevated and low relative energies of bestrophin-1 subunit dimerization could both be harmful to the normal structure and function of the CaCC.

We managed to establish a statistically significant QSPR model correlating BDSI with the modulus of computed relative dimerization energy:

$${\text{BDSI}} = 0.24935 \cdot \left| {\Delta \Delta {\text{E}}_{ \dim } } \right| + 6.56527$$

(6)

(for details see Fig. 7) and linking structural biology considerations to the clinical findings. This validated QSPR model also enables prediction of clinical symptoms (predicted mean visual acuity averaged over both eyes at a given patient age, PMVA or logMAR_ave) of other similar BEST1 variants and estimation of the pathogenicity of BEST1 variants based on |ΔΔE_dim| computed for a given variant by molecular mechanics:

$${\text{PMVA}} = {\text{logMAR}}_{\text{ave}} = \frac{2.6}{200} \cdot \left[ {0.24935 \cdot \left| {\Delta \Delta {\text{E}}_{ \dim } } \right| + 6.56527} \right] \cdot {\text{e}}^{{ \frac{\text{A}}{100}}}$$

(7)

where A is the age of the patient (set at 40 years in our predictions, Table 2). In the logMAR scale (0 to 1.3 or more) higher values mean worse visual function, so that BEST1 variants leading to predicted logMAR_ave values higher than 0.5 (corresponding to |ΔΔE_dim| ≥ 77 kcal mol⁻¹) can be classified as potentially pathogenic. Calculated logMAR_ave values lower than 0.5 are regarded as less likely to be pathogenic (Table 2). The threshold 0.5 is based on the WHO criteria for low vision using the logMAR scale [31].

Fig. 7

Plot of regression equation of Best Disease Severity Index (BDSI) versus modulus of relative energy of dimerization |ΔΔE_dim| of variant BEST1 subunits: BDSI = 0.24935 · |ΔΔE_dim| + 6.56527, Eq. (6). BDSI was derived from clinically evaluated BCVA (best-corrected visual acuity in logMAR scale) and age of 27 probands, Eq. (5), Table 2; ΔΔE_dim was computed by molecular mechanics from the 3D model of human CaCC. Statistical parameters of the correlation: number of available data points: n = 20, squared correlation coefficient: R² = 0.854, leave-one-out cross-validated correlation coefficient: ${\text{R}}^{ 2}_{\text{xv}}$ = 0.827, Fischer F-test: F = 76.164, significance of correlation: α > 95%, number of outliers removed: n_o = 5

Finally, out of the six new BEST1 variants reported here, four residue replacements [p.(Val9Gly), p.(Asn179Asp), p.(Glu292Gln) and p.(Asn296Lys)] seem to affect Ca²⁺ ion binding and CaCC activation (Table 2).

Discussion

Computational modelling of the tertiary structure of bestrophin-1 and quaternary structure of CaCC made it possible to link the detected variants to possible molecular mechanisms that may be involved in channel function impairment. In relation to the location of individual amino acid replacements in the 3D structure of bestrophin-1, the variants may affect channel shape, stability, activation, gating, ion selectivity, ion throughput, and probably other features as well. Our computational structural biology approach was therefore based on characteristics involving protein structure and function which are more complex than similarity/diversity of variant amino acid residues.

It is complicated to quantify the clinical features of probands with Best vitelliform macular dystrophy. Best disease stages of the left and right eye as diagnosed in our probands were not correlated with BCVA in the logMAR scale (Table 1). We, therefore, selected BCVA and designed our own age-adjusted quantitative Best’s disease severity index (BDSI), see Eq. 5. The BDSI index showed a good correlation with the computed modulus of bestrophin-1 subunit dimerization energy |ΔΔE_dim| (Fig. 7) which is assumed to reflect the potential harmful effect of individual missense variants on BEST1 CaCC structure and function. During regression analysis, it was necessary to remove five outlier points (probands P13, P24, P25, P27 and P31) identified by the leave-one-out cross-validation algorithm. The corresponding point variants (p.(Tyr29Cys), p.(Arg218Ser), p.(Ile232Asn), p.(Ala243Thr) and p.(Glu292Gln) are mainly residue replacements leading to low-to-medium positive and negative ΔΔE_dim values, which are thus expected to be linked to small changes in CaCC structure and function. We can therefore assume that the QSPR model and the proposed computational tool for BEST1 variant pathogenicity prediction apply to conservative as well as radical variants. Although many intrinsic factors are known to co-determine the biological consequences of gene diversity and the resulting impairment of health or well-being [32], we were able to establish a robust but simple QSPR model linking computed relative energies of human bestrophin-1 protein with the symptoms of BVMD diagnosed in a cohort of patients. To illustrate application of the model we compare the PMVA of two probands, P4 and P17, of similar age bearing different BEST1 variants and displaying different clinical symptoms (Table 3). As we can see, patient P17 has greater visual impairment (higher logMAR_ave and mean BDSI values). Structural biology assessment of his BEST1 variant predicts larger alterations of stability and structure of assembled bestrophin-1 subunits (larger |ΔΔE_dim|), possibly detrimental to CaCC function. Consequently, predicted mean visual acuity (PMVA) estimated for P17 at age 40 years is larger than the threshold of 0.5 (PMVA = 0.59) and the corresponding BEST1 variant p.(Ser108Arg) can thus be considered likely pathogenic.

Table 3 Comparison of two BEST1 variants

Further calibration of this computational pathogenicity estimate on a much larger cohort of patients and wider range of amino acid variants is still needed to set the borderline between low and high pathogenicity and confirm the reliability of the resulting clinical outcome predictions.

Conclusions

Similar computational biology approaches based on molecular mechanics energies or other structural properties of native and variant proteins and their interactions suggest the possibility of exploring the molecular basis of genetic diseases and perhaps predicting likely clinical outcomes of genetic variants. Worth to note, the predicted pathogenicity descriptor PMVA can also be used for estimating the development of patient’s BCVA over the years by substituting selected age (A) into Eq. (7).

We also reported the results of genetic testing in 57 Italian BVMD patients, including 36 subjects with a BEST1 variant. Compared with other reports in the literature for different populations, our results indicate that there are no specific ethnic variants [2], while segregation study confirmed the variable penetrance and expressivity of the disease in patients with the same BEST1 variant, even in the same family [33, 34]. Six new variants are also reported, thus broadening the BEST1 variant spectrum.