Molecular evolution and functional divergence of the bestrophin protein family

Abstract

Background

Mutations in human bestrophin 1 are associated with at least three autosomal-dominant macular dystrophies including Best disease, adult onset vitelliform macular dystrophy and autosomal dominant vitreo-retinochoroidopathy. The protein is integral to the membrane and is likely involved in Ca²⁺-dependent transport of chloride ions across cellular membranes. Bestrophin 1 together with its three homologues forms a phylogenetically highly conserved family of proteins.

Results

A bioinformatics study was performed to investigate the phylogenetic relationship among the bestrophin family members and to statistically evaluate sequence conservation and functional divergence. Phylogenetic tree assembly with all available eukaryotic bestrophin sequences suggests gene duplication events in the lineage leading to the vertebrates. A common N-terminal topology which includes four highly conserved transmembrane domains is shared by the members of the four paralogous groups of vertebrate bestrophins and has been constrained by purifying selection. Pairwise comparison shows that altered functional constraints have occurred at specific amino acid positions after phylogenetic diversification of the paralogues. Most notably, significant functional divergence was found between bestrophin 4 and the other family members, as well as between bestrophin 2 and bestrophin 3. Site-specific profiles were established by posterior probability analysis revealing significantly divergent clusters mainly in two hydrophilic loops and a region immediately adjacent to the last predicted transmembrane domain. Strikingly, codons 279 and 347 of human bestrophin 4 reveal high divergence when compared to the paralogous positions strongly indicating the functional importance of these residues for the bestrophin 4 protein. None of the functionally divergent amino acids were found to reside within obvious sequences patterns or motifs.

Conclusion

Our study highlights the molecular evolution of the bestrophin family of transmembrane proteins and indicates amino acid residues likely relevant for distinct functional properties of the paralogues. These findings may provide a starting point for further experimental verifications.

Background

The bestrophins are a phylogenetically conserved family of integral membrane proteins initially identified in Caenorhabditis elegans [1]. Homologous sequences are found in animals, fungi, and prokaryotes, but not in protozoans or plants [2]. Conservation is mainly restricted to the N-terminal 350–400 amino acids with an invariant motif arginine-phenylalanine-proline (RFP) of unknown functional properties.

