Effects of active site mutations in haemoglobin I from Lucina pectinata: a molecular dynamic study

Haemoglobin I from Lucina pectinata is a monomeric protein consisting of 142 amino acids. Its active site contains a peculiar arrangement of phenylalanine residues (PheB10, PheCD1 and PheE11) and a distal Gln at position E7. Active site mutations at positions B10, E7 and E11 were performed in deoxy haemoglobin I (HbI), followed by 10 ns molecular dynamic simulations. The results showed that the mutations induced changes in domains far from the active site producing more flexible structures than the native HbI. Distance analyses revealed that the heme pocket amino acids at positions E7 and B10 are extremely sensitive to any heme pocket residue mutation. The high flexibility observed by the E7 position suggests an important role in the ligand binding kinetics in ferrous HbI, while both positions play a major role in the ligand stabilisation processes. Furthermore, our results showed that E11Phe plays a pivotal role in protein stability.


Introduction
Haemoglobins and myoglobins have been under intense scientific scrutiny because of their essential role in manyi mportant biological systems, and their widespread presence in all organism kingdoms. Although, there are manye fforts focused on the understandingo fl igand selectivity and stability, manya spects of the ligand binding andu nbinding processes observed in these hemeproteins are still under debate. Fore xample, it is well known that in most haemoglobins and myoglobins the ligand binding reaction consistso fa tl east two well separated kinetic processes: bimolecular process and geminate process [1]. Research has also provided some insight into how the haemeoglobin hemep ocket microenvironment electronically affects ligand interaction and the effects it has in the ligand binding -unbinding kinetics. It was initially suggested, however, that the sizeo ft he distal residue was the main factor controlling ligandbound structuresand ligand binding kinetics, but today it is knownt hat electrostatic forces play am ajor role. Despite these observations, acombination of factors such as heme pocket steric constrain, hydrogen bonding, local polarity, proximal and globally structurale ffects are still being proposed in order to explain ligand binding-unbinding kinetics andh ence, their selection ands tability in vertebrates and non-vertebrates haemoglobins [2].
In this respect, non-vertebrate haemoglobins have been increasingly studied because they can provide information relevanttothe evolutionofboth structural and functional aspects of the globin family [3]. These haemoglobins occur in widely anatomical sites among invertebrates such as the cytoplasm, red blood cells and body fluids,a nd are believedt ob ind and transportl igands other than O 2 . They exhibit much wider variationsi nt heir primary and quaternary structure [3]; however,t he tertiarys tructure surrounding the hemeprostheticgroup [Mb fold] is highly conserved among vertebrate as well as non-vertebrate haemoglobins [4]. The Mb fold consists of six to eight ahelicalsegments connectedbyshort b -turns or loops [5,6], which form at hree-on-three helicals andwich, where the heme group binds to the protein through ap roximal histidine.Ithas been proposed that ligand binding kinetics of thesen on-vertebrate haemoglobins ares trongly influenced by the structure of the heme cavity, particularly the size and polarity of residues occupying the distal portion which exerts steric and dielectric effects [5,7].
However, investigations into vertebrate myoglobins and haemoglobin have suggested that ligand binding and unbinding kinetics are not only affected by the sizea nd polarity of the hemep ocket residuesb ut also by other residues away from theh eme. Theseh emep ocket arrangements provide al igand stabilisation mechanism with oxygen through ah ydrogen bond with histidine [8], since most of vertebrate myoglobins and haemoglobins have ahistidine and aleucine in the distal positions E7 and B10, respectively.I nv ariousn on-vertebrate haemoglobins, the E7 and B10 positions are usually occupied by glutamine and tyrosine resulting in at ight cage for oxygen, which exhibits higher binding affinities relative to vertebrate myoglobins [3]. Ah ydrogen bond formation betweeno xygen and the TyrB10 and GlnE7ist hought to be responsible fort he ligand stabilisationi nt hese haemoglobins. Several studies have been performed with vertebrate myoglobins [9,10] to achieve the high oxygen affinities observedi nn on-vertebrate haemoglobins. Double and triple mutants of spermw hale myoglobins have been studied (ArgCD3Asp, HisE7Val, ThrE10Arg [9]; LeuB10Tyr, HisE7Gln,T hrE10Arg [10]) showing correct alterations in the active site, but these mutations did not reproduce the ligand binding kinetics properties of wild type non-vertebrate myoglobins and haemoglobins. The difference in tertiary structure between mammalian myoglobins and non-vertebrate globins is proposed as an explanation for the low degreeofsuccess in these studies, due to the differencesi nt he primary sequence conservation (and no phylogenetic relationship) [3]. Ishikawa et al. [1]d emonstratedt hatt he ligand binding-unbinding kinetics are not only affected by the amino acids in the hemepocket moiety, but by residuesout of the active site centre.
In order to evaluatel igand selection and stability in hemeproteins, we have performed molecular dynamic simulations using an invertebrate haemoglobin from Lucinap ectinata as am odel system in viewo fits unique properties. This clam contains three different haemoglobins: haemoglobin I( HbI),h aemoglobin II (HbII)a nd haemoglobin III (HbIII). Each haemoglobin has adifferent biological function, as well as different physical-chemical properties. Fore xample, HbI is am onomeric proteino f 142 amino acids residues [11],which exhibits the classical Mb fold [12].The distal residuesPheB10, PheCD1, GlnE7 and PheE11 comprise the HbI heme pocket [13],which is considered an atural occurring mutant haemoglobin [14] compared to mammalian haemoglobins. These amino acid residuesform an array of nearly parallel aromatic residues near the heme. In its ferric state, HbI has an extraordinary affinity for hydrogen sulfide( H 2 S) and is also capable of binding and forming stable complexes with nitric oxide, azide (N 3 )and cyanide in this oxidation state. Deoxy HbI has one of the fastestcombination rates with oxygen [O 2 ] among globins [11]. Theligand binding properties of HbI may be influenced by the lowdistal pocket polarity and the aromatic nature resulting from this array of phenylalanyl and aromatic residues. This particular arrangemento f residuesa tt he distal ligand-binding site is unusualf or haemoglobin and has not been observed before in the globin family [13]. Hence, acomparative analysis between of the native HbI and various HbI mutated systemsbased on computational techniques will yield relevant information concerning the relatives tructuralc hangeso fH bI and howt hese changesm ight affect ligand binding properties and selectivity.Specifically,the study intends to provide amolecular description of the effect of active site mutationsinthe deoxy species,which preludes the ligand migration and binding in the haemoglobin. Ther esults presented here suggest that PheE11 may be involved in structuralstability,aswell as ligand diffusion through the proteinwhile both GlnE7 and PheB10 are responsible for ligand binding and stabilityint he HbI heme moiety.
Allsimulations were performed usingthe GROMACS simulation package version 3.2.1 [19]. The GROMOS96 43a1 forcefield [20] was used to model the intramolecular proteini nteractions and the intermolecular interactions betweenthe protein and water molecules. The Single Point Charge (SPC) water modelw as used to describet he solvent molecules [21].The HbI molecule was placed in a cubicb ox with periodicb oundaryc onditions in all directionsand solvation was performed usingthe standard procedure of the GROMACS package.Aminimum distanceb etweent he protein and its periodic image was set to 2nmt op revent the interactions between them. Initially the energy of each system was minimised with the steepestd escent algorithm and 600 ps of position restrain dynamics. Thefi nal configuration of this procedure was used as starting point for the production MD runs. The simulations were carried out in aconstant temperature, pressurea nd number of molecules ensemble with at emperature of 298 K. AB erendsen thermostat with a coupling constant of 0.1 ps was used to maintain constant temperature, while the pressure was isotropically maintained constant by coupling it to aB erendsen barostat at 1atm, [22]. Long electrostatic interactions were calculated using the particle mesh Ewald algorithm [23] and at win range cut-off of 0.9/1.4 nm was applied for van der Waals interactions. Neighbour lists were utilised and updated everyfi fth integration step. All proteina nd water bond lengths were constrained usingt he LINCS and SETTLE algorithm, respectively [24,25]. Thes imulations were run for 10 ns,w hich was enough time to obtain proper convergence of the computed properties.The time stepfor all the simulations was set to 2fs. Initially, the data analysiswas performed usingdifferent average structures with variablen umbers of frames. However, it was found that the best average structure was obtained whenthe last 2nsofthe simulation wereused. All the trajectory analyses were performed using GROMACS, while the molecular graphic images were generated using visual molecular dynamicss oftware [26]. Root mean square deviation (RMSD)a nalyses were performed to evaluate proper convergenceo ft he system. The RMSDw as evaluated during the wholetrajectory by comparing each frame with the initial configuration, and the average values of RMSD were calculated for the last 2nsofthe trajectory. Figure2 shows the evolution of the C a atoms RMSDv alues with respect to the crystal structure for GlnE7His (panelA)and PheB10Leu/GlnE7His( panelB )s ystems,w here the convergenceo ft he systemsc an be observed owing the RMSDs mall fluctuation during the last 2ns [ 27]. Similar results were obtainedinthe othermutated systems.
To investigate the overall structural behaviour of the native and mutated HbI systems, the molecular dynamic  simulations were analysed based on the radius of gyration (Rg), hydrophobic/hydrophilic/totals olvent accessibility surface area (SASA), secondary structure evolution (SSE), distances analysis and root mean squarefl uctuations (RMSF).The RMSDanalysis was alsoused to describe the flexibility of the systems. Thel ast 2nso ft he simulations were used to create as tructure representing the average position of all the atoms. This structure was used to calculate Rg, SSE, SASA, RMSF and distances analyses of the systems. For each of the structural properties mentioned asingle valuefor each property was computedand reported as the average value. These structuralp roperties were analysed for the whole protein (all the atoms) and specific proteinsecondary structure domains (Ca atoms), that were select based on the changesinthe average structure of the native HbI after the 10 ns simulation. Time line diagrams were generated to determinet he time evolution of the proteinsecondary structure. Figure 3shows the averaged structure of native HbI after a 10 ns simulation, where the amino acids described by Torres-Mercado et al. [16] as comprising the molecule helical domains are labelled and highlighted to simplify the structural analysis. These helical domains are represented as follows: the residues in the segments HA (3 -18-blue), HB (20 -35-r ed), HC (36 -42-y ellow), HD (52 -57orange), HE (58 -77-tan), HF (88 -98-green), HG (103 -120 -c yan) and HH (121-140 -p urple). The random coil or turn segments can be observed in Figure 3c oloured in silver, as well as some residues that after the simulation became part of helices. The secondary structure domains assigned to the systems in the SSE structural analysis represent the conformation that the amino acid residues predominantly adoptedd uringt he last 2nso ft he simulation. As mentioned above, after the 10 ns simulation of the native HbI several residues initially comprised in random coil or b -turn conformations became part of helical domains;s pecifically,t he residue sequences 43 -48, 78-79, 82-87a nd 99-100. Also, it was observed that a 3-10-helix [117 -121] changed to form aregular a -helix.

