Three-dimensional Modelling of the Voltage-gated Sodium Ion Channel from Anopheles gambiae Reveals Spatial Clustering of Evolutionarily Conserved Acidic Residues at the Extracellular Sites.

BACKGROUND: The eukaryotic voltage-gated sodium channel(e-Nav) is a large asymmetric transmembrane protein with important functions concerning neurological function. No structure has been resolved at high resolution for this protein.
METHODS: A homology model of the transmembrane and extracellular regions of an Anopheles gambiae para-like channel with emphasis on the pore entrance has been constructed, based upon the templates provided by a prokaryotic sodium channel and a potassium two-pore channel. The latter provides a template for the extracellular regions, which are located above the entrance to the pore, which is likely to open at a side of a dome formed by these loops.
RESULTS: A model created with this arrangement shows a structure similar to low-resolution cryoelectron microscope images of a related structure. The pore entrance also shows favorable electrostatic interface.
CONCLUSION: Residues responsible for the negative charge around the pore have been traced in phylogeny to highlight their importance. This model is intended for the study of pore-blocking toxins. Copyright© Bentham Science Publishers; For any queries, please email at epub@benthamscience.org.

INTRODUCTION

The eukaryotic voltage-gated sodium channel (e-Nav) is an essential component involved in the conduction of impulses in nerve tissue. The formation of action potential has been observed in in vitro patch-clamp experiments across membranes obtained from neurons. The Hodgkin-Huxley model [1-3] used electrical components to model electrical response across such a membrane. Besides membrane capacitance and voltage difference, voltage and time-dependent conductances were needed. These correspond to voltage-gated ion channels, consisted of two independent parts, i.e. potassium and sodium ion conductances respectively. Biochemical studies subsequently showed that two separate channel proteins, the eukaryotic voltage-gated potassium channel (Kv) and sodium channel (e-Nav) were responsible.

Many biochemical studies are devoted to understand the aspects of the structure of the Nav channel [4, 5] voltage sensing and gating [6, 7]. The structure of the simpler eukaryotic Kv channel [8] allowed detailed studies on the e-Kv channel. The e-Nav consists of two subunits named alpha and beta. The alpha subunit is the main channel with the transmembrane region. Its secondary structure has been determined and shown to be similar to the e-Kv channel [8]. Prokaryotic Nav channel (p-Nav) structures were also determined [9, 10].

Both p-Nav and e-Kv channels form a homotetramer, with subunits around ~300-350 amino-acids (AA), with very similar conserved structure. The e-Nav and closely related eukaryotic voltage-gated calcium channels (e-Cav) [11] consist of one ~2000AA chain, with four duplicated and significantly diverged transmembrane domains (D1 to D4) corresponding to subunits of the homotetrameric p-Nav. Large protein regions in the extracellular and intracellular side between these repeats account for the significantly longer sequence than the total length of the p-Nav tetramer.

Atomic resolution structural determination of eukaryotic Nav channel is extremely challenging by virtue of its size, membrane-bound state and asymmetry. It falls out of scope for current x-ray crystallography and NMR techniques. Newer cryo-EM tomography methods are in their infancy and are better at resolving membrane associated molecules which exhibit symmetry [12], though these limitations are being challenged [13]. Low resolution 19Å cryo-EM density maps [14] of the structure of Nav from Electrophorus neurons gives an approximate idea of the goal [15], while other studies add to the data [16]. The density maps show a bell-shaped molecule, with entrance to the channel present at the side of a large dome extending into the extracellular side, which is likely to be rich in helical structures [15]. The shape of the transmembrane region is, however, similar to the p-Nav channel. The importance of the Nav channel leads to the possibility that its structure may be determined in the near future. Currently however, homology modelling may be used to obtain a possible model.

