Heat-resistant cytosolic malate dehydrogenases (cMDHs) of thermophilic intertidal snails (genus Echinolittorina): protein underpinnings of tolerance to body temperatures reaching 55°C

The physiological and biochemical determinants of organisms' environmental optima and tolerance ranges are an important focus of research in environmental physiology and global change biology (Somero et al., 2017). Efforts to identify adaptive variation in physiological and biochemical systems not only can help to account for current distribution patterns of species, but also can be helpful in predicting the consequences of global change. A major focal point in studies of biochemical adaptation to the environment has been on temperature effects on proteins (Somero, 2010; Fields et al., 2015). Adaptive variation in functional properties, for example, ligand binding, as indexed by the Michaelis–Menten constant (K_M), and structural stability, has been documented in many comparative studies of orthologous proteins from differently thermally adapted species (e.g. Fields et al., 2015). Studies that have extended these lines of analysis to elucidate the underlying molecular mechanisms of adaptation have shown that a relatively small number of amino acid substitutions – as little as a single substitution in some cases (Dong and Somero, 2009) – can account for adaptive differences in function and stability. These adaptive amino acid substitutions lead to alterations in weak bonding interactions (hydrogen bonds, ion-pair interactions and hydrophobic interactions) that alter stabilities of proteins (Scandurra et al., 1998; Xiao and Honig, 1999; Fitter and Haber-Pohlmeier, 2004; Dong and Somero, 2009). An important conclusion from these analyses is that proteins are under selection to maintain a balance between stability and flexibility (Kohen et al., 1999; Fields, 2001; Somero, 2010; Fields et al., 2015), which leads to corresponding states of protein structural stability at normal cellular temperatures of different species (Somero et al., 2017). Whereas this general stability–flexibility balance has been repeatedly demonstrated, the precise regions of the protein for which this balance is most important remain to be clearly elucidated. Identifying these regions was the primary objective of the present study. To attain our goal of identifying the regions of the protein where temperature-adaptive structural changes may be most important, we employed a method not heretofore used in comparative studies of temperature adaptation of enzymes: molecular dynamics (MD) simulations.

The study system we have chosen comprises two species of heat-tolerant intertidal snails, the congeneric littorine species Echinolittorina malaccana (Philippi 1847) and E. radiata (Souleyet 1852) that occur at different heights in the intertidal zone. The use of differently thermally adapted congeners allows a strong ‘signal-to-noise’ ratio in comparative analyses of proteins: only minor differences in amino acid sequence are commonly present among orthologs, which allows adaptively important structural changes to be readily identified. Our focus on intertidal species reflects the importance of thermal gradients in these habitats in governing vertical distribution patterns of organisms; congeneric species found at different heights commonly display adaptive differences in their physiological and biochemical properties (Somero et al., 2017). The rocky intertidal zone has been a common study system for examining the relationships between environmental characteristics, notably temperature, and adaptive variation among species. This environment is characterized by large differences in temperature owing to the tidal cycle, vertical position on rocky surfaces and hour-by-hour variation in physical conditions, e.g. solar heating (Helmuth et al., 2002; Walters et al., 2012). Many intertidal organisms encounter intensive thermal stresses, especially species with limited abilities to regulate body temperature behaviorally (Harley, 2011; Walters et al., 2012; Seabra et al., 2016). Previous studies have shown that many intertidal organisms currently live close to their upper thermal limits and thus are likely to be significantly challenged by future increases in temperature (Stillman, 2002; Somero, 2012; Dong et al., 2015; Marshall et al., 2015). Thus, the rocky intertidal zone is an ideal model system for investigating the relationships between conditions of thermal exposure and adaptation, as well as for elucidating the susceptibility of organisms to the impact of climate change (Mislan et al., 2014).

