The structural determinants of the Asp-rich domain were modeled by using the X-ray structure of TSP-1 T3₆ repeat [36] as a template, because it exhibits the largest number of Ca²⁺-binding residues conserved across bestrophins family (Figure 1). The number of Ca²⁺ ions bound to the Asp-rich domain is unknown. Typically, as observed from other Ca²⁺-binding domains in the PDB database, each Ca²⁺ ion is usually coordinated by three Asp/Glu residues [38]. Since the template has five Ca²⁺ ions bound and 33% of negatively charged residues more than the Asp-rich domain, it is plausible to expect that the last would bind five or less Ca²⁺ ions. Therefore we considered models including five (M5), four (M4), and three (M3′, M3″) Ca²⁺ ions (Figure 2).

In all molecular dynamics (MD) simulations, two ions (Ca₁ and Ca₂ in Figure 3) remained bound. Ca₁ was bound to D302, D304, E306 and D312, Ca₂ to D303, D304, E306, D312, Q316, S318 and D323 (Table S1). Instead, the other ions were partially or completely solvated. The decrease in the number of bound Ca²⁺ ions with respect to TSP-1 T3₆ can be caused by a combination of factors: first of all, the lower number of negatively charged residues in Asp-rich domain with respect to T3₆ (Figure 1), but also the limitations associated with the force field (see Text S1).

MD simulations showed that residues D302, D304, E306, D312 and D323 (Table S1) play a major role for Ca²⁺ binding. In order to identify other Ca²⁺-binding residues, we used an enhanced sampling method, i.e. metadynamics [39]. In this approach, the MD simulation is altered by biasing the selected collective variables (CVs) with a history-dependent potential. The bias, constructed as a sum of Gaussians centered along the trajectory of the system, reduces the free-energy barriers and enhances the rate of transitions, in the present case the binding and unbinding of Ca²⁺ ions.

The starting structure for the metadynamics simulation was selected among our energy-minimized models. We chose M3′ model, because it kept the largest number of bound Ca²⁺ ions (three) during the MD simulations for as long as 25 ns (Figure 3). The selected CVs measured the extent of binding of Ca₂ and Ca₄, respectively, since Ca₁ was always bound to a well-defined binding pocket in all previous simulations. The limitations of the force field along with the large number of conformational degrees of freedom of the system prevented us from calculating the free-energy profile associated to Ca²⁺ binding. Therefore, the simulation was only used to sample stable configurations, different from those studied in the MD.

Throughout the metadynamics run, Ca₁ kept its coordination shell as in the previous MD simulations (Table S1 and Figure S1,A). Due to the bias potential, Ca₂ and Ca₄ instead turned out to bind and unbind several times (Figure S1,B) exploring different configurations (Figure S2). Ca₂ was either bound to D301, D302, D303, D312 and D323 (Figure 4) or it was solvated. Analogously, the third Ca²⁺ ion (Ca₄ in Figure 4) was bound to D301, D302, D303, D304, E306, D312 and D323 in some of the stable conformations. Therefore, our metadynamics simulations confirmed that the residues identified by MD can play a role for Ca²⁺ binding. In addition, they further suggested a possible role of D301 and D303 for such binding.

We attempted to identify mutations that could have an effect on Ca²⁺ binding. We constructed structural models in which putative key Asp and Glu residues among those identified in the above section were mutated to alanine. These are D301A, D302A, D303A, D304A, E306A, D312A and the double mutant E306A+D323A (see Table S2). Adducts with five Ca²⁺ ions underwent 20 ns of MD.

At the end of MD simulations, D301A bound only one Ca²⁺ ion (Table S2), while its 3-D structure remained similar to that of WT.

The Asp-rich domain bearing D304A mutation unfolded, possibly because D304 is the central residue of the domain and therefore responsible for its structural stability; all Ca²⁺ ions were partially or fully solvated (Table S2). Since the structural information of the channel is lacking, it is difficult to establish if the unfolding of D304A would occur in the real molecular surrounding, that is, when the Asp-rich domain is attached to the remaining bestrophin part.

Mutants D302A, D303A, E306A, and D312A resembled the WT as they bound at least two Ca²⁺ ions (Table S2). In fact, other residues contributed to Ca²⁺ coordination replacing the mutated ones.

Based on these results we concluded that D301A and D304A mutations have the largest effect on Ca²⁺ binding, while D302A, D303A, E306A, and D312A affect Ca²⁺ binding to a lesser extent.

