Thermal history and gape of individual Mytilus californianus correlate with oxidative damage and thermoprotective osmolytes

Across broad spatial and temporal scales, adaptive variation in organismal defenses against (and susceptibility to) environmental stress generally correlates with mean environmental conditions (Vernberg, 1962; Somero, 2005; Somero et al., 2017). At the other extreme, micro-scale environmental variation, for instance among individuals within a population, potentially contributes to larger-scale patterns of physiology, ecology and evolution (Bolnick et al., 2011; Pruitt and Ferrari, 2011; Pelletier and Garant, 2012; Farine et al., 2015; Dowd et al., 2015; Jimenez et al., 2015; Lathlean et al., 2016). Each individual performs based on its unique environmental experience, yet knowledge of the effects of micro-scale spatial and temporal environmental variation on individual animals' physiology remains sparse (but see Helmuth and Hofmann, 2001; McGaughran et al., 2010). Such knowledge may be pivotal to forecasting the biological effects of global change, as it will help clarify the roles of adult physiological acclimatization, developmental physiological plasticity and physiological adaptation in coping with environmental changes.

Temperature variation, in particular, affects nearly all aspects of ectotherm physiology, from macromolecular structure and function to reproductive performance and survival. Within species, correlations between mean habitat temperature and thermal stress tolerance have been documented between populations, but few studies have overcome the logistical constraints that hamper attempts to study individual thermal physiology in the field (see, for example, Miller et al., 2015).

Temperature variation plays an especially prominent role in the rocky intertidal zone, the dynamic interface between marine and terrestrial environments. In these habitats, temperature varies across latitudinal gradients, vertical tidal heights and also within sites across micro-scales on the order of centimeters (Bingham et al., 2011; Denny et al., 2011). In fact, variation in body temperature among individual Mytilus californianus mussels within a single bed can exceed large-scale latitudinal differences (Helmuth et al., 2006; Denny et al., 2011). This persistent micro-scale variation contributes to among-individual physiological variation (e.g. in antioxidant capacities), particularly when the environment varies around a stressful mean (Jimenez et al., 2015). This effect presumably acts via physiological plasticity as each individual attempts to ‘match’ its physiology to its unique micro-environment. However, the relationship between an individual's current physiological status and its past environmental experience remains elusive for most organisms.

We take advantage of novel individual monitoring technology, the sessile nature of many rocky intertidal zone invertebrates and the high degree of spatial thermal variation in this habitat to assess the physiological responses of individual sea mussels (Mytilus californianus Conrad 1837) to micro-scale environmental variation. We combine monitoring of body temperature and valve gaping behavior in the field with several biochemical measures of thermal defenses and macromolecular damage. Following collection from wave-exposed (cool) and wave-protected (warm) origin habitats, individuals were outplanted to three different intertidal sites (tidepool, low and high), and body temperature and valve gaping behavior were continuously recorded for 23 days. Oxidative damage, antioxidant capacities and tissue contents of organic osmolytes – some of which are known thermoprotectants and/or antioxidants – were measured to characterize individual physiological profiles. We hypothesized (1) that higher body temperatures would correspond with increased oxidative damage, and (2) that this increased damage would be counteracted physiologically via increased antioxidant capacities and modulation of organic osmolytes. Overall, our findings do implicate recent thermal history as a predominant driver of an individual mussel's current physiological status.

For a complete description of the MusselTracker system, see Miller and Dowd (2017). Briefly, the system consisted of a custom-built circuit board, real-time clock, micro-SD memory card and ports to attach sensor packages for each individual mussel. These sensors included a calibrated thermocouple that was inserted through a small hole in the shell of each individual mussel to record its mantle cavity temperature (i.e. body temperature, resolution=0.25°C) and a Hall effect magnetic sensor mounted near the apex of the shell that was used to monitor valve gaping behavior (Wilson et al., 2005; Dowd and Somero, 2013). The selected placement of the thermocouple ensured that no tissue other than mantle or gonad was compromised. Data were collected from each individual at 1 Hz. We do not believe the MusselTracker system caused additional stress to the monitored individuals, because (1) the attached sensors are similar to the abundant epifauna covering M. californianus shells, such as barnacles and limpets (Paine, 1976), (2) gaping behavior and byssal thread formation of instrumented individuals followed expected patterns and (3) levels of oxidative DNA damage and catalase activity were comparable to those of uninstrumented mussels collected directly from the field at a similar tidal height (Dowd et al., 2013; L.U.G. and W.W.D., unpublished data).