For studies of marine invertebrates, cytosolic malate dehydrogenase (cMDH; EC 1.1.1.37, l-malate: NAD⁺ oxidoreductase) has been an especially important experimental system for elucidating adaptive patterns among orthologous proteins and analyzing how differences in amino acid sequence contribute to adaptive variation in function and structural stability (Dahlhoff and Somero, 1993; Fields et al., 2006; Dong and Somero, 2009). Functional studies have focused strongly on ligand binding, as indexed by the K_M for the cofactor reduced nicotinamide adenine dinucleotide (NADH) (K_M^NADH). Ligand binding is highly sensitive to temperature, and at temperatures near the upper thermal tolerance limits of a species, sharp increases in K_M^NADH are generally observed. These rapid decreases in binding ability have been interpreted as arising from temperature-induced changes in enzyme conformation that alter the geometry of the ligand-binding sites (Somero et al., 2017). Orthologous enzymes of warm-adapted species not only can sustain appropriate ligand binding activities at higher temperatures, but also exhibit a relatively high resistance to thermal denaturation of overall structure. Thus, it is generally accepted that one of the most common strategies to enhance protein thermal stability is to increase the rigidity of the protein (Giver et al., 1998; Van den Burg et al., 1998). However, how the overall stability of a protein affects the thermal sensitivity of the conformation of ligand binding sites remains to be fully understood. The predicted differences in net free energy of protein stabilization resulting from the minor variation in sequences among closely related congeners are small. The effects of minor variations in sequence on stability may not be great enough to allow a priori predictions of whether an amino acid substitution at a particular site will lead to a predictable change in thermal stability (Fields et al., 2015).

A path forward for this type of analysis might involve computer-based molecular dynamics methods that can simulate the rapid movements of amino acid side-chains over time (Phillips et al., 2005). For example, the amount of variation in the position of a residue or side-chain, as indexed by the root mean square deviation (RMSD) in its position, could provide important insights into the sensitivity of localized regions of a protein, such as the binding site, to changes in temperature.

Here, we employ MD analysis in conjunction with more conventional measurements of temperature sensitivity of ligand binding and overall structural stability in studies of cMDH orthologs of littorine gastropods of the genus Echinolittorina, which are among the most heat-tolerant animals known. Echinolittorina snails have the highest heat tolerance known for intertidal species (Marshall et al., 2011; Marshall and McQuaid, 2011; Wang et al., 2014), making these species excellent subjects for studies of adaptation to high temperature. Littorine gastropods are among the dominant macrobenthos on the high intertidal shore throughout temperate and tropical oceans (Reid et al., 2006, 2007). They experience prolonged emersion periods during neap tides, as well as extreme environmental temperatures over 50°C in the case of Echinolittorina species (Williams and Morritt, 1995; Uglow and Williams, 2001; Marshall et al., 2010). Snails of the genus Echinolittorina are distributed widely on the high intertidal shore and exhibit high upper thermal limits (>55°C) and a series of behavioral, morphological, physiological and biochemical responses to cope with thermal stresses (Marshall et al., 2011; Wang et al., 2016; Li, 2012; Wong, 2013). Echinolittorina malaccana is a common species in the splash zone and high intertidal zone in southern China and Southeast Asia, reaching higher tidal levels than any other sympatric congeners (Reid et al., 2006; Marshall and Chua, 2012). The mean body temperature of E. malaccana was recorded up to a range between 41 and 44°C in August/September 2009 on artificial seawalls at Pantai Tungku in Brunei Darussalam (4°32′N, 114°43′E) (Marshall et al., 2010). The congener species E. radiata is distributed along the entire Chinese mainland coast and has a vertical distribution that lies below the range of E. malaccana (Ohgaki, 1985; Williams, 1994; Ma, 2004; Reid et al., 2006). The high thermal tolerance of the species has been demonstrated in several types of studies. In studies of whole-organism and cardiac thermal tolerance, no mortality was observed after 24 h exposure at temperatures between 42 and 46°C for E. malaccana, and the values of Arrhenius breakpoint temperature for heart rate of E. malaccana and E. radiata were 46.5–48.2°C (Marshall et al., 2011; Li, 2012). For E. malaccana, the temperature at which 50% of a population died (LT₅₀) was 56.3–57.1°C (Wang et al., 2016); E. radiata showed a slightly lower upper thermal limit of 55.5°C (Li, 2012).

Here, we show that one of the underpinnings of this high thermal tolerance is the extreme heat stability of the congeners' cMDHs. Further, we show that only two differences in amino acid sequence may explain the inter-congener differences observed in ligand binding, protein denaturation kinetics and molecular flexibility of localized regions that are critical to binding and catalysis.