In order to verify this, the experiments of the WT and mutant mBest2 were performed. The putative Ca²⁺-binding domain of mBest2 has almost the same sequence as that of hBest1 (Figure 1). To investigate whether the WT structures are affected by the differences, MD simulations were performed also on mBest2. As a result of MD simulations both models showed very similar behavior, with Ca₁ and Ca₂ that remained stably bound, while Ca₃ and Ca₅ were solvated. Ca₄ bound more peptide O-donors in mBest2 than in hBest1. However, this difference should not be overestimated, since even starting from the same initial system (e.g. hBest1), running two independent MD simulations of the same length could result in slightly different events.

The changes of Ca²⁺ coordination in the aforementioned mutations may affect the amplitude of bestrophin Ca²⁺-activated Cl⁻ currents, which can be experimentally measured. This was done here for mBest2. WT and mutants currents were compared after their heterologous expression in the HEK-293 cell line (Figure 5).

To determine if mutant proteins are delivered to the plasma membrane as WT, we used immunostaining with antibodies against mBest2 [11]. The WT protein showed increased localization in plasma membrane and/or in close sub-membrane space in a significant portion of transfected cells, although some fraction of protein was localized also intracellulary (Figure 5A). These results are similar to data previously published for hBest2 [8] and xBest2 [40]. The localization of the proteins was further investigated by co-staining with fluorescently labeled wheat germ agglutinin (WGA), a selective marker of glycosylated surface-expressed proteins. Co-localization of both signals, which indicates membrane expression of the protein, was observed. The distribution of mBest2-related signal and its co-localization with WGA in mutants were very similar to WT, indicating that the proteins were expressed and delivered to the plasma membrane (Figure 5). Consequently, changes of the measured current were caused by a modification of the channel functionality and not by the inability of the mutant proteins to reach the plasma membrane.

We then measured currents in HEK-293 cells transfected with mBest2 WT or mutants both in nominally 0 free Ca²⁺ level and in 22 µM free Ca²⁺ intracellular concentration. Currents were measured in the whole-cell voltage-clamp configuration applying voltage-steps of 20 mV from −100 mV to +100 mV from a holding potential of 0 mV. The current values at +60 mV were normalized to the capacitance of each cell and averaged.

The WT current in nominally 0 Ca²⁺ was significantly smaller than that in 22 µM Ca²⁺ concentration (Figures 6–7). This shows that Ca²⁺ ions activated a current in agreement with previous results [5], [11]. Most of the Ca²⁺-induced current was reversibly blocked by the Cl⁻ channel blocker SITS indicating the current is indeed carried by Cl⁻.

In contrast, D304A and D301A produced negligibly small currents in both Ca²⁺ levels (Figure 7). Two other mutants, E306A and D312A, produced currents: whilst the amplitude of the current of E306A was similar to that of WT (Figure 7), D312A showed a somehow unexpected result. It produced a substantial amount of current in the nominally 0 Ca²⁺ concentration (Figure 7). In 22 µM Ca²⁺ concentration, the current significantly increased and was larger than that of WT. The current was also blocked by SITS in both Ca²⁺ levels, indicating that it is a Cl⁻ current (data not shown). Therefore, the D312A mutation not only allowed a Ca²⁺-activated Cl⁻ current, but it also modified channel gating increasing its current also in the absence of Ca²⁺.

Finally, the double mutant D304A+E306A abolished the current, similarly to the single mutant D304A (Figure 7).

Experiments showed that WT, E306A and D312A mutants allowed a substantial Ca²⁺-activated current, whereas D301A and D304A mutants produced a negligibly small Ca²⁺-activated current. These results were fully consistent with our theoretical predictions.

The sequences of the human TSP-1 T3_5-7 repeats (PDB code: 1UX6) [36] were aligned with the Ca²⁺ bowl from human and mouse BK_Ca channels [26]–[30], [41], [42] and the Asp-rich domains of hBest1 [1], [2], mBest1 [2], hBest2 [43], mBest2 [44], hBest3 [43], mBest3 [45], hBest4 [43], ceBest1 [4], dvBest1 [46], dmBest2 [47] and dmBest4 [47] (Swiss-Prot entries and notations are reported in the legend of Figure 1). The alignment was initially performed with ClustalW [48] and then slightly refined to avoid apolar residues in Ca²⁺ binding positions. Sequence similarity was determined considering amino acid pairs having positive values in BLOSUM62 [49] substitution matrices. The homology models of hBest1 and mBest2 Asp-rich domains from W287 to D323 were built with the Modeller 6v2 program [50], [51], using the crystal structure of the TSP-1 T3₆ motif as a template (PDB code: 1UX6) [36]. In the case of hBest1, models were generated where five to three Ca²⁺ ions present in the template were kept. These are: (i) M5, with five Ca²⁺ ions coordinated. (ii) M4 is composed by four Ca²⁺ ions which are Ca₁, Ca₂, Ca₄ and Ca₅. (iii) M3′ and M3″ are including three Ca²⁺ ions, Ca₁, Ca₂, Ca₄ or Ca₁, Ca₂, Ca₅, respectively.

