Variation in developmental temperature alters adulthood plasticity of thermal tolerance in Tigriopus californicus

As the earth warms, organisms are increasingly impacted by the effects of high environmental temperatures (e.g. Wiens, 2016; Cohen et al., 2018; Pinsky et al., 2019). Indeed, the geographic range limits of many species have already shifted as a result of anthropogenic climate change, and in general these shifts have been towards regions with cooler temperatures (e.g. Parmesan and Yohe, 2003; Perry et al., 2005; Chen et al., 2011). The extent to which these effects have occurred, and will continue to occur, depends largely on the adaptive and plastic capacities of organisms to adjust key physiological traits, such as thermal tolerance limits (especially in aquatic ectotherms; Sunday et al., 2012; Pinsky et al., 2019), in response to increased temperatures (e.g. Crain et al., 2008; Somero, 2010; Bay et al., 2017; Kellermann and van Heerwaarden, 2019). In particular, given that rapid phenotypic changes are necessary because of high rates of environmental change (e.g. Barrett and Hendry, 2012; Fox et al., 2019), phenotypic plasticity may play a critical role in the resilience of populations and species to the effects of climate change (Merilä and Hendry, 2014; Seebacher et al., 2015; Donelson et al., 2019; Morley et al., 2019).

Phenotypic plasticity occurs across life stages and generations (e.g. Kelly et al., 2011; Schulte et al., 2011; Beaman et al., 2016; Burggren, 2015). For example, temperatures experienced during development or adulthood often have irreversible or reversible effects on physiological traits (e.g. Schulte et al., 2011; Beaman et al., 2016), and multi- or trans-generational effects of thermal variation are commonly observed (e.g. Crill et al., 1996; Massamba-N'Siala et al., 2014; Zizzari and Ellers, 2014; Donelson et al., 2018). Thus, physiological phenotypes have the potential to be shaped by effects of plasticity across different life stages. However, compared with the effects of adaptation on phenotypic plasticity (e.g. Crispo, 2007; Hendry, 2015; Donelson et al., 2019; Kelly, 2019), effects of plasticity at one life stage on the expression of plasticity at another life stage have received relatively little attention (Beaman et al., 2016). That said, developmental conditions are known to alter the adulthood plasticity of several traits (reviewed in Beaman et al., 2016). For instance, adult plasticity of swimming performance and metabolic rate depends on developmental environment in mosquitofish (Gambusia holbrooki; Seebacher et al., 2014; Seebacher and Grigaltchik, 2015). Yet, despite the likely biogeographic importance of thermal tolerance limits (Sunday et al., 2012), and many published examples of thermal tolerance limit plasticity in ectothermic organisms as a result of developmental or adulthood temperatures (e.g. Stillman and Somero, 2000; Ford and Beitinger, 2005; Fangue et al., 2006; Angiletta, 2009; Overgaard et al., 2011; Cooper et al., 2012a,b; Tepolt and Somero, 2014; Jakobs et al., 2015; Troia et al., 2015; Kingsolver et al., 2016; Pereira et al., 2017; Diamond et al., 2018; Mueller et al., 2019; Yanar et al., 2019), relatively few studies have assessed the potential for developmental temperatures to shape the phenotypic plasticity of upper thermal tolerance in adults (although see Schaefer and Ryan, 2006; Kellermann et al., 2017; Kellermann and Sgrò, 2018). Here, we examined these effects, and their potential mechanistic basis in populations of the intertidal copepod Tigriopus californicus.

Tigriopus californicus are small (∼1.2 mm) harpacticoid copepods with short generation times (3–4 weeks) that inhabit supralittoral tidepools along the west coast of North America from Baja California, Mexico, to southern Alaska, USA. Populations of this species occur on rocky outcrops isolated by sandy beaches, conditions that result in very low gene flow and high levels of genetic divergence among populations (Burton and Lee, 1994; Burton, 1997, 1998; Edmands, 2001; Peterson et al., 2013; Pereira et al., 2016; Barreto et al., 2018). Although much of this divergence is likely a result of small effective population sizes and genetic drift acting on selectively neutral variation, signatures of directional selection have been detected across the transcriptome (Pereira et al., 2016). This suggests that at least a portion of the genetic differentiation among populations is likely adaptive. Moreover, several common-garden studies in laboratory-raised individuals have demonstrated differences in upper and lower thermal tolerance limits that are consistent with local thermal adaptation in response to the latitudinal temperature gradient across the species range (Willett, 2010; Kelly et al., 2012; Wallace et al., 2014; Pereira et al., 2014, 2017; Leong et al., 2018; Willett and Son, 2018; Foley et al., 2019). This variation among populations has also been associated with genetically based differences in the function and regulation of heat shock protein genes (Schoville et al., 2012; Barreto et al., 2015; Tangwancharoen et al., 2018) and in the maintenance of mitochondrial ATP synthesis rates at high temperatures (Harada et al., 2019). Few studies have examined temperature-mediated phenotypic plasticity in these traits in T. californicus. However, elevated developmental temperature is known to increase upper thermal tolerance regardless of population (Kelly et al., 2012, 2017; Pereira et al., 2017), and adult plasticity in this trait is thought to be limited [although only relatively short acclimation periods have been examined (e.g. 1 day); Pereira et al., 2017]. Taken together with short generation times and ease of laboratory culture, these observations make T. calfornicus an ideal study system in which to investigate the effects of developmental temperature on adulthood plasticity in an aquatic ectotherm.

In the present study, we used laboratory-raised T. californicus to test two hypotheses: (1) variation in developmental temperatures changes the expression of phenotypic plasticity of upper thermal tolerance in adults, and (2) the physiological mechanisms involved in local thermal adaptation among populations are also involved in thermal limit plasticity. We first expanded our previous study (Harada et al., 2019) that put forward methods to estimate upper thermal tolerance with critical thermal maximum (CT_max) measurements in this species. We then used this method to facilitate experiments examining the effects of developmental temperature on the plasticity of this proxy for upper thermal tolerance in adults of four Californian populations of T. californicus. Finally, we assessed the effects of developmental and adulthood temperatures on mechanisms involved in local thermal adaptation in this species: the thermal performance curve of ATP synthesis rate, and the mRNA expression levels of heat shock protein genes and mitochondrial-encoded genes following acute heat stress.

Adult copepods [Tigriopus californicus (Baker 1912)] were collected from supralittoral tidepools across 10 locations along the west coast of North America, which spanned ∼21.5° of latitude: San Roque, Mexico (SR), La Bufadora, Mexico (BF), San Diego, California (SD), Bird Rock, California (BR), Abalone Cove, California (AB), Estero Bay, California (EB), San Simeon, California (SS), Santa Cruz, California (SC), Pescadero, California (PE), and Pacific Crest, Canada (PC) (Table S1, Fig. S1A,B). Collected animals were transported to Scripps Institution of Oceanography (La Jolla, CA, USA) in 1 litre plastic bottles containing seawater obtained from the same tidepools. The collection for each location was divided across several laboratory cultures, which were maintained at 20°C, 36 ppt and 12 h:12 h light:dark photoperiod using filtered seawater and deionized water to adjust salinity as necessary. Laboratory cultures were maintained for at least two generations (∼2 months) prior to experiments. During laboratory acclimations and experimental treatments, copepods consumed natural algal growth within the cultures, as well as a mixture of ground spirulina (Salt Creek, Inc., South Salt Lake City, UT, USA) and TetraMin Tropical Flakes (Spectrum Brands Pet LLC, Blacksburg, VA, USA) that was added approximately once per week.

Upper thermal tolerance was estimated by critical thermal maximum (CT_max) trials using loss of locomotor performance as the assay endpoint (as in Harada et al., 2019). In brief, 16 adult copepods of each population (eight females and eight males for all populations; divided across five trials) were transferred to 10 cm Petri dishes containing filtered seawater (20°C and 36 ppt) with no food overnight. In the morning, copepods were individually transferred into 0.2 ml strip tubes with 100 μl of water from the Petri dish. Tubes were left uncapped, and were placed in an Applied Biosystems SimpliAmp™ Thermal Cycler (Thermo Fisher Scientific, Waltham, MA, USA). After 5 min at 20°C, the temperature was increased using the AutoDelta function at rates of 0.1°C per 20 s from 20 to 32°C, and 0.1°C per minute from 32°C to the temperature at which the last individual in the trial lost locomotor performance. Loss of performance (i.e. knockdown) was monitored by cycling 40 μl of water in each tube with a pipettor. Typically, this procedure results in erratic swimming behaviour in T. californicus; however, at extremely high temperatures, this swimming response ceases, and copepods passively sink to the bottom of their tube. Endpoints were determined when an individual did not respond to three sequential tests with the pipettor, and CT_max was recorded as the temperature at which the endpoint was observed. After CT_max was determined, copepods were returned to 10 cm Petri dishes with 20°C filtered seawater (36 ppt) for recovery, and survivorship assessed 8 h after the trials was >90%.

To assess variation in the phenotypic plasticity of upper thermal tolerance in adult copepods as a result of differences in developmental temperature, gravid SD and BR females with mature (i.e. red) egg sacs were removed from laboratory cultures (24 females from four cultures per population). Egg sacs were dissected from the females (which synchronizes hatching), placed individually in wells of 6-well plates containing filtered seawater (20°C and 36 ppt), and allowed to hatch overnight. In the morning, nauplii (i.e. larvae) from each egg sac were counted, and split evenly across six treatments in 10 cm Petri dishes (Fig. 1). Three treatments were developed at 20°C for 14 days, and three treatments were developed at 25°C for 10 days. Preliminary trials with the SD population determined that these developmental times were those required for the majority of individuals to reach adulthood (and to observe the first gravid female) at each temperature, suggesting that the temperature coefficient (Q₁₀) of developmental rate equals ∼2 in this species. At the end of the developmental periods, the developmental treatments at each temperature were transferred to one of three adult acclimations: 20°C for 14 days, 25°C for 10 days or 25°C for 14 days. These lengths of acclimation were chosen to allow 2 weeks of acclimation at 20°C, and comparisons between equivalent acclimations using either absolute time or physiologically adjusted time at 25°C (assuming a continued Q₁₀ of ∼2 for life history traits). On the days that the adult acclimations were completed, CT_max was determined as described above for 16 copepods from each treatment and population. Note that individuals used in the tolerance trials were transferred to fresh filtered seawater without food at their acclimation conditions on the evening before the end of the acclimation treatments (i.e. the day before trials), and CT_max trials started from the acclimation temperatures in all cases.

To examine the potential for local thermal adaptation of the effects of developmental temperature on adult thermal tolerance plasticity, we performed a second experiment beginning with gravid females from the SC and PE populations. This experiment was conducted as described above for the SD and BR populations; however, the 25°C–10 day adult acclimation treatments were excluded, meaning egg sacs for each population were split across four treatments in total (Fig. 1). Again, CT_max was determined for 16 copepods for all treatments except the PE 25°C development and 25°C adulthood treatment, for which n=15.

To examine plasticity of the thermal performance curve for ATP synthesis rate, we compared two temperature treatments: 20°C for both development and adulthood versus 25°C for both development and adulthood (Fig. 1). Gravid SD and BR females carrying mature egg sacs (60 per population) were transferred from laboratory cultures to 10 cm Petri dishes (six per population) containing ∼60 ml of filtered seawater (20°C and 36 ppt) with food. The egg sacs from the majority of these females (6–10 per plate) hatched overnight. All females were removed in the morning; egg sacs that were still carried by females were dissected free and returned to their respective dishes. Dissected egg sacs hatched within 3 h, and once all egg sacs had hatched, three Petri dishes for each population were transferred to 25°C. As described above, juveniles in the 20 or 25°C dishes developed for 14 or 10 days, respectively, and adult acclimations at the two temperatures were also 14 or 10 days, respectively.

ATP synthesis rates were measured at 20, 25, 30, 33, 35 and 37°C using procedures similar to those of Harada et al. (2019). On the day before the end of the adult acclimation treatments, groups of 32 copepods (six groups per population×treatment) were held at their acclimation temperatures in 10 cm Petri dishes with filtered seawater (36 ppt) and no food overnight. In the morning, the groups of copepods were rinsed with 200 μl homogenization buffer (400 mmol l⁻¹ sucrose, 100 mmol l⁻¹ KCl, 70 mmol l⁻¹ Hepes, 6 mmol l⁻¹ EGTA, 3 mmol l⁻¹ EDTA, 1% w/v BSA, pH 7.6), which had been chilled on ice. Each group was transferred to a 2 ml glass Teflon homogenizer, and homogenized in 800 μl of fresh buffer. Following homogenization, mitochondria were isolated by differential centrifugation in 1.5 ml microcentrifuge tubes (Eppendorf, Hamburg, Germany). First, the tubes were centrifuged at 4°C and 1000 g for 5 min, and supernatants were transferred to fresh tubes. Second, the new tubes were centrifuged at 4°C and 11,000 g for 10 min. The resulting supernatants were discarded, and mitochondrial pellets were resuspended in 205 μl assay buffer (560 mmol l⁻¹ sucrose, 100 mmol l⁻¹ KCl, 70 mmol l⁻¹ Hepes, 10 mmol l⁻¹ KH₂PO₄, pH 7.6). Isolated mitochondria were divided into eight 25 μl aliquots: six for synthesis reactions (one per temperature), one for initial ATP concentration determination, and one for measuring DNA content, which was used to normalize ATP synthesis rate. DNA content was assayed with Invitrogen™ Quant-iT™ PicoGreen™ dsDNA reagent following the manufacturer's protocols (Thermo Fisher Scientific).

ATP synthesis reactions were conducted in 0.2 μl strip tubes, and were initiated by adding 5 μl of a saturating substrate cocktail (final assay substrate concentrations 5 mmol l⁻¹ pyruvate, 2 mmol l⁻¹ malate, 10 mmol l⁻¹ succinate and 1 mmol l⁻¹ ADP in assay buffer), resulting in electron donation to both complex I and complex II (CI+II) of the electron transport system. Following substrate addition, tubes were immediately transferred to an Applied Biosystems SimpliAmp™ Thermal Cycler and incubated at the desired assay temperatures for 10 min. At the end of the reactions, 25 μl of each assay was added to 25 μl of CellTiter-Glo (Promega, Madison, WI, USA), which stops ATP synthesis and is used for ATP quantification. To determine initial ATP concentrations in the assays, one aliquot of each mitochondrial isolation was added to CellTiter-Glo immediately following substrate addition. All assays were held in the dark at room temperature for 10 min after addition of CellTiter-Glo, and then luminescence was determined with a Fluoroskan Ascent^® FL (Thermo Fisher Scientific). ATP concentrations were calculated by comparison to a prepared standard curve (5 nmol l⁻¹ to 10 μmol l⁻¹ in assay buffer), and synthesis rates at each temperature were determined by subtracting the initial ATP concentration from the final ATP concentrations at each temperature for each mitochondrial isolation.

Variation in gene expression following heat shock was assessed for the same treatments as those used to examine ATP synthesis rates: 20°C for both development and adulthood, and 25°C for both development and adulthood (Fig. 1). Again, offspring from 60 SD and 60 BR females were divided between these treatments (repeated as described above). In the evening prior to the last day of the adult acclimations, groups of 15 copepods (18 groups per population×treatment) were transferred to 15 ml Falcon™ conical tubes (Thermo Fisher Scientific) containing 10 ml of filtered seawater (36 ppt) with no food at the acclimation temperature of the copepods. In the morning, tubes were transferred to water baths held at 35 or 36°C for 1 h (six per population×treatment at each temperature), and then returned to the acclimation temperature of the copepods for 1 h as in Barreto et al. (2015). The remaining six tubes (per population×treatment) were handled in the same manner, but were kept at the acclimation temperature of the copepods for the entire 2 h. At the end of all heat shock trials, copepods were frozen at −80°C until RNA isolation.

Briefly, RNA was isolated using TRI Reagent^® (Sigma-Aldrich, St Louis, MO, USA) with half-volume reactions according to the manufacturer's instructions. RNA pellets were resuspended in 12 μl of Invitrogen™ UltraPure™ DNase/RNase-Free Distilled Water (Thermo Fisher Scientific) and were incubated at 56°C for 5 min. Isolations were treated with DNase using Invitrogen™ TURBO DNA-free™ Kits (Thermo Fisher Scientific) following the supplied protocols, and RNA concentration was determined with a NanoDrop spectrophotometer (Thermo Fisher Scientific). RNA integrity was confirmed by gel electrophoresis using two high-concentration samples. A total of 100–150 ng of total RNA was used to synthesize cDNA for each sample with Applied Biosystems™ High-capacity RNA-to-cDNA™ Kits (Thermo Fisher Scientific) according to the manufacturer's instructions, and the resulting cDNA samples were normalized to 2 nginput RNA μl⁻¹.

The mRNA expression levels of heat shock protein beta 1 (hspb1), heat shock protein 70 (hsp70), mitochondrial-encoded ATP synthase membrane subunit 6 (mt-atp6) and glyceraldehyde 3-phosphate dehydrogenase (gapdh) were assessed by quantitative real-time polymerase chain reaction (qRT-PCR). Primers for hspb1, hsp70 and gapdh for the SD population were obtained from Barreto et al. (2015). If necessary, owing to single nucleotide polymorphisms between the populations, equivalent primers were designed for the BR population using a population-specific reference genome (Barreto et al., 2018). Primers for mt-atp6 for each population were designed using population-specific mitochondrial genomes (DQ913891; Burton et al., 2007; Barreto et al., 2018). All primer sequences are listed in Table 1. Fifteen-microlitre qRT-PCR reactions were prepared in duplicate with 4 μl cDNA, 5 pmol of each primer and 7.5 μl iTaq Universal SYBR Green Supermix (Bio-Rad Laboratories, Inc., Hercules, CA, USA). All reactions were conducted in an AriaMx Real-time PCR System (Agilent Technologies, Inc., Santa Clara, CA, USA) with the following protocol: 95°C for 2 min, then 95°C for 10 s followed by 58°C for 30 s for 40 cycles. The presence of a single amplicon was confirmed by a melting curve analysis after each reaction. Samples for each population were quantified relative to population-specific five-point standard curves that were included on all reaction plates, and were prepared by serial dilution (1× to 1/625×) of a high-concentration heat shock sample from each population. Transcript levels of hspb1, hsp70 and mt-atp6 were then expressed relative to those of gapdh, which has been confirmed to be an appropriate housekeeping gene for heat shock studies in T. californicus (Schoville et al., 2012; Barreto et al., 2015; Harada and Burton, 2019). Final qRT-PCR sample sizes for the majority of our treatments and genes were n=6; however, for some groups, n=4 or 5 owing to a combination of insufficient RNA for cDNA synthesis, failed reactions as a result of extremely low hspb1 expression levels under control conditions, or insufficient cDNA to assess all genes (one instance resulting in no estimate for mt-atp6). As a result, the final sample sizes for all qRT-PCR data are presented in detail in Table S2.

All analyses were performed with R v3.4.0 (https://www.r-project.org/) with α=0.05. Latitudinal variation in CT_max across all populations failed to satisfy the assumptions for parametric statistics even after log transformation. Thus, differences among populations were assessed by Kruskal–Wallis ANOVA followed by Nemenyi post hoc tests. Potential effects of sex on CT_max were assessed by Wilcoxon rank-sum tests within each population. In contrast, variation in CT_max associated with developmental and adulthood temperatures met the assumptions for parametric tests in both the SD and BR, and the SC and PE experiments. These data were assessed by general linear models followed by ANOVAs with population, developmental temperature and adult acclimation treatment as factors. Post hoc comparisons among groups were performed with Tukey tests. After log transformation to meet assumptions of normality and homogeneity of variance, variation in ATP synthesis rate was assessed with a mixed-effect linear model followed by ANOVA with fixed effects of population, acclimation treatment and assay temperature, and a random effect of mitochondrial isolation. All interactions between factors were not significant in the initial model (P≥0.12 for all), and were removed from the final model for this test. Planned pairwise post hoc comparisons were conducted between assay temperatures within each population×acclimation treatment, between acclimation treatments within each population×assay temperature, and between populations within each acclimation treatment×assay temperature with Student's t-tests (84 comparisons). The resulting P-values were corrected for multiple tests with the Benjamini–Hochberg method (Benjamini and Hochberg, 1995). As an alternative method to examine variation in the thermal sensitivity of ATP synthesis rate, rates were normalized to the 25°C rate within each mitochondrial isolation. Variation in normalized ATP synthesis rate was assessed at assay temperatures from 30 to 37°C with similar methods to those described above for the unnormalized rates (although log transformation and removal of interactions among model factors were not required). Finally, mRNA expression data were all log transformed to meet necessary assumptions, and then differences among groups were examined by two-way ANOVAs with acclimation treatment and heat shock exposure as factors followed by post hoc Tukey tests. Note that comparisons of gene expression between populations were not made, because the expression levels of each population were quantified relative to standard curves that were population specific and for some genes the qRT-PCR primers were population specific as well (Table 1). Full ANOVA tables for all models tested have been uploaded to the Dryad Digital Repository (https://doi.org/10.5061/dryad.np5hqbzp0; Healy et al., 2019).

CT_max demonstrated significant variation among populations (P<2.2×10⁻¹⁶; Fig. 2), and there were no differences in this trait between females and males in the present study (P≥0.19 within all populations). Overall, CT_max increased from northern to southern populations (Fig. 2), suggesting that variation in this thermal limit generally parallels previously published latitudinal variation in lethal temperatures among populations (Willett, 2010; Kelly et al., 2012; Pereira et al., 2014; 2017; Leong et al., 2018; Willett and Son, 2018; Foley et al., 2019). This makes CT_max ideal for examining effects of phenotypic plasticity on upper thermal tolerance in experiments utilizing designs where offspring from individual egg sacs are divided among treatments.

Development at 20 or 25°C resulted in variation in CT_max among SD and BR copepods that was significantly affected by a three-way interaction among population, developmental temperature and adult acclimation treatment (P=0.03). This effect was resolved by post hoc tests, although there was little evidence for differential effects between the populations (Fig. 3A). For all SD and BR copepods that developed at 20°C, CT_max values were similar regardless of adult acclimation temperature or time (range of means±s.e.m.: 37.9±0.1 to 38.3±0.1°C; P=0.054 for 20°C-developed SD acclimated to 20°C for 14 days versus 25°C for 10 d as adults, and P≥0.22 for all other comparisons), suggesting that in 20°C-developed T. californicus there is no upper thermal limit plasticity in adults. In contrast, relative to development at 20°C, development at 25°C resulted in significant increases in CT_max for both SD and BR when adults were also acclimated at 25°C (range of means±s.e.m.: 38.5±0.1 to 39.1±0.1°C; P≤0.04 for all comparisons within populations) and these effects were similar at 10 and 14 days of acclimation (P≥0.16 between times for both populations). However, if 25°C-developed SD and BR copepods were acclimated at 20°C as adults, there was a significant loss of tolerance (i.e. decrease in CT_max) in both populations (means±s.e.m.: 38.1±0.1 and 37.7±0.1°C for SD and BR, respectively; P<0.01 for all within-population comparisons).

To explore these effects in populations known to be locally adapted to lower temperatures than the SD and BR populations, we conducted a second experiment examining plasticity of CT_max in the SC and PE populations. In SC and PE copepods, variation in adult CT_max was unaffected by the three-way interaction among population, developmental temperature and adult acclimation treatment (P=0.60). Moreover, there were no significant effects of interactions between population and either developmental temperature (P=0.83) or adult acclimation treatment (P=0.23), or of population alone (P=0.20). However, the interactive effect of developmental temperature and adult acclimation treatment significantly affected CT_max (P=9.9×10⁻⁵). Again, post hoc comparisons resolved this effect, and suggested similar patterns of variation among treatments for SC and PE as those described above for SD and BR (Fig. 3B). There was no significant variation in CT_max as a result of adult acclimation temperature in 20°C-developed copepods from either SC or PE (means±s.e.m.: 36.4±0.1 and 36.8±0.1°C for SC, and 36.8±0.1 and 37.1±0.1°C for PE; P≥0.39 between adult temperatures for both populations). In contrast, development at 25°C resulted in significantly higher CT_max in both SC and PE copepods, but only when adults were acclimated at 25°C (means±s.e.m.: 38.0±0.2 and 37.6±0.1°C for SC and PE, respectively; P≤0.04 for all comparisons within populations). If 25°C-developed copepods from either population were acclimated to 20°C as adults, there was a significant decrease in CT_max (means±s.e.m.: 36.7±0.2 and 36.7±0.01 for SC and PE, respectively; P<0.01 for both within-population comparisons).

Maintaining T. californicus at 20 or 25°C for both development and acclimation as adults resulted in significant variation in the thermal performance curve for CI+II ATP synthesis rates in isolated mitochondria (Table 2). Specifically, these curves were affected by population (P=4.4×10⁻⁵), development and adult acclimation temperature (P<2.2×10⁻¹⁶), and assay temperature (P<2.2×10⁻¹⁶). In both developmental and acclimation treatments, synthesis rates initially increased and then decreased with increasing assay temperatures in SD and BR (Table 2), and post hoc tests found no evidence for differences between the SD and BR copepods within each treatment×assay temperature combination (adjusted P≥0.11 for all). Yet, across all assay temperatures in both populations, 25°C development and adult acclimation resulted in higher ATP synthesis rates compared with those measured following development and acclimation at 20°C (adjusted P≤0.046 for all). This overall vertical shift in the thermal performance curve has the potential to mask variation between these treatments in the extent to which ATP synthesis rates are maintained during acute exposures to high temperatures. However, in SD copepods, rates of ATP synthesis first significantly declined with temperature between assay temperatures of 30 and 33°C for the 20°C treatment (adjusted P=7.3×10⁻⁴), whereas for the 25°C treatment the first decrease occurred between assay temperatures of 33 and 35°C (adjusted P=1.7×10⁻³). This suggests that copepods developed and acclimated as adults at warmer temperatures maintained synthesis rates at higher temperatures, at least in this population.

As an alternative approach to examine maintenance of ATP synthesis at high temperatures, we normalized synthesis rates across temperatures to those measured at 25°C for each mitochondrial isolation, which allows comparisons of the proportional changes in synthesis rate with assay temperature (Fig. 4). Note that this normalization could also reasonably be done to the rates measured at 20°C, but the results would be similar (Fig. 4). After normalization, proportional changes in rates of ATP synthesis were affected by a three-way interaction among population, temperature of development and adult acclimation, and assay temperature (P=0.02). Post hoc comparisons revealed similar patterns of variation among assay temperatures as those detected for the unnormalized rates (as would be expected), and when assayed at 37°C, 20°C-developed and -acclimated SD copepods maintained higher synthesis rates than 20°C-developed and -acclimated BR copepods (adjusted P=0.03). Additionally, at high assay temperatures, both populations maintained greater rates of ATP synthesis following development and adult acclimation at 25°C than following development and adult acclimation at 20°C (adjusted P≤5.5×10⁻³ for 33 to 37°C in SD and adjusted P=0.02 for 37°C in BR).

The mRNA expression levels of both heat shock proteins examined in the present study (hspb1 and hsp70) demonstrated similar effects of heat shock and developmental and adult acclimation temperature regardless of population (SD or BR; Fig. 5). In all cases, gene expression was affected by a significant interaction between the heat shock treatment, and the temperature of development and adult acclimation (P=0.04 for hsp70 in SD, and P≤5.4×10⁻³ for all others). In general, copepods developed and acclimated as adults at 25°C expressed higher levels of hspb1 and hsp70 than copepods developed and acclimated at 20°C (particularly after heat shock), although these patterns were not always resolved by post hoc tests (Fig. 5).

Given the potential role of mitochondrial performance in determining upper thermal tolerance (e.g. Harada et al., 2019) and a previous demonstration of decreased mitochondrial-encoded mRNA levels following heat shock in T. californicus (Schoville et al., 2012), we also examined variation in the expression of mt-atp6. In the present study, there were interactive effects of heat shock treatment and developmental and adult acclimation temperature on the mRNA expression of mt-atp6 in both the SD and BR populations (P≤0.02; Fig. 6). In SD copepods, regardless of the temperature of development and adult acclimation, mt-atp6 levels were similar in the control and 35°C heat shock treatments (P=1.00 for both). Within both of these treatments, expression levels were significantly higher in copepods developed and acclimated as adults at 25°C than in those developed and acclimated at 20°C (P<0.001 for both). In contrast, in SD copepods that had been developed and acclimated as adults at 20°C, heat shock at 36°C increased mt-atp6 expression (P≤0.02), whereas in those that had been developed and acclimated at 25°C, the same exposure decreased mt-atp6 expression (P<0.001 for both). As a result, there was no effect of the temperature experienced throughout development and adulthood on mt-atp6 mRNA levels in the 36°C heat shock treatment (P=0.88). In BR copepods, variation in mt-atp6 expression demonstrated somewhat different patterns than those observed for SD copepods. For both developmental and adult acclimation temperatures, there were trends for decreasing mt-atp6 mRNA levels with increasing heat shock temperatures, but these patterns were only resolved in post hoc tests in BR copepods developed and acclimated as adults at 25°C between the control treatment and the 35 and 36°C heat shock treatments (P<0.001 for both; P≥0.06 for all others). However, mt-atp6 expression was greater in BR copepods developed and acclimated as adults at 25°C than at 20°C in the control treatment and in the 35 or 36°C heat shock treatments (P≤0.02).

The results presented here provide experimental support for both of our proposed hypotheses. First, variation in developmental temperature resulted in differences in the plasticity of upper thermal limits in adult T. californicus. Regardless of population, 25°C-developed copepods demonstrated clear plasticity of CT_max between adulthood temperatures of 20 and 25°C. In contrast, there was no evidence of plasticity of this trait in adults that had developed at 20°C. Second, differences in developmental and adulthood acclimation temperatures were associated with plastic changes in two physiological mechanisms that are thought to contribute to the basis of local adaptation of upper thermal tolerance in this species. Furthermore, these effects were consistent with the differences in CT_max between these developmental and adult acclimation treatments. Therefore, our data suggest that adaptive processes may have the potential to shape the effects of developmental temperatures on the plasticity of thermal tolerance owing to shared underlying mechanisms, despite similar patterns of plasticity observed in the four locally adapted populations of T. californicus examined in the present study.

In general, dynamic and static thermal tolerance assays (i.e. gradual ramping exposures to high temperatures and abrupt exposures to a constant high temperature) resolve similar patterns of variation among experimental groups or treatments (e.g. Ford and Beitinger, 2005; Jørgensen et al., 2019), and our previous study suggested this was also the case among three Californian populations of T. californicus (distributed across ∼3° latitude; Harada et al., 2019). In the present study, there was a clear pattern of CT_max variation among populations that was consistent with substantial latitudinal thermal adaptation of upper thermal tolerance in T. californicus, as has been suggested in studies using static assays (Willett, 2010; Kelly et al., 2012; Pereira et al., 2014, 2017; Leong et al., 2018; Willett and Son, 2018; Foley et al., 2019). Overall, CT_max increased from northern to southern populations, although there was somewhat limited statistical resolution among populations, likely owing to a combination of relatively small differences in several comparisons and nonparametric post hoc tests. In general, thermal tolerance limits are expected to decline approximately linearly with latitude (Sunday et al., 2011), whereas our results suggest this is not the case in this species (Fig. S1C). This may reflect the relatively gradual latitudinal thermal gradient at higher compared with lower latitudes across the species range (particularly in the summer; Fig. S1B), and consistent with this possibility, both Leong et al. (2018) and Pereira et al. (2017) demonstrated approximately linear changes in upper thermal tolerance with differences in habitat air temperatures among T. californicus populations.

Despite this clear signature of local adaptation associated with the latitudinal thermal gradient along the west coast of California, interpreting the results of the present study in the context of T. californicus habitat temperatures is challenging. Splashpool temperatures vary substantially throughout the day (e.g. Harada and Burton, 2019), and the intertidal is a ‘mosaic’ habitat in which local thermal conditions may not necessarily reflect expected patterns with latitudinal variation in sea surface or air temperatures (Helmuth et al., 2002, 2006; Sanford and Kelly, 2011). For instance, variation in the daily timing of tidal cycles with latitude can paradoxically result in higher temperature exposures at northern compared to southern latitudes (Kuo and Sanford, 2009), although this is less likely to be relevant for supralittoral tidepools. There are limited published temperature data for T. californicus tidepools; however, at least comparing the SD and SC populations, overall summer temperatures tend to be warmer for the more southern population (i.e. SD; Leong et al., 2018), which is consistent with the difference in CT_max between these populations. Both average and maximum temperatures may be important for local adaptation, and it is likely that for upper thermal limits, maximum temperatures are more relevant (e.g. Somero, 2005). Yet, with the limited available data there does not appear to be a tight relationship between maximum temperatures and CT_max for SD and SC copepods (Leong et al., 2018). In part, this may be a consequence of CT_max only representing a proxy of the lethal thermal limit, which is typically justified as an ‘ecological death’ associated with an inability to escape predation or harmful conditions owing to loss of locomotor performance (e.g. Beitinger et al., 2000), and the extent to which this justification is relevant for T. californicus is unknown. That said, there is clearly local adaptation of CT_max in T. californicus, and comparing our results with those of Pereira et al. (2017) suggests that there is at least reasonable concordance between variation in CT_max and variation in lethal temperatures across populations.

Several other factors may influence comparisons of CT_max and splashpool temperatures, and additional temperature recordings will be necessary to examine this relationship in a comprehensive manner. For our data, these comparisons are likely dependent on the acclimation temperatures in the present study. Temperatures of 20 and 25°C are not atypical average weekly or monthly summer tidepool temperatures for the SD and SC populations (e.g. Leong et al., 2018), but it is unclear whether average conditions control field acclimatization (e.g. Fangue et al., 2011), and variable or cycling temperatures may alter acclimation responses compared with constant conditions (e.g. Paaijmans et al., 2013). Moreover, the thermal ramping rates in our CT_max trials are likely faster than natural rates of temperature increase in T. californicus tidepools in most cases, which may indicate that our CT_max values underestimate upper thermal tolerance under habitat conditions (Harada and Burton, 2019). Taken together, results in T. californicus to date consistently suggest that latitudinal temperature variation plays an influential role in inter-population variation in upper thermal tolerance, but the roles of local-scale differences in temperature and of habitat variability in determining upper thermal limits are yet to be fully resolved.

One distinct result of the present study was the lack of variation in CT_max between the sexes (see Fig. S1C). There is an overall consensus that, in comparison to males, female T. californicus are more tolerant of stressful conditions for a wide range of abiotic factors, including temperature (Willett, 2010; Willett and Son, 2018; Foley et al., 2019). Insufficient statistical power due to nonparametric tests, and relatively low sex-specific sample sizes (n=8), may explain the lack of sex effects in the present study. However, two other studies have also failed to detect differences in upper thermal tolerance between the sexes (Pereira et al., 2014, 2017). Regardless, CT_max was not statistically affected by sex in our study, and as a result we did not consider variation between females and males further here.

Across the SD, BR, SC and PE populations of T. californicus, we consistently observed variation in CT_max plasticity in adults as a result of temperatures experienced during development. After development at 25°C, CT_max was higher in copepods acclimated to 25°C in adulthood than in copepods acclimated to 20°C, whereas 20°C-developed copepods displayed no difference in upper thermal limits between the adult acclimation treatments. These patterns could be the result of differences in reversible adult plasticity owing to developmental temperatures, or of temperature-dependent reversibility of developmental plasticity, but in either case this phenotypic variation is consistent with an interactive effect between developmental and adulthood temperatures. To our knowledge, this is the first demonstration of interactive effects between temperatures in development and in later stages of life on thermal limit plasticity in a marine ectotherm. Alternatively, these patterns could be potential consequences of differences in developmental survival between 20 and 25°C as we did not directly monitor survivorship in this study; however, Pereira et al. (2017) and Harada et al. (2019) observed little evidence of differential survival at these temperatures across most T. californicus populations. Furthermore, previous studies have also detected interactive effects of developmental and adulthood temperatures on upper thermal tolerance in zebrafish (Danio rerio; Schaefer and Ryan, 2006) and fruit flies (Drosophila spp.; Kellermann et al., 2017; Kellermann and Sgrò, 2018), which, in combination with the results of the present study, suggest there is mounting evidence that these effects may be common for this trait.

In T. californicus, the effects of developmental temperature on adulthood plasticity of CT_max were relatively large, as plasticity of CT_max was completely absent in adults after development at 20°C, whereas after development at 25°C adulthood plasticity was observed in all populations. Moreover, the adulthood acclimation response ratio in 25°C-developed copepods (i.e. ΔCT_max °C⁻¹) was typical for aquatic ectotherms (∼0.2; Gunderson and Stillman, 2015). Similarly, variation in developmental temperature results in presence–absence differences in adulthood plasticity in D. melanogaster (Kellermann et al., 2017), although in Drosophila spp., cooler developmental temperatures tend to increase plasticity in adults (Kellermann et al., 2017; Kellermann and Sgrò, 2018), whereas our results suggest warmer developmental temperatures increase adult plasticity in T. californicus. The loss of adult CT_max plasticity at an only moderately reduced developmental temperature in T. californicus is potentially surprising given the prevalence of at least a modest capacity for acclimation of this trait across many species (e.g. Gunderson and Stillman, 2015). Acclimation to constant conditions, rather than cycling thermal regimes that more closely resemble natural tidepool conditions, may influence both this lack of plasticity and the effects of developmental temperature on adulthood plasticity in the present study. In addition, we examined only a relatively small range of adulthood temperatures (20–25°C), which may contribute to the lack of observed plasticity, and the extent to which plasticity may alter CT_max in 20°C-developed copepods over a wider range of adult acclimation temperatures remains an open question.

The short generation times of T. californicus and Drosophila spp. likely increase the concordance between developmental and adulthood temperatures (particularly in the habitats of T. californicus). This may contribute to the effects of developmental temperature on adulthood plasticity of thermal tolerance in this species, because developmental effects are expected to be stronger if conditions in development are predictive of those experienced as adults (Cooper et al., 2010, 2012a,b; Nettle and Bateson, 2015; Beaman et al., 2016). Consistent with this possibility, in the comparatively long-lived zebrafish, Schaefer and Ryan (2006) observed only subtle shifts in CT_max plasticity in adults as a result of differences in developmental temperatures. Although the interactive effect of developmental and adulthood temperatures with respect to patterns of CT_max plasticity was relatively strong in T. californicus, the maximum difference in CT_max among treatments was approximately 1°C. Regardless, our data demonstrate that variation in developmental temperatures can have substantial effects on the adulthood plasticity of upper thermal tolerance in aquatic ectotherms.

The possibility of interactions between adaptive processes and phenotypic plasticity is well established (e.g. Crispo, 2007; Hendry, 2016; Donelson et al., 2019; Kelly, 2019), and thus there is also the potential for local thermal adaptation to shape effects of developmental temperature on the plasticity of upper thermal tolerance in adults. Furthermore, if heat hardening is used as a metric of adulthood plasticity, there is some evidence for adaptive variation in these effects among Drosophila spp. (Kellermann and Sgrò, 2018). However, when this possibility was assessed with adult acclimations in temperate and tropical D. melanogaster, developmental effects on adulthood plasticity were similar among populations (Kellermann et al., 2017). Similarly, although development at 20 or 25°C was associated with variation in CT_max plasticity in adults from all of the T. californicus populations examined in the present study, there was no variation in this effect among populations. Indeed, in 25°C-developed copepods, the average CT_max difference between adult acclimations of 20 and 25°C were remarkably similar across the four populations (SD: 1.0°C, BR: 0.8°C, SC: 1.3°C and PE: 0.9°C). Thus, our results suggest that effects of developmental temperature on adulthood plasticity of CT_max have not been altered substantially by local thermal adaptation among these populations of T. californicus. However, we also found that physiological mechanisms that putatively underlie latitudinal variation in upper thermal tolerance in this species (e.g. Schoville et al., 2012; Harada et al., 2019) show patterns of variation in response to developmental and adulthood temperatures that parallel variation in CT_max.

Harada et al. (2019) demonstrated that during acute exposures to elevated temperatures, the temperatures at which maximal ATP synthesis rates first decline are correlated with CT_max across the SD, AB and SC populations of T. californicus. Consistent with this relationship, several studies have demonstrated loss of ATP synthesis capacity in heart mitochondria at temperatures that are approximately equal to or immediately below the upper thermal limits in species of fish (Iftikar and Hickey, 2013; Christen et al., 2018; O'Brien et al., 2018). Temperature-mediated plasticity of mitochondrial functions is also often observed in ectothermic species (e.g. Guderley, 2004; Seebacher et al., 2010; Chung and Schulte, 2015; Chung et al., 2017a,b, 2018; Bryant et al., 2018), and our data suggest that this is the case in T. californicus. In the SD and BR populations, copepods that were developed and acclimated as adults at 25°C had greater ATP synthesis rates than those that were developed and acclimated at 20°C, which was consistent with higher expression levels of mt-atp6 in these 25°C treatments than these 20°C treatments under control (i.e. non-heat shocked) conditions. Additionally, developmental and adult acclimation temperatures of 25°C compared with 20°C resulted in greater maintenance of ATP synthesis rates at high temperatures, which is consistent with difference in CT_max between these treatments. In comparison to Harada et al. (2019), the thermal performance curves observed in our study were horizontally shifted to moderately lower temperatures, and were remarkably flat with maximum Q₁₀ values of approximately 1.4–1.5 across treatments. Although relatively thermally insensitive physiological rates have been observed previously in T. californicus (e.g. Scheffler et al., 2019), there is clearly an unknown source of variation in these ATP synthesis curves among studies. It is possible that culturing conditions could contribute to this variation, as Harada et al. (2019) examined copepods taken directly from stock cultures (i.e. 400 ml beakers), whereas in the present study we raised groups of copepods specifically for these measurements in 10 cm Petri dishes. Associated with this difference, these studies likely also differ somewhat in densities of copepods, algal growth and, potentially, oxygen levels under holding conditions, which have the potential to result in plastic variation in mitochondrial performance between the studies. Regardless, with the exception of temperature, our 20 and 25°C treatments were held under equivalent conditions, and therefore the difference in high-temperature maintenance of ATP synthesis rates between these treatments is likely robust to any variation in thermal performance curve estimates across studies.

Schoville et al. (2012) examined genetically determined differences in the transcriptomic response to acute heat stress between the SD and SC populations of T. californicus. Both the strongest response and largest difference between the two populations was the extent to which heat shock protein genes were induced following heat stress. Particularly for hspb1 and hsp70, heat shock protein mRNA expression was increased to much higher levels in the warm-adapted SD population than in the relatively cold-adapted SC population. As heat shock proteins are molecular chaperones that mitigate the negative effects of high temperature associated with damaged and unfolded proteins (Hochachka and Somero, 2002), these transcriptomic patterns suggest that differences in heat shock protein expression may contribute to the difference in upper thermal tolerance between the SD and SC populations. The evidence for a correlation between large inductions of heat shock protein expression and increased tolerance of high temperatures is somewhat mixed among genes and species (e.g. Healy et al., 2010; Gleason and Burton, 2015); however, studies in fruit flies (D. melanogaster) and marine snails (Chlorostoma funebralis) generally support a positive relationship between these two traits (e.g. Bettencourt et al., 2008; Tomanek, 2008). In T. californicus, Tangwancharoen et al. (2018) demonstrated putatively adaptive functional variation associated with sequence differences among populations in both the regulatory and coding regions of hspb1, and Barreto et al. (2015) utilized RNA interference to show that knockdown of transcripts for this gene directly decreases survivorship following acute thermal stress. Therefore, the increased inductions of hspb1 and hsp70 we observed in SD and BR copepods developed and acclimated as adults at 25°C compared with those developed and acclimated at 20°C are likely beneficial effects of plasticity, and are consistent with CT_max differences between these treatments. Part of this variation in heat shock protein expression may be associated with the differences in recovery temperatures in our study, as each developmental and adult acclimation treatment was recovered at its acclimation temperature. However, in all cases, the fold differences in expression between the 20 and 25°C developmental and adulthood treatments following heat shock (2.0–6.4) are greater than would be expected as a result of thermodynamic effects on transcription rates alone (1.4–1.7 given an expected Q₁₀ of 2–3, which is consistent with thermal sensitivities of transcriptional elongation rates; e.g. van Breukelen and Martin, 2002).

Taken together, our results indicate that both the extents to which ATP synthesis rates are maintained and heat shock proteins are induced at high temperatures are elevated in T. californicus that are developed and acclimated at 25°C compared with those that are developed and acclimated at 20°C. The acute temperature exposures used to assess these mechanisms here matched those of previous studies in this species (Schoville et al., 2012; Barreto et al., 2015; Harada et al., 2019), but these exposures were notably different than the thermal ramp experienced during our CT_max trials, which may affect comparisons among these traits (Harada and Burton, 2019). Despite this, the patterns of variation in ATP synthesis rates and heat shock protein expression observed here were consistent with differences in CT_max between the copepods that were developed and acclimated as adults at 25°C and the copepods that were developed and acclimated at 20°C. This suggests that maintenance of ATP synthesis rates and induction of heat shock proteins likely contribute to the basis for plasticity of upper thermal tolerance associated with developmental and adulthood temperatures in this species. Yet, the extent to which the effects of developmental temperature, specifically, on plasticity in adults can be attributed to these mechanisms requires additional research, as these traits were not assayed in 25°C-developed copepods that were transferred to 20°C in adulthood in the present study. However, our data suggest these mechanisms may play a role in plastic effects associated with developmental temperatures in general.

The effects of environmental change on organisms ultimately depend on the capacities to modulate key physiological traits to facilitate performance and persistence. Phenotypic plasticity, adaptation and interactions between these two processes all play important roles in these responses (e.g. Kellermann and van Heerwaarden, 2019). Here, we show that temperatures experienced in development also shape the adulthood plasticity of upper thermal limits in the intertidal copepod T. californicus. These effects may be particularly relevant for aquatic ectotherms as thermal tolerance limits likely underlie geographic range limits in many of these species (Sunday et al., 2012; Pinsky et al., 2019). Our results highlight that beneficial effects of developmental plasticity with respect to environmental change have the potential to be overestimated if considered without accounting for temperature variation in adulthood. Additionally, the data presented here suggest that the physiological mechanisms that may underlie these effects (e.g. shifts in the thermal performance curve for ATP synthesis and the regulation of heat shock genes) are, at least to some extent, shared with the mechanisms associated with local thermal adaptation in T. californicus. This mechanistic overlap indicates the potential for interactions among local adaptation and plasticity at difference life stages to shape variation in upper thermal tolerance in ectothermic organisms.