Adult mussels (n=30, mean valve length 66.8±3.3 mm, range 60.8–72.8 mm) were collected from a wave-splashed ‘exposed’ and a wave-sheltered ‘protected’ site, situated 24 m apart, within a single population at Hopkins Marine Station (HMS), Pacific Grove, California (36.6203°N, 121.9042°W). Although seawater temperatures at the two sites are the same (2015 range 9.4°C–21.0°C; HMS Marine Life Observatory, http://mlo.stanford.edu/sst.htm), individuals at the exposed site experience cooler body temperatures due to more frequent wave-splash during low tide periods compared with the protected site (Denny et al., 2011). The abiotic conditions and physiological profiles of mussels at these two sites have been extensively characterized (Denny et al., 2011; Dowd et al., 2013; Jimenez et al., 2015). Mussels from the two origin sites were allocated among three outplant sites (Table 1). Individuals were held under common garden conditions in the laboratory for 7 days after attachment of the MusselTracker sensors. During this time they formed byssal thread attachments to the acrylic plates on which the systems were mounted. Non-instrumented live mussels (40–70 mm in size) were packed between the focal animals to mimic natural, dense mussel beds. MusselTracker plates were bolted to the rock substrate at each of three locations (high intertidal, low intertidal and a continuously submerged tidepool), which have distinct thermal and tidal profiles (Table 1). Upon deployment, each plate was covered with a loose-fitting plastic mesh (5 mm openings) for 2 days to prevent immediate dislodgement of animals by wave action.

Following 23 days in the field (15 July–6 August 2015), individuals from all three outplant sites were collected during low tide. Immediately prior to collection, gape recordings indicated that individuals at all sites were closed and, therefore, not feeding. Gill tissue was dissected from each individual and immediately frozen in liquid nitrogen. Samples were stored at −80°C until analysis. The mussel gill constitutes a significant fraction of total tissue mass (allowing subsampling for multiple assays) and has been the subject of numerous physiological studies (e.g. Lockwood et al., 2010; Tomanek and Zuzow, 2010; Jimenez et al., 2015; Jimenez et al., 2016). It is the principal site of gas exchange (and, thus, is a highly aerobic tissue that directly experiences ambient fluctuations in oxygen saturation), and it serves as the suspension-feeding apparatus. Therefore, we expect biochemical measures in the gill to reflect organismal consequences of environmental variation.

The data from each MusselTracker plate were downloaded from the micro-SD cards and concatenated into continuous time series. Those series were manually edited to remove periods when the plates were being maintained, when a sensor had failed, or following instances of mortality. Following this quality control, the MusselTracker data set included nearly continuous records of both body temperature and valve gape for each of 21 individuals. Note that we retained data for additional mussels for which only one of these measures was available for later analyses.

All body temperature data were corrected based on individual thermocouple calibrations performed under a range of temperatures (5–45°C) in a laboratory water bath. Corrected temperatures were rounded to the nearest 0.25°C, the nominal resolution of the thermocouples. The mean daily maximum temperature, overall maximum temperature, cumulative degree-hours above 25°C (Denny et al., 2011) and additional metrics (see below) were calculated for each mussel.

Valve gaping behavior determines the availability of oxygen (Jansen et al., 2009; Nicastro et al., 2010), which is necessary for sustained aerobic metabolism. However, high oxygen levels and/or cycles of hypoxia and reoxygenation can also lead to increased amounts of oxidative damage sustained by an individual (Li and Jackson, 2002). Hall effect gape data were processed to create a scaled output for each individual. The data were filtered with a first order Butterworth filter (10-s window) to smooth the signal. After smoothing, the lower 1st percentile and the upper 99th percentile (the maximum observed change from baseline for that individual) values were identified and used to scale the data into relative opening values between 0 and 100%. For full details, see Miller and Dowd (2017). For cases in which magnets or Hall effect sensors were dislodged during deployment in the field, a new zero point was determined immediately after reattachment. For each individual we calculated the proportion of time spent open during the field monitoring period. A gaping threshold of 20% of the maximum gape was established from empirical cumulative distribution functions of each individual’s gape data (Miller and Dowd, 2017). Specifically, these functions indicated an inflection point at roughly 20% of maximum gape, with mussels spending a high proportion of time both above and below 20% open but very little time around this value. In addition, M. californianus of the same size as those used in this study generally started pumping water between 15% and 20% opening in a laboratory tank, so 20% represents an opening value at which we may generally assume that valves are open enough to sustain water flow (L.P.M., personal observations). Lastly, previous studies have shown that a gaping threshold of 20% represents a 95% probability of valves being closed in bivalve mollusks (Jou et al., 2013; Ballesta-Artero et al., 2017). Thus, we considered a mussel to be ‘open’ and to have access to food and oxygen when the gape was greater than 20% of the maximum. Gaping behavior was only observed in mussels during high tide; they did not gape at low tide when exposed to air (Miller and Dowd, 2017). Interestingly, even constantly submerged tidepool individuals closed up as a group during low tide.

The oxidative DNA adduct 8-hydroxy-2′-deoxyguanosine (8-OHdG), one of the most prevalent types of oxidative damage to DNA, can serve as an oxidative stress biomarker in a wide range of organisms including mammals and aquatic invertebrates (Liu et al., 2004; Li et al., 2005; Valavanidis et al., 2006; Halliwell and Gutteridge, 2007; Lister et al., 2015). Heat stress increases the amount of 8-OHdG (Huang et al., 2012, 2015), and this oxidative DNA damage can result in mutagenesis (Kasai, 1997; Halliwell, 2000). We used a competitive, enzyme-linked immunosorbent assay (ELISA) for quantitative detection of 8-OHdG (Japan Institute for the Control of Aging, Nikken Seil Co., Ltd, Fukuroi City, Japan). DNA was extracted from approximately 25 mg of gill tissue from each individual with a NucleoSpin Tissue Kit (Macherey-Nagel, Düren, Germany), digested with nuclease P1 to fragment the DNA into single nucleotides, and treated with alkaline phosphatase to convert single nucleotides to single nucleosides. After this sample pretreatment, the 8-OHdG ELISA kit was used according to the manufacturer’s instructions, loading 16.5 µg of digested DNA for each individual. Each individual was run in duplicate.

Lipid peroxidation occurs along cellular membranes via a chain reaction of peroxyl radicals under oxidative stress. Lipid peroxidation is expected to increase following exposure to high temperatures (Ando et al., 1997; Deschaseaux et al., 2010). It can compromise membrane structure and function, alter ion transport and inhibit metabolic processes (Halliwell and Gutteridge, 1990; Girotti, 1998; Catalá, 2009). To quantify this oxidative damage to lipids, we performed a microplate-based version of the ferrous oxidation of xylenol orange (FOX) assay (Wolff, 1994; Hermeslima and Storey, 1995; Gay and Gebicki, 2003) as in Jimenez et al. (2015). Roughly 40 mg of gill tissue from each individual was used.

We measured antioxidant capacity against peroxyl radicals (involved in the chain reaction leading to formation of LOOH) and catalase enzymatic activity against hydrogen peroxide. Gill tissue (0.1 g) was diced and put in a 1:9 (mass:volume) solution of 75 mmol l⁻¹ PBS (pH 7.0). Tissues were homogenized in a bead homogenizer for two 2-min rounds at 50 Hz and placed on ice for 2 min in between rounds. Homogenates were then centrifuged for 30 min at 17,900 g in 4°C, and supernatants were frozen at −80°C until use in assays.

Although this response is not universal (see Jimenez et al., 2016), antioxidant capacity against peroxyl radicals might be expected to increase following heat stress in order to prevent LOOH damage. A microplate-based version of the competitive ORAC assay (Cao et al., 1999; Prior and Cao, 1999) was used to estimate antioxidant capacities against peroxyl radicals according to the methods in Jimenez et al. (2016).

Catalase enzyme activity is predicted to increase under heat stress (Steare and Yellon, 1994; Almeida et al., 2015). Catalase activity was measured as described in Jimenez et al. (2016) by monitoring the rate of decrease in hydrogen peroxide absorbance at 240 nm with a temperature-controlled spectrophotometer (Beers and Sizer, 1952). Samples were run in triplicate; replicates were averaged to give a single measure for each individual.

Under heat stress conditions, the tissue contents of some osmolytes can increase (Li et al., 2009; Kumar et al., 2014; Wang et al., 2014; Ghaedi and Andrew, 2016), and these select osmolytes can act as thermoprotectants (Somero and Yancey, 1997; Yancey, 2005). Contents of organic osmolytes were determined from the mussels from the field sites. In a parallel laboratory study, osmolytes of common-gardened animals were quantified before, during and shortly after exposure to a single, acute heat stress event to test for any short-term changes in osmolyte contents.

Osmolytes were isolated from mussel gill using the extraction and high-performance liquid chromatography (HPLC) quantification method of Wolff et al. (1989). Gill samples (mean mass=150 mg) were thawed on ice, each side of each hemibranch segment was blotted twice on a laboratory wipe and the sample was weighed. Samples were homogenized in perchloric acid in a bead homogenizer for 3 min at 20 Hz plus 1 min at 30 Hz, then centrifuged to pellet proteins. The supernatants were neutralized with 2 mol l⁻¹ KOH and run through C18 hydrophobic resin cartridges to remove lipids. Five solute compounds were quantified by HPLC: alanine (Ala), glycine betaine (GlyBet), glucose (Glc), glycine (Gly) and taurine (Tau). Standards for each of these compounds were run in triplicate. Note that the contents reported in mmol kg⁻¹ are considerably lower (by a factor of approximately four) than the actual intracellular concentrations (mmol l⁻¹), owing to dilution with extracellular fluids that are generally low in organic osmolytes (Yancey, 2005). Because the tissue processing method was uniform across samples, this dilution does not impact the statistical interpretations.

Mussels 60–75 mm long (n=6 per time point) were collected from a wave-exposed site (MLLW +0.9 m) at HMS and held under common-garden conditions in a flow-through seawater system for 2 weeks (Dowd and Somero, 2013). During this time, mussels were continuously immersed at an average temperature of 14°C and fed every other day with a mixture of algal species (Shellfish Diet 1800©, Reed Mariculture, San Jose, CA, USA). ‘Baseline’ group mussels remained immersed in the flow-through system during the 4-h treatment period. ‘Heated’ and ‘4 h post-heating’ group mussels were emersed and exposed to heat stress using an electric coil heater and a Micro-Infinity temperature controller (Newport Electronics, Santa Ana, CA, USA). Temperature was ramped at +7°C h⁻¹ from 15°C to 36°C to simulate solar heating on a hot day during low tide (Tomanek and Somero, 1999; Lockwood et al., 2010; Denny et al., 2011). The mussels were held at 36°C for one additional hour before the heated group was dissected. The post-heating mussels were reimmersed in the flow-through seawater system for 4 h prior to dissection, in order to detect any osmolyte peaks that might occur during recovery from heat stress, as has been observed for Hsp70 (Gracey et al., 2008). Gill tissue was removed from the mussels, frozen in liquid nitrogen, and stored at −80°C for later extraction. Osmolytes were isolated from these gill tissue samples and quantified via HPLC using the techniques described above.

For our statistical analyses, each individual mussel was the unit of experimentation (we have individual data records for each mussel, not one sensor per acrylic plate). In addition, because a single plate was placed at each outplant site, we cannot make site-level generalizations to other tidepool, low or high intertidal sites.

For each of the nine physiological variables measured in the MusselTracker-instrumented animals, we used an automated, exhaustive model selection approach to select the terms that provided the most informative general linear model (GLM) using the glmulti package in R (Calcagno and Mazancourt, 2010). The full model included origin (exposed or protected), outplant location (high, low and tidepool) and their interaction as factors. Individual metrics of thermal history and valve gape were included as continuous predictors. Both of these metrics were included as predictors because they are not perfectly correlated in the field; for example, valve closure during aerial exposure does not always coincide with warm temperatures (Miller and Dowd, 2017). Several MusselTracker thermocouples were damaged prior to the end of the outplant period, preventing use of the complete 23-day dataset to generate metrics of thermal history. Therefore, to enhance the available sample size, we used individuals' maximum body temperatures on 20 July 2015, a warm day on which 26 of 30 thermocouples were still logging data, as our metric of thermal history. This value was very strongly positively correlated with other metrics of cumulative thermal stress (e.g. Pearson's r=0.89, 0.90 and 0.76 with mean daily maximum temperature, overall maximum temperature and degree-hours above 25°C, respectively; P<0.001 for each). For the metric of gaping behavior, we used the proportion of time with the valves spent in an open position (>20% of maximum) that would allow for gas exchange, feeding and/or waste excretion. As for thermal history, other metrics of gaping behavior were strongly positively correlated with the chosen metric (e.g. Pearson's r=0.89 with the proportion of time greater than 40% of maximum gape). Lastly, based on graphical exploration of the data, we included in the full model interactions between outplant location and thermal history and between outplant location and gape to allow for site-specific patterns. We also visually examined residual plots using the car package in R (Fox and Weisberg, 2011) to confirm that there was no evidence of heteroscedasticity. Lipid hydroperoxide data (LOOH) were Box–Cox transformed to satisfy the assumption of normality; all other variables were untransformed.

The 95% confidence set of best-ranked models (the set of models whose cumulative Akaike's information criterion weights add up to 0.95; Burnham and Anderson, 2002) for each physiological variable was chosen from 50 possible combinations of terms based on a modified version of Akaike's information criterion adjusted for small sample sizes (AICc) (Hurvich and Tsai, 1989). From this confidence set, the best statistical model with the lowest AICc was selected. Type III ANOVA tests were run for the best model for each physiological variable to determine statistical significance for each predictor variable, with α=0.05. Note that these P-values do not indicate whether this model with the lowest AICc value is statistically the best model; rather, for each model, these P-values indicate whether each particular predictor variable significantly influences the given physiological variable.

Because none of the best models had strong Akaike weights >0.9, conservative model coefficients for each physiological variable were derived via model averaging across the 95% confidence set. In addition, to further assess the relative importance of each predictor variable for each of the nine physiological variables, predictor weights (w) were also calculated (by summing the Akaike weights for each model in the confidence set in which that variable appears) based on all models in each variable's 95% confidence set (Burnham and Anderson, 2002). The predictor with the largest predictor weight is deemed to be the most important.

For each physiological variable, added-variable plots were constructed using the car package in R (Fox and Weisberg, 2011) to illustrate the effect of the predictor of interest while controlling for all other terms in the best statistical model (Wang, 1985; Cook and Weisberg, 1999). These plots compare two sets of residuals. The x-axis represents residuals resulting from regression of the focal predictor variable (e.g. thermal history) against the remaining predictor variables, and the y-axis represents residuals computed from regression of the response variable (e.g. DNA damage) against the predictor variables but omitting the focal predictor variable. A statistically significant slope indicates a significant contribution of the focal predictor in the full GLM, and the slope corresponds to that predictor's coefficient β in the statistical model.

We used principal components analysis (PCA) in the FactoMineR package in R (Le et al., 2008) to assess global patterns of individuals' physiological profiles. The nine measures of physiological status for each individual were used as inputs to this analysis. For each of the first three principal components, which explained a total of 64.9% of the overall variance in the physiological dataset, we determined the categorical and continuous predictors that significantly correlated with that principal component using the dimdesc function of FactoMineR (Husson et al., 2010).

Our data reveal mussels' physiological responses to sublethal environmental variation in the field. There were no mortalities attributed to abiotic conditions, despite some individuals reaching body temperatures that approached the LT₅₀ for this species, 38.2°C (Denny et al., 2011). However, three mussels at the high site were consumed by black oystercatchers, Haematopus bachmani, during low tides.

Mussels at the high intertidal site experienced considerably warmer maximum body temperatures and slightly cooler minimum body temperatures than individuals at either the low or the tidepool site (Table 1). We also observed a significant degree of inter-individual variation in daily maximum body temperatures within each outplant plate (Table 1; Table S1). The high site exhibited the most variation, with a mean range of individuals' daily maximum body temperatures of 7.0°C; on one day the difference between the warmest and coolest mussels within the high site exceeded 14°C (Table 1). The maximum temperature on 20 July, our metric for thermal history, ranged from 17.25°C for an individual at the low site to 35.75°C for a mussel at the high site.

Mussels at the high outplant site spent a significantly lower proportion of time gaping (mean=18.56%) compared with the tidepool (mean=60.98%) and low (mean=61.42%) outplant sites (one-way ANOVA, P<0.001; Tukey HSD post hoc comparisons against high site, P<0.001 for both). Across all three sites, the proportion of time spent gaping ranged from 12.10% in an individual at the high site to 71.89% for a mussel at the low site. Within-site, inter-individual variation in the proportion of time spent gaping was most pronounced at the tidepool and low outplant sites, with values ranging from ∼45 to 70% in each (Table S1).

The 95% confidence set for 8-OHdG DNA damage included 17 candidate models, and 11 of these included the maximum temperature on 20 July (i.e. thermal history) as a predictor (Table S2). The best model included origin site (exposed<protected), outplant site (tidepool and low<high), their interaction, thermal history and the interaction between outplant site and thermal history. All predictors except origin site had significant effects in this best model (Table 2). The variable with the highest predictor weight was outplant site (w=0.89; Table 2), and the eight mussels at the high intertidal site included the seven highest measures of oxidative DNA damage (Fig. S1A). Across all three sites, 8-OHdG was positively correlated with thermal history (averaged model coefficient±confidence interval: β=0.030±0.077 ng ml⁻¹ °C^–1 on 20 July; F_2,11=8.4, P=0.015). This relationship was site-specific, as indicated by the significant interaction between outplant site and thermal history and the strong positive slope among mussels at the high intertidal site (Table 2, Fig. 1A). Interestingly, all three individuals from the wave-exposed origin site that were deployed in the tidepool had significantly lower values of 8-OHdG (t-test, P=0.03, d.f.=3.70; Fig. S1A). This is likely a result of differences in gaping behavior; these individuals were open for a smaller percentage of time compared with wave-protected origin mussels in the tidepool (Table S1).

There were nine candidate models in the 95% confidence set (Table S2). The best statistical model for catalase activity included a positive effect of thermal history (β=9.17±13.33 U °C^–1 on 20 July; F_1,18=6.7, P=0.020) and a significant effect of origin site (protected>exposed; β_prot=147.36±85.79 U; F_1,18=1.9, P=0.002, d.f.=1; Table 2, Figs 1B, 2; Fig. S1C). Similarly, the two predictors with the highest predictor weight were origin site (w=0.914) and thermal history (w=0.75; Table 2).

Mussels' individual histories significantly influenced the contents of three of the five gill osmolytes. The best model for taurine content included thermal history and gape, although only the effect of thermal history was statistically significant (β=0.88±1.01 mmol kg⁻¹ °C^–1 on 20 July; F_1,16=12.3, P=0.003; Table 2, Fig. 1C; Fig. S1D). Thermal history was included in 11 of the 14 candidate models in the 95% confidence set, and thermal history was also the predictor variable with the highest predictor weight (w=0.853; Table 2). Gape was included in nine of the 14 candidate models (Table S2), and gape had the second highest predictor weight (w=0.625, β=28.03±66.71 mmol kg⁻¹; Table 2).

There were 16 candidate models in the confidence set for glycine betaine content (Table S2). The best model for glycine betaine content included positive effects of both thermal history (β=0.97±1.67 mmol kg⁻¹ °C^–1 on 20 July; F_1,16=10.0, P=0.006) and gape (β=44.03±73.79 mmol kg⁻¹; F_1,16=8.7, P=0.010; Table 2, Fig. 1D,E; Fig. S1E,G). Gape (w=0.785) and thermal history (w=0.685) also had the highest predictor weights (Table 2).

Of the 13 candidate models for glycine content, the best model only included a positive effect of gape (β=4.31±8.80 mmol kg⁻¹; F_1,17=14.8, P=0.001; Table 2, Fig. 1F; Fig. S1H). Gape also had the highest predictor weight (w=0.601; Table 2). Although the model in which glycine was regressed only on thermal history was statistically significant (β=−0.018±0.15 mmol kg⁻¹ °C^–1 on 20 July; F_1,17=9.4, P=0.007; Fig. S1F), the AICc value was more than three units larger for this model, and the predictor weight for thermal history was also much lower (w=0.215; Table 2).

There were 15 models in the 95% confidence set for lipid damage (Table S2). Although the best model indicated no influence of any of the predictor variables (Table 2), 10 of the 15 models in the candidate set included thermal history (β=2.810E–3±1.479E–3 mol l⁻¹ °C^–1 on 20 July), and the maximum temperature on 20 July also had the highest predictor weight (w=0.622; Fig. S1B; Table 2).

Model selection for the osmolyte alanine generated more ambiguity than in the six variables above. Out of the 15 candidate models (Table S2), the best statistical model for alanine content included outplant site only (low and tidepool<high, β_low=−0.42±1.37 mmol kg⁻¹, β_tidepool=−0.77±2.04 mmol kg⁻¹; F_2,16=8.7, P=0.003), and outplant site also had the highest predictor weight (w=0.432) (Table 2). However, three other variables (origin site, maximum temperature on 20 July and gape) had similar predictor weights (w_origin=0.242, w_temp=0.343, w_gape=0.327). Thermal history had a positive effect on alanine content and gape had a negative effect.

The best model for antioxidant capacities against peroxyl radicals (ORAC, eight candidate models) and gill glucose contents (nine candidate models) was the null model (Table 2; Table S2). Origin site had the highest predictor weight for both variables (w_ORAC=0.192, w_Glc=0.379; Table 2), although these values are very low in comparison to the other physiological variables.

The PCA distinguished the high site from the low and tidepool sites along the first principal component (PC1) dimension (ANOVA, F_2,22=36.8, P<0.001; Fig. S2); there was no significant effect of origin site along PC1 (F_1,23=0.75, P=0.397) or PC2 (F_1,23=1.41, P=0.247). PC1 explained 27.3% of the total variation in the physiological profiles among individual mussels (Fig. S2). The physiological variables that were significantly positively correlated with scores on PC1 were 8-OHdG, taurine and alanine; glycine content was strongly negatively correlated with PC1 (Fig. 3; Fig. S2). Individual thermal history also was significantly positively correlated with mussels' scores on PC1 (Pearson's r=0.83; Figs 3 and 4). In contrast, gape was significantly negatively correlated with mussels' scores on PC1 (Pearson's r=−0.65, P<0.001). Neither thermal history nor gape was significantly correlated with individual scores on PC2, which explained 23.4% of the overall physiological variation (Fig. S2). Contents of glycine betaine and taurine, as well as catalase activity and LOOH lipid damage, were positively correlated with scores on PC2. Alanine and glucose contents were negatively correlated with scores on PC2 (Fig. 3; Fig. S2).