Models of the Asp-rich domain of hBest1 were inserted in a box of edges 49, 50 and 53 Å filled with ∼4,100 water molecules [52]. Cl⁻ (three in M5, one in M4) or Na⁺ (one in M3′ and M3″) counterions were added to ensure systems charge neutrality. The force fields for the protein, ions and water molecules were the Amber parm99 [53], [54] and TIP3P [52] force fields, respectively.

For M5, energy minimization of water and counterions was performed with position restraints on the Asp-rich domain and Ca²⁺ ions, with a force constant of value 25.0 kcal/mol Å⁻². After the first 300 cycles of steepest descent minimization, conjugate gradient method was used until 50,000 cycles. Then steepest descent minimization of the system without restraints was performed for another 150 cycles, following by conjugate gradient minimization until 50,000 cycles. After that, constrained solute and Ca²⁺ ions underwent 60 ps of linear heating MD from 0 K till 300 K. Restraints on Asp-rich domain and Ca²⁺ ions corresponded to a force constant of 5.0 kcal/mol Å⁻². A time step of 2 fs was used. Finally, the system underwent NPT MD simulations for 75 ps at 1 fs time step. The system was maintained at the reference pressure and temperature of 1 bar and 300 K, respectively, by coupling the system to a Berendsen thermostat and barostat [55]. A cut-off of 9 Å was used for calculating non-bonded interactions.

Finally, unrestrained MD simulations were carried out for an overall length of 70 ns. A time step of 2 fs was used. All bond lengths were kept fixed applying the LINCS algorithm [56]. Temperature (300 K) and pressure (1 bar) were kept constant as aforementioned. Periodic boundary conditions were applied treating long-range electrostatic interactions with the particle-mesh Ewald technique using a short-range cutoff of 10 Å [57]. The same cutoff was used for the van der Waals interactions. Pair lists were updated every 20 steps.

Models with four or three Ca²⁺ ions were built starting from the MD snapshot of M5 at 9 ns (after the initial relaxation) by manually removing one or two ions (see above). In M4, M3′ and M3″ models Ca₃ was always removed since it is immediately solvated during the simulation of M5 (Figure 3). M4, M3′ and M3″ underwent then MD simulations for 20 ns (M4) and 35 ns (M3′, M3″), respectively, using exactly the same protocol as that used for the M5 model.

The models of the Asp-rich domain mutants (D301A, D302A, D303A, D304A, E306A, D312A, E306A+D323A) were built, based on the M5 model, using the program Molden [58]. They underwent the same protocol of minimization and relaxation of the WT followed by 20-ns MD runs.

The model of mBest2 Asp-rich domain underwent the same protocol of minimization and relaxation of hBest1, followed by a 60-ns MD run.

All the systems underwent minimization using Amber 8.0 program [59] and then MD simulations using the Gromacs 3.2.1 program [60], [61]. Trajectory analysis was performed by using the VMD 1.8.6 program [62].

