Temperature adaptations of the thermophilic snail Echinolittorina malaccana: insights from metabolomic analysis

The fact that temperature has pervasive effects on biological structures and rates of physiological activity underlies the importance of elucidating the fundamental mechanistic determinants of thermal optima and thermal limits, which play crucial roles in determining species distribution patterns and sensitivities to anthropogenic climate change (Somero, 2012; Somero et al., 2017). Although much progress has been made in clarifying the bases of biochemical and physiological adaptation (Angilletta, 2009; Cossins and Bowler, 1987; Somero et al., 2017), particularly at the level of individual proteins (Dong et al., 2018; Liao et al., 2017; Somero, 2004), metabolic pathways (Noda-Garcia et al., 2018; Götze et al., 2020; Xu et al., 2020), and membrane structure and function (e.g. homeoviscous and homeophasic adaption; Hazel, 1995; Logue et al., 2000; Somero et al., 2017), there remain many unanswered questions about the mechanisms that set the thermal limits for life. Among these are questions about the types of adaptive responses available to species that withstand extremes of high temperature. Studies of these species could shed important light on the types of mechanisms that will be important for ectotherms as they attempt to cope with a warming world. Thus, issues of adaptability to temperature are becoming increasingly important as human-induced climate change is increasing the global mean temperature and the frequency of extreme temperature events, leading to changes in biodiversity and biogeographic distributions (Deutsch et al., 2020; Dong et al., 2015; Parmesan, 2006; Parmesan and Yohe, 2003).

Here, we focus on the thermal physiology of one of the most heat-tolerant eukaryotes, the rocky intertidal periwinkle snail Echinolittorina malaccana (Marshall et al., 2011; Liao et al., 2017). The upper lethal temperature of this species is near 55°C, so the species can serve as an appropriate study organism for investigating mechanisms of extreme thermal tolerance in animals from rocky intertidal environments. The rocky intertidal ecosystem is one of the most vulnerable ecosystems on Earth to climate change (Helmuth et al., 2006). In response to the extreme thermal stress on the shore, intertidal organisms have developed adaptive responses in terms of behavior, physiology and protein evolution, and have become a ‘model’ system for exploring thermal adaptive mechanisms and for studying the impacts of climate change on different levels of biological organization, from the cellular to the biogeographic scale (Helmuth et al., 2006; Kuo and Sanford, 2009; Menge et al., 2008; Mieszkowska et al., 2005; Ng et al., 2017; Somero et al., 2017).

The present study builds on earlier work with rocky intertidal invertebrates that has revealed the importance of cardiac physiology and metabolic reorganization in setting thermal tolerance limits (Stenseng et al., 2005; Dong et al., 2017). These two phenomena are closely linked: sustaining cardiac activity is essential for allowing energy metabolism [adenosine triphosphate (ATP) generation] to occur through efficient aerobic processes (Marshall et al., 2011; Marshall and McQuaid, 2011). When high temperatures compromise cardiac function, it may be necessary to shift ATP production to less efficient anaerobic pathways and greatly reduce rates of ATP turnover (metabolic depression) (Anestis et al., 2010; Han et al., 2017). Metabolic depression is adaptive in that it allows extended survival at thermal extremes by conserving energy stores while allowing adequate provision of ATP for essential functions (Marshall et al., 2011; Marshall and McQuaid, 2011; Pörtner, 2001; Sokolova et al., 2012). The shift in metabolic balance that occurs with entry into metabolic depression is revealed by changes in metabolic end-products that are produced. For example, in Ostrea edulis, levels of succinate, malate and fumarate (biomarkers of anaerobic metabolism) were significantly elevated under heat stress (Eymann et al., 2020).

Sustaining the functions of organs and the metabolic pathways that power their activities necessitates adequate stability in cellular proteins and membranes, two types of biochemical structures that are intrinsically highly sensitive to thermal perturbation (Somero et al., 2017). Membrane stability is modulated largely through adjustments in lipid composition, including acyl chain saturation and head group composition (Hazel, 1995). Stabilization of the cellular proteome under thermal stress is achieved by three primary mechanisms: evolutionary adjustment of intrinsic protein stability (Liao et al., 2017, 2019), reversal of heat-induced damage to proteins through activities of molecular chaperones known as heat-shock proteins (Feder and Hofmann, 1999; Sørensen et al., 2003), and upregulation of protein-stabilizing organic solutes (osmolytes) that can offset thermally induced denaturation (Yancey et al., 1982; Yancey, 2005; Somero et al., 2017). The role of protein-stabilizing organic osmolytes is well understood in thermophilic bacteria and archaea (Somero et al., 2017), but remains largely unexplored in marine invertebrates. These species are iso-osmotic with seawater, and a major share of the osmotically active solutes in their cells are small organic molecules, such as free amino acids and their metabolic breakdown products (e.g. taurine) and methylammonium compounds [e.g. trimethylamine-N-oxide (TMAO) and glycine betaine] (Gleason et al., 2017; Li et al., 2019; Sokolov and Sokolova, 2019; Zhang et al., 2017). Amino acids play multiple roles as building blocks of proteins and substrates for ATP generation, as well as osmolytes (Somero and Bowlus, 1983; Venter et al., 2018; Yancey, 2005). Changes in the amino acid pool during heating may induce negative consequences for protein synthesis and metabolic homeostasis (Haider et al., 2020). Thus, the composition and utilization of the free amino acid pool may play multiple important roles in marine invertebrates in the face of thermal stress.