Overallstructural analysis
The RMSDs howed no significant fluctuations around its average value during the last nanosecond of the simulations. Hence, as explained in the literature [27] oscillation aroundaconstant valueo ft he RMSD implies that the system has reached stable or metastables tates. The global RMSD analysiss howed that the systems mutated with valine (GlnE7Val and PheE11Val), including thed oublem utated system PheB10Leu/PheE11Val, exhibited ah igher degree of flexibilityc ompared to the others ystemsw ith RMSD values of 2.5, 2.9 and 2.8 A˚, respectively. The RMSD analysis by secondary structure revealed that helicesA-D and the b -turn EF (Figure 3) were more flexible in the mutated systemst han in the native protein. As will be discussed later,n os ignificant secondary structure changeswere observed in the HA and HB domains: hence the high RMSD observed is assigned to aw hole helicald omain displacement (translation). On the otherh and, the overall Rg and SASA analyses showed only slight differencesbetweenthe native HbI and the mutated systems. Specifically,t he native HbI had a lowerR gc omparedt ot he mutateds ystems,w ith differences not exceeding 5%. In all cases, when the native structure was compared to the mutated systems, the total hydrophilic SASA valuesi ncreased between2and 5%. Helices A-Dshowed increments in their hydrophilic SASAa nd ad ecrease in theirh ydrophobic SASA. The randomcoil or b -turn domains showed small changes in their hydrophobic/hydrophilic SASA,except for the CD segment, which exhibited at endency to decrease its hydrophobic SASA. Hence, we suggest that the increment in SASA corresponds primarily to an increase in the exposure of the proteinhydrophilic portionstothe solvent. In general, single or doublem utations induced as light increasei nt he SASA,p articularlyf or thes ystems involving valine mutationsa tp osition E11.
To further investigates econdary structural changes upon protein mutation,RMSF analyses were performed in native HbI, as well as the mutated systems. Figure 4shows the difference betweent he mutant HbI RMSF and the native RMSFf or the single mutated (a) and double mutated (b) systems. Interestingly, only the mobility of the chains in the 30 -60region showed asignificant increment upon mutation. The PheE11Val system againexhibited the higher degree of side chains flexibility, which extends to residues1 0-57. However, not all these residuesw ere involved in the unfolding process observed in this system. Similar results were obtainedw ith GlnE7Asn, which showed high flexibility in the side chains 42-48, but exhibited secondary structure changes involving only residues4 7-48. Hence, we suggestt hat the side chain rearrangementsi nduced by them utations are not necessarily related to ad irect change in the secondary structure of the residue,but to structure displacements that change ther esidue-residuea nd/orr esidue-solvent interactions. Table 2s ummarises the average distance between the C a of the mutated residue and the C a of all the residuesinthe affected proteinsegment exhibitingamajor change in the secondary structure domains (Table 1). For example, thes egment 43 -48w as affected by the PheB10Val; thus,w em easured the distancef rom the C a of Val28 to the C a for all the residuesc omprising the segment [Phe43,Ser44, Gly45, Leu46, Phe47 and Ser48]. After thesed istances were measured an average value was calculated, which is reportedi nT able 2. Our results showed that ar esidue change in the heme pocket moiety affected notonly the environment closetothe residue, but also the secondary structures ( Table 1) as far as 16.4 Aå way from the altered position. However, the mutations appearn ot to affect significantly the secondary structure domains close to the heme pocket.T he GlnE7Asn and GlnE7His induced secondary structuralc hangesa sm uch as 11.7 and 11.9 A˚away from the active site positions, respectively. Similarly, the PheB10Val and PheB10Leu affected the secondary structure as far as 10.9 and 11.1 A˚, respectively. In general,t he most dramatic effect was observed when the PheE11w as changed, leading to secondary structuralc hanges1 5-16.4 A˚away from the E11 position. The PheB10Leu/GlnE7His and PheB10Leu/PheE11Val systemsa ffected segments on an average distancefrom the B10 and E7 positions of 10.5 to 14.5 and 13.1 to 15.3 A˚,respectively. Interestingly, in the GlnE7His system its range of effects increased in the PheB10/Gl-nE7His system, while in the case of the PheE11Val the range of effect decreasedi nt he PheB10Leu/PheE11Val. Hence, it cannot be establishedh ow thee ffecto f individual alterations are combined in ad oublem utated system,because the GlnE7His and PheE11Val do nothave the same effect on residuesint he studied systems.
To understand the displacement effect of the mutations in the affected segments, the distances betweent he simulated native HbI molecule and the mutated systems were measured. Figure 5s hows as uperposition of the simulated native HbI (light gray) and the mutated protein (dark gray) where the displacements mentioned can be easily observed. The spheres represent the C a of specific residuesi nt he primary sequence; i.e. the C a of Gly45 in the native and mutated system. The distancemeasurement was performed always betweent he same residue in both molecules. The distanceanalyses revealed that the systems mutated with Asn,L eu and His showed changess maller than 3Å ,thus exhibiting no significant displacement from the native HbI. On the other hand, the valine mutations showed changes up to 9.1 A˚relative to the native structure.