Specimens of E. malaccana and E. radiata (body length of ∼0.5 cm) were collected from a rocky shore habitat in Xiamen, China (24°28′N, 118°05′E), in May 2016 and transported live to the State Key Laboratory of Marine Environmental Science, Xiamen University (Xiamen, China) within 2 h. All individuals were held at 16°C in aquaria with no emersion for at least 48 h before euthanization.

A total of approximately 1 g of foot muscle was dissected from animals held on an ice-cooled tray; approximately 100 specimens of each species were need to obtain a muscle sample of this size. Muscle was homogenized in 10 volumes of ice-cold potassium phosphate buffer (50 mmol l^–1, pH 6.8 at 4°C) with a rotor-stator homogenizer (FLUKO, Shanghai, China) in three 10 s bursts at maximum speed (35,000 g). The homogenate was centrifuged at 12,000 g for 30 min and the supernatant was removed and heated at 60.5°C for 3 min. The mitochondrial isoform of MDH (mMDH) is more thermally labile than the cytosolic paralog and was fully eliminated by this heat treatment (Fig. S1). cMDH in the heated supernatant was precipitated by the addition of ammonium sulfate to give a final concentration of 45% saturation. The sample was stirred at 4°C for 15 min on a mechanical stirrer (IKA, Guangzhou, China), held on ice for 30 min and then centrifuged at 12,000 g for 30 min. The supernatant was decanted and brought to 80% saturation with ammonium sulfate. After a further 15 min of stirring and 30 min of standing on ice, the sample was centrifuged again at 12,000 g for 30 min. The supernatant was discarded and the pellet was resuspended in a potassium phosphate buffer (50 mmol l^–1) and dialyzed against this buffer overnight at 4°C.

The kinetics of cMDH orthologs were measured by a modified method based on previous studies (Fields et al., 2006; Dong and Somero, 2009). Values of K_M^NADH were determined at five temperatures (20, 25, 30, 35 and 40°C). Activities were determined using a UV-1800 spectrophotometer (Shimadzu, Kyoto, Japan) equipped with a temperature-controlled cell attached to a water bath (TIANHENG, Ningbo, China). Temperature of the cuvette was maintained to within ±0.1°C of designated temperatures. An imidazole chloride buffer (200 mmol l^–1, pH 7.0 at 20°C) was used to ensure that the pH level in the experimental system corresponded to the intracellular value over the range of measurement temperatures (Yancey and Somero, 1978). For each K_M^NADH determination, five concentrations of NADH were used: 20, 40, 75, 100 and 150 μmol l^–1. The co-substrate oxaloacetic acid (OAA) was present at a starting concentration of 200 μmol l^–1. K_M^NADH values were calculated from initial velocity measurements at different concentrations of NADH using a non-linear least-square fit to the Michaelis–Menten equation with Prism 5.0 software (GraphPad Software, La Jolla, CA, USA).

Optimal temperatures for cMDH activities were determined by a method similar to that described above. The determinations were conducted at six temperatures (25, 30, 35, 40, 45 and 50°C) in a reaction mixture containing 200 mmol l^–1 imidazole-HCl (pH 7.0 at 20°C), 150 μmol l^–1 NADH and 200 μmol l^–1 OAA. Enzyme activity was normalized against the highest activity for each ortholog. The relationship between enzyme activity and temperature was analyzed by a second-order polynomial regression analysis model.

Thermal stabilities of cMDH orthologs were characterized as the residual activities following incubation at different temperatures for different times (Fields et al., 2006; Dong and Somero, 2009). After overnight dialysis, enzymes were diluted to equivalent activity and incubated at 47.5, 52.5 and 57.5°C. Samples of each enzyme were transferred to ice at t=0, 10, 20, 30, 45 and 60 min and activity was determined (n=5) at 25°C. Residual activity was defined as the ratio between the mean activity at time t and the mean activity at time 0.