No significant differences were detected among the three laboratory heat stress time points for any of the five measured osmolytes (Glc P=0.29, Tau P=0.11, Ala P=0.56, Gly Bet P=0.50, Gly P=0.90, d.f.=2 for all, one-way ANOVA; Table S3). For parameter estimates and standard errors, see Table S3.

We utilized novel MusselTracker technology to record in situ body temperature and gaping behavior of M. californianus individuals placed at high intertidal, low intertidal and tidepool sites for 23 days. When combined with biochemical measurements on the same individuals, these data indicate substantial environmentally mediated physiological variation among individual mussels. In support of our hypotheses, we found that higher body temperatures significantly correlate with increased oxidative DNA damage, as well as with catalase antioxidant capacity and taurine and glycine betaine osmolyte contents. Notably, these correlations emerged despite conducting the experiments in an uncontrolled field setting. Overall, inter-individual variation in recent thermal history appears to be a significant driver of inter-individual physiological variation within mussel beds.

This study took advantage of the fact that adult M. californianus are sessile, overcoming logistical constraints that have prevented development of similar datasets in other species. Previous studies, including several in mussels, have used biomimetic devices (Denny et al., 2011; Jimenez et al., 2015), measured live animals at only a single time point (Elvin and Gonor, 1979), or measured the temperature of the surrounding environment to approximate body temperature (Bay and Palumbi, 2014). Thus, these studies do not adequately capture the thermal experience of individuals. To date, we are not aware of any other studies that have continuously monitored field body temperatures of the same individuals that were sampled for subsequent biochemical analyses.