Distal cavity fluctuation
To evaluate the effect of the changesi nt he heme pocket moiety, the distances betweent he heme pocket residues and the hemeiron were measured for each mutated system. Figure 6p lots the distanceb etweenh eme irona nd the heme pocket amino acidr esiduesi np ositions E7,E 11, B10 and CD1. Positions E11a nd CD1 did not change dramatically upon active site mutation except in the systemsinvolvingthe E11 position. Theaverage distance fluctuations were 3.4 and 1.1 A˚for the E11p osition, and 1.5 and 2Å for the CD1 position in the PheE11Val and PheB10Leu/PheE11Val, respectively. Displacement of CD1Phe has been observed in various ferrous myoglobin systemsd ue to as imultaneous movement of the distal residue E7 and the porphyrin macrocyclei nt he presence of an externall igand [31,32]. It has been shown recently [34] that upon ligand photodissociation the doming of the heme towards ad eoxy configuration, as well as the relaxation of the residue at position E7 provokea displacement of the CD1Phe away from the heme. Therefore, since no ligand was included in ourcalculations and the heme macrocycle was already in ad eoxy 'out of plane' configuration, movement of CD1Phe was not expected here. Nevertheless,r ecent molecular dynamic simulations of ferric unligated HbI showed the flexibility of CD1Phe to be influenced by the opening of the E7Gln gate [33].Their results suggestedthat whenthe E7Glnwas in its closed conformation CD1Phe moves towards the Fe III atom. However, as Figure6shows, movement of the residuesatthe E7 position away or towardthe ferrous iron did not provokeas ubstantial displacement of the CD1Phe.A sm entioned above, movement of CD1Phe in myoglobin is induced in part by the configuration of the  porphyrin group. Thus, differences in the electronic structure of unligated ferrous and ferric HbI can induce variation in their heme chromophore spin states and hence, their heme configurations, which may account for the variation of CD1Phe flexibility in both systems. Position B10a nd E7 on theo ther hand,w erev ery sensitivet oa ctives item utations.T he mutationsP he-B10Leu,G lnE7Hisa nd GlnE7Val inducedd isplacements exceeding2 .5 A˚compared to native HbI. ThePheB10Tyr, PheB10Vala nd PheE11Vali nduced displacementsi nt he B10positionof1.1,1.2 and3.0 A˚,respectively, compared to native HbI. Hence, theser esultsr evealedt hatt he E7 andB 10 positionsw eret he most affected upon mutation no matter them utated position.T hisb ehaviour couldb e observed in Figure 6w here theC a of residues at position B10and E7 show afluctuationrange of , 4.4Å ,while the E11and CD1exhibited alimited rangeofdisplacements for most of thesystems.Interestingly,these positionshavebeen identifieda sp laying an importantr olei nl igandb inding ands tabilisationprocesses.Rizzi et al. [ 34]suggestedthat ferric HbIs tabilisest he H 2 Sl igandb ymeans of hydrogen bondingwithE7Gln andthatthisresidue exhibitedahigher degree of flexibility in comparison to theE7His commonly foundinmammalian haemoglobins.Thisobservation ledto thep roposalt hatt he flexibility of E7Glni nf erricH bI facilitatest he exit of waterm olecules locatedi nt he heme pocket,a llowingh ydrogens ulphide bindingt oo ccur. Indeed,recentmolecular dynamicsimulations of ferric HbI in itsunligated stateconfirmthisinterpretation [33]. In this study, it wasobservedthatE7Gln interchanges an open and closed configuration within 200psa llowingt he entrance andstabilisation of theH 2 Sligand. As Figure 6shows,the mutationso fd istalr esiduesa tp ositions B10a nd E11 inducedl arge displacementso fG lnE7,s uggesting ah igh degree of flexibility fort hisr esidue in thed eoxy species. Thus,t hisp ositionm ight play an importantr ole in ligand bindinginthe deoxyHbI species, as well as in theferricone. Theseresults areconsistentwiththe hypothesis drawnfrom infrared experimental data,i nw hich afl exible E7Glnw as observed forb othf errica nd ferrousH bI complexes [ 35]. Moreover,i na ccordancew itho ur moleculard ynamic results, it is also suggestedthatnot only E7Glnwas involved in ligand reactivity in theHbI ferroussystems,but that Phe in theB 10 position wasa lsod irectlyi nvolved. This interpretation is furtherconfirmedbyanalysingthe binding kinetics of oxygen with variousHbI mutatedsystems [35]. Fori nstance, ourw orkd emonstrates that replacemento f E7Glnb ye itherA sn or Valp rovokedadisplacement of theser esiduesa sw ella st he residuei nt he B10p osition away from theh emei ronc entre, producinga no pendistal cavity that mayfacilitatethe entrance of ligands. Indeed,the associationr atec onstantsf or theseH bI mutateds ystems increasesfrom190 £ 10 6 M 2 1 s 2 1 forthe native protein, to 230 £ 10 6 and4 90 £ 10 6 M 2 1 s 2 1 fort he GlnE7Asn and GlnE7Val mutants, respectively,i ndicatingt hatt he open cavity observed in ourstudy allows O 2 to reactreadily with HbIh emei ronc entre. Ourr esults also suggestt hat substitution of B10Phe by Leuinduces adistaldisplacement similartothatobservedinGlnE7Asn, whilereplacement by apolar residuelikeTyr in that B10positionhas theopposite Figure 6. Distance between heme iron and heme pocket residues at positions E7, E11, B10 and CD1 for each mutated system. Molecular Simulation 723 effect.A sF igure6shows, in theP heB10Leu mutation an open cavity is produced,w hile in theP heB10Tyr ac losed heme distal site is formed duet ot he movement of the B10Tyr towardst he heme iron centre.O xygenb inding kinetics of theseHbI B10mutants indicate that forthe Leu system,t he associationc onstantr emains practicallyt he same as theG lnE7Asn, whilea3 0-fold decrease was observed for theT yr mutant (from 190 £ 10 6 to 6.8 £ 10 6 M 2 1 s 2 1 ). Hence, theset heoretical ande xperimentalr esults suggestas imilar pathway forl igand movement into thed istalc avityi nt he GlnE7Asn and PheB10Leum utants,a nd ab arrier forl igandb inding into theh eme iron centre fort he PheB10Tyrs ystemd ue to its movement towardsf errous centre.T aken together,t hese resultsd emonstratet he direct implicationo fE 7Gln and B10Phe in ligand bindinga nd stabilityi nt he ferrous HbIsystem.

Conclusion
In thiswork, we presented molecular dynamic simulations of native HbI and several HbI mutated systems solvatedby water.T he results were analysed based on the following properties: the SSE to characterised domain changes, RMSD, RMSF, Rg and SASA.D istancesb etweent he hemeiron to the heme pocket residues, as well as between the heme pocket residuest ot he residuesc omprising affected proteind omains upon mutation were measured to characterise the deoxy HbI systems. The active site mutationsd id nota ffect significantly the secondary structure of domains close to the heme pocket, such as helixes Band E, but had asignificant effect on domains far away from the active site. The mostaffected domains were the Ca nd Dh elicess uffering partial or total unfolding upon certainm utations,e nhancedi nt he mutations including position E11. In addition, RMSFa nalysis indicates that the PheE11Val system exhibits higher degreeo fs ide chains flexibility and segmentsd isplacements, which may affect significantly ligand diffusion through the proteinm atrixa si th as been observed before in humanhaemoglobin. These findingsstronglysuggesta direct role of residue at the E11p osition in the overall structural stability and hence, ligand internal movements in the native HbI protein.