A 60-ns long metadynamics simulation was performed by biasing two collective variables (CVs) with a history-dependent potential. These are the coordination numbers of Ca₂ and Ca₄ to Asp and Glu residues of the Asp-rich domain. They were defined as the number of contacts shorter than 4.3 Å among one Ca²⁺ ion and all carboxylate carbons of the studied domain:where i runs over Ca₂ and Ca₄, j runs over all carboxylate carbons, R_ij is the distance between the two atoms, R ₀ = 4.3 Å is a cutoff distance chosen in order to discriminate between bonded and non-bonded Ca²⁺ ions (this cutoff distance ensures that Ca²⁺ ion is considered bonded when it is at most 3 Å away from carboxylate oxygen atoms (see e.g. [38]), and n_ij is a function switching smoothly between 1 (bonded Ca²⁺ ion) and 0 (non-bonded Ca²⁺ ion). Similar CVs have been employed in a previous study of Ca²⁺-binding polymers [63]. The 2-D Gaussians forming the metadynamics bias potential have width σ₁ = σ₂ = 0.2 (along S₁ and S₂ respectively) and a height of 0.2 kcal/mol. A new Gaussian is added along the trajectory every 1 ps. Metadynamics is implemented in a modified version of the Gromacs-3.3.1 package [39]. The code has been kindly provided by Alessandro Laio.

The clone of mBest2 in pCMV-Sport6 mammalian expression plasmid was obtained from the RZPD (Berlin, Germany) collection (Deutsches Ressourcenzentrum für Genomforschung). Site-specific mutations of mBest2 were made using a PCR-based site-directed mutagenesis kit (Quiagen). Mutations were confirmed by DNA sequencing.

mBest2 and its mutants were transfected into HEK-293 cells along with a vector that expressed EGFP (pEGFP; Invitrogen) at a 1∶8 ratio using Fugene-6 transfection reagent (Roche). 2 days after transfection, cells were dissociated and replated. Transfected cells were identified by EGFP fluorescence and used for patch clamp experiments within 3 days after transfection.

To test successful expression of proteins, as well as their spatial distribution in HEK-293 cells, we used co-staining of antibodies against mBest2 and labeled wheat germ agglutinin (WGA) as a marker of glycosylated surface-expressed proteins – a method previously used to show localization of proteins in cell membrane and/or close submembrane space [40], [64].

Transfected cells (without the EGFP expressing vector) were fixed with 4% paraphormaldehyde for 20 min, washed, incubated with 5 µg/ml WGA conjugated to Alexa Fluor 488 (Invitrogen) for 20–30 min, incubated for 30 min in blocking solution containing 1% BSA (Sigma), 0.1% Triton X-100 (Sigma) and 1% FCS (Sigma), then for 2 hours with an antibody against mBest2 [11] (diluted in 1∶500 ratio with blocking solution), washed extensively and incubated with a A594-conjugated anti-rabbit secondary antibody (Invitrogen, diluted 1∶500 in blocking solution) for 1 hour, washed and mounted with Vectashield (Vector Laboratories, Burlingame, CA). Immunoreactivity was analyzed by confocal microscopy with a Leica TCS SP2 confocal microscope. All incubations were made at room temperature; all substances were diluted in PBS, unless otherwise indicated.

Currents were measured with an Axopatch 200B patch-clamp amplifier (Molecular Devices, CA, USA) in the whole-cell voltage-clamp mode. The amplifier was controlled via a Digidata 1440A (Molecular Devices, CA, USA).

Patch pipettes were made using borosilicate capillaries (WPI, Sarasota, Florida, USA) and pulled with a Narishige PP83 puller (Narishige, Tokyo, Japan), using a double stage pull. The pipette resistance was 2–4 MΩ. Data were sampled at 20 kHz and low-pass filtered at 10 kHz. Acquisition and storage of data were performed with the PClamp 10 software (Axon Instruments, CA, USA). Data analysis and figures were made with the Clampfit 10 software (Axon Instruments, CA, USA) and Igor 6.0 software (Wavemetrics, Lake Oswego, OR, USA).

The holding potential was 0 mV and a voltage protocol from −100 mV to +100 mV, with 20 mV steps of 200 ms duration was applied (starting 2 minutes after reaching the whole-cell configuration). Cell capacitance was estimated from a 10 mV step test protocol in the presence of the extracellular Cl⁻ current blocker SITS. Current densities were calculated normalizing the current recorded at +60 mV to the capacitance of each cell. Data are shown as mean±standard error of the mean. N is the number of cells. Statistical significance was determined using un-paired t-tests. P values<0.05 were considered significant.

The standard extracellular solution contained (in mM): 140 NaCl, 5 KCl, 2 CaCl₂, 1 MgCl₂, and 10 HEPES, 10 Glucose, pH = 7.4. The pipette solution contained (in mM): 130 CsCl, 2 MgCl₂, 10 HEPES, 5 EGTA without added Ca²⁺ for the nominally 0 Ca²⁺ solution, or with 5 mM CaCl₂ for the 22 µM Ca²⁺ free pipette solutions (free Ca²⁺ was calculated with the program WinMAXC, C. Patton, Stanford University, Palo Alto, CA, USA). pH was adjusted to 7.3 by adding HCl or CsOH. The osmolarity was adjusted with sucrose to 310 or 300 mOsm for the extracellular and intracellular solution, respectively.

To block Cl⁻ currents, the extracellular blocker 4-acetamido-4′-isothiocyanato-stilben-2,2′-disulfonate (SITS) was directly dissolved in the extracellular solution at 2 mM. All chemicals were purchased from Sigma. All experiments were carried out at room temperature (20–22°C).