Total RNA was purified from foot muscle (n=5) using Trizol reagent (Invitrogen, Carlsbad, CA, USA). Reverse transcriptase (RT) reactions were performed using PrimeScript™ RT reagent kits (TaKaRa, Dalian, China). Partial sequences of the 5′ ends of cDNA of cMDH in E. malaccana (EmcMDH) were obtained based on a de novo transcriptome sequencing (Wang et al., 2014). 3′ gene-specific primers (GSP-Ma-1 and GSP-Ma-2; Table S1) were designed based on the partial sequences above. The full-length cDNA of EmcMDH was obtained using the rapid amplification of cDNA ends (RACE) protocol with the primers above. PCR (94°C for 3 min, followed by 35 cycles of 94°C 30 s, 57°C 1 min and 72°C 2 min, with a final 10 min extension at 72°C) was used to amplify the 3′ ends of the cDNA. The target product showing a single band of expected size was visualized in 1.2% agarose gel and sent to a commercial company for sequencing (Invitrogen Biotechnology Co., Shanghai, China). Based on the open reading frame (ORF) sequences of EmcMDH, two pairs of primers (MaRa-F1, MaRa-F2, MaRa-R1 and MaRa-R2) were designed to amplify the full-length cDNA coding region of cMDH in E. radiata (ErcMDH) using PCR as above (Table S1). The cDNA sequences were assembled using DNAMAN software (Lynnon Biosoft, San Ramon, CA, USA) and the deduced amino acid sequences were aligned using the ClustalX2 algorithm (Larkin et al., 2007).

Three-dimensional models were constructed based on a template structure file in PDB format of porcine (Sus scrofa) cMDH (5MDH, Chain-A) (Chapman et al., 1999), with Swiss-Model software (Biasini et al., 2014). The sequence identities between 5MDH-A and EmcMDH and ErcMDH were 70.91 and 70.61%, respectively. The computations of hydrogen bond content and distance between atoms were processed with Swiss-PdbViewer 4.1.0 (Guex et al., 2009) and the visualizations of three-dimensional structures were processed with Visual Molecular Dynamics (VMD) program (Humphrey et al., 1996). The verification of protein structure was performed using PROCHECK (Laskowski et al., 1993), Verify3D (Bowie et al., 1991) and ERRAT (Macarthur et al., 1994) procedures.

MD simulation covers length scales ranging from the atomic level to entire microstructures, and produces a complete trajectory history of different protein properties. These include thermodynamic features such as the equation of state and dynamical properties such as structural characteristics, of which the radial distribution function is the simplest (Rapaport, 1999).

MD simulations of cMDH orthologs were assigned at 42 and 57°C for 5 ns. All simulations were performed in the isobaric–isothermal (NPT) ensemble, using the NAMD simulation software (Phillips et al., 2005) with a CHARMM force field (MacKerell et al., 1998). Computed three-dimensional structures of EmcMDH and ErcMDH constructed above were used as the starting models of the simulations. Each starting structure was solvated in an independent water box to create simulation conditions that more closely resembled the cellular environment. The size of the water box was created with a layer of water 5 Å in each direction from the atom with the largest coordinate in that direction. Local interaction distance common to both electrostatic and van der Waals calculations (cutoff) and the distance between pairs for inclusion in pair lists (pairlistdist) were 12 and 16 Å, respectively. The SHAKE algorithm (Ryckaert et al., 1977) was used to constrain bond length with a time step of 2 fs. Pressure and temperature were maintained by the Langevin piston and Langevin dynamic method at 1 bar and the set temperature, respectively. Long-range interaction was applied using the particle mesh Ewald method. After energy minimization using the steepest descent algorithm for 5000 steps, each system was assigned velocity from the Maxwell–Boltzmann distribution in triplicate. Trajectories of the structures were collected every 0.002 ns.

where r_i represents the position at time i, r₀ represents the reference value, r represents the average value of the root mean square fluctuation and N represents the number of atoms.

Data were analyzed using SPSS statistical software (Version 11, IBM, Armonk, NY, USA). The difference in biochemical characteristics between EmcMDH and ErcMDH was analyzed using an independent sample t-test. Differences in K_M^NADH and residual activity at different temperatures were analyzed using one-way ANOVA followed by a Duncan post hoc multiple range test. To test for significant differences in the relationship of K_M^NADH versus temperature and temporal changes of residual enzymatic activity, the slopes of regressions were analyzed using an analysis of covariance (ANCOVA) followed by a least significant difference (LSD) multiple comparisons test (P=0.05).

The K_M^NADH values of both cMDH orthologs increased with increasing measurement temperature (Fig. 1). Relative to ErcMDH, EmcMDH had a significantly smaller increase in K_M^NADH with increasing measurement temperature at 35 and 40°C (F_8,7.34=0.10, P=0.002; F_8,7.62=0.12, P=0.002, respectively). The slope of K_M^NADH versus temperature for EmcMDH was significantly lower than that for the ErcMDH ortholog (F_2,46=52.06, P<0.001; Fig. 1, inset).