The MusselTracker system also measured individuals' valve gaping behavior in the field. Previous studies investigating differences in Mytilus gaping behavior were primarily conducted under controlled laboratory conditions or with small sample sizes for short periods of time in the field (Wilson et al., 2005; Anestis et al., 2007; Dowd and Somero, 2013; Olabarria et al., 2016). Therefore, our work provides novel insight into how gaping behavior over the course of multiple tidal cycles in the field affects inter-individual variation in physiology.

Univariate statistical analyses indicated that inter-individual variation in thermal history significantly affected four of the nine physiological variables measured. For all four of these variables, the best model included a positive correlation with thermal history. These included levels of oxidative DNA damage, as well as activity of the antioxidant enzyme catalase and content of the osmolytes taurine and glycine betaine.

Our biochemical analyses focused on oxidative damage and antioxidant capacities (Jimenez et al., 2015, 2016). Although reactive oxygen species (ROS) are a normal by-product of aerobic metabolism (Abele et al., 2002; Buttemer et al., 2010), elevated temperatures during low tide and/or cycles of hypoxia and reoxygenation owing to shell valve movements can cause enhanced production of these damaging intermediates (Abele et al., 2002; de Oliveira et al., 2005; Heise et al., 2006; Tomanek and Zuzow, 2010; Dowd and Somero, 2013; Rivera-Ingraham et al., 2013). (Because M. californianus close their valves during hot low tides in the field, to definitively disentangle the effects of temperature and valve movements on ROS, future manipulative laboratory studies will need to be performed that expose individuals to varying temperatures while keeping valve opening constant.) Overall, when the level of ROS being produced exceeds the capacity of antioxidant defenses, oxidative stress occurs in the form of damage to macromolecules such as proteins, lipids and DNA (Halliwell and Gutteridge, 2007). In the present study, oxidative damage to DNA was dependent on individual thermal history during the outplant period. Although results for oxidative lipid damage were more ambiguous, 10 of the 15 models in the 95% confidence set included an effect of body temperature. This potential relationship between lipid damage and body temperature is supported by previous work, which suggested lipid peroxidation in M. californianus rises rapidly with acute thermal stress and then returns to baseline levels within 24 h (Jimenez et al., 2016).