Humans and mammals have greater diversity in sodium ion channels [11, 17, 18] compared to insects. Humans have around 10 paralogs of the sodium channel, with nomenclatures Nav1.1 to Nav1.9 and Nax [19]. These Nav nomenclatures are used mainly for mammalian homologs. They have been the target of drug research concerned with neurological disorders such as epilepsy. They are also a target for neurotoxins, both peptides and compounds, natural and artificial. Insecticides, such as pyrethroids, are lipophilic and target the channel from the lipid layer, preventing it from deactivating [20]. Pore-blocking toxins such as alpha-scorpion toxins [21] which are peptides, or toxins such as tetrodotoxin [TTX] [22] which are small molecules act by deactivating the channel. Different kinds of toxins and their actions on Nav channels have been classified [23] and is an active area of research. The extracellular loop regions play a major role in peptide toxin binding [24]. These may have evolved in an arms race between toxin producing organisms and target ion channels which evolved immunity to them: the evolution of immunity to toxins in garter snakes is one example [25].

This study concerns modelling the pore, transmembrane scaffold and the extracellular regions of an insect e-Nav channel – the malarial mosquito Anopheles gambiae. This model may then be subjected to theoretical studies such as molecular dynamics and docking to test scenarios and explain behaviour [26]. The human counterparts may also be modelled from this. The scope of this study currently does not consider aspects such as voltage gating or conductivity, though they are planned. Insect Nav channels have been extensively studied for the purpose of insecticide design against acquired resistance [28].

Structural Overview

The voltage-gated sodium ion channel ‘para’ from mosquito Anopheles gambiae is represented by Uniprot [30] protein sequence A5I843_ANOGA [28] with a sequence length of 2139 residues. The sequence shows homology (~81.5% identity) with the annotated sequence SCNA_DROME corresponding to the para [31] channel of closely-related Diptera model organism Drosophila melanogaster, so named because knockout of this gene caused paralysis in the fly. There is no other known e-Nav paralog in Drosophila.

The closest homologues are the voltage-gated eukaryotic channel proteins with resolved complete structures, including both pore and voltage-sensing subdomains, are homo-tetrameric e-Kv channels represented by pdb: 3LUT [32], pdb:2A79 [33] and pdb:2R9R [34]. The p-Nav structures available are NavAb pdb:3RVY [9] from Acrobacter butzleri and its calcium-conducting mutant pdb:4MS2 [35] differing only by three residues near the selectivity filter region, NavRh pdb:4DXW [10] from Rickettsiales sp. HIMB114 representing a closed form, and the cryo-EM structure pdb:4BGN [36] from Caldalkalibacillus thermarum.

The e-Cav and e-Nav channels have four duplicated transmembrane domains each homologous to the e-Kv and p-Nav sequence fused into a single chain and are therefore pseudotetramers. Each domain consists of transmembrane helices S1-S6 and other structures in between, such as the intermembrane Pore helix and the filter loop. Substantial segments of the protein between these transmembrane domains may form supersecondary structures present on the cytosolic side, of currently unknown function. Only the cytosolic C-terminal region has been determined to have a calcium-binding (pdb:4JPZ) [27], but this is irrelevant to the current study.

The four domains (see Fig. ) are D1:133–463(330), D2:731–1059(328), D3:1287–1596(309) and D4:1605–1855(255) and intermediate regions I1:464–730(266), I2:1060–1286(226), I3:1596–1605(9) and CT:1856–2139(283). The I3 section is functionally important and called the ‘Na channel gate’ or IFM loop. The I1 and I2 regions are called the DUF3451 superfamily and the Na-trans_assoc (transport associated) subfamily domain respectively, which are both conserved across e-Nav but with unknown function and structure. Each of the duplicated domain (D1..D4) have diverged significantly from each other (Table ), and the core channel is composed of four related but very different units. The low sequence identity of <20% however is compensated by high structural similarity and compatibility to the channel of interest.

For each of these repeats, four helices (S1..S4) form the voltage sensing subdomain, with helix S4 characterized by alternating R (arginine) and sometimes K (lysine) residues whose protonation/deprotonation can affect the shape of the pore. The pore itself is defined by the pore subdomain, consisting of helices S5, S6, the pore helix and the filter region. The connecting helix, labelled ‘Conn’ transmits the force necessary for changes in orientation in the voltage-sensing subdomain to the pore subdomain.