Enzyme activity was normalized against the activity at 35 and 40°C (the temperature of the highest activity for each ortholog) for E. malaccana and E. radiata, respectively. Similar patterns in activities were observed between EmcMDH and ErcMDH orthologs at different temperatures (Fig. 2). The optimal temperature for both EmcMDH and ErcMDH activities was 36°C. At temperatures between 24.5 and 47.2°C, the normalized activities remained at least 60% of maximum.

Both EmcMDH and ErcMDH had high thermal resistance (Fig. 3). After 60 min incubation at 47.5°C, EmcMDH retained 95.0% of initial activity, whereas that of the ErcMDH ortholog was 92.6%. After 60 min incubation at 52.5°C, EmcMDH still possessed 65.5% residual activity, compared with 59.9% for the ErcMDH ortholog. After incubation at 57.5°C for 60 min, ErcMDH was almost fully inactivated, whereas the residual activity of the EmcMDH ortholog was 11.7%. The slopes of the regressions of residual activities at different times of incubation versus temperature of ErcMDHs were significantly steeper than those of EmcMDH orthologs at all treatment temperatures (F_2,75=30.90, P<0.001 for 47.5°C; F_2,73=340.22, P<0.001 for 52.5°C; F_2,55=980.83, P<0.001 for 57.5°C).

The sequences of the two cMDH cDNAs shared high identity. The 999 bp coding regions of genes for EmcMDH (thermally tolerant) and ErcMDH (thermally sensitive) had 42 synonymous and four nonsynonymous substitutions (Fig. S2), which resulted in two amino acid differences, those at residues 48 and 114 (Fig. 4). At those sites, EmcMDH had glycine (Gly-48 and Gly-114), whereas ErcMDH had serine (Ser-48 and Ser-114). The accuracies of the molecular modeling for EmcMDH and ErcMDH were verified (Table S2). The two amino acid substitutions occurred in two alpha-helices, αC and αDE (Fig. 5). The Ser-114 in ErcMDH is predicted to form an additional hydrogen bond, that between O Lysine-110 (Lys) and OG Ser-114, relative to EmcMDH (Fig. 6, Table 1, bottom line).

RMSD values were determined to evaluate overall protein flexibility at different temperatures, as reflected by fluctuation in backbone atom position; the higher the RMSD value, the more flexible the structure is overall. The MD simulations of the cMDH orthologs at 42 and 57°C reached equilibrium RMSD values in approximately 1.0 and 0.7 ns, respectively (Fig. 7). The equilibrium (1–5 ns) RMSD values of EmcMDH and ErcMDH at 42°C were approximately equal to 2.1 Å. At 57°C, the RMSD value of the more heat tolerant EmcMDH remained essentially constant (1.99 Å), whereas that of the ErcMDH ortholog increased up to 2.44 Å (Table 2). The equilibrium RMSD value of ErcMDH at 57°C was significantly higher than that of EmcMDH at 57°C (F_4,2.81=1.69, P=0.001), indicating a more flexible overall structure for ErcMDH.

Conformational stability for individual residues was indexed by measuring RMSF values; a higher RMSF value indicates a greater flexibility for a given residue. The overall equilibrium (1–5 ns) RMSF value of the more heat-stable EmcMDH remained near 1.70 Å at 42 and 57°C, whereas it increased from 1.68 to 1.88 Å for the more heat-sensitive ErcMDH ortholog (Table 2), consistent with a less stable, more flexible structure for ErcMDH. Consistent with this observation, at 57°C the equilibrium RMSF value of ErcMDH was significantly higher than that of the EmcMDH ortholog (F_4,3.48=1.17, P<0.001; Table 2).