Patterns of osmolyte accumulation and catalase activity allow us to draw tentative conclusions regarding the kinetics of physiological responses to thermal history. First, the laboratory time-series experiment revealed that none of the five measured osmolyte contents changed after a single, acute heat stress. This suggests that organic osmolyte contents gradually adjust to the environmental context, rather than responding to a single low tide event. In addition, measurements of catalase enzyme activity in this study, together with previous field and laboratory results (Dowd et al., 2013), suggest that catalase activity is initially determined by the site in which the animal settles (origin site effect), although we cannot conclude yet whether such differences reflect effects of developmental environment on physiology or genetically predetermined variation (see below). What is clear is that subsequent thermal experience and/or food availability superimposes adjustments on these baseline differences in catalase over multiple tidal cycles. Although catalase activity is plastic, antioxidant capacity against peroxyl radicals seems not to be. This result from animals in the field agrees with previous work in a controlled laboratory setting, which detected differences among mussel tissues in their peroxyl radical scavenging capacities, yet found no evidence for an increase in peroxyl ORAC activity with a history of emersion at high temperature (Jimenez et al., 2016).

Although our data help to clarify the relationships between thermal history and physiological status, additional studies are required to determine the relative importance of the various metrics that can be used to quantify an individual's thermal experience [e.g. mean daily maximum temperature, degree-hours (i.e. thermal time) above some threshold, total number of warm days, overall maximum body temperature]. Because all thermal history metrics were strongly correlated, our experimental design did not allow us to determine which of these possible indices of thermal stress drives physiological responses. Future work must manipulate different thermal history metrics in controlled conditions to provide further insight.