In the present study, we performed a suite of field and laboratory experiments designed to reveal new aspects of the thermal tolerance mechanisms of E. malaccana. We sought to perform an integrated analysis that shed light on the types of changes in the metabolome that could play roles in protein and membrane stabilization, and in reorganization of ATP-generating pathways during increased temperature. To this end, we determined the field operative temperature and cardiac performance curves of E. malaccana, measured concentrations of metabolites under different temperatures by using untargeted metabolomics, and quantified the content of free amino acids and compatible osmolytes (TMAO, glycine betaine, creatinine, choline and carnitine) with targeted metabolomics. This study provides novel insights into the mechanisms used by this extremely thermophilic mollusc to tolerate some of the highest temperatures animals can withstand.

Specimens of Echinolittorinamalaccana (Philippi 1847) (body length 8±2 mm) were collected in Xiamen, Fujian, China (24°33′N, 118°09′E) in May 2017 and June 2019. All snails were transported back to the laboratory at Xiamen University within 2 h. They were then randomly allocated to plastic containers (16×11×6 cm) at fixed densities (∼30 individuals in each container) and acclimated at 25°C for ∼60 days, a period likely to be sufficient to eliminate effects of past thermal history (Stenseng et al., 2005). Snails were fed with the green alga Chlorella vulgaris twice each week.