The relative conformational flexibility and its temperature dependence became more complicated when local fluctuation in structure was considered. Local changes in flexibility can in some cases be linked to the conformational status of regions of the enzyme known to be involved in substrate binding and catalysis (Fig. 8). The more heat-tolerant EmcMDH ortholog largely preserved its geometry at the two different simulative temperatures, as shown by a relatively consistent degree of fluctuation in most regions of the protein at both temperatures. In contrast, ErcMDH showed a more pronounced temperature dependence. Two regions of the sequence of known functional significance showed differences between the two orthologs. One region was located at residues 90‒100, which comprises the substrate binding site and the active site loop. The latter region of the sequence collapses during substrate binding to form the catalytic vacuole, the functional configuration of the enzyme that permits catalysis. At 42°C, the RMSF value of ErcMDH at this region was notably lower than that of the EmcMDH ortholog, a difference that correlates with a weaker substrate binding affinity in the more heat-sensitive ErcMDH ortholog. This inter-ortholog difference suggests that higher conformational flexibility of this region may facilitate stronger binding, albeit the mechanistic basis of this conjectured effect remains to be examined. Another flexible region with a significantly larger RMSF value compared with the overall equilibrium value occurs at residues 230‒245, a region involved in conversion of substrates to products. This region exhibited similar movement in the two orthologs at 42°C. However, the RMSF values of this region in ErcMDH increased greatly when the simulation temperature increased to 57°C, relative to the change observed for EmcMDH. Finally, the RMSF value of residues 300–310, which occur near the C terminus of the protein, showed a larger increase with rising temperature in ErcMDH compared with EmcMDH.

Biochemical adaptation of proteins to temperature has been documented for a variety of proteins across a diversity of ectothermic taxa (Tattersall et al., 2012; Fields et al., 2015; Lim et al., 2016). A key factor that underlies these temperature-adaptive differences is the need for an appropriate degree of structural stability at normal cellular temperatures (Somero et al., 2017). Thus, proteins are under selection to maintain a balance between stability and flexibility. Stability must be sufficient to ensure that the protein does not unfold and potentially aggregate under normal physiological conditions. In contrast, flexibility must be sufficient to allow the conformational changes necessary for catalysis or regulation to be feasible at cellular temperatures (Fields, 2001; Somero, 2010; Fields et al., 2015). As inhabitants of the high intertidal zone on tropical rocky shores, Echinolittorina snails attain body temperatures higher than almost all other eukaryotes, so proteins of these species may provide important windows into how eukaryotes withstand extreme heat stress and reach the appropriate stability/flexibility balance in their structures.

One of the kinetic properties of enzymes that is an index of thermal stability of function is the Michaelis–Menten constant for substrates and cofactors (Somero et al., 2017). Values of K_M^NADH of most cMDH orthologs in marine invertebrates are maintained within a relatively narrow range at physiological temperatures, but generally increase rapidly above a certain high temperature (Dahlhoff and Somero, 1993; Fields et al., 2006; Dong and Somero, 2009). The K_M^NADH values of the two cMDH orthologs from Echinolittorina congeners at their habitat temperatures are similar to those measured in mesophilic invertebrates at their lower habitat (body) temperatures (Dahlhoff and Somero, 1993; Fields et al., 2006; Dong and Somero, 2009).

Protein stability of the two orthologs was high, as indicated by the thermal optima of activity (Fig. 2) and resistance to heat denaturation (Fig. 3). Activities of EmcMDH and ErcMDH were maintained at high levels (≥60% of maximal activity) at temperatures between 25 and 45°C. The optimal temperature for both EmcMDH and ErcMDH activities was near 36°C. The thermal stabilities of cMDH orthologs from Echinolittorina congeners were notably higher than those of other intertidal species (Fields et al., 2006; Dong and Somero, 2009). The residual activities of Echinolittorina cMDHs were approximately 95% after 60 min incubation at 47.5°C, whereas the residual activities of mussel (genus Mytilus) cMDHs were only 10% after 60 min incubation at 42.5°C (Fields et al., 2006). The high residual activities of cMDH orthologs provide a partial explanation for the high thermal tolerance of the genus Echinolittorina. Even when incubated at temperatures above the Arrhenius breakpoint temperature (52.5°C), the Echinolittorina cMDH orthologs still retained substantial residual activities (>60% for 60 min). Thus, the present study shows that the high thermal stabilities of structure and ligand binding (K_M^NADH) of cMDHs of Echinolittorina congeners are likely to contribute to the species' abilities to thrive in their extreme thermal environment.