Lastly, while this study strongly suggests that variation in temperature stress plays a significant role in generating physiological variation, previous work has shown that food availability also affects antioxidant enzyme capacities such as catalase in M. californianus (Dowd et al., 2013). Food availability to individual mussels could not be controlled in our field study. In addition, responses to hypoxia such as metabolic suppression, expression of heat shock proteins and/or increased antioxidant capacity that are associated with tight-valve closure physiology could generate resistance to extremes of another variable such as temperature (Shick et al., 1988; Anestis et al., 2010; Nogueira et al., 2017). Therefore, to confirm that variation in thermal history is indeed the primary driver of the observed physiological differences, further work should manipulate thermal history against a backdrop of constant food availability and gaping status.

The oxidative DNA damage marker 8-OHdG appears to be an indicator of accumulated macromolecular damage from thermal stress in M. californianus. Notably, 8-OHdG's positive relationship with thermal history was most pronounced at the high, warm intertidal site, in accord with the recent proposal that variation around a stressful mean temperature plays a significant role in eliciting physiological differences among individuals (Jimenez et al., 2015). Accumulation of 8-OHdG has been observed with increasing age (Gruber et al., 2015) and with exposure to toxicants (Steinert, 1999) in other marine bivalve molluscs. Here, any age-related effects on 8-OHdG levels, such as those that might exist between origin sites due to variation in size-at-age (Connor and Robles, 2015), were controlled for by including origin as a predictor in the model selection procedure. To our knowledge, no other studies have investigated the relationship between oxidative DNA damage and temperature stress in bivalves. Our data indicate that 8-OHdG is a useful metric for examining the effects of cumulative thermal stress on wild populations.