The first human bestrophin cloned, bestrophin 1, is encoded by the vitelliforme macular dystrophy type 2 (VMD2) gene on chromosome 11q13 and was shown to be associated with Best macular dystrophy (BMD, OMIM #153700) also known as Best disease [3, 4]. Subsequently, mutations in this gene were also found to cause adult onset vitelliform macular dystrophy (AVMD, OMIM #608161) and autosomal dominant vitreo-retinochoroidopathy (ADVIRC, OMIM #193220). The currently known 106 disease-causing mutations [5] are associated with a dominant pattern of inheritance with the majority being missense mutations located in four clusters near the RFP motif and the predicted transmembrane domains (TMDs) [6].

Based on immunocytochemical studies in macaque, porcine and human eyes, bestrophin 1 was shown to localize to the basolateral plasma membrane of the retinal pigment epithelium (RPE) [7, 8]. In addition, bestrophin 1 is broadly expressed in other epithelia including the intestine and the lung [9], whereas bestrophin 2 expression appears confined to the olfactory epithelium [10]. On the functional level, whole cell patch clamp experiments suggest that the bestrophins act as Ca²⁺-dependent transporters of chloride ions across epithelial borders [11–16]. Functional proteins likely exist as multimeric complexes eventually by forming homo- or heteromeric associations between bestrophin family members.

Uncertainty exists as to the specific function of bestrophin 1. As an alternative to its activity as a chloride channel, bestrophin 1 may act as an accessory protein in the regulation of voltage-gated calcium channels [17]. In addition, bestrophin 1 could be involved in the volume sensitivity of RPE cells during phagocytosis of photoreceptor outer segments [11]. Consequently, distinct functional aspects of bestrophin 1 may contribute to RPE dysfunction and thus may explain the variable phenotypes associated with a mutant protein.

Phylogenetic studies of protein families can be a valuable tool to determine conserved but also divergent regions, potentially leading to functional predictions [18]. In this study we elucidated the evolutionary history of the bestrophins and identified structural and putative functional motifs of the bestrophin protein family by a comprehensive bioinformatics/phylogenetic approach. This has led us to predict distinct amino acid residues that may be of importance in the functional divergence of the bestrophin paralogues.

Results

Phylogenetic analysis of the bestrophin family

We first retrieved the available bestrophin sequences from the currently sequenced genomes. Querying major databases and unfinished genomes with the full-length amino acid sequences from the four human bestrophin paralogues identified 173 homologous proteins in vertebrates (mammals, birds, amphibians and fishes), urochordates (sea squirt), and invertebrates (insects, nematodes) (see Additional file 1). While C. elegans reveals a large number of bestrophins (n = 26), most other organisms harbour three or four family members. In each of the recently sequenced genomes of the urochordata Ciona inestinalis and Ciona savigny, the closest relatives of the craniates, we identified only a single bestrophin sequence strongly suggesting that gene duplication events may have occurred in the lineage leading to the vertebrates.

Our further study focused on vertebrate bestrophins. After exclusion of unfinished and partial protein sequences, a phylogenetic tree with 53 selected bestrophins (Table 1) was constructed by the neighbour-joining (NJ) method based on a gamma corrected Jones-Taylor-Thornton (JTT) distance matrix [19, 20] (Figure 1). Accordingly, the vertebrate bestrophins can unambiguously be separated into distinct clusters. With the urochordata as outgroup, the vertebrate bestrophins form four separate monophyletic groups, all with high bootstrap support indicating that the formation of the paralogous subfamilies occurred before the divergence of individual species (Figure 1). The phylogenetic branches of bestrophins 1 and 3 separated considerably earlier in evolution than bestrophin subgroups 2 and 4, partially explaining the level of sequence conservation between each of the two subfamilies. The high level of sequence identity within a subfamily suggests evolutionarily conserved functions. Overall, the data indicate that the bestrophin subfamilies have evolved by gene duplication events, in good agreement with previous findings demonstrating that large-scale gene duplications have occurred during chordate evolution [21, 22].

Table 1 Vertebrate bestrophin homologues (n = 53) used for phylogenetic analysis

Figure 1

Phylogenetic tree of 53 homologues of vertebrate bestrophins. The phylogenetic tree of the bestrophin protein family was inferred by the neighbor-joining method (NJ) method applying gamma corrected distances. Urochordata (Ciona intestinalis [Ci] and Ciona savigny [Cs]) were used as out-group. The support for each phylogenetic group was tested with 1000 bootstrap pseudoreplicates. Protein sequences are given in Table 1.

N-terminal topology of the vertebrate bestrophins

We next analyzed whether the strong sequence conservation between the N-terminal regions of the vertebrate bestrophins also predicts a common topology. Putative transmembrane domains and hydrophobic regions were modelled in silico (Figure 2). Hydropathy plotting of 53 bestrophins (orthologues and paralogues, Table 1) revealed six hydrophobic regions including four putative transmembrane domains with hydropathy values reaching or above the threshold of 1.6 (Figure 2). These four hydrophobic peaks are strongly conserved in the vertebrate bestrophins suggesting membrane insertion of all bestrophin family members similar to bestrophin 1 [23].

Figure 2

Average hydropathy plot of 53 homologues of vertebrate bestrophins. Hydropathy plot was generated from 53 protein sequences as given in Table 1. A threshold value of 1.6 is marked by the dotted line. Positive values indicate hydrophobic regions. Putative transmembrane domains were determined by the TOPPRED II algorithm and are depicted as gray and white boxes on the top.

Variable selective pressures among amino acid sites

To analyze positive or negative selection of specific amino acid regions within the full-length protein sequences of mammalian bestrophins, substitution rate ratios of non-synonymous (dN) versus synonymous (dS) mutations (dN/dS or ω) were calculated for vertebrate bestrophins as given in Table 1. A dN/dS value <1 is indicative of purifying selection acting against amino acid changes, whereas dN/dS values >1 indicate an excess of amino acid changes suggesting adaptive evolution [24, 25]. Amino acids in a protein sequence are expected to be under different selective pressure and to have different underlying dN/dS ratios. Table 2 shows that the ratio values for the vertebrate bestrophins are substantially lower than 1 suggesting that the N-termini of mammalian bestrophins within each subgroup were under strong purifying selection pressure. In order to test for positive selection at individual amino acid codons, the site-specific models implemented in codeml were used. Likelihood rate tests (LRTs) were performed between model M0 (one ratio) and M3 (discrete), M1a (nearly neutral) and M2a (positive selection), and M7 (beta) and M8 (beta and ω) (Table 2). The selection model (M2a) does not suggest presence of positively selected sites (P = 1). M0 was rejected when compared to M3 (P < 0.001), although no positively selected sites were detected. This could be explained by the fact that the majority of the protein is subjected to constant purifying selection, while a few sites undergo positive selection [26]. From the three models which allow for selection to be tested (M2a, M3 and M8), only M8 was significantly favoured over M7 (P < 0.001) and detected 44 sites at the 95% level (Table 2). This provides evidence for the Darwinian selection in vertebrate bestrophins.

Table 2 Likelihood values and parameter estimates for the vertebrate bestrophin genes

Analysis of functional divergence

To further investigate whether amino acid substitutions in the highly conserved N-termini of the bestrophins could have caused adaptive functional diversification, amino acid residues 1–367 from 53 vertebrate bestrophins (Table 1) were used for posterior analysis by the DIVERGE program algorithms [27, 28]. Pairwise comparisons of paralogous bestrophins from subfamilies 1 to 4 were carried out and the rate of amino acid evolution at each sequence position was estimated. This identified residues potentially responsible for functional divergence. The results are given as a coefficient of functional divergence (θ) with standard errors and significance levels (Table 3). Significant functional divergence was detected between bestrophin 2 and 3 (P_corr = 6.3 × 10^-9) as well as between bestrophin 4 and the other bestrophin subgroups (P_corr ≤ 8.9 × 10^-8).

Table 3 Functional divergence estimated in 53 vertebrate bestrophin paralogues

Next, single amino acid residues responsible for the functional divergence were predicted based on site-specific profiles (Figure 3A) in combination with suitable cut-off-values derived from the posterior probability of each comparison (Table 3). Residues predicted to be functionally divergent in bestrophin, were mapped onto topology models of human bestrophins 2 and 4, respectively (Figure 3b). The predicted functional sites are not equally distributed throughout the respective bestrophin, but instead are clustered in the hydrophilic loops between predicted transmembrane domains 1 and 2 as well as 5 and 6 and immediately C-terminal to transmembrane domain 6 (Figure 3B, see Additional file 2, and Table 4). Despite the high global sequence identity of mammalian bestrophins 2 and 3, functionally divergent amino acids were also identified between these bestrophins. For example, at position 264 (numbered according to bestrophin 1) the tyrosine residue is completely conserved within bestrophin 1, 2 and 3 paralogues but highly divergent in bestrophin 4 paralogues. Similarly, residue 331 (numbered according to bestrophin 1) is highly divergent in bestrophin 1, 2 and 3 paralogues but strictly conserved in bestrophins 4 (Figure 3B, see Additional file 2, and Table 4). We thus speculate that these two sites may be of functional importance for the respective bestrophin subfamilies. In addition, 11 out of 44 positively selected sites detected by LRT in codeml were also found to be functionally divergent between bestrophin paralogues (marked by asterisk in Table 4).

Figure 3

Site specific profiles for evolutionary rate changes in the vertebrate bestrophin protein family. A, The posterior probabilities of functional divergence for vertebrate bestrophins 1 to 4 were obtained with Diverge [27]. Individual cut-off values for each comparison are marked with red horizontal lines. B, Residues with predicted functional divergence between bestrophin subfamilies are mapped onto the membrane topology models of bestrophin 4 (top) and bestrophin 2 (bottom). Divergent amino acids are shown and listed in Table 4. In the model drawing, residues excluded from the analysis are depicted in gray.

Table 4 Amino acid changes associated with the divergence of human bestrophins

We finally scanned the predicted divergent amino acids and surrounding sequences in bestrophins 2 and 4 for the presence of domains and functional sequence patterns with MotifScan and InterProScan. Identity to known protein motifs, except for several short sequence tags including putative phosphorylation sites, could not be detected.

Discussion

Our phylogenetic analysis suggests that the bestrophins originated by duplication and divergence of a common protein at the base of the eukaryotic tree, some 700 Myr ago. Analysis of the bestrophins from a phylogenetic perspective may provide the basis for understanding the functional diversity within this conserved protein family. Both small-scale and large-scale gene duplications are known to contribute to the complexity of eukaryotic organisms [29, 30]. The presence of bestrophins in prokaryotes [2] and their occurrence in phylogenetically distant eukaryotic species such as insects and mammals highlight their general functional importance.

High levels of sequence homology were generally found between the N-terminal regions of the bestrophins, while the C-termini differ substantially particularly between paralogues. At the N-terminal sequences, vertebrate bestrophins show nearly identical hydrophobicity plots with four major and two minor hydrophobic peaks, indicating that they share a highly similar membrane topology [16, 23]. The hydrophobic regions with high sequence conservation also seem to contribute to oligomerization [15, 31] and the second putative transmembrane domain is thought to be especially important for channel function [13, 14, 32].

Our data and the study by Hagen et al. [2] demonstrate that the vertebrate bestrophins cluster in organismal groups, i.e. the homologues are more related to each other from different species than to their species-specific paralogues. Whenever high quality genome sequences are available (e.g. human, chimp, dog, rat, and fugu), four bestrophin paralogues can be identified. This strongly argues for a high phylogenetic conservation of the bestrophin subfamilies. Missing paralogues in some species may be due to currently incomplete genomic sequences. One exception may be X. laevis, for which a high quality genome sequence has already been released but still this species reveals only three copies of bestrophins, interestingly all highly similar to bestrophin 2. However, as a tetraploid species, X. laevis appears to have undergone unique evolutionary changes. Another exception may be the murine bestrophin 4 (gene name: VMD2-L2) which our group has previously shown to represent a functionless pseudogene [33].

In the vertebrate phylogeny, bestrophins 1 and 3 have separated and diverged at a higher rate than bestrophins 2 and 4. To estimate the selection forces behind this we calculated the substitution rate ratios for the four subgroups. Calculating dN/dS ratios across the entire length of the N-terminal amino acid sequences (aa1-367), we could not, however, find an indication for positive selection. In contrast, the low dN/dS values imply that a strong purifying selection may lead to the high sequence conservation observed. Examples such as the chaperonins [24] show that positive selection reaches dN/dS values well above 1 for a high selective pressure to fix amino acid substitutions in evolution. Nevertheless, we detected positive selection at several amino acid sites located predominantly in the hydrophilic regions of bestrophin. Testing several models only one (M7vsM8) suggested positively selected sites with high posterior probabilities. Another comparison (M0 vs M3), although significant failed to identify positively selected residues with high posterior probability. This points to insufficient information from the alignments or, alternatively, may be explained by the fact that purifying selection in conserved regions could have masked single diversification signals.

Posterior probability analysis for pairwise comparisons of bestrophin paralogues identified significant functional divergence between bestrophin 2 and bestrophin 3 as well as between bestrophin 4 and its remaining paralogues. Generally, the specific amino acid substitutions in the bestrophin paralogues suggest an optimization of protein function, and/or subfunctionalization of the gene duplicates. Eleven out of the 44 positively selected sites were also found to be functionally divergent between bestrophins. The most prominent finding of the posterior probability analysis is that bestrophin 4 diverged significantly during evolution in a number of amino acids from its paralogues. As shown previously [34], bestrophin 4 can be activated by free Ca²⁺ on the cytoplasmic side although the kinetics of the activation/deactivation is much slower than for typical Ca²⁺-activated chloride currents [35]. Furthermore, heterologously expressed human bestrophin 4 shows very slow voltage-dependant current relaxations, in contrast to the lack of voltage-dependent current relaxations by human bestrophin 1 and 2. These distinctive properties of bestrophin 4 could be attributable to the additional 15 unique amino acids within the extracellular loop between putative transmembrane domain 5 and 6. This loop contains three positively charged lysine residues (K266, K269 and K272) within a proline rich region (P267, P273, P277 and P279). It can be speculated that the function of the divergent P279 residue could be to maintain structural or conformational stability of the loop region.

Human bestrophin 1 and 3 have considerably longer amino acid sequences compared to bestrophin 2 and 4 (585 and 668 vs 509 and 473 aa). Human and mouse bestrophin 3, heterologously expressed in HEK-293 cells, produce no chloride currents in a physiological range [16, 36], in contrast for example to mouse bestrophin 2 [13, 14]. Such a functional divergence has been explained to be mediated by an auto-inhibitory (AI) domain composed of seven critical residues ₃₅₆IPSFLGS₃₆₂ present in mouse bestrophin 3 [36]. Interestingly 2 out of 5 amino acid residues implicated in the functional divergence between bestrophin 2 and 3 are located in close proximity to this highly conserved region (Q354 and D366 in bestrophin 2, respectively).

Conclusion

This study addresses the evolutionary history of the bestrophin protein family. The precise functional properties of the four family members are still unclear but may comprise activities as Ca²⁺-activated chloride channels or aspects of Ca²⁺ channel regulation. Functional divergence between the paralogous bestrophins suggests that the bestrophin family members have evolved different functional properties after gene duplication events which have occurred early in the vertebrate lineage. This is most evident for bestrophin 4 where a number of experimental studies provide further support for the in silico data. Our study also demonstrates that amino acids critical for functional divergence are located in regions of the proteins that are likely accessible to soluble ligands, supporting the biochemical and physiological significance of this prediction.

Methods

Data collection

PSI-BLAST and TBLASTN searches with protein sequences of the four human bestrophins were performed in protein databases or unfinished genome sequencing projects at NCBI [37], ENSEMBL [38], the Sanger Institute [39], UCSC Genome Bioinformatics Group [40], and the Joint Genome Institute [41]. Proteins identified by the BLAST search algorithms were considered as potential homologues when amino acid identity was above 25% over a stretch of = 200 amino acids. After removal of expressed sequence tags, splice variants and redundant sequences, the initial data set included 173 distinct sequences from 43 species (see Additional file 1).

Sequence alignment and phylogenetic tree reconstruction

Amino acid sequences were aligned with CLUSTAL W (version 1.82) [42] and alignments were refined manually in Genedoc [43]. Incomplete sequences, and highly divergent regions or gaps resulting in uncertain alignments were excluded from the analysis. The final data set includes a total of 53 sequences from 21 vertebrate species (Table 1). A total of 367 N-terminal amino acids were aligned for further studies. Accurate nucleotide alignments were obtained with PAL2NAL [44], a program which constructs multiple codon alignments from matching protein sequences. ProtTest v1.4 [45], implementing the Akaike Information criterion (AIC) was used to estimate the most appropriate model of amino acid substitution for tree building analyses. The best fit model of protein evolution for the bestrophin protein family according to ProtTest corresponds to a JTT+I+G+F model [19].

Tree reconstructions were done by the Neighbor-joining method (NJ) [20] from the protein alignment done in the MEGA v3.1 software package, with the gamma distribution model implemented to account for heterogeneity among sites [46], and rooted with the bestrophins from Urohordata. The shape parameter of the gamma distribution (α) was estimated using baseml from the PAMLv4.0 [47] to be α = 0.51. Support for each phylogenetic group was tested using 1,000 bootstrap pseudoreplicates. Tree topology assessed by maximum parsimony (using Mega v3.1), was substantial similar with the NJ tree.

Topological structure prediction

Topological structure prediction was done with the TOPPRED II software [48], which is based on the Kyte and Doolittle algorithm [49]. The average hydrophobicity values of putative transmembrane domains of 20–23 amino acid residues were calculated according the Eisenberg scale [50]. An average hydropathy plot of 53 bestrophin-related protein sequences was generated by the PEPWINDOWALL program with a window of 19 amino acids [49].

Positive selection assessment

DNA sequences and related multiple protein sequence alignments were submitted to the PAL2NAL web server [51] which converts a multiple sequence alignment of proteins and the corresponding DNA sequences into a codon alignment. The resulting codon alignments and NJ tree were used in the program codeml from the PAMLv4.0 software package [47] to calculate the dN/dS (or ω) ratio for each site and to test different evolutionary models. The site specific models recommended by Anisimova et al. [26] were tested: Model M0 (one ratio), M1a (nearly neutral), M2a (positive selection), M3 (discrete), M7 (beta) and M8 (beta+ ω). Model M0 assumed a constant ω-ratio, while in models M1a and M2a ω-ratio is estimated from the date (0 < ω₀ < 1) while ω₁ = 1 is fixed. M7 and M8 assume a β-distribution for the ω-value between 0 and 1. Models M2a, M3, and M8 allow the occurrence of positively selected sites (ω >1). Subsequent likelihood rate comparisons of M0 with M3, M1a with M2a, and M7 with M8, respectively, were performed to test which model fits the data significantly better. Twice the difference in log likelihood between the models is compared with a chi-square distribution with n degrees of freedom, n being the difference between the numbers of parameters of the two models.

Functional divergence and detection of amino acids critical for altered functional constraints

Bestrophin sequence duplication events were tested for type I functional divergence [52] based on the method by Gu et al [27]. The analysis was carried out with Diverge (version 2.0) [28]. This method is based on maximum likelihood procedures to estimate significant changes in the rate of evolution after the emergence of two paralogous sequences. Type I sites represent amino acid residues conserved in one subfamily but highly variable in another, implying that these residues have been subjected to different functional constraints. A set of 53 protein sequences was included in the study (Table 1, see Additional file 2). Due to gaps and a shorter length of bestrophin 2 and 4, a total of 15 residues from human bestrophin 4 (codons 264–278) and one (codon 356) from human bestrophin 2 were excluded from the analysis. Consequently, each bestrophin paralogue was restricted to 367 amino acid residues. A new NJ tree was constructed within Diverge with Poisson distance and re-rooted. The coefficient of functional divergence (θ) and the posterior probability for the functional divergence were calculated for each position in the alignment. To detect amino acid residues reflecting functional divergence, bestrophin subfamilies were pair-wise compared to each other. The cut-off value for the posterior probability was determined by consecutively eliminating the highest scoring residues from the alignment until the coefficient of functional divergence dropped to zero.