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

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 2+ 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 2+ clasps) and three binding sites for Cl − for each subunit in the channel pore [9,10].
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 2+ clasps [9,10], Fig. 2. 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.

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.

Fig. 2
Crystal structure of the "Ca 2+ 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 2+ 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 For the identification of variants, the Retina International database (http://www.retin a-inter natio nal.org), the Human Gene Mutation Database Professional version 2017.2 (https ://porta l.bioba se-inter natio nal.com/hgmd/ pro/), the Exome Variant Server (http://evs.gs.washi ngton .edu/EVS/) and dbSNP (https ://www.ncbi.nlm.nih. gov/snp/) were also consulted.

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 (Uni-ProtKB 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/embos s_needl e/). 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 2+ 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 energyminimization 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 2+ 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]: 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 2+ 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: where E tot {Z} aq and E sol {Ca 2+ } 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 2+ ion.
(1)  [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 2 , 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.

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).
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).

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 2+ 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.
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%): The summation in Eq.      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: (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: 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 −1 ) 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].  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. 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. 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 2 = 0.854, leave-one-out cross-validated correlation coefficient: R 2 xv = 0.827, Fischer F-test: F = 76.164, significance of correlation: α > 95%, number of outliers removed: n o = 5