Temperature sensors (iButtons, DS1922L, Maxim Integrated, WI, USA) were deployed on the rock surface of the high intertidal zone where E. malaccana lives to monitor habitat temperatures at the sampling site. iButtons were attached to the rock using Z-Spar Splash Zone Compound (Kopper's Company, PA, USA). Temperatures were measured for 3 months (from July to September 2014) at a recording frequency of once every 60 min. Because the iButtons used for this long-term temperature monitoring were lost after 2014, we had to use the temperature data collected in 2014 as test data of habitat temperature. Although this caused a time difference between the collection of environmental (in 2014) and biological (in 2017 or 2019) data, the environmental data collected in 2014 can represent the general pattern of the thermal regime that the experimental E. malaccana faced.

Collections of snails and measurements of heart rates were performed in 2017. After acclimation at 25°C for ∼60 days, snails were placed in air-filled glass containers (N=30, one individual per container) covered with foil containing holes for gas exchange. The snails were then heated from the acclimation temperature (25°C) to the temperature at which the heartbeat stopped [flatline temperature (FLT)], at a rate of 0.1°C min⁻¹ in a water bath (Grant TXF 200, Grant, UK). During heating, heart rate was measured in air using a non-invasive method (Burnett et al., 2013; Dong and Williams, 2011) utilizing infrared radiation. The infrared reflective sensors were adhered to the shells of snails directly over the heart with Blu Tack adhesive (Bostik, Staffordshire, UK), and the infrared signal (variations in the light-dependent current) was amplified (AMP03, Newshift, Leiria, Portugal), filtered, smoothed and recorded with a PowerLab converter (ADInstruments, March-Hugstetten, Germany). Heart rate data were analyzed using LabChart (v.7.2) software. Arrhenius breakpoint temperatures (ABTs), the temperature at which further heating leads to a sharp fall in heart rate, were calculated by a linear regression using Origin version 9.0 (OriginLab, MA, USA). Cardiac performance curves were generated using an exponentially modified Gaussian (EMG) analysis for each individual snail and using an additive mix modeling analysis (GAMM) for all individuals (Dong et al., 2017). Thus, GAMM was used to describe the thermal sensitivities of the heart rate of a population. An additive model permitted a non-linear response to temperature, without requiring us to specify a function (Angilletta et al., 2013; Zuur et al., 2009). EMG was used to describe thermal performance curves, and the complex structure of EMG should provide a better fit to non-linear data (Angilletta, 2006). EMG and GAMM were carried out in R (v. 3.6) with the packages mgcv (Wood, 2004) and nlme (Pinheiro et al., 2014).

The snails used for enzyme activity assays were sampled in 2019 from the same site used to collect specimens for the heart rate study. The activities of two important enzymes in regulating the glycolytic pathway, hexokinase (HK) and pyruvate kinase (PK), were measured. HK initiates pathway activity by providing phosphorylated glucose for catabolic enzymes. The metabolic switch from aerobic to anaerobic pathways of ATP generation is controlled by the reversible phosphorylation of PK and the allosteric regulation that affects its activity (Li et al., 2017; Meng et al., 2018; Moullac et al., 2007). Samples (N=3) used for enzyme analysis were heated according to the protocol described below (‘Metabolomic profiling’ section) for metabolomic analyses but were not cooled back to 25°C. For the enzyme activity assays, 80±20 mg of foot muscle from 10 specimens was pooled together as one sample (Table 1). Tissue of each sample was homogenized in nine volumes of extraction buffer (BC0540, BC0740, Solarbio, Beijing, China). The homogenate was then centrifuged (8000 g, 10 min, 4°C) to obtain the supernatant fraction that was used to detect the activity of HK and PK using commercial assay kits (BC0540, BC0740, Solarbio) according to the manufacturer's protocol. The enzyme activity was determined at 25°C by recording the change in absorbance at 340 nm due to oxidation of NADH, using a spectrophotometer (UV-1800, SHIMADZU, Japan). One unit of enzyme activity corresponded to the oxidation of 1 μmol NADH (340 nm) min⁻¹.

Specimens were collected and metabolite analyses performed in 2019. Both targeted and untargeted metabolomics were applied. Untargeted metabolomics can be used for global metabolome analysis of all the detectable metabolites. Targeted metabolomics, which focuses on a limited number of metabolites, can quantify the metabolite precisely by employing authentic analytical standards for the metabolites of interest. In the untargeted metabolomics analysis, the peak area of each metabolite is used as the parameter in the statistical analysis to define the concentration differences between different biological samples being measured. This type of analysis yields what is termed a relative quantification. In the targeted metabolomics analysis, the internal standards were used to construct calibration curves, which were further used to quantify the content of each metabolite; this approach is referred to as absolute quantification (Griffiths et al., 2010; Naz et al., 2014).

After acclimation at 25°C, 80 samples (foot tissue from eight snails was combined for a single sample, for a total of 640 snails) were used for metabolomic analyses (Table 1). The control samples were kept at 25°C (C25) and the other samples were heated in air from 25°C to the designated treatment temperatures (37°C, H37; 45°C, H45; and 52°C, H52) at a rate of 0.1°C min⁻¹ using glass chambers held in a water bath (Grant TXF 200, Grant, UK). After reaching each designated temperature, samples in the H37, H45 and H52 treatment conditions were cooled to 25°C at a rate of 0.1°C min⁻¹ and then maintained in 25°C seawater for 2 h. Recovery was used to optimize detection of the changes in metabolite content (Malmendal et al., 2006; Teets et al., 2012). After recovery for 2 h, each sample was frozen in liquid nitrogen and then thawed for homogenization. After grinding mixed samples (eight snails were combined for a single sample for analysis). All homogenates were stored at −80°C until metabolite extraction.

To evaluate the global changes in metabolites during heat stress, untargeted metabolomics was performed to compare the differences in content of all detectable metabolites among the different temperature treatments. For the control treatment and each heat treatment, 20 replicates for each treatment were used for untargeted metabolomics analysis. The extraction protocol was as described in Wang et al. (2016). Each sample for the control treatment and each heat treatment was placed into a solution [acetonitrile:methanol:water=2:2:1, 2 μmol l⁻¹ internal standard; 75-05-8 (acetonitrile), 67-56-1(methanol), 103616-89-3 (IS,2-chloro-l-phenylalanine), CNW Technologies, Shanghai, China] for metabolite extraction. After a 30 s vortex, the samples were homogenized at 35 Hz for 4 min and sonicated for 5 min in an ice-water bath. The homogenization and sonication cycles were repeated twice, then the sample was incubated at −40°C for 1 h and centrifuged for 15 min at 9600 g and 4°C. The supernatant was dried in a vacuum concentrator at 37°C. Then, the dried samples were reconstituted in 50% acetonitrile by sonication on ice for 10 min. The preparation was then centrifuged for 15 min at 16,250 g and 4°C, and the supernatant was collected and transferred to a fresh glass vial for further analysis. Ultra-high-performance liquid chromatography (UHPLC) separation was carried out using a UHPLC system (1290 II Infinity series, Agilent, USA), equipped with a UPLC BEH Amide Column (2.1×100 mm, 1.7 μm, Waters). The column temperature was 25°C and the auto-sampler temperature was 4°C. A mass spectrometer (TripleTOF 6600, AB Sciex, USA) was used for its ability to acquire tandem mass spectrometry (MS/MS) spectra on an information-dependent basis during a liquid chromatography (LC)/mass spectrometry (MS) experiment. In this mode, the acquisition software (Analyst TF 1.7, AB Sciex, USA) continuously evaluates the full scan survey MS data as it collects data and triggers the acquisition of MS/MS spectra depending on preselected criteria.

To evaluate the changes in organic osmolytes during heat stress, targeted metabolomics was performed to determine the content of 25 free amino acids (1-methyl-l-histidine, 3-methyl-l-histidine, 4-aminobutyric acid, 4-hydroxyproline, 5-hydroxylysine, beta-alanine, l-asparagine, l-aspartic acid, l-citrulline, l-methionine, l-ornithine, l-phenylalanine, l-proline, l-threonine, l-tryptophan, l-tyrosine, l-valine, glycine, l-alanine, l-arginine, l-glutamic acid, l-glutamine, l-histidine, l-lysine and l-serine), TMAO, glycine betaine, creatinine, choline and carnitine. For the control group and each heat treatment group, foot tissue from eight specimens was combined for a single sample for analysis, leading to a total of 10 replicates for each treatment for the analysis of free amino acids, and four replicates for each treatment for the analysis of TMAO, glycine betaine, creatinine, choline and carnitine (Table 1). Sample for analysis of free amino acids was placed into the solution (acetonitrile:methanol:water=2:2:1, 1 μmol l⁻¹ internal standard), and sample for analysis of TMAO, glycine betaine, creatinine, choline and carnitine was placed into the solution [acetonitrile:water=4:1, 0.1% formic acid (64-18-6, CNW Technologies), 200 nmol l⁻¹ internal standard], for metabolite extraction and further analysis. The contents were measured by UHPLC-MS/MS analysis (1290 Infinity II series UHPLC System and 6460 Triple Quadrupole Mass Spectrometer, Agilent) and calculated with reference to calibration curves. Working standard solutions were individually prepared by dissolving or diluting each standard substance to give a final concentration of 10 μmol l⁻¹. A series of calibration standard solutions was then prepared by stepwise dilution of this standard solution. The isotopic compounds corresponding to each target compound are used as the internal standard (IS) at a concentration of 200 nmol l⁻¹. Calibration solutions were subjected to UHPLC-multiple reaction monitoring (MRM)-MS/MS analysis. In calibration curves, y is the ratio of peak areas for analyte/IS, and x is the concentration (nmol l⁻¹) of the analyte. Least squares method was used for the regression fitting, and 1/× weighting was applied in the curve fitting as it provided highest accuracy and correlation coefficient (R²).

The differences in enzyme activity and the concentrations of free amino acids, TMAO, glycine betaine, sarcosine, glycerol or inositol among the C25 (the control treatment), H37, H45 and H52 treatment groups were analyzed with one-way ANOVA using SPSS Statistics 22 (IBM, Armonk, NY, USA).

To analyze untargeted metabolomic data, MS raw data (.wiff) files were converted to the mzXML format by ProteoWizard and were processed by R package xcms (version 3.6). The process includes peak deconvolution, alignment and integration. The minFraction (parament of peak alignment) and cutoff were set as 0.5 and 0.3, respectively, and resolution was set as 20,000. The metabolites were identified according to HMDB_ALL_IN_BIOTREEMSMSBANK (BioTree Biotechnology, Shanghai, China). The mass error tolerances for MS1 and MS2 matches were set to ±25 ppm and ±35 ppm, respectively. Compounds with a MS2 spectral score of >0.6 were further analyzed (Liu et al., 2020 preprint). The peak width (the value corresponding to the sample number) and retention time (RT) reproducibility (rtmed in metadata) were shown in metadata, and the metadata and MS spectra have been submitted to the China National Microbiology Data Center (accession number NMDCX0000105). MS2 spectral similarity scores were set to range from 0 to 1, and metabolites with a score higher than 0.6 were further analyzed (Liu et al., 2020 preprint).

The MetaboAnalyst 4.0 website (https://www.metaboanalyst.ca/) was used to perform statistical analysis (normalization, multivariate, univariate and cluster analysis) and pathway analysis (Figs S1 and S2). Normalization by reference sample (Dieterle et al., 2006), generalized logarithm transformation and auto scaling were carried out to normalize raw data. In pathway analysis, the negative log-transformed P-values (–log P) were used to explain the degree of difference in metabolic pathways between different temperature treatments.

Considering that some metabolites containing easily protonated functional groups (such as acylcarnitines with a trimethylamine group) can only be detected in the positive ionization mode, whereas some compounds with easily deprotonated functional groups (such as fatty acids) are mainly detected in the negative ionization mode, duplicate analysis was conducted in positive and negative modes of all samples – one with each ionization mode – to increase metabolite coverage (Arbulu et al., 2015; Calderón-Santiago et al., 2016). LC-MS/MS analysis was conducted in both positive and negative modes to detect metabolites that have different polar groups. Compounds identified in positive and negative modes were combined to obtain the overall set. For compounds identified in both positive and negative modes, those with higher MS2 scores were retained if they displayed consistent changing trends, whereas those with contradictory trends were removed. One-way ANOVA and partial least squares-discriminate analysis (PLS-DA) were performed to determine metabolites for which content significantly differed among the four treatments [adjusted P<0.05, variable importance for the projection (VIP) score >1].

To understand the relationship between thermotolerance and certain metabolites and metabolic pathways, we compared our results with literature data on critical temperature [the temperature that is lethal for one-half of the study population (LT₅₀) and FLT], induction temperature (the temperature that induces increases in metabolites) and content of amino acids of common molluscs, including Littorina littorea (Melatunan et al., 2011; Newell et al., 1971), Littorina saxatilis (Sokolova et al., 2000a,b), Theodoxus fluviatilis (Wiesenthal et al., 2019), Haliotis fulgens (Tripp-Valdez et al., 2019), Crassostrea gigas (Ghaffari et al., 2019; Haider et al., 2020; Moreira et al., 2017; Moullac et al., 2007; Pazos et al., 1996; Yang et al., 2016), Ostrea edulis (Eymann et al., 2020), Mytilus galloprovincialis (Digilio et al., 2016; Han and Dong, 2020; Yao and Somero, 2012), Mytilus californianus (Yao and Somero, 2012), Perna canaliculus (Dunphy et al., 2018) and Laternula elliptica (Clark et al., 2017; Peck et al., 2002).

Habitat temperatures of E. malaccana were monitored at a rocky intertidal site typical for the species during the hottest season from July to September 2014 in Xiamen. During the period of measurement, the highest temperature was 53.03°C, and the mean habitat temperature was 32.07±5.73°C (Fig. 1).

Considerable variation in heart rate (beats min⁻¹), but not in the pattern of the response of heart rate to rising temperature, was found among individuals (Fig. 1). The cardiac thermal performance curve generated by incorporating data from all individuals was bimodal in response to rising temperature (Fig. 1, solid red line). Heart rate increased from 28.85±7.95 beats min⁻¹ (mean±s.d.) at 25°C to 42.00±13.10 beats min⁻¹ at ∼37°C, reaching the first peak. Heart rate then slowly decreased from 37°C to 45°C (35.20±12.68 beats min⁻¹). As the temperature increased above 45°C, heart rate again increased with rising temperature, attaining a maximum rate (67.35±15.87 beats min⁻¹) at ∼52°C. The calculated ABT and FLT were 51.24±1.92°C and 54.73±1.55°C, respectively.

Considering cardiac performance in the context of the relative frequencies of field temperatures the snail would likely experience, the temperature range over which a rise in temperature led to the first peak in heart rate, 25°C to 37°C, comprised 79.76% of the field temperatures. The relative frequency of temperatures between 37°C and 45°C, when heart rate slowly decreased, was 14.53%. The relative frequency of temperatures between 45°C and 52°C was 5.57%, and a temperature above 52°C, which was higher than the ABT of E. malaccana, occurred only once.

The mean HK activities of muscle tissue of E. malaccana were significantly different among the treatment groups (one-way ANOVA, F_3,8=60.71, P<0.001) (Fig. 2A, Table S1). The HK activity in H45 (11.33±0.30 U, mean±s.d.) was higher than that in C25 (P<0.001), H37 (P<0.001) and H52 (P<0.001), and the activity in H52 (7.50±0.28 U) was higher than that in C25 (P=0.025) and H37 (P=0.007). There was no significant difference in HK activity between C25 (6.26±0.67 U) and H37 (5.88±0.44 U) (P>0.05).

The mean PK enzyme activities of muscle tissue in E. malaccana were significantly different among the four heat-stress treatments (one-way ANOVA, F_3,8=6.394, P=0.0161) (Fig. 2B, Table S1). The PK activity was significantly lower in H37 (6.00±1.76 U, mean±s.d.) than that in C25 (10.35±0.70 U) (P=0.039), H45 (12.59±2.11 U) (P=0.006) and H52 (12.76±2.08 U) (P=0.005). There was no significant difference in PK activity among C25, H45 and H52 (P>0.05).

A hierarchical clustering analysis combining positive-mode and negative-mode measurements was performed to reveal the similarities and differences of the samples according to the relative values of metabolites under different experimental conditions. The results of cluster analysis showed that the top 25 metabolites that were significantly different among different treatments (P<0.05) could be grouped into three clusters (Fig. 3A). The metabolites in Cluster 1 accumulated with increasing temperature, and the levels of metabolites in Cluster 2 decreased with increasing temperature. The levels of metabolites in Cluster 3 first increased and then decreased as the temperature increased. The top 25 metabolites that were significantly different among the four treatments (P<0.05) were related to glycerophospholipid metabolism, the tricarboxylic acid (TCA) cycle, amino acid metabolism, sugar metabolism and purine metabolism, and these metabolites were interconnected: most metabolites in Cluster 2 or Cluster 3 can generate metabolites in Cluster 1 (Fig. 3B).

PLS-DA was used to establish a relationship model between the relative value of metabolites and experimental conditions, in order to allow the distinguishing of differences among the different treatments and to predict the classification of samples. PLS-DA projection results showed that samples in this study could be separated into four groups in positive mode (R²=0.86, Q²=0.57) and negative mode (R²=0.89, Q²=0.63) (Fig. S3). The biomarkers were identified by one-way ANOVA and PLS-DA (adjusted P<0.05, VIP score >1). According to the functions of these metabolites, these biomarkers were divided into four groups (Fig. 4).

The metabolites in Group 1 were glucose 1-phosphate (Fig. 4A) and glucose 6-phosphate (Fig. 4B). The concentration of glucose 1-phosphate in H37 was significantly lower than that in C25 (P=0.011), H45 (P=0.021) and H52 (P<0.001), and its concentration in H52 was significantly higher than that in C25 (P=0.048), H37 (P<0.001) and H45 (P=0.022). The concentration of glucose 6-phosphate in H52 was significantly higher than that in C25 (P<0.001), H37 (P<0.001) and H45 (P=0.008).

Group 2 metabolites comprised lactate (Fig. 4C) and succinate (Fig. 4D), which are related to anaerobic metabolism. The content of lactate significantly increased at 37°C (P=0.002), 45°C (P=0.049) and 52°C (P<0.001), and the content in H52 was significantly higher than that in H45 (P=0.039). The content of succinate also significantly increased at 37°C (P=0.01) and 52°C (P<0.001), and the content in H52 was significantly higher than that in H45 (P=0.019).

Group 3 contained sn-glycerol-3-phosphoethanolamine (GPE; Fig. 4E) and CDP-ethanolamine (CDP-EA; Fig. 4F), which are components of cell membrane lipids. The content of GPE increased with rising temperature and was significantly increased at 45°C (P=0.044) and 52°C (P<0.001); the content in H52 was significantly higher than that in H45 (P=0.002). The content of CDP-EA in H52 was significantly higher than that in H37 (P=0.022).

Group 4 contained S-methyl-5′-thioadenosine (MTA), GDP-l-fucose, glutathione disulfide (GSSG), and the dipeptides Tyr–Met, Tyr–Ala and Tyr–Asp (Fig. 4G–L), which are related to the cellular stress response. The content of MTA was significantly higher in H52 compared with that in the other three treatments (P<0.001). The content of GDP-l-fucose significantly increased at 45°C (P=0.004) and 52°C (P<0.001), and the content in H52 was significantly higher than that in H45 (P=0.004). In addition, the content in H37 was lower than that in C25, but not significantly (P>0.05). The content of GSSG was significantly higher in H52 compared with that in C25 (P=0.021), H37 (P=0.036) and H45 (P=0.045). The content of Tyr–Met decreased with increasing temperature and significantly decreased at 45°C (P=0.02) and 52°C (P<0.001); there was no significant difference between the content in H37 and H45 (P>0.05). The trend for Tyr–Ala content was the same as that for Tyr–Met: the content significantly decreased at 45°C and 52°C (P<0.001). The Tyr–Asp content also decreased with increasing temperature and was significantly lower in H45 (P=0.002) and H52 (P<0.001) compared with that in C25. Fig. 4M–O show the relationships among LT₅₀, FLT and the induction temperatures (temperature at which significant increase in the concentrations of metabolites in Groups 1, 2 and 3 occurs) for E. malaccana and other marine invertebrates. The temperatures for E. malaccana were all higher than those in the other molluscs surveyed.

These differentially abundant metabolites (i.e. biomarkers) among treatments were used to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway analysis. The results showed that, compared with 25°C, there were 11, 4 and 25 significantly enriched pathways (P<0.05) at 37°C, 45°C and 52°C, respectively (Fig. 5A–C). Four pathways, including the citrate (TCA) cycle, tryptophan metabolism, pyruvate metabolism and lysine degradation, were significantly influenced throughout the heating process (Fig. 5C). There were 12 pathways enriched by differentially abundant metabolites between H37 and H45, 22 pathways between H37 and H52, and 21 pathways between H45 and H52 (P<0.05) (Fig. 5D–F). The pathway related to alanine, aspartate and glutamate metabolism was sensitive to temperature increases and was the pathway most significantly affected in H37 (–log P=8.27) and H52 (–log P=6.04) compared with C25, but was not significantly affected in H45 (–log P=2.48) (Fig. 5A–C). The pathway related to cysteine and methionine metabolism was the most significantly affected pathway in H52 compared with C25 (–log P=19.374, Fig. 5C) and H45 (–log P=24.067, Fig. 5F). The pathway related to amino sugar and nucleotide sugar metabolism was significantly affected in H52 compared with H37 (–log P=18.944). Most metabolic pathways of E. malaccana remained at steady state under moderate thermal stress (37–45°C) (P>0.05), but they were significantly affected at extreme high temperature (52°C) (P<0.05, –log P>2.996).

A total of 25 amino acids and derivatives were identified in the muscle tissue of E. malaccana, and the concentrations of 14 amino acids were above 0.1 mmol kg wet weight⁻¹. The levels of glutamate and serine generally increased with increasing temperature (Fig. 6A,B). The content of glutamate in C25 was significantly lower than that in H45 (P=0.038) and H52 (P<0.001), but there was no significant difference in glutamate content among the H37, H45 and H52 treatment groups (P>0.05). The concentration of serine in H52 was significantly higher than that in C25 (P<0.001), H37 (P<0.001) and H45 (P=0.005), but there was no significant difference in serine concentration among the C25, H37 and H45 groups (P>0.05). The content of most free amino acids was maintained at stable levels under heat stress (Fig. 6C–I).

TMAO and creatinine were not detected in the muscle tissue of E. malaccana, but other members of the methylamine family – including choline, glycine betaine and carnitine – were identified and showed significant differences among the four heat-stress treatments (P<0.05) (Fig. 6J–L). The concentration of betaine significantly increased with increasing temperature (P<0.05) and reached a maximum in H52. The concentration of choline also increased with increasing temperature, and the concentration in C25 was significantly lower than that in H37 (P=0.04), H45 (P=0.008) and H52 (P=0.008). The concentration of carnitine was significantly higher at 52°C compared with that at 25°C (P=0.005). Compared with other metabolites, betaine was present at a high concentration (>40 mmol kg wet weight⁻¹) in C25, and its concentration was above 90 mmol kg wet weight⁻¹ under heat stress.

The thermophilic snail E. malaccana employs a variety of mechanisms to adapt to the extremely high temperatures of its habitat, which may approach 60°C. We found that a number of physiological responses, including initiation of anaerobic metabolism and metabolic depression, accumulation of metabolites protecting cells from damage by free radicals, increases in glycerophospholipid metabolism for cell membrane stabilization, and elevation of compatible osmolytes and amino acids for protein structure stability, were induced at thermal stress of different intensities. Together, these physiological adjustments combine to help E. malaccana survive at the thermal extremes of its habitat (Fig. 7).

Maintaining an adequate supply of energy to support the cellular stress response as well as normal metabolic requirements is essential for surviving at high temperature. Moreover, high temperature and hypoxia commonly co-occur in intertidal molluscs (Deutsch et al., 2020), and, under these combined stresses, animals either have to maintain a low aerobic metabolic rate or provide energy via anaerobic metabolism (Han et al., 2017; Hui et al., 2020; Marshall et al., 2011; Meng et al., 2018). According to the cardiac performance curve of E. malaccana, metabolic depression occurs from ∼37°C to ∼45°C, yet snails can survive for a period when the ambient temperature exceeds the cardiac ABT (∼52°C). These results indicate that metabolic depression enables the snail to maintain a relatively low metabolic rate for surviving intermediate high temperatures. This capacity is an important physiological adaptation against thermal stress for the species inhabiting the high intertidal zone (Marshall et al., 2011; Hui et al., 2020). These results also imply that anaerobic metabolism should be an important alternative or supplementary source of energy supply for coping with the elevated energy requirements in the face of thermal stress.

The metabolomics results suggest that anaerobic metabolism at intermediate thermal stress (from 37°C to 45°C) plays an important role in energy supply. The content of the anaerobic end products lactate and succinate increased at 37°C and 52°C. The accumulation of products of anaerobic metabolism under hypoxic stress, which can be induced by high temperature, high partial pressure of CO₂ (pCO₂) or anoxia, is common in intertidal organisms. For example, in the periwinkle L.littorea, the levels of lactate and succinate increase when snails are exposed to elevated temperature or pCO₂ (Kluytmans et al., 1977, 1978; Melatunan et al., 2011; Zurburg and Kluytmans, 1980). During anoxia, the end products of anaerobic metabolism, lactate and tauropine, accumulate in the abalone Haliotis midae (O'Omolo et al., 2003). In summary, an increase in anaerobic metabolism can be a physiologically adaptive response to short-term hypoxia or anoxia, ensuring adequate energy supply during periods of overall metabolic depression.

The changes in the activities of HK and PK in the glycolysis pathway are important for energy production and thus for organisms to survive under stress (Faiz-ur-Rahman et al., 1974; Kubota and Ashihara, 1990; Mutuku and Nose, 2012). Glucose 1-phosphate and glucose 6-phosphate, intermediate metabolites of glycolysis, were also reported to accumulate in rice for surviving under stress (Mutuku and Nose, 2012). In the present study, the activity of HK increased at 45°C, and glucose 6-phosphate, as a product of HK, accumulated at 52°C. The PK activity significantly decreased at 37°C, possibly as a result of phosphorylation, which reduces the activity of PK (Moullac et al., 2007). The depression of PK activity could be related to increased anaerobic metabolism under moderate stress. The finding that the content of metabolites (glucose 1-phosphate and glucose 6-phosphate) and activities of HK and PK in E. malaccana decreased at moderate thermal stress and then increased at extremely high temperature could indicate that glycolysis is inhibited under mild stress (37–45°C) and then activated under extreme thermal stress (45–52°C). These results further confirm that metabolic depression in response to moderate thermal stress is an important physiological adaptation to thermal stress on the rocky shore.

Alanine and asparagine play important roles in protecting organisms from hypoxic damage. High levels of alanine could help organisms to cope with the challenge of anoxia (Gleason et al., 2017). Echinolittorinamalaccana contains a high level of alanine, which might be associated with its anaerobic ability. Asparagine is derived from aspartate, and aspartate is made by transamination of oxaloacetate, which is an intermediate product of the TCA cycle. KEGG metabolic pathway analysis results showed that the pathway related to alanine, aspartate and glutamate metabolism in E. malaccana was the most sensitive pathway to high temperature. Compared with other species, the higher levels of alanine and asparagine in E. malaccana under non-stressful conditions could be important for enhancing the anaerobic capacity (Haider et al., 2020; Wiesenthal et al., 2019). It is worth noting that aspartate decreased with the accumulation of alanine under hypoxic stress, which would maintain the total amino acid pool at a constant level (Baginski and Pierce, 1978; Zurburg and Kluytmans, 1980). It is plausible that the low content of aspartate in E. malaccana is for maintaining the amino acid pool balance when the alanine content is high under non-stressful conditions.

The cellular stress response is activated in response to stress-induced macromolecular damage, as may occur under extreme temperature and concentrations of reactive oxygen species (ROS) (Kültz, 2020). Previous studies showed that peptides may have various physiological modulatory and regulatory functions, and some dipeptides can protect cells from damage by free radicals. Tyr-containing dipeptides have the highest radical-scavenging activities among dipeptides (Kati et al., 2006; Zhang et al., 2016; Zheng et al., 2016). In E. malaccana, Tyr–Met, Tyr–Ala and Tyr–Asp, all of which are antioxidative dipeptides, and glutathione disulfide (GSSG) were identified. The content of GSSG increased only at 52°C, but dipeptides decreased at 37°C. The decrease in antioxidants and the accumulation of oxidative stress markers indicated that antioxidative dipeptides can scavenge free radicals and ROS, to maintain cellular homeostasis and protect organisms from oxidative damage under mild thermal stress, although dipeptides seem less likely to prevent oxidative damage at the higher temperatures studied, 45°C and 52°C.

GDP-l-fucose is a sugar nucleotide and the activated nucleotide sugar form of l-fucose. l-Fucose is a constituent of many structural polysaccharides and glycoproteins in various organisms and plays important metabolic roles in complex carbohydrates and glycoproteins. It can also induce antioxidant activity against oxidative stress (Kandasamy et al., 2014; Vermeer et al., 2001). The accumulation of GDP-l-fucose at 45°C and 52°C indicated that the antioxidant protection of E. malaccana may change as thermal stress intensifies; reliance on dipeptides is reduced and other systems like those induced by GDP-l-fucose assume greater importance.

MTA is a metabolite of methionine. Previous reports have confirmed that MTA can inhibit protein methylation, either through direct interaction with protein methyltransferases or indirectly through the inactivation of S-adenosylhomocysteine hydrolase (Avila et al., 2004; Williams-Ashman et al., 1982). As an intermediary metabolite, MTA plays regulatory functions in the cellular response to stress and influences cellular proliferation, differentiation and apoptosis (Douglas and Shang, 1977). Coincidentally, the content of methionine, and the pathway related to methionine, which is closely associated with MTA in E. malaccana, was enhanced significantly under extreme stress. The high level of MTA at 52°C implies that MTA may regulate cell function and protect organisms from methylation damage under extreme thermal stress.

The enhancement of glycerophospholipid metabolism under extreme heat stress in E. malaccana could be important in stabilizing cellular membranes. Glycerophospholipids (mainly phosphatidylethanolamine and phosphatidylcholine) are important components of cell membranes (Nguyen et al., 2020; Overgaard et al., 2008), and have close relations with the permeability of the cell membrane and the activity of ATPases (Nemova et al., 2017). The regulation of glycerophospholipid content under thermal stress may be an adaptive response to warmer temperature (Bolstridge et al., 2010; Hunsberger et al., 2014). In the present study, the levels of GPE and CDP-EA generated by diethanolamine in glycerophospholipid metabolism increased. When E. malaccana is exposed to extremely high temperature, elevated levels of GPE and CDP-EA could lead to stabilization of the cell membrane. The finding that glycerophospholipid metabolism in E. malaccana was only affected at 52°C, not at lower temperatures, suggests that cell membranes of E. malaccana may remain stable under moderate thermal stress, but can be damaged at extreme temperature.

Serine can function as a modifying group in the synthesis of phosphatidylserine (PtdSer). As an important membrane phospholipid, PtdSer can regulate the function and state of proteins in the membrane and plays an important role in cell metabolism (Pérez-Lara et al., 2016). In the present study, serine accumulated significantly with increasing temperature, and its accumulation may help in regulating cellular metabolism and maintaining the structural stability of cell membranes (Vance, 2008; Vance and Tasseva, 2013).

The concentrations of the compatible osmolytes glycine betaine, choline and carnitine increased with increasing temperature, which could offset the effects of high temperature on protein stability. Methylamines such as TMAO, betaine and carnitine can enhance protein stability by stabilizing conformation and subunit–subunit interactions, leading to favorable effects on enzyme activity and thereby helping to ensure the provision of sufficient cellular energy metabolism under stress (Peluso et al., 2001; Yancey, 2005). The content of protective osmolytes such as betaine and taurine shows adaptive change under heat stress in marine invertebrates; for example, the levels of betaine and taurine increase in the intertidal mussel M.californianus at high temperature (Gleason et al., 2017). In the present study, betaine was present at high concentrations (40–100 mmol kg wet weight⁻¹) compared with those in the other intertidal molluscs surveyed (Gleason et al., 2017; Zhang et al., 2017), and the content of betaine, choline and carnitine significantly increased with increasing temperature, which could potentially protect the protein structure from denaturation in high temperatures.

Previous studies have found that elevated concentrations of free amino acids under stress might play important roles in helping organisms survive in harsh environments (Haider et al., 2020; Wiesenthal et al., 2019; Zemanová et al., 2017). In the present study, in E. malaccana, under non-stressful conditions, alanine, asparagine, methionine, glutamate, arginine, glutamine and serine content was higher than that in other molluscs (Fig. S4, Table S2). However, the levels of most free amino acids did not significantly change under thermal stress, indicating that they might be accumulated by this species in a stress-preparatory or preventive protection mechanism. The stress-preparatory and protective effects of free amino acids might involve several different pathways. For example, alanine and asparagine are involved in anaerobic metabolism (Gleason et al., 2017; Haider et al., 2020; Wiesenthal et al., 2019), and serine is important for stabilizing membrane structure (Vance, 2008; Vance and Tasseva, 2013). Protein degradation increases and may induce the accumulation of amino acids under stress (Hildebrandt, 2018). However, E. malaccana amino acid levels maintain a steady state under heat stress, which suggests that protein degradation is minimal. Stabilization of the proteome of E. malaccana may derive from the high intrinsic thermal stability of its proteins (Dong et al., 2018; Liao et al., 2019), as well as stabilization from thermoprotective osmolytes.

The periwinkle snail E. malaccana experiences periods of intensive thermal stress during the summer. However, even during the hottest season (from July to September 2014) in Xiamen, E. malaccana generally encounters what are likely to be relatively non-stressful temperatures, as indexed by the physiological and biochemical studies reported here. Temperatures between 25°C to 37°C occur with the highest frequency (79.76%). From 37°C to 45°C (the relative frequency of temperature is 14.53%), E. malaccana suffered what we define as moderate thermal stress: over this temperature range, snails showed metabolic depression for energy saving as found in tropical shore molluscs (Marshall et al., 2011; Marshall and McQuaid, 2011). When temperature increased from 45°C to 52°C (5.57% of field temperatures), snails suffered what we define as extreme thermal stress: heart rate increased and a number of metabolites related to coping with extreme thermal stress were produced. Overall, most of the field temperatures fall into a range in which thermal stress as defined here will not occur, and, when temperatures rise into the range in which moderate or extreme stress occurs, E. malaccana appears to have a battery of adaptive responses available to cope with thermal challenges.

We combined analyses of field temperature data, cardiac responses to heating and metabolomics to elucidate how the thermophilic snail E. malaccana survives in an extremely stressful thermal environment. Results showed that temperatures over a range at which cardiac activity alters its thermal response are associated with changes in the metabolome that are consistent with a shift to anaerobic metabolism and the production of metabolites that are protective against oxidative stress, methylation reactions, and protein and membrane destabilization.

The authors thank Prof. G. N. Somero for comments, discussion and revisions to the manuscript.