Heat stress also causes other forms of DNA damage in marine mussels, yet the genotoxic effects of abiotic stress have received relatively little attention. For example, acute heat stress induced both single- and double-stranded breaks in DNA in laboratory studies of Mytilus congeners (Yao and Somero, 2012). The physiological costs, such as those related to DNA repair (Yao and Somero, 2013), and possible mutagenic consequences of heat-induced DNA damage merit further study.

Although a handful of studies have demonstrated that osmolytes accumulate under laboratory heat stress conditions in wheat, yeast and aphids (Li et al., 2009; Kumar et al., 2014; Wang et al., 2014; Ghaedi and Andrew, 2016), to our knowledge this is the first study showing that the content of an osmolyte molecule with known thermoprotective properties is dynamically regulated in response to thermal stress in nature. Individual thermal history was significantly and positively correlated with taurine content. Taurine also had the highest overall content out of the gill osmolytes measured, as has been observed in the osmotic profiles of other shallow-dwelling benthic molluscs (Yancey, 2005). These results are consistent with a role for taurine as an endogenous thermoprotectant in mussels.

Accumulation of taurine could provide several biochemical benefits. This sulfonic ß-amino acid is a protein-structure stabilizer that, at least for the model enzyme lysozyme, was found to be the most effective thermoprotectant among the common osmolytes (Arakawa and Timasheff, 1985). Thermoprotectants are ‘extrinsic’ stabilizers that act to maintain appropriate marginal stability of proteins and exert similar stabilizing effects across the proteome (Somero and Yancey, 1997; Yancey, 2005; Somero et al., 2017). Taurine's positive relationship with thermal history could also be related to stabilization of phospholipid membranes via its interactions with embedded proteins (Huxtable, 1992). Heat stress over the tidal cycle drives homeoviscous remodeling of cellular phospholipid membrane composition in M. californianus (Williams and Somero, 1996); taurine might further stabilize these remodeled membranes. Moreover, taurine may act indirectly as an antioxidant, possibly by regulating components of the mitochondrial electron transport chain (Jong et al., 2012). For instance, taurine decreased oxidative damage and increased total DNA recovery after damage in calf thymus DNA (Messina and Dawson, 2000). Alternatively, the correlation of taurine content and individual thermal history could relate to taurine's role as a prominent osmoregulatory agent (Lange, 1963). Heat stress has been documented to impair osmoregulatory ability (Tang et al., 2014), and the maintenance of active osmoregulation has been hypothesized to contribute to thermal tolerance (Jian and Huang, 2001; Xu et al., 2013). Further mechanistic work is required to document the relative importance of these various putative benefits of taurine accumulation under episodic exposure to thermal stress.

For tissue contents of taurine to rise, there must be an increase in its retention and/or production. The observed accumulation of taurine with heat stress could be driven by changes in the abundance of the transmembrane taurine transporter (TAUT) (Huxtable, 1992; Hosoi et al., 2005, 2007; Lin et al., 2016). The expression of TAUT mRNA increases under moderate heat stress in a deep-sea mussel (Nakamura-Kusakabe et al., 2016), and transcription of TAUT is induced by heat shock factor 1 in mice (Jung et al., 2013). Following thermal stress, taurine production may increase as irreversibly damaged proteins are degraded into their constituent amino acids, which may be reused or further modified. Taurine is derived from cysteine in a temperature-dependent fashion in mammals and fungi (Soboleva et al., 2004; Kumar et al., 1983). If mussels accumulate thermoprotective breakdown products such as taurine after episodes of proteotoxic heat stress, then a de facto acclimatization mechanism is in action.

