1. Introduction

A big societal challenge is to increase the yields of crop plants without detrimental impact on the environment. Among the major crops, seed crops are one of the main sources for human sustenance (http://www.fao.org/docrep/007/y5609e/y5609e02.htm). In this respect, plants like the common bean representing the genus Phaseolus are very attractive as they are able to fix atmospheric nitrogen. The genus Phaseolus contains five domesticated species with distinct traits: common bean (Phaseolus vulgaris L.), runner bean (Phaseolus coccineus L.), lima bean (Phaseolus. lunatus L.), tepary bean (Phaseolus. acutifolius A. Gray), and year bean (Phaseolus. polyanthus Greenman). P. coccineus is closely related to P. vulgaris and it is the most important grain legume consumed worldwide for its edible seeds and pods ^[1,2].

The seed development can be divided into three phases: (1) the embryo biogenesis, differentiation and histodifferentiation stage, which occurs with a slow increase in biomass rate, (2) the maturation and filling stage featuring an intense metabolism and a huge biomass increase and (3) the desiccation stage. Recently, a proteomic study performed for the stages, from late embryogenesis to seed desiccation, has highlighted the dynamics of the protein profile during the development of P. vulgaris seeds ^[3]. An increased proportion of proteins involved in the general metabolism including the tricarboxylic acid cycle and glycolysis as well as DNA and RNA metabolisms characterizes the first stage. The active synthesis of storage macromolecules occurs during the second stage while in the last one, the energy transformation is slow down and the synthesis of proteins related to the oxidation/reduction processes increases. Early physiological studies on P. lunatus seeds show the correlation between the initial rate of respiration and seed germination ^[4]. The seed germination starts with a relatively fast growth, which implies a high requirement of catabolism based on large storage macromolecules, mainly proteins and starch, as well as anabolism, i.e. de novo synthesis of proteins, lipids and sugar polymers as well as other metabolites required for the growth and development of the plant.

The intensive respiratory and metabolic activity observed during both seed development and germination requires an efficient exchange of metabolites between the mitochondrial matrix and the cytoplasm. This exchange involves transport of solutes through the mitochondrial outer and inner membranes (MOM and MIM, respectively), Several specific transporters that belong to the mitochondrial carrier family (MCF) ^[5], potassium channels ^[6,7,8] and ABC transporters ^[9,10] enable to transport solutes through MIM. The voltage-dependent anion-selective channel (VDAC) is the major pathway for metabolite (e.g., ATP) and ion transport through MOM. As such, it may have central role in the regulation of mitochondrial and cell metabolism and thereby in the seed development and germination. Despite the fact that VDAC is a predominant proteins of MOM, there is growing evidence indicating its extra-mitochondrial localization ^[11]. For instance, three out of the four VDAC isoforms of Arabidopsis thaliana, i.e. AtVDAC1, AtVDAC2 and AtVDAC3, may be targeted to both mitochondria and the plasma membrane (PM), and the latter one is the most abundant in PM ^[12,13]. The PM localization was also shown for VDAC1 in mammalian cells ^[14,15]. The function of plant VDAC in PM is so far unknown. For mouse cells, it has been suggested that VDAC1 could function as a redox enzyme ^[16] although, the mechanism underpinning this function remains to be elucidated.

The first evidence of plant VDAC dates back from the sixties of the 20^th century when electron microscopy studies on P. aureus mitochondria showed the presence in MOM of large (2.5-3.0 nm in diameter) "pits" ^[17]. Then, it has been shown that the "pits" co-occur at a high surface density (2 × 10⁴ protein/µm²) with a 30 kDa protein ^[18] that is permeable to a polymer of about 5 kDa ^[19]. The protein was then purified from mitochondria of various plants and its typical VDAC features were detected by electrophysiology ^[11]. Moreover, a protein, named PcVDAC32, in reference to its apparent molecular mass, was purified from P. coccineus seeds and was used for both structural and functional studies ^[11]. Importantly, P. vulgaris genome and transcriptome sequence availability ^[20] enables to get insight into specific features of the different VDAC isoforms in this genus. As the main grain legume consumed worldwide, the common bean (P. vulgaris) is generally considered as a model for food legumes. In parallel to the genetic improvement made to increase their agricultural and food quality, efforts should be undertaken to understand their physiology at a molecular level.

Mitochondria are the place of ATP required for cell metabolism and contribute to plant response to biotic and abiotic stresses. Genetic modifications that improve some plant trait should also preserved or enhanced the mitochondrial function notably to avoid adverse consequence in cell bioenergetics and in plant resistance to environmental clues. VDAC plays a key role in the exchange of ions and metabolites (e.g., ATP, ADP, NADH, succinate, malate, tRNA) between the cytosol and mitochondria ^[11]. Importantly, there is a growing number of data indicating a role of VDAC in stress-induced regulated cell death in plant cells. The molecular mechanisms underpinning the phenomenon are still elusive. As VDAC plays a central role at the cytoplasm-mitochondrial interface, it is of prime importance to identify its various isoforms that are expressed in different plants, to study their expression, to decipher their structure-function relationship and to understand regulating mechanisms. In this framework, we provide here the current state of our knowledge of Phaseolus VDAC functional and structural properties hoping that it will pave the way for future research.

---

2. Functional Features

Two VDAC isoforms were purified from seeds of Phaseolus coccineus var. streamlined and named PcVDAC32 and PcVDAC31 at the basis of their apparent molecular mass estimated by SDS-PAGE ^[21]. It has been shown that PcVDAC32 is three times more abundant than PcVDAC31 but they have the same channel activity ^[21] and overall secondary structure ^[22]. Moreover, in two-dimensional gels PcVDAC32 and PcVDAC31 isoforms migrate as two polypeptides with pI of 7.6 and 8.8, respectively that permits their separation by chromatofocusing. In addition, the two isoforms display a different proteolytic pattern upon limited digestion with the protease V8 from Staphylococcus aureus that suggests a significant difference in their amino acid sequence ^[21]. Because of its greater abundance and because of the possibility to purify PcVDAC32, this isoform was chosen for functional studies.

The functional properties of VDAC are usually studied by electrophysiology after reconstitution of one or several proteins in a planar lipid bilayer for single channel or multichannel analysis, respectively ^[23]. Channels are usually characterized by their conductance, selectivity and gating. The conductance reflects the rate of ion flow through a channel. In a salt concentration close to that found in vivo (0.1 M) PcVDAC32 open state has a conductance of ~0.4 nS (Figure 1). The conductance increases linearly with KCl concentration (Figure 2), suggesting the absence of binding site for K⁺ and Cl^- ions.

For historical reason the conductance of VDAC channels is usually recorded in 1 M KCl. At this concentration, the open state conductance is nearly equal for various studied plant VDACs (Table 1). The open state selectivity of VDAC is weak. Nonetheless, its permeability towards inorganic anions is higher than that of inorganic cations (Table 1). This is consistent with its transport function since most metabolites flowing through VDAC are negatively charged (e.g., ATP⁴, ADP^3-, succinate^2-). VDAC is referred to as voltage-dependent because it usually switches from the fully open state to subconductance states in the presence of a difference in membrane electric potential (membrane potential |V| ≥ 20 mV, where |V| denotes the absolute value). A typical transition of a single channel from the open state to subconductance states is shown in (Figure 3A). It occurs at both positive and negative membrane potentials. In multichannel experiments the voltage-dependence features a symmetrical bell shape (Figure 3B). It can be quantified by plotting the change in the relative conductance of the membrane as a function of the applied voltage ^[24,25] and using the Boltzmann distribution: ln (P_o(V)/(1 - P_o(V)) = ± ne(V - V_h)/kT where P_o(V) is the probability of occurrence of the open state, V_h is the voltage at which half of the channels are in their open state, n is a measure of the steepness of the voltage dependence, V is the voltage applied across the membrane, e is the elementary charge, k is the Boltzmann constant and T is the absolute temperature. P_o(V) is calculated from the relative change in membrane conductance: P_o(V) = (G(V) -G_min)/(G_max - G_min), where G_max and G_min are the maximal and the minimal conductance, respectively. G_max is obtained at low applied voltages (|V|< 20 mV) when the channels are predominately in the fully open state and G_min is calculated at high applied voltages (|V| > 50mV) when the channels switch to subconductance states. The values of n and V_h calculated for PcVDAC32 are 2.5 and 27 mV, respectively ^[24]. When the reconstituted VDAC undergo transition from the open state to subconductance states, the change in the membrane conductance is about 30% (Figure 3B). The magnitude of the conductance change depends on the membrane lipid composition and salt concentration ^[24].

The voltage-mediated change in PcVDAC32 conductance is usually correlated to a reverse selectivity, i.e. the channel becomes cation selective in the subconductance states (up to P_K/P_Cl ~ 40 in the presence of a fivefold KCl gradient) ^[21]. The fact that PcVDAC32 adopts switchable states with opposite selectivity suggests that, as VDAC of other organisms, it might control the metabolite exchanges through MOM. In support of this assumption, it is worth noting that Neurospora crassa VDAC is permeable to ATP in the open state but not after switching to a subconductance state ^[26].

A large number of VDAC proteins purified from mitochondria of organisms representing different phylogenetic groups share the same fundamental properties in respect to the conductance, selectivity and voltage dependence. Thus, despite the low similarity of their amino acid sequence ^[11], their basic functional properties are conserved in the course of evolution although this does not exclude the existence of specific VDAC regulation in cells of given organisms. For instance, PcVDAC32 selectivity is modulated by stigmasterol but not by cholesterol, the sole mammalian membrane sterol, and the PcVDAC32 voltage-dependence is inhibited at low ionic strengths that prevail in vivo ^[24], which is not observed for Neurospora crassa VDAC ^[26,32,33].

---

3. Structural Features

The three-dimensional (3-D) structures of mouse and human VDAC1 isoforms (mVDAC1 and hVDAC1, respectively) as well as of zebrafish VDAC2 isoform were determined at an atomic level using NMR, X-ray crystallography, or a combination of both the methods ^{[34,35,36,37]}. These isoforms form a large β-barrel-shaped pore made of 19 β-strands and one N-terminal segment folded into α-helix that is aligned almost parallel to the plane of the membrane at mid-height of the pore. Thus, the N-terminal segment creates a constriction in the central region of the pore. The results of solid-state NMR study confirms this orientation of the helix ^[38] and the propensity of the N-terminus to adopt an helical fold was shown by experimental and computational approaches ^[39,40]. It should be worth mentioning that these atomic resolution structures has been questioned ^[41,42] because they are not consistent with previously proposed empirical functional models. However, there are theoretical and experimental clues of their biological relevance ^{[11,43,44,45,46,47,48]}.

As mentioned above, across different phylogenetic groups, VDAC amino acid sequences have a very low identity (≤5%) when considered as a whole but a relatively higher (≥25%) pairwiseidentity ^[11] suggesting that the protein sequences display common features. Firstly, they contain long stretches of alternating hydrophobic and hydrophilic amino acid residues typical of β-strands forming the b-barrel structure and a region with amino acid content characteristic for α-helix at N-terminus. Secondly, two conserved amino acid motifs have been detected close to C-terminus, namely the eukaryotic porin motif (PROSITE entry: PS00558) ^[49] and the β-motif regarded as essential for eukaryotic β-barrel protein import and assembly ^[50,51,52]. In addition, a motif specific to "green plants" (Viridiplantae) and located in N-terminal region has been identified ^[49]. Thus, VDAC proteins of different organisms might share a similar structure. For proper functioning of a hydrophilic pore formed by a transmembrane β-barrel protein a strict conservation of amino acids is not essential because the conservation of a proper arrangement of hydrophobic and hydrophilic residues should be suitable to preserve the β-barrel fold.

Spectroscopic studies indicate that VDAC proteins have the same secondary structure. Far-UV circular dichroism spectra of fungal, plant and human VDAC, solubilized in detergent, display negative extremum at 218 nm and positive peak at 195 nm, which are typical of β-pleatedsheet ^{[53,54,55,56,57]}. The secondary structure of PcVDAC32, reconstituted in a planar lipid bilayer, was also investigated using the attenuated total reflection-Fourier transform infrared spectroscopy (ATR-FTIR). The results strongly indicate the occurrence of α-helix and antiparallel β-strands that are tilted by an angle of approximately 45 degrees and a proportion of amino acid residues involved in β-strand formation which is consistent with the three-dimensional structure data of mammalian VDAC1 ^[11,22]. Thus, comparative modelling based on mVDAC1 crystallographic and NMR studies was performed to build 3-D model of PcVDAC32 ^[11]. Normal mode analysis suggests that this 3-D model features a higher flexibility of β-strands 1 to 6 ^[11], in agreement with results of NMR studies ^[34,58], crystallographic B factor analysis and computational calculations ^[58] performed on the mammalian VDAC1. The electrostatic free energy profiles computed for the PcVDAC32 model and mVDAC1 crystal structure are similar and show that the transport of Cl^- is favored relative to that of K⁺. Finally, for the PcVDAC32 model the distribution of the hydrophobic and charged amino acid residues are akin to those of mVDAC1 crystal structure. Notwithstanding these similarities, the structural model of PcVDAC32 must be treated as a workinghypothesis to guide experiments and it must be tested against experimental results.

Substitutions of charged amino acid residues in the yeast Saccharomyces cerevisiae VDAC1 impact its selectivity indicating the major contribution of the electrostatic interactions to the VDAC selectivity ^[59]. The selectivity of N. crassa VDAC and PcVDAC32 are salt solution concentration-dependent ^[60,61]. This is a strong experimental evidence supporting the electrostatic contribution to the selectivity. A macroscopic one-compartment electrodiffusion model containing the effective fixed charge and the mobile ions diffusing inside the pore describes correctly PcVDAC32selectivity ^[61]. The conductance and selectivity values computed from Brownian Dynamics (BD) simulations performed for the PcVDAC32 3-D model agree with the experimental data ^[61]. Interestingly, the concentration dependence of both selectivity and conductance measured for PcVDAC32 is fairly well describe by the BD simulation performed for mVDAC1 structure. Thus, the data suggest a conserved transport mechanism of ion transport.

---

4. Phaseolus VDAC Isoforms

Both P. vulgaris and P. coccineus belong to the same clade, which is relevant for comparative genomic and transcriptomic studies. We used genomic and transcriptomic data obtained for Andean landrace P. vulgaris ^[20] available at http://www.Phytozome.net/ to get further insight into P. vulgaris VDACs.

We found seven sequences described as P. vulgaris VDACs but not verified previously. Compared to other organisms, plants have a larger number of isoforms that reflects their molecular complexity. Whereas the origin of this complexity is still elusive there are some clues suggesting that this diversity might be correlated to more specific protein-protein interactions and isoform-specific functions ^[62]. Plants are sessile organisms that must face various biotic and abiotic stresses. Therefore, we might hypothesize that selective pressures associated with plant stress defenses have led to some VDAC isoform specificities after gene duplication. Moreover, a differential expression pattern of VDAC isoforms has been shown in different plant tissues ^[11], which suggests that another reason for functional specificities might be related to tissues specificities.

The theoretical physicochemical properties of the encoded protein sequences (named here PvVDAC1-PvVDAC7), supplemented by those of PcVDAC32, are summarized in Table 2.

The sequences of PvVDAC6 and PvVDAC7 are significantly shorter than the other ones. Moreover, the net charge of PvVDAC6 is negative whereas PvVDAC1-4, PvVDAC7 and PcVDAC32 have a positive net charge, which is consistent with anionic selectivity of VDAC in the open state. In agreement with experimental results ^[21], the predicted proteins with exception of PvVDAC6, have basic pI. The theoretical pI for PcVDAC32 is similar to the experimental value(pI = 8.8) ^[21]. The Gravy score (a measure of the average hydrophobicity) for the studied proteins is negative that is a typical feature of hydrophilic proteins. This is also observed for bacterial porins and β-barrel proteins of the chloroplast envelop ^[64]. It is due to the low content of hydrophobic residues, which are required in transmembrane β-barrel proteins to interact with membrane lipids (one out of two in transmembrane β-strands).

The phylogenetic tree generated for the predicted PvVDAC proteins contains three main clusters (Figure 4). The pairwise distances separating the sequences are given in Table 3. PcVDAC32 and PvVDAC4 are located very closely to each other and their sequences display 98% identity (99% similarity) whereas PvVDAC4 and PvVDAC5 sequences share 74% identity (88% similarity). In the most distant group, PvVDAC6 and PvVDAC7 sequences have 55% identity (74% similarity). When compared with PcVDAC32, both PvVDAC6 and PvVDAC7 display 36% identity (54% and 50% similarity, respectively). This denotes high amino acid sequence variability of putative PvVDAC proteins. Remarkably, a good correlation between the PvVDAC sequence clustering (Figure 4) and the corresponding transcript expression (Figure 5) seems to exist. PvVDAC5 and PvVDAC4 are highly expressed and the later display the highest expression. This reminds PcVDAC31 and PcVDAC32, the latter being three time more abundant ^[21]. It appears that organisms have at least one abundant canonical isoform of VDAC that shares similar electrophysiological properties and sequence similarities with its cognate isoform from mitochondria of other species. For instance, this is the case for the mammalian VDAC1, yeast VDAC1 and plant PcVDAC32. The high amino acid identity/similarity between PvVDAC4 and PcVDAC32 suggests that they have similar structure and consequently function. Furthermore, we might assume a close relationship between PcVDAC31 and PvVDAC5 because their N-terminal regions comprising thirty amino acids [VKGPGLYTDIGKRTRDLLFKDYQNDHKFTI] are identical and their pI are close to neutrality ^[21]. Other predicted PvVDAC proteins have a low level of their respective transcript expression. Interestingly, despite shorter sequences predicted for PcVDAC6-7, their transcript expression occurs, which suggests that they might have a function.

Multiple sequence alignment of VDAC sequences from different organisms (Figure S1) confirms the presence of VDAC conserved motifs ^[49]. As mentioned in § 3 (Structural features), two conserved amino acid motifs, namely the β motif ^[50] and the eukaryotic porin motif (EPM) are common for mitochondrial β-barrel proteins including VDAC proteins. The latter is the sorting signal of mitochondrial β-barrel proteins. These two motifs are conserved in all the predicted PvVDAC proteins and in PcVDAC32 (Figure 6). Additionally, the Plant motif conserved in "green plants" ^[11,49] is also conserved in the studied sequences although in the case of PvVDAC6 and PvVDAC7 this motif is truncated at half of its length (Figure 7). Importantly, full sequences sharing an identity (≥80%) with PvVDAC6 or PvVDAC7 were also found for other genera, namely Vigna angularis, Glycine max and Cajanus cajan (Figure S2). They all belong to the same plant tribe (Phaseoleae) of the Faboideae, a subfamily of the Fabaceae, which explains the high identity score. As these sequences were found in different plant genera by independent laboratories, we can reasonably assume that these shorter sequences did not arise from a sequencing error.

Inside the EPM motif, the conserved pentapeptide VKA[K/R]V sequence has been proven to be essential for the assembly of yeast VDAC into MOM ^[66]. The EPM located at the C-terminus comprises 23 residues (Figure 6). In the sequence, the fifth residue is critical for targeting of Arabidopsis thaliana VDAC isoforms ^[67]. The position is occupied by residues P223, H223, P223 and Q221 in AtVDAC1, AtVDAC2, AtVDAC3 and AtVDAC4, respectively. Only proteins harboring proline at this position are exclusively targeted toward mitochondria ^[67]. With the exception of PvVDAC3, the other PvVDAC proteins contain proline at this critical position (Figure 8). In PcVDAC3, the proline residue is substituted for serine. This is a conservative substitution so we might expect that all the bean isoforms are targeted toward MOM.

We have already mentioned here that PcVDAC32 is structurally and functionally similar to mammalian VDAC1. To assess whether predicted PvVDAC proteins might fold like mammalian VDAC1 we applied homology modelling, which is a method of choice to generate 3-D model of a protein of known amino acid sequence. For this purpose we used the SWISS-MODEL available at http://swissmodel.expasy.org/ ^[71]. The 3-D structure model was generated using mVDAC1 X-ray structure (3emn) as a template (Figure 9). The backbone of the model generated with PvVDAC1-5 and PcVDAC32 sequences is consistent with that of mVDAC1 3-D structure. However, both PvVDAC6 and PvVDAC7 models are at odds with mVDAC1 structure developed by crystallography.

The missing amino acid residues in PvVDAC6 and PvVDAC7 start at position 25 of their sequences viz. about the beginning of the first transmembrane β-strands indicating that the conformation of these two isoforms may significantly differ from the experimental 3-D structures. We resorted to the SPARK^X fold recognition method ^[72] to generate their putative folding patterns (Figure 10A). The results indicate that four β-strands at N-terminus are disrupted in both PvVDAC6 and PvVDAC7. Additionally, for PvVDAC7, the first predicted β-strand is too short (three predicted residues) to form a transmembrane segment. Taking into account the 19 β-strands proposed the mammalian VDAC1, this suggests that the maximal number of β-strands in PvVDAC6 and PvVDAC7 structure might be 15 and 14, respectively. The number of residues between the first and last β-strands (7 and 9 residues, respectively) in PvVDAC6 or between the second and the last β-strands (8 and 9 residues, respectively) in PvVDAC7 should be sufficient to close properly the barrel. We noticed that the position of the N-terminal α-helix in both isoforms perfectly matched that of the mVDAC1 3-D structure (Figure 10B). Moreover, the three key residues involved in the main energy barrier for the inorganic ions and metabolites permeation are conserved and point toward the diffusion pore. Thus, this suggests that the helix might have a similar function in all PvVDAC isoforms.

Knowing that genes encoding PvVDAC6 and PvVDAC7 are expressed in plants (at least at the mRNA level) and that their predicted amino acid sequences include motifs conserved in mitochondrial β-barrel proteins (including a partial green plant VDAC motif) we can reasonably assume that the protein might adopt a β-barrel fold. The lower number of β-strands forming the pore will produced a barrel with a wall-to-wall diameter shorter (~2.5 nm) than the canonical VDAC (3.5 nm) ^[46]. We can thus anticipate that PvVDAC6 and PvVDAC7 might have a significantly different function than other PvVDACs. The issue of how these two proteins are folded and whether they have channel activity similar to VDAC is now open.

---

5. Conclusions

Our current knowledge of structure-function relationships in the case of plant VDAC proteins is mainly based on data obtained for PcVDAC32 from P. coccineus. However, although two PcVDAC isoforms are known to be expressed in P. coccineus seeds, little is known about the number of plant VDAC isoforms and their role in plant cell physiology. In this respect, seven putative VDAC encoding sequences termed here PvVDAC1-7 are found in the genome of P. vulgaris. PvVDAC4 has numerous properties similar to those of PcVDAC32, which suggests that they have similar structure and function. PvVDAC5 seems to be similar to PcVDAC31 and might have similar electrophysiological properties. PvVDAC1-3 sequences share structural similarities with the canonical VDAC found in plants, algae, fungi and animals. We can therefore assume that they belong to the VDAC family. The high expression of PvVDAC4 and PvVDAC5 transcripts suggests that the encoded proteins might be purified from mitochondria. The other ones would probably require the use of an efficient heterologous expression model that preserve the native protein properties. The sequences termed as PvVDAC6 and PvVDAC7 are shorter, and could form a β-barrel but with a smaller number of β-strands. Therefore, their annotation as members of the VDAC family should be subject to caution. This overview of Phaseolus vulgaris VDAC putative proteins highlights new information on their structure and function. Notably, their structural and/or functional similarities with VDAC of other plants pave the way for future research aiming at deciphering their role in plant physiology.

---

Acknowledgments

This work was supported by FRS-FNRS/PAN JOINT RESEARCH PROJECT. HK and FH are the project coordinators for Poland and Belgium, respectively. F.H. is a Research Director from the FRS-FNRS (Belgium).

---

Conflicts of Interest

All authors declare no conflicts of interest in this paper.

---