To further explore the relationship between structure and function of these two cMDH orthologs, we examined the dynamic behavior of structure at different temperatures using an MD simulation method (Phillips et al., 2005). Low RMSD values, i.e. values below approximately 3‒4 Å, indicate a high structural rigidity as manifested by slight amounts of fluctuation in conformation at the measurement (simulation) temperature. Previous studies of proteins of thermophilic organisms have shown that these species' proteins are characterized by high thermal stability and low levels of structural fluctuation in MD simulation studies (Li et al., 2005). For example, the equilibrium RMSD value of thermophilic isopropylmalate dehydrogenase from Thermus thermophilus was 3.5 Å at 64°C (Sharma and Sastry, 2015); the thermophilic β-fructosidase from Thermotoga maritima reached the stable state with a value of 1.97 Å at 80°C (Mazola et al., 2015). The hyperthermophilic amylase from Bacillus licheniformis is thermostable at very high temperature; the native state of the protein and enzymatic activity is strongly conserved at temperature below 90°C, whereas 50% of the molecules are unfolded at 103°C (Fitter et al., 2001; Fitter and Haber-Pohlmeier, 2004; Fitter, 2005). In concert with these observations, RMSD values increased from approximately 3 Å at 27°C to ∼10 Å at 127°C, an indication of complete unfolding of the protein (Søgaard et al., 1993; Haki and Rakshit, 2003; Lu, 2009). Additional linkages between protein stability and RMSD values has been found in studies that employed site-directed mutagenesis. The thermal stability of an already heat-stable enzyme, isopropylmalate dehydrogenase, from the thermophilic bacterium T. thermophilus was further increased by introducing mutations that reduced the enzyme's conformational fluctuations (Sharma and Sastry, 2015).

In the present study, the low RMSD values (2.09–2.44 Å; Table 2) of cMDH orthologs of the Echinolittorina congeners are consistent with the observed thermostability of the proteins (Fig. 3). At a simulation temperature of 42°C, which lies near typical body temperatures seen in summer, low RMSD values were found for both orthologs, consistent with their observed thermostability of function. Whereas both cMDH orthologs were thermally stable at typical body temperatures, a significant change in the RMSD value of the ErcMDH ortholog was observed when the enzyme was examined at 57°C, a potentially lethal temperature for this species. In comparison, even at the highest MD simulation temperature, the rigidity of the more heat-stable EmcMDH ortholog was maintained at a high level, with a mean RMSD value of 1.99 Å. A similar phenomenon of high rigidity (2.7 and 2.8 Å for temperatures of 27 and 100°C, respectively) was observed in hyperthermophilic rubredoxin from an archaeon (Pyrococcus furiosus); intrinsically low RMSD values were approximately constant over the whole range of simulation temperatures used, 27 and 100°C (Grottesi et al., 2002).

The structural flexibility revealed by RMSF and its temperature dependence revealed additional details about the proteins' structural states when fluctuations in local regions of the sequence were considered. For example, the RMSF value of the C-terminal region of ErcMDH increased at 57°C. Similar high flexibilities of N- and/or C-terminal regions have been observed in several proteins examined at different MD simulative temperatures (Grottesi et al., 2002; Tang and Liu, 2007; Skjaerven et al., 2011). One relevant stabilizing feature discovered in thermophilic 3-isopropylmalate dehydrogenase from T. thermophilus was the shortening of N- and C-terminal regions (Wallon et al., 1997). However, this type of putative adaptive change was not observed in the Echinolittorina orthologs, which have the same number of amino acid residues in their sequences. The functional significance of the enhanced flexibility of the C-terminal region of the ErcMDH ortholog is not known.

Two other relatively flexible regions in the two cMDH orthologs, those comprising residues 90–100 and residues 230–245, contain the residues involved in ligand binding and catalysis (Chapman et al., 1999). Binding involves residues 92–98, the active site loop comprises residues 91‒101, and the region important for subunit–subunit interactions comprises residues 230, 238, 243, 244 and 248 (Chapman et al., 1999; Dong and Somero, 2009). These are sites where evolutionary trade-off between activity and stability may be especially important (Tattersall et al., 2012; Zeiske et al., 2016). A potential general pattern in temperature adaptation of enzymes was proposed that conjectured that the sensitivity of an enzyme to temperature and, thereby, the intensity of selection for temperature adaptation of structure depends on the degree of conformational change that the enzyme undergoes during catalysis (Lockwood and Somero, 2012). Thus, selection would be expected to favor conservation of local flexibility in these regions of the protein across species (corresponding states of structure) and, for a given species, across its normal range of cellular temperatures. The differences in temperature effects on RMSF values between the Echinolittorina orthologs, notably in residues 230‒245, may be a reflection of this pattern of selection. A similar phenomenon was observed in ribonuclease HI enzymes from thermophilic T. thermophilus and mesophilic Escherichia coli. In these enzymes, the regions containing many basic and aromatic residues important for catalysis were the most flexible, as shown by relatively higher RMSF values (2–5 Å) compared with other regions of the proteins (Tang and Liu, 2007). In summary, adjustments in the flexibility of regions in and around the active sites of enzymes are proposed to favor conservation of adequate flexibility of structure to allow effective function at low temperatures, yet ensure adequate conformational stability at higher temperatures.