The filter residues, shown by circles representing ions in Fig. (, have a characteristic conserved DEKA tetrad motif or filter contributed by D1 to D4 respectively. In addition, the p-loop region extends into the extracellular part in domain D1 and D3. The D3 and D4 domains however have extracellular extension labelled q-loop by analogy with the p-loop, and, in the case of D4, a very short p-loop. D2 has a longer C-terminal side region (represented by additional helix), and a short p-loop, longer than that of D4. The transmembrane regions consist of four repeats, which have diverged from each other but still retain basic common structure.

Fig. (1)

Structure of Nav channel. The NT (N-terminal), I1, I2 I3 and CT (C-terminal) regions are represented as bins because their secondary structure is unknown and irrelevant to this study. The structure of the CT region however has been determined for human Nav channels, and is shown to have a calcium binding function, interacting with calmodulin and other factors. The circles indicate position of ion near the filter loop, which contribute to the DEKA tetrad.

Additional cytosolic side I1...I3 and CT regions of the protein are not modelled, and the transmembrane region is treated as a heterotetramer, and accessory subunits ignored. Around 30% of the structure is glycosylated [37] by N-linked glycans with negative charge, such as sialic acid. Removal of these glycans does not impair function or folding [37]. However, this adds a constraint to any model, as predicted glycosylated Asparagine residues must be exposed or surface-accessible for glycosylation.

Homology modelling has been used to model the protein [20] restricted to domains D2 and partially domain D3, as this region is known to interact with the lipophilic DDT (dichlorodiphenyltrichloroethane) and pyrethroid pesticides such as pyrethrin I. The goal of these studies were to understand the cause of immunity gained by insects to these insecticides. Other attempts at modelling these channels are seen on ModBase [38] and elsewhere [39], concentrate on one or another singular domain. A recent attempt at homology modelling looks at the pore region exclusively [40] and has been used for studying the permeability of the channel. The extracellular loop region has been ignored in this model, treated as loops, and the helices S5 and S6, the pore helix and filter region are the main focus of the study. An attempt to model the extracellular loops in the eel electrophorus Nav channel used a postulate that the structure may be an alpha-turn-beta structure [37, 41], but this hypothesis has not been verified. The low sequence identity to other membrane proteins of known structure and lack of secondary structure predicted by commonly used secondary structure prediction methods such as PSIPred [42] or JPred [43], adds to the challenge of modelling these portions of the protein.

The human two-pore potassium channel pdb:3UKM [44] is a dimer representing a channel protein with extracellular regions. Each subunit of the dimer may be split or cleaved into two parts, each correspond to parts of a repeat domain of the voltage gated channel. This part representing the core helices S5, pore helix and filter region, and helix S6. Only the first part has an extracellular extension, corresponding to that of D2 and D4.

Methods

Full-length Alignment and Phylogeny

The target putative sequence of sodium channel protein A5I843_ANOGA [28] from Anopheles gambiae obtained from UNIPROT [30] was aligned first with other known sodium channel proteins. The Diptera para channel protein SCNA_DROME [14] from Drosophila melanogaster was included, as it is well-annotated. Two mammalian sodium channel proteins, one from human and one from rat were included in the alignment.

The detailed annotations in SCNA_DROME were used to determine boundaries of other sequences by alignment. The four transmembrane domain repeat boundaries D1, D2, D3 and D4 were marked and used to delineate similar boundaries on the other sequences. Transmembrane helices as well as the intermembrane helices which denote S1-S6 and the pore helix respectively were marked. The characteristic filter sequences are aligned.

This alignment with around 1500–2000+ residues was used as a reference, and not for modelling purposes, as there is no structure of a eukaryotic sodium or calcium channel resolved yet. The JalView [45] interface was used, and the alignment programs MUSCLE [46], CLUSTALW_Realign [47] and MAFFT [48] were used via this interface, to obtain the alignment. This alignment was extended by addition of total 242 e-Nav sequences obtained from UNIPROT [30], sampling deuterostomes such as amniotes, fishes, tunicates etc. and protostomes such as annelids, molluscs and specifically insects (see Supplementary file 1). RAxML [49] was used with the maximum likelihood method to obtain phylogenetic tree with bootstrap 100 (Supplementary files 2 and 3).

Per-repeat alignments

Each repeat D1 - D4 were aligned to each subunit of the tetrameric template 4MS2 and e-Kv sequences. Repeats D1 and D3 were aligned to type A subunit of the template 3UKM and repeats D2 and D4 to type B subunit of the template 3UKM to yield a heterotetramer (3UKM_mod). Flexible structural superposition was used to obtain the initial sequence alignment using MATT [50] and the alignment was adjusted manually by visual inspection.

Loop alignments

The extracellular p-loops of domains D1..D4 were separately aligned. It was observed that these domains were either conserved in length (D2, D4) or variable (D1) but have common conserved segments. These conserved segments were used to define putative secondary structure of the extracellular loops.

The conserved segments may represent helices, and the extracellular region may be similar to the structure seen in 3UKM. The loop alignments for D1, D3 and D4 are shown in Figs. (, and ).

Fig. (3)

Alignment of the D3 p-loop. It shows the conserved cysteines, and highlights important acidic residues.

Fig. (4)

Alignment of the D4 p-loop. This loop forms the entrance to the pore in the model and has conserved DD or ND motif. Asparagine (N) residues serve as attachment points for glycans like Sialic Acid which are negatively charged and enhance the functioning of the channel.

Modelling the domains

The domain alignments were concatenated into a longer alignment. A short linker (I3) between D3 and D4 was also included with its corresponding secondary structure predictions. The program MODELLER [51] version 9.15 was used. The whole tetrameric structures of templates 4MS2 and 3UKM_mod were used as templates. While pdb:3RVY represents a sodium channel and pdb:4MS2 a calcium mutant, there is very little root mean square difference (RMSD) between the two. Ions are not present in pdb:3RVY. Calcium ions present in pdb:4MS2were used as a structural reference.

The extracellular loops in D1 and D3 differ in size, unlike the template 3UKM with which they align. The loops from D1 seem to contribute three helices and that from D3 contributes one helix. The structure obtained is similar to the template, which has two helices from each subunit. One opening is left open to the channel, while the other side is blocked by the extra loops connecting the helices. Trial and error adjustments lead to the putative model shown in Fig. .

The secondary structure predictions and disulphide connections were incorporated in the MODELLER script as restraint patches. Ions were treated as ‘block residues’ and distances from backbone atoms in interacting residues to stationary ions were added as restraint to aid modelling the filter region. To reproduce the putative arrangement shown in Fig. (, the 3UKM_mod was included again, with most parts unaligned and outside the main alignment (Fig. ). Only the the second helix of the second subunit (h4 in Fig. ) extracellular region was aligned with the third helix position of D1 (For complete alignment, see Supplementary file 4).

Fig. (5)

Proposed secondary structure configuration of the extracellular region, fitted to the template provided by the human two-pore potassium channel structure pdb:3UKM. The proposed structure takes into account the length of D1 and D3 p-loops and their conservation in alignment. At least 4 cysteine residues are present, with the KC sequence common for D1 and D3 p-loops at around the same place, implying symmetry. Sequence-based structure prediction does not predict such a structure, due to absence of structural representatives in known structures. The D2 p-loop provides the negative electrostatics necessary at the pore by the presence of aspartates, with D1 and D3 providing additional aspartate and glutamates. Attachment points for N-linked glycans are also present around the pore entrance.

The presence of several cysteines in the extracellular loops suggest that many disulphide bonds combinations are possible. Two such combinations were tested, one involving C207 and one involving C216 from D1 in the bond with C758 on D3. 20 models for each set were generated and ranked. Further refinement was done to understand the possibilities. Top-ranked models which showed unnatural knots and entanglements between the subunits were discarded. The current model discussed involves a bond involving C207. A conserved cysteine is also present on the D4 q-loop.

The resultant model shows possibility of more refinement and adjustment, such as the possibility of an anti-parallel beta sheet formed by the q-loop of D4 and these extra loops, as shown in Fig. (. I-TASSER [52] predicted secondary structure predict the presence of beta strands which may form these. However, at this stage, the focus is on the opening to the pore. The resultant extracellular region in top view is shown in Fig. (.

Fig. (7)

Configuration of the modelled extracellular regions with lower areas clipped. The ions which were included as block residues are seen as spheres in center. This is comparable to the schematic in Fig. 5. The domains are labelled D1 to D4, with italic labels denoting the pore subdomain and bold normal text denoting the voltage sensing subdomain. The helix positions are labeled h1 to h4.

Simulation

A short simulation was carried out to check the stability of the model in a lipid-embedded environment. GROMACS 5.1 [53] was used. The membed [54] protocol was used to embed it in a dipalmitoylphosphatidylcholine(DPPC) pre-equilibrated and extended bilayer provided by Slipid [55]. The position of membrane was calculated using the PPM server [56, 57] and VMD [58] was used to align the bilayer to the protein. The Slipid parameters in conjunction with AMBER FF99SB [59] forcefield modified for Gromacs (ffamber) [60, 61]. Simulation was carried out for 10ns after a 1ns equilibration. ‘Block’ residues in the MODELLER generated model were replaced by Na+ and corresponding 4 neutralizing Cl- ions added. A total of 375468 particles, with 1101 protein residues, 1062 DPPC residues and 73146 water molecules was run on a GPGPU-based system. The RMSD and radius of gyration plots are provided in supplementary files 5 and 6 respectively.

Results

Whole sequence phylogeny analysis using eukaryotic channel Nav sequences shows that insect Nav sequences like Anopheles gambiae para-like and Drosophila para channel cluster close to the SCN4A sequence of human sodium channels.

The p-loops of the domains D1-D4 show variability in length in D1. Fig. shows the alignment of the p-loop of the D1 domain, and some flanking transmembrane regions. Anopheles gambiae (ANOGA) sequence, the Drosophila homolog para (SCNA_DROME), human representatives, one rat representative and a channel from the electric eel Electrophorus have been included. The latter is performed because of the availability of its cryo-EM structure [15].

The sequence of SCN4A and the insect channels is particularly different. The SCN4A human has repeats of the sequence -WYGNDS- with point mutations. The SCN4A rat sequence is similar with one such repeat conspicuously absent (see Fig. ). These additional sequence lengths could contain more secondary structures or disordered regions or glycan binding sites that may prevent certain toxins from binding. The model in Fig. can account for these, as the additions would be in the loop regions, conserving the overall scaffold. The Electrophorus channel sequence however has similar length to the insect channel sequence.

The model chosen here gave a favorable DOPE score, i.e. the lowest score model without showing structural abnormalities such as interchain knots. It also showed a well-formed pore entrance, similar to that predicted in the cryo-em data [15]. The cryo-EM data suggests a dome-shaped structure, with a pore entrance to one of the sides, not directly open vertically as in most resolved voltage-gated channel structures. Our model postulates a possible arrange-ment of secondary structures forming this dome. The charge distribution around the entrance to the pore is, as expected in cation channels, with an overall negative charge. There are also prominent pockets with attachment points for ASN N-linked glycans. While the ASN residues themselves may contribute a net positive charge, negatively charged glycans like sialic acid can provide a conducive electrostatic environment for attracting cations to the mouth of the channel (Fig. ).

The residues responsible for providing this negative charge were traced in an alignment for the four domain p-loops. In particular, the domains D1, D3 and D4 are mainly responsible for contributing this charge. Five such residues have been noted, numbered GLU 218 (350 as in UNIPROT numbering), GLU 762 (1490), ASP 1056 (1784), ASP 1057 (1785) and ASP 1092 (1820) in the resultant model. The numbers in brackets are sequence numbers on single-chain Uniprot sequence A5I843_ANOGA. The first is present on D1, second on D3, the third and fourth on the floor of the entrance on D4 p-loop and fifth on the D4 q-loop as seen in Fig. (.

Fig. (9)

The ribbon view of structure shown in Fig. (). The acidic residues surrounding the mouth of the pore are visible.

Conserved acidic residues are also present on D2. However, loop of D2 does not appear to be spatially proximate to the pore entrance in this model. There is, however, a possibility that the side chains of these residues may line the inner part of the pore within the “dome” as seen in the cryo-em image.

Residue 218 (350) (highlighted in Fig. ), along with a neighboring acidic aspartate residue, is present in insect ion channels, but not present in other ion channels. A toxin specific to binding these residues may be specific to insect channels.

The residue 762 (1490) is marked and is not conserved even in insects (Fig. ). The consensus sequence residue for these representative sequences is Lys. Residues 1056(1784), 1057(1785) are highly conserved on the D4 p-loop as seen in Fig. (-mostly as a sequence of two residues DD or sometimes

ND. Residue 1092 (1820) is present on the D4 q-loop. The D4 q-loop also has a highly conserved cysteine. The per-residue conservation and overall alignment with possible templates is shown in Supplementary File 9.

Minimization and simulation showed that the secondary structure of the extracellular loops remain stable. Minimization was done after placing the protein in the lipid bilayer. The ProSA–web interface [62] was used to assess the equilibrated model (Supplementary 7,8). It showed slightly lower quality in the transmembrane regions, which is expected as ProSA assumes a soluble molecule. However, the extracellular regions have a good score, with an overall Z-score of -4.58. A single bend near residue K755 in the D3 p-loop was the only part with unfavourable energy. Sodium ions move in and out of the channel during the simulation from the extracellular side, showing free movement in the absence of a voltage or ion concentration gradient. In the presence of a gradient, the ions may move into the channel. The negative electrostatic areas near the channel mouth prevent the ion from going too far and allow it to return.

Conclusions

This study shows a possible model of the sodium ion channel from an insect Anopheles gambiae. Anopheles is the vector for Plasmodium falciparum which causes malaria. The model shows a clearly formed pore with favorable electrostatics. The model mimics the two-pore channel in its extracellular regions. The two-pore channel and the Nav channel may be share homology and relationship of the extracellular regions of the pore channel may be similar to the e-Nav channels with their p-Nav homologs. While the basic architecture is similar, compensatory reshuffling where helix h4 (Fig. ) has been incorporated into D1 from D3 allows the pore helix to have a singular opening.

The favourable electrostatics seen near the entrance to the pore is contributed by acidic residues such as aspartates and glutamates (D, E) which surround the entrance. Some of these aspartates or glutamates are conserved, while some are very specific to insects. The DE motif in Anopheles Nav D1 p-loop is not present in mammalian channels, which have a proline in its place.

Previous knowledge on three-dimensional information of such transporters were limited to either low-resolution or not the entire transmembrane domain. To our knowledge, this is the first modelling attempt of the full transmembrane domain of sodium channel of Anopheles gambie, including all the extracellular loop regions. The current model has satisfied several structural compatibility and stability criteria. There is, however, scope for improvement and refinement.

The model will be employed for future studies. A large number of chemical compounds are known to interact with this protein, particularly pyrethroids, DDT and other insecticides [20]. A preliminary blind docking using Autodock Vina 1.1.2 [29], utilizing the lipophilic Pyrethrin-1 as ligand, showed the docking site in the transmembrane region between domains D2 and D3, similar to the case described by O’Reilly et al. [20]. The mode is shown in Supplementary file 10.

The current focus is to understand the specificity of pore-blocking toxins towards sodium channels. The interacting residues of the extracellular loops with these toxins and their conservation in phylogenetic clades may help us understand toxin specificity to an insect such as Anopheles gambiae and related clades, while not affecting other insects or mammals.

Blind docking using the soluble tetrodotoxin as ligand showed that this ligand is at least capable of entering through the extracellular loop region to interact with the filter residues as shown by Lukacs et al. [40]. However, it is still some distance away. This mode is shown in Supplementary Fig. ().

Molecular dynamics simulations may further show interactions with the extracellular loops before the toxin moves into the pore. Apart from the filter residues, these regions may play an important role in toxin-channel interaction and cannot be ignored.

In future, such modelling studies could enable the detailed examination of human sodium channels, which are active targets for drugs treating epilepsy and other neurological disorders [63].

Table 1

Sequence identity of the four duplicated domains in A51843_ANOGA.

	D1	D2	D3	D4
D1	100
D2	23.84	100
D3	22.09	24.27	100
D4	20.57	29.26	24.14	100
3RVY	19.8	19.9	20.9	24.3
KCNA2_rat	16.31	15.91	14.9	21.59