The best model for glycine betaine content included positive effects of both thermal history and gaping, complicating interpretation of this osmolyte's mechanistic significance. High contents of glycine betaine, together with its activity as a weak to moderate thermoprotectant (Gopal and Ahluwalia, 1993; Somero and Yancey, 1997; Auton et al., 2006), suggest it could fulfill a thermoprotectant role similar to that of taurine. For example, addition of glycine betaine to cultured mammalian kidney cells suppresses HSP70 elevation during heat shocks (Sheikh-Hamad et al., 1994). Indeed, glycine betaine and taurine contents were positively correlated across individual mussels (Pearson's r=0.80, P<0.001, d.f.=23). However, unlike taurine, glycine betaine content was also correlated positively with gape, and it did not load significantly on PC1.

Glycine was distinct in the PCA analysis as the only osmolyte whose content loaded negatively on PC1, corroborating its positive univariate relationship with gape. Given the inverse relationship between glycine and thermal history on PC1, it is plausible that glycine decreases in thermally stressed individuals at the high site to offset accumulation of taurine and maintain osmotic balance. However, the weak negative correlation between taurine and glycine contents across individuals was insignificant (Pearson's r=−0.33, P=0.104, d.f.=23). Glycine is a prominent osmotic effector in other bivalves, but its kinetics in response to stress events vary markedly across taxa (Bishop et al., 1994; Paynter et al., 1995). In the cases of both glycine betaine and glycine, production and clearance kinetics need to be determined at higher temporal resolution over the course of events characterized by thermal stress, oxygen limitation or both.

The statistical results were the most ambiguous for the osmolyte alanine. The best univariate model for alanine content included outplant site only, although three other variables had roughly equal predictor variable weights. In the PCA, alanine content was positively correlated with PC1 scores and negatively correlated with PC2 scores; thermal history was positively correlated and gape history was negatively correlated with PC1. Though alanine can be a weak protein stabilizer, the contents estimated in M. californianus gill tissue (2.7–7.6 mmol kg⁻¹) were probably not high enough to impact general protein stability. The changes in alanine content more likely relate to anaerobic metabolism during emersion. Alanine is made by transamination of the glycolytic end-product pyruvate, and both compounds accumulate in mussels under anoxic/emersion conditions (Isani et al., 1995; Connor and Gracey, 2012).

Overall, alanine's known correlation with anaerobic metabolism could indicate that variation in its content is more linked to vertical location in the intertidal zone (here, outplant site) and to corresponding constraints on gaping behavior, rather than to body temperature. Specifically, mussels at the high intertidal site used in this study experienced longer periods of emersion and spent a significantly lower proportion of time with the valves open; however, we cannot make any generalizations to other high intertidal sites because our study was not replicated at the level of outplant site. The links among gaping behavior, reliance of mussels on anaerobic metabolism during emersion (Bayne et al., 1976) when they also are prone to experience thermal stress, and accumulation of alanine in the tissues are potentially complex. Ultimately, further manipulative experiments are needed in order to clearly distinguish the effects on alanine contents of heat stress during emersion and variation in oxygen availability related to gaping behavior.

Mussels in this study were collected from two origin sites (wave-exposed and wave-protected), but there was little evidence for origin-specific biochemical patterns. Of the nine physiological variables measured, only one, catalase enzyme activity, was unambiguously affected by origin site. Similarly, origin site was uncorrelated with overall physiological profiles in the multivariate analysis. Although individuals at the wave-exposed site experienced more intense wave action and cooler body temperatures during their post-settlement life compared with mussels from the wave-protected site (Denny et al., 2011), these factors appear to be largely superseded by recent experience. Adult phenotypic plasticity, particularly in the antioxidant system, is one plausible mechanism for these observations. Similarly, adult mussels from the same two origin sites exhibited comparable degrees of mean plasticity of antioxidant capacities in another recent study (Jimenez et al., 2015). Interestingly, the results of that study did implicate origin as a factor contributing to inter-individual variation in physiological plasticity. Specifically, the magnitude of variation around the mean antioxidant status changed more among mussels from the protected site upon manipulation of the magnitude of variation in body temperature (Jimenez et al., 2015). Furthermore, M. californianus individuals from these two origin sites differed in their mean metabolic rates even after common gardening (protected<exposed). Overall, these current and previous findings suggest that a complex set of interacting factors determines an individual’s current level of defense against environmental stress. Characterizing the relative influences of genetic constraint, developmental environment and adult physiological plasticity in response to recent history represents a major challenge for environmental physiologists, particularly in the context of climate change. Future work must clarify the time scales over which environmental histories influence physiology. Tools such as the MusselTracker system provide an outstanding opportunity to achieve this goal, for example by correlating repeated measures of individuals' biochemical status with measures of their in situ experience.

Individual-level monitoring provides novel insight into environmental and behavioral drivers of physiological variation within populations. Mussels' physiological profiles seemed to respond in a plastic fashion to recent experience in the field; however, repeated samples of the same individuals in different environmental contexts are needed to demonstrate conclusively that compensatory phenotypic plasticity is occurring. To our knowledge this is the first study to demonstrate a link between individual thermal history in the field and the tissue content of organic osmolytes with thermoprotective properties. Furthermore, variation in behavior among individuals also contributes to physiological variation, particularly when a behavior such as valve gaping so intimately relates to physiological function (in this case, access to oxygen for aerobic respiration). Going forward, more work examining inter-individual physiological and behavioral variation is needed in field settings in order to make accurate predictions regarding physiological, ecological and evolutionary consequences of a changing environment.

Shaina Alves provided field and laboratory support. W.W.D. would like to thank the faculty and staff of Hopkins Marine Station, especially Mark Denny and Steve Palumbi, for hosting the summer field season. Two anonymous reviewers and an expert statistician provided valuable comments on the manuscript.