Because the vertical distribution of E. radiata extends only up to the lower part of the range of E. malaccana (Mak and Williams, 1999; Li, 2012), E. radiata experiences lower temperatures and, correspondingly, has a slightly lower thermal tolerance than E. malaccana (Marshall et al., 2011; Li, 2012). The kinetic and structural characteristics of the congeners' cMDHs reflect these differences in thermal exposure and tolerance. The cMDH of E. radiata has higher K_M^NADH values at temperatures of 35°C and above, and more rapid loss of activity during heating than the ortholog of its congener. These observations can be interpreted in part by the results of the MD simulations, which showed that ErcMDH had a significantly greater overall flexibility, as indexed by the RMSD value, at higher temperature than EmcMDH. The difference in temperature adaptation between EmcMDH and ErcMDH orthologs is consistent with that seen in dehydrogenase enzymes in other molluscan species that occur at different heights in the intertidal zone and/or at different latitudes (Dong and Somero, 2009; Lockwood and Somero, 2012; Fields et al., 2015).

Previous studies have demonstrated that minor amounts of amino acid substitution can be sufficient to shift the thermal characteristics of enzymes (Fields et al., 2006; Dong and Somero, 2009). For example, differences in temperature adaptation of K_M^NADH and structural stability between cMDH orthologs in the congeneric limpets Lottia digitalis and L. austrodigitalis result from a single amino acid substitution, a serine (L. austrodigitalis) for a glycine (L. digitalis), which enabled additional hydrogen bonds to form in the subunit of the warm-adapted species (L. austrodigitalis) (Dong and Somero, 2009). However, this pattern of amino acid substitution does not hold for the relationship between structural stability and thermal tolerance in the Echinolittorina cMDHs. Here, in the more heat-tolerant EmcMDH, glycine residues replaced serines at the two sites where sequence differences occurred. This finding indicates that the effects of a single substitution cannot be fully interpreted by examining only localized weak bonding involving the residue in question (here, residue 114), but must include a broader analysis of changes in structure at other regions of the protein. MD simulation approaches may provide an important means for such a holistic analysis. The MD simulation analyses show that the more heat-resistant EmcMDH ortholog has, in fact, a higher intrinsic stability than the ortholog of its congener, even though one might predict a priori that the substitution at site 114 would favor higher stability in ErcMDH. The RMSD and RMSF values of EmcMDH were highly preserved during simulation at different temperatures, whereas the values for more heat-sensitive ErcMDH increased with temperature. The difference in structural thermal sensitivity between the EmcMDH and ErcMDH orthologs is consistent with the temperature adaptation observed in β-fructosidase of T. maritima (Mazola et al., 2015). The active site geometry of this thermophilic β-fructosidase was strongly conserved during simulations at different temperatures, whereas an increase in RMSD value was observed for a mesophilic ortholog at elevated temperature.

The high thermal stabilities of function and structure of the cMDHs of Echinolittorina snails may contribute to the survival of these species at temperatures near the upper thermal limits of eukaryotes. The adaptive variation in stability of ligand binding and structure is founded on only two differences in amino acid sequence: substitution of glycine for serine at two positions in the cMDH of the more warm-adapted species. Despite the fact that the glycine to serine substitutions lead to a predicted reduced hydrogen bonding near the substitution site (114) in E. malaccana, overall protein stability indexed by heat denaturation in vitro and predicted by in silico MD simulations is higher in the ortholog of this species. The differences in conformation arising from the glycine/serine substitutions likely have important effects on local structural flexibility that stabilize overall protein structure and ensure the enzyme's abilities to bind ligands at active sites and carry out catalysis.