Introduction

Global change already has a major impact on marine ecosystems through changes in environmental variables such as water temperature, salinity, pH, and oxygen content (Wiltshire et al. 2010; Doney et al. 2012; Poloczanska et al. 2013; Boersma et al. 2016; Molinos et al. 2016; Boyd et al. 2018). However, patterns of change vary spatially depending on regional scale habitat characteristics. For instance, how anthropogenic or other types of climate changes drive environmental conditions will differ between open ocean, estuarine, and nearshore habitats; coastal areas in particular are likely to experience comparatively drastic changes in environmental conditions (reviews Hofmann et al. 2011; Duarte et al. 2013). These areas are increasingly affected by rising surface temperatures and changes in river runoff (Gräwe et al. 2013; Hiddink et al. 2015; Robins et al. 2016; Reusch et al. 2018).

As in most cases with two or more climate-driven environmental variables (= drivers: Crain et al. 2008; Piggott et al. 2015; Orr et al. 2020; Tekin et al. 2020), the combination of lower salinity in a scenario of increased temperature is likely to elicit complex interactive responses. The interplay between drivers can evoke many non-independent effects, where one driver can enhance (synergistic effect) or weaken (antagonistic effect) the impact of the other driver beyond the additive effect expected from each driver in isolation (e.g. see reviews Orr et al. 2020; Tekin et al. 2020). In coastal species, tolerance to lower salinity can be influenced by temperature, leading to either synergistic or antagonistic effects depending on the species (e.g. Janas and Spicer 2008; Nasrolahi et al. 2012; Spitzner et al. 2019; Torres et al. 2021a).

A critical point is that survival and physiological responses to temperature and salinity may vary among populations, especially for those species with a wide spatial distribution. In situations where the magnitude of variation is important, single population studies may not fully characterise the responses of a given species. There is an increasing body of work showing that responses to climate drivers vary considerably at the intraspecific level (Carter et al. 2013; Durrant et al. 2013; Applebaum et al. 2014; Spitzner et al. 2019). Intraspecific variation in responses to suboptimal conditions arise from phenotypic plasticity and evolutionary adaptation (Chevin et al. 2010; Reusch 2014; Sinclair et al. 2013; Boyd et al. 2018), and local genetic adaptation is likely to be expressed as variation among populations in the responses (De Villemereuil et al. 2016). Irrespective of the mechanism driving phenotypic variation, recent studies have further pointed out the necessity to incorporate intraspecific trait variation into models of community structure in response to climate change (Bolnick et al. 2011; Violle et al. 2014; Moran et al. 2016). Intraspecific variation in responses can drive portfolio effects (Schindler et al. 2015; Price et al. 2021), whereby diversity in the types of responses buffers a species from environmental variation. The type of effect will depend on the species. For example, in species with dispersive stages, when such variation occurs among populations, a spatial portfolio effect should occur if dispersion ensures sufficient connectivity among populations, because populations with more resilient genotypes would subsidise other populations (Casaucao et al. 2021). In species of commercial interest (such as salmon), portfolio effects can ensure food security, but climate change, and anthropogenic activities can lead to portfolio simplifications (Price et al. 2021). On the other hand, in the case of invasive species, the diversity of responses underpinning the portfolio effect should theoretically enhance invasion potential across environmental gradients; the portfolio effect might be manifested as sustained invasion in spite of environmental variation. Along this line, quantifying intraspecific variations, but also the current portfolio, is useful to identify populations with higher invasion potential or more resilient or sensitive to climate change.

Within each population, developmental, growth, and other physiological rates must be maintained in a scenario of increased temperatures, either through interpopulation connectivity or local adaptive mechanisms. These rates are key to determine fitness-related traits such as body size, and the amount of nutritional reserves needed to survive through development under unfavourable conditions, e.g. food limitation. Salinity and temperature are known to govern physiological rates and survival in many marine organisms (Pörtner 2001; Anger 2003; Somero 2005; Sokolowa et al. 2012). Understanding how such rates are driven by temperature and salinity is important for the development of mechanistic physiological based models predicting species’ responses. In most metazoans, the physiological processes associated with movement, feeding, growth, and development rely on oxygen consumption (Hochachka and Somero 1984; Burggren and Roberts 1991). Physiological rates can also be viewed as a series of linked traits governing body mass at critical life history events, such as metamorphosis (Werner 1988; Hentschel 1999; Gotthard 2001; Heyland et al. 2011; Torres and Giménez 2020). Body size, or reserves at metamorphosis, are indicators of post-metamorphic fitness because they explain variations in post-metamorphic survival (Jarrett 2003; Giménez et al. 2004; Pechenik 2006; Tremblay et al. 2007; Torres et al. 2016). Variations in body mass or reserves at metamorphosis may not only reflect passive stress responses to a given condition, but instead constitute a part of an adaptive compensatory mechanism (Gotthard and Nylin 1995). For instance, increased developmental time in response to suboptimal food conditions may constitute a compensatory response to maintain body mass at metamorphosis despite reduced growth rates. In this example, extended development may minimise the fitness costs associated with metamorphosing with a reduced body size, at expenses of (presumably lower) costs associated with delaying metamorphosis. We know little about such multiple trait responses to salinity and temperature in early life history stages, when growth and development are critical. In ectotherms, increased temperatures accelerate developmental rates, but the effects on growth may vary. Moreover, increased temperature in organisms exposed to optimal salinity results in increased respiration rates (Ikeda 1970; Anger 2001; Saborowski et al. 2002) up to a point where oxygen demands cannot be matched by appropriate increase in oxygen supply (Pörtner 2001, 2012; González-Ortegón et al. 2012). At such levels, the metabolic demands for oxygen exceed the organism’s ability to distribute oxygen to tissues, and organismal performance declines (Pörtner and Knust 2007; Giomi and Pörtner 2013; Deutsch et al. 2015).

We also know that responses to decreased salinity can vary among species depending on the degree of euryhalinity, with increases in respiration rates observed in euryhaline species and decreases in stenohaline species, or at salinities beyond the tolerance range for euryhaline species (Kinne 1971; Anger 2001, 2003). Osmoregulation involves ion regulation (performed by the Na+-K+ ATPase) at expenses of ATP; thus, high energy demands (and hence respiration/oxygen consumption rates) are expected in euryhaline species when exposed to extreme salinities. In stenohaline or weak euryhaline species, low salinity depresses metabolism (Johns 1982; Yagi et al. 1990; Anger et al. 1998), most likely through perturbations of cellular level processes and eventually damage. In such species, low salinity may reduce growth rates (Torres et al. 2011) and slow development (Anger 1991, Nurdiani and Zeng 2007; Torres et al. 2008, 2020, Šargač et al. 2021, see Anger 2003 for review of earlier literature). We know indeed much less about how developmental, growth and oxygen consumption rates in larval stages respond to the combined action of lower salinity and increased temperature. This multiple driver scenario may reduce the critical temperature at which metabolic depression occurs (Yagi et al. 1990) or may depress metabolic rates at all temperatures. Overall, the interactive nature of multiple driver effects and the fact that oxygen consumption rates may both increase or decrease in response to suboptimal conditions would diminish our trust in extrapolation based on current information.

Here, we quantify the effects of reduced salinity and increased temperature on the oxygen consumption, developmental, and growth rates of larvae of the shore crab Carcinus maenas obtained from three European populations (Cádiz: Cádiz Bay, Helgoland: German Bight, Kerteminde: Baltic Sea), located along its native range (Fig. 1a; Leignel et al. 2014; Young and Elliott 2020). C. maenas is native to Europe but it is considered a global invader (Roman and Palumbi 2004; Young and Elliott 2020) and is listed in the top 100 worst invaders (Lowe et al. 2000). The life history of C. maenas is characterised by a complex life cycle (Spitzner et al. 2018); females are benthic (= bottom dwelling), carry eggs and release larvae in the water column. Larvae drift in the sea and develop through four pelagic stages (zoeae), followed by a settling stage (megalopa) that colonises the parental habitat (sheltered intertidal zones). We were interested in determining how responses of different performance rates shape the integrated phenotype in terms of body mass, carbon and nitrogen content. Additionally, we investigated whether such responses were consistent among populations located as far away as the Cádiz Bay and the German Bight or the Baltic Sea, as a way to test if a general response pattern emerges. The population of Cádiz, near the equatorial side of the species’ home distribution range, is exposed to higher temperatures than those of the North and Baltic Seas. The population of the North Sea is located towards the centre of the latitudinal distribution which includes also the Norwegian coast. Those differences may contribute to the among-population variation in larval responses to temperature and salinity. This study is based on (1) unpublished information on oxygen consumption rates for the three populations, (2) unpublished information on survival, growth and developmental rates from a population from Cádiz, and (3) published information (Šargač et al. 2021) of similar rates from the Helgoland and Kerteminde populations (presented as grey sections in the figures) for comparison. We know that larvae of C. maenas can exhibit an antagonistic response by which the negative effect of lower salinity (range 20–32 PSU) on survival and growth is mitigated at higher temperatures (range 15–24 °C). This response has been called “Thermal Mitigation of Low Salinity stress” (abbreviated as TMLS, Spitzner et al. 2019) and has been found in other coastal crustaceans (e.g. Janas and Spicer 2008; Nasrolahi et al. 2012). We also know that for survival and development, the magnitude of TMLS can vary within and among populations (Spitzner et al. 2019; Torres et al. 2020; Šargač et al. 2021). However, we do not know how larvae respond in local populations experiencing higher temperatures such as the one of Cádiz Bay.

Fig. 1
figure 1

Carcinus maenas. Experimental setup and design. a Native range and collection sites: berried females with eggs were collected from three different populations (black stars) across the native European range (blue area). Symbols (♀) correspond to the number of females that successfully produced larvae from each population. b Larval rearing: freshly hatched larvae from each population and female were reared in a factorial design consisting of eight combinations of temperature, T (15, 18, 21, and 24 °C) and salinity, S (25 and 32.5 PSU) in a common garden setup. Each combination of population × female × T × S had five replicate bowls (60 mL) with ten randomly assigned larvae in each. c Respiration measurements design: once the larvae moulted to the megalopa stage, they were randomly taken into 4 mL vials and oxygen consumption rates were measured. Six megalopae were taken from each combination of population × female × T × S, fasted for 2 h and oxygen consumption was measured for each megalopa individually. Additionally, six replicate vials without animals were placed as a control

Carcinus maenas is representative of a coastal species with important intraspecific variations in the responses to lower salinity and higher temperature at the regional scale. Because pelagic dispersal is key for population connectivity (Armsworth 2002; Cowen and Sponaugle 2009), larval physiology is crucial for a potential portfolio effect by recovering local populations from mass mortality (Giménez et al. 2020). More in general, marine larvae are particularly vulnerable to environmental changes (Pandori and Sorte 2019) and larval physiological quality can drive recruitment and connectivity (Tremblay et al. 2007; Shima and Swearer 2009). Larval supply and settlement is central for long-term persistence of marine benthic populations (Morgan 2020).

The question is whether populations located at the southern distribution thermal edge, such as Cádiz, show similar responses to those experiencing coastal marine waters further north (Helgoland: Spitzner et al. 2019; Irish Sea: Torres et al. 2020) or brackish environments (Kerteminde: Šargač et al. 2021). It is also unknown whether oxygen consumption rates respond to salinity and temperature consistently among populations, regardless of distance and habitat type. We quantified oxygen consumption rates at the megalopa stage for which we have previous information on the effect of temperature on such rates (Dawirs 1983), but not about the combined effect of temperature and salinity nor about whether responses vary among populations. One can derive several hypothetical scenarios in responses of oxygen consumption rates with that of survival and other physiological rates. One may find consistent responses across all rates such that among population differences are manifested as a global type of response. Alternatively, differences between populations in the rates may not match, which would indicate differences in the sensitivity to the drivers across different physiological processes.

Materials and methods

Female collection, transport, and maintenance

Carcinus maenas berried females (with early embryos, carapace length = 2.55–5.10 cm) were collected during their reproductive period in the same year from three European locations (Fig. 1a): Helgoland (North Sea, Germany, coordinates: 54° 10′ 49.176′′ N, 7° 53′ 20.198′′ E, salinity = 30–33 PSU); Cádiz (Gulf of Cádiz, Spain, coordinates: 36° 33′ 34.5′′ N 6° 12′ 23.1′′ W; salinity = 34–35 PSU) and Kerteminde (Baltic Sea, Denmark, coordinates: 55° 26′ 59.9′′ N, 10° 39′ 40.1′′ E; salinity = 16–20 PSU). Females from the Cádiz population were hand-collected in February–March from the intertidal muddy shore of the San Pedro River next to the Institute of Marine Science of Andalusia (Cádiz, Spain). Females from Kerteminde were caught in April–May using traps (Torres et al. 2021b) in the subtidal fjord next to the Marine Biological Research Centre, University of Southern Denmark (Kerteminde, Denmark). Prior to transport to the Marine Biological Station at the Alfred-Wegener-Institute (Helgoland, Germany), berried females were kept in aquaria with a flow through system of natural water with gentle aeration reflecting natural salinity and temperature conditions on the site (35–37 PSU and 17 ± 3 °C in Cádiz; 16–20 PSU and 10 ± 1 °C in Kerteminde). In preparation for transport, females were placed individually in 1 L plastic containers, half-filled with strongly aerated water from the site of origin, and a wet towel. All containers were transported in a Coleman® cooler box to ensure a constant temperature of the natural local conditions during transport (ca. 24–36 h). Females from the Helgoland population were collected manually during April–June from the intertidal rocky shore and directly transported to the Marine Biological Station. For those females, we simulated (for 24–36 h) the transport conditions experienced by the animals from the other two populations. After arrival to the laboratory (or the simulation of transport stress), females from each population were placed individually in 5 L aquaria filled with UV-treated filtered natural seawater from Helgoland (32.5 ± 1 PSU) and gentle aeration. Animals were placed at the same temperature as the time of collection which corresponds to their natural habitat of origin (i.e., 10, 15, and 16 °C for Kerteminde, Cádiz, and Helgoland, respectively), then temperature was gradually increased (0.5 °C per day) until the same common embryonic temperature (18 °C) was reached for all populations. Females were fed twice a week (frozen shrimp Crangon crangon) and water was changed daily.

Experimental design

Freshly hatched larvae originated from each female were randomly assigned to eight combinations of four temperatures (15, 18, 21, and 24 °C), and two salinities (25 and 32.5 PSU) in a factorial design (Fig. 1b). Larvae were reared in groups of ten individuals in five replicate bowls (60 mL) per treatment combination (Fig. 1b). This design was repeated with larvae of each different female separately. A total of five females from each population were used in the experiments, except for those originated from Kerteminde, where only larvae from three females were available to be included in the experiments (Fig. 1a). Note that for the Kerteminde population, measurements at the lowest temperature (15 °C) and salinity (25 PSU) were not possible due to a high sensitivity of early stages to these stressful conditions resulting in failure to develop to the megalopa. Overall, 5200 larvae were cultured during this experiment (i.e., 10 larvae × 5 replicates × 2 salinities × 4 temperatures × 13 females of origin). The design contains two separate three-way interaction terms (temperature × salinity × female and temperature × salinity × population) with “female” being nested within population.

Experiments were conducted in temperature-controlled rooms (± 0.5 °C) with a fixed 12–12 h light–dark cycle. Aerated natural seawater from Helgoland, previously UV-treated and filtered (filter pore size = 0.2 µm), was used in all the experiments. Diluted seawater was prepared by mixing natural seawater with appropriate amounts of tap water to reach the appropriate salinity to the nearest ± 0.1 PSU (salinometer: WTW Cond 3110 SET1) and stored in 60 L tanks at each of the experimental temperatures. Salinity was monitored, daily in the storage tanks, and before and after water changes in the experimental bowls to ensure that it remained constant during the experiments. During the daily water change, larvae were fed ad libitum with daily freshly hatched Artemia sp. (Great Salt Lake Artemia, Sanders, USA); moults and dead larvae were reported and discarded.

Oxygen consumption measurements

All oxygen consumption measurements were performed at the same temperatures and salinities used for larval rearing: the same four experimental temperatures in temperature-controlled rooms (15, 18, 21 and 24 °C) and the same salinities (25–32.5 PSU). Prior to the oxygen consumption measurements, recently moulted megalopae (< 24 h after moult) from each treatment were starved for 2 h to avoid confounding effects of excretion on the oxygen consumption measurements. Oxygen consumption measurements were conducted in 4 mL vials integrated with an optically isolated oxygen sensor type PSt5 at its bottom (PreSens, Regensburg, Germany). Each experimental vial was first filled with O2-saturated water (i.e., water was bubbled until the moment of filling the vials for oxygen consumption measurements) at the rearing temperature (15, 18, 21 or 24 °C). Immediately after, one megalopa per replicate (six replicates for each treatment) was placed in the measurement vial, and subsequently the vial was closed underwater to avoid air bubbles (Fig. 1c). The controls for the oxygen consumption measurements consisted of six additional vials at each salinity and temperature without animals. Vials with samples and control vials were placed in a 24-well plate and placed on a SDR SensorDish® Reader (PreSens, Regensburg, Germany). This system consists of a 24-channel reader of oxygen luminescence quenching and function as the incubation chamber, providing a high-quality measurement without oxygen consumption or gas exchange between the environment and the vial (Warkentin et al. 2007). The PreSens system was calibrated with a two-point calibration (0 and 100% oxygen saturation) before the experiment for each water condition (temperature and salinity) according to the manufacturer’s protocol. The 24-channel reader and vials were placed on a rocking platform shaker (IKA Rocker 2D digital, Staufen, Germany) at 40 revolutions per minute (rpm) to avoid oxygen stratification within the vials. As primary producers can increase the oxygen content in the vials, an opaque black plastic box was placed to cover the vials to block the light during measurement. Oxygen consumption was measured for ca. 4 h, and a longer period in the colder temperatures to compensate for slower oxygen consumption rates. Readings were recorded every 15 s as concentration of O2 (mg L−1). In cases when unfavourable conditions had led to unsuccessful development to the megalopa, measurements were made only when there were at least three replicates. After completion of measurements, animals were inspected to confirm they were alive, and sampled for the biomass quantification.

Biomass, survival and development

In order to explain the patterns observed during oxygen consumption measurements and its potential consequences, we included already published data for survival, duration of development, and growth of the Kerteminde and Helgoland populations (Šargač et al. 2021). In this study, we added additional analyses including the new data of the Cádiz population together with the already published data set of the other two populations. Note that already published data is shaded in grey in figures (for significant differences between Kerteminde and Helgoland see Šargač et al. 2021), for easier comparison between the populations.

After the oxygen consumption measurements were concluded, megalopae were sampled individually for the determination of biomass (i.e., dry mass and elemental composition) following Šargač et al. (2021). Megalopae were transferred directly from the vials, gently rinsed 10 s with distilled water, shortly blotted dry and placed into tin cartridges and stored at − 20 °C for further analyses. Afterwards, samples were freeze-dried for 48 h (Christ Alpha 1–4 freeze drier) and dry mass was determined using a microbalance (Sartorius SC2, to the nearest 0.1 µg). Carbon and nitrogen contents were determined with an Elemental analyser (vario MICRO cube CHNS analyser, Elementar Analysensysteme). Cumulative survival was calculated as proportion of larvae surviving from hatching to the megalopae according to Šargač et al. (2021). In order to avoid situations of log(0) values, we rescaled the proportions (p) with the formula p′ = [p (N−1) + 0.5]/N, where N is the number of larvae per glass (= 10 larvae) at the start of the experiment. Developmental time was determined as the time required from hatching to reach the megalopa. Analyses for biomass, survival, and duration of development were performed for all three populations and an additional analysis of the Cádiz population alone.

Data analysis

The experimental design consisted of three fixed factors (population: P; temperature: T, and salinity: S) and female of origin (F) as a random factor nested in population. Mixed model analyses were first conducted through backward model selection (Zuur et al. 2009) using the nlme package (Pinheiro et al. 2019). For mixed models, the function lme was used, and for fixed models gls. Model selection was performed in two steps: we first compared models for random terms after restricted maximum likelihood (REML), then after refitting the model with the best random structure, we checked the fixed terms with maximum likelihood (ML). Corrected Akaike information criteria (AICc) and likelihood ratio tests (LRT) were used to compare models and the Tukey HSD test was used for post-hoc comparisons. The model with the lowest AICc was always selected for further analysis. In cases when the difference between models was ΔAICc < 3 and the more complex model had lower AICc, we used LRT to test which model was more suitable: (1) when models differed significantly (p < 0.05), the model with lower AICc was chosen; (2) when there was no significant difference between the models, we chose the model with the lowest number of parameters.

The oxygen consumption rate in each vial was calculated by linear regression, using the part of the curve where oxygen concentration vs. time showed a linear trend. Individual oxygen consumption rates (VO2) were calculated for each population as the difference in oxygen consumed per hour between control vials (without animals) and vials containing megalopae. Oxygen consumption rates were first calculated as averages of replicates for each treatment, female of origin and population. Then, we determined the average for all females in each treatment and population. Oxygen consumption rates were presented as µg of oxygen consumed per individual per hour (µg O2 ind−1 h−1). Dry mass (DM), carbon (C) and nitrogen (N) content were calculated for post-moult megalopae, and to calculate weight-specific oxygen consumption rates (QO2) we divided oxygen consumption rates by the dry mass (µg) of each individual animal that was used in the measurement. Additional statistical analyses, such as means, standard deviations and standard errors, were computed in R (version 4.0.2) through RStudio via “plyr” package.

The full mixed model for oxygen consumption rates (Table 1) had three factorial fixed components: population, salinity and temperature (~ P × S × T). The replicate units (n = 6) were the vials used to measure oxygen consumption. Replicates for oxygen consumption did not fully coincide with the replicates for larval rearing because we had to ensure that oxygen consumption measures were taken with larvae with the same moulting age. Because some larvae inevitably die, and due to limitations in the number of parallel measurements to be taken per day, it was not possible to assign larvae from different rearing glasses to separate oxygen consumption vials. Importantly, because of the quick developmental changes occurring within the moult cycle (Anger 2001), confounding due to larval age is more relevant than potential confounding associated to the replicate glasses; hence, we focused on minimising the former. In order to reduce the latter source, we focused on variation in the responses among larvae from different females and kept “female of origin” (F) in the random term of all selected models. Female of origin was entered in the full model in interaction with fixed factors with terms coded as “random = list (F = pdDiag(~ S × T))”; we also checked that the full covariance matrix (pdSymm) did not result in model convergence. For more simple models, the random component was also tested with the full as well as the diagonal matrices, e.g. as “random = list (F = pdDiag(~ S))” or “random = list (F = pdDiag(~ T))”. The simplest model (for the random term) was coded as “random =  ~ 1|F”. Due to insufficient moulted megalopae from the Kerteminde population at 15 °C and 25 PSU, the following two analyses were performed: (1) for all three populations between all the treatments except 15 °C; (2) additional analysis between Cádiz and Helgoland containing 15 °C (Table 1).

Table 1 Carcinus maenas. AICc values for the model selection on random terms were carried out through restricted maximum likelihood (REML) fitting

For biomass and elemental composition (DM, C and N), and duration of development and instant growth rates of the megalopae, a mixed model was performed (Table S2) as explained for oxygen consumption rates. We compared populations with a three-way fixed interaction (~ P × S × T). First, we analysed all three populations excluding the groups exposed to 15 °C, and then we compared Cádiz and Helgoland with all eight treatments. Female of origin was kept as a random factor throughout all the steps coded in a model as “random =  ~ 1|F”. Survival analyses were performed for all eight treatments among the three populations (Table S1). Additional ANOVAs were conducted for the Cádiz population for each parameter (survival, biomass, duration of development, and growth) to check the significant interactions between the factors (salinity × temperature), for later comparison between the three populations.

Results

We first report the rates of oxygen consumption for the three populations. Then, we describe the responses for the Cádiz population for survival, development, biomass, and growth, and compare it with those of Helgoland and Kerteminde, already published in Šargač et al. (2021).

Oxygen consumption rates

The best statistical model for oxygen consumption rates per individual (VO2) did not retain any interaction between the tested terms (Table 1), while oxygen consumption rates increased with temperature (T). Temperature was the main driver for oxygen consumption rates per individual, where a higher temperature generally led to higher rates in all three populations (Fig. 2). The population effect was only retained in the model comparing Helgoland and Cádiz, whereby the latter showed on average slightly higher VO2 values.

Fig. 2
figure 2

Carcinus maenas. Oxygen consumption rates per individual (VO2) for megalopae. Comparison between three populations measured at two larval salinities (25 PSU, blue and 32.5 PSU, green) and four temperatures (15, 18, 21 and 24 °C). Symbols represent each combination of factors per population. Data shown as average individual oxygen consumption rates (µg O2 ind−1 h−1) ± SE among larvae produced by different females (n = 5 or 3)

Weight-specific oxygen consumption rates (QO2) increased with temperature for all populations and reached higher values for the Cádiz population as compared to Helgoland and Kerteminde (Fig. 3). The best statistical model retained the interaction between temperature and population (Table 1) when the three populations were compared (without 15 °C treatments), and as main factors for the comparison of Cádiz and Helgoland (with all treatments). Salinity was not retained either as main factor or in interaction, in any of the statistical models.

Fig. 3
figure 3

Carcinus maenas. Weight-specific oxygen consumption rates (QO2) per individual for megalopae. Comparison between three populations measured at two larval salinities (25 PSU, blue and 32.5 PSU, green) and four temperatures (15, 18, 21 and 24 °C). Symbols represent each combination of factors per population. Data shown as average weight specific oxygen consumption rates (µg O2 ind−1 h−1 µg DM−1) ± SE among larvae produced by different females

Populations differed in the sensitivity of QO2 to temperature with Helgoland being less sensitive than the other two. For instance, differences in average QO2 (= ΔQ) between Helgoland and Kerteminde were larger (ΔQ = 2.0 × 10–3, p = 0.01) at 24 °C than at 18 °C (ΔQ = 0.29 × 10–3) and 21 °C (ΔQ = − 1.1 × 10–3) where no significant differences were found (p > 0.60, all Tukey test: data by salinity and female of origin pooled; 15 °C not considered). In addition, differences between Cádiz and Kerteminde were twice as large at 18 °C (ΔQ = − 2.7 × 10–3, p = 0.00011) and 21 °C (ΔQ = − 2.6 × 10–3, p < 10–5) than at 24 °C (ΔQ = − 1.2 × 10–3) where no significant difference was found (p = 0.47). A further check, through a model using temperature as a continuous linear predictor (all temperatures included), also retained the interaction of population and temperature (ΔAIC = 2, likelihood ratio test: p < 0.05 vs additive model), and computed ca. twice higher regression slopes in Cádiz (1.17 × 10–3) and Kerteminde (1.56 × 10–3) than on Helgoland (0.67 × 10–3).

Survival

Larval survival for the Cádiz population showed a pattern consistent with TMLS; detrimental effects of the low salinity treatment (= 25 PSU) were not detected at the higher tested temperatures). In addition, survival increased with temperature, with higher survival at the highest temperature (24 °C) in comparison to Kerteminde and Helgoland (Fig. 4). This made the Cádiz population slightly different from the one of Helgoland and Kerteminde, where in 32.5 PSU there was no increase of survival at 21–24 °C, as compared to 15 or 18 °C. In addition, the response of the Cádiz (and Helgoland) populations differed from that of Kerteminde where survival was depressed at 21 and 24 °C at 32.5 PSU. Best statistical models for survival retained female as a random factor (Table S1). Most of the 2-way interactions were retained (S × T and P × T for raw data; and all 2-way interactions for logarithmic data), indicating interactive effects of population, salinity and temperature. However, at the lowest temperatures, the survival was consistently lower at the low salinity treatment (p < 0.05 for 15 °C and 18 °C). In general, significant differences (p < 0.05) in survival among salinities indicated lower survival at lower temperatures (15 and 18 °C) for all populations.

Fig. 4
figure 4

Carcinus maenas. Cumulative average survival to megalopa. Comparison between three populations reared at two larval salinities (25 PSU, blue and 32.5 PSU, green) and four temperatures (15, 18, 21 and 24 °C). Symbols represent each combination of factors per population. Data shown as mean values ± SE among larvae produced by different females (n = 5 or 3). Asterisks represent significant differences between salinities for larvae of the Cádiz population for each temperature. Already published data are presented in grey sections (for significant differences in those populations see Šargač et al. 2021)

Biomass and elemental composition

The dry mass, carbon and nitrogen content of megalopae from the Cádiz population did not show clear differences among temperature treatments and were consistently lower at lower salinity across the whole temperature range (Fig. 5). Contrary to what was observed for the Helgoland population, there was no evidence of mitigation at high temperatures (21–24 ºC) in the Cádiz population. Instead, the response pattern to the combination of temperature and salinity resembled the one found for Kerteminde. The megalopae from the Cádiz population showed the lowest dry mass, carbon and nitrogen content. In addition, initial dry mass of freshly hatched larvae (Fig. S1) showed no significant differences between Cádiz and Kerteminde (9.16 µg ind−1 for Cádiz; 8.12 µg ind−1 for Kerteminde), while Helgoland larvae had significantly higher dry mass at hatching in comparison to the other two populations (10.7 µg ind−1; p < 0.05). In all biomass analyses, female of origin was kept as a random factor (Table S2) highlighting important variations among larvae from different females. The best models kept some 2-way interactions (P × T for DM; S × T for C and N). Analyses for Cádiz and Helgoland between all eight treatments showed no significant interaction for DM and significant 2-way interaction (S × T) for C and N. For all cases, the additive terms were important and were retained in the best statistical models.

Fig. 5
figure 5

Carcinus maenas. Dry mass, carbon and nitrogen content (µg ind−1) of megalopae among three populations reared at two larval salinities (25 PSU, blue and 32.5 PSU, green) and four temperatures (15, 18, 21 and 24 °C). Symbols as in Fig. 4

Duration of development and growth

Duration of development to megalopa showed similar patterns in Cádiz to the other two populations. Analyses comparing all three populations, as well as a comparison for Cádiz vs Helgoland for all eight treatments, showed important interactive effects of salinity and temperature (S × T) and no variation among populations (Table S2). In the Cádiz population, higher temperatures mitigated negative effects of lower salinity at 15 °C (p < 0.05), following the pattern observed in Helgoland (Fig. 6).

Fig. 6
figure 6

Carcinus maenas. Duration of larval development (days) from hatching to megalopa among three populations reared at two larval salinities (32.5 PSU, green and 25 PSU, blue) and four temperatures (15, 18, 21 and 24 °C). Symbols represent each combination of factors per population. Data shown as averages of developmental duration ± SE among larvae produced by different females. Asterisks represent significant differences between salinities for the Cádiz population in different temperatures. Already published data is presented in grey sections (for significant differences in those populations see Šargač et al. 2021)

Best model for instantaneous growth rates of megalopae (Table S2) kept the 3-way interaction (P × T × S) and showed significant differences among the populations and treatments. There was an effect of salinity in the Cádiz population, where lower salinity significantly caused a decrease in growth at 15, 18 and 21 °C (Fig. 7, p < 0.05). In comparison to the other two populations, megalopae from the Cádiz population had significantly (p < 0.05) lower growth at 24 °C at both salinities, and at 21 °C when exposed to 25 PSU. In addition, the Kerteminde population had significantly (p < 0.05) faster growth at 18 °C in both salinities, compared to Cádiz and Helgoland.

Fig. 7
figure 7

Carcinus maenas. Instantaneous growth rates (dry mass) from hatching to megalopa among three populations reared at two larval salinities (32.5 PSU, green and 25 PSU, blue) and four temperatures (15, 18, 21 and 24 °C). Symbols represent each combination of factors per population. Data shown as average individual growth (µg ind−1 day−1) ± SE among larvae produced by different females. Asterisks represent significant differences between salinities for the Cádiz population in different temperatures. Already published data is presented in grey sections (for significant differences in those populations see Šargač et al. 2021)

Discussion

We studied larval performance of the shore crab C. maenas from three European populations and found important variations among descriptors of performance in response to temperature and salinity. Such variation is, therefore, part of the diversity (or portfolio) of responses to temperature and salinity. We found population differences in (1) the overall response to temperature and salinity and (2) the response to high temperature, with Cádiz showing increased tolerance to high temperature compared to Helgoland and Kerteminde. These two general findings are discussed below in separate sections, followed by the implications of the portfolio of responses to climate change.

Interactive responses to temperature and salinity

Mitigation effects of lower salinity at higher temperatures did not occur consistently across all rates and on body mass. TMLS was present in survival in all three populations, with a strong effect in Helgoland, and a weaker pattern in Cádiz and Kerteminde. TMLS has been found in both the North and the Irish Sea (Spitzner et al. 2019; Torres et al. 2020), i.e., in areas that are characterised by moderately high salinities (30–35 PSU), but it is much weaker at Kerteminde especially at salinity 20 PSU (Šargač et al. 2021) as compared to that of Helgoland and Cádiz (see Figure S2 for comparison). The weak performance of larvae from Kerteminde at low salinities is striking because that population experiences also low salinities (around 15 PSU). Our findings confirm that TMLS is also present in the southern population (Cádiz) coming from the warmer area. This indicates that TMLS is a general pattern affecting survival in C. maenas from the north to the south of their native distribution, irrelevant of their natural osmotic or thermal conditions. For Cádiz, however, the negative effects of low salinity on growth and body mass were not mitigated at high temperature (unlike Helgoland).

For Cádiz, mitigation effects at the higher tested temperatures (21, 24 °C) were found for survival and developmental time, but not for growth rate and body mass. Hence, for Cádiz and in consistence with previous analyses (De Laender 2018; Torres and Giménez 2020), responses at a given level of organisation are not necessarily found at other levels. Because of non-linearity and compensatory effects, the patterns found in physiological responses (e.g. growth rates) do not necessarily appear in those of fitness (e.g. larval survival). Compensatory responses for the Cádiz population may be understood here as a hierarchy of processes (Knowlton 1974), operating at lower salinity but contingent on the temperature conditions. At high temperatures, survival in response to lower salinity was prioritised at expenses of growth; however, at low temperatures reductions in growth were not associated with consistently high survival. Reductions of body mass observed at any temperature may, however, have consequences for fitness at the post-metamorphic stages, which were not quantified here. For instance, reduced body size or mass has been associated with low survival rates in many invertebrates (Giménez et al. 2004; Pechenik 2006; Torres et al. 2016) and increased vulnerability to cannibalism and predation in C. maenas (Moksnes 2004).

We were interested in exploring links among metabolic, developmental, and growth rates because such rates are expected to drive body size at metamorphosis, which is a fitness-related trait (Werner 1988; Hentschel 1999; Gotthard 2001; Heyland et al. 2011; Torres and Giménez 2020). For Helgoland, growth and developmental rates were depressed at lower salinity only when larvae were exposed to the lowest tested temperature, but not when they were exposed to higher temperatures. This is striking because larvae from the Helgoland population are not expected to experience temperatures higher than 18 °C, except perhaps in the Wadden Sea at the peak of summer. One would therefore hypothesise for such a population, that low salinity would lead to metabolic depression at 15 °C or lower temperatures that are expected in the field. However, oxygen consumption in the megalopae responded consistently to temperature, but not to salinity. In addition, the pattern found for Helgoland did not differ much from those found in Cádiz and Kerteminde, where high temperatures did not mitigate the effect of low salinity on growth and body mass. One may argue that, at least for Helgoland and Kerteminde, there was some level of reduction of VO2 driven by lower salinity in larvae exposed to the lowest temperature, but that the variability in the responses was too large for an interaction between temperature and salinity to be retained in a statistical model. However, when oxygen consumption rates were expressed on a per dry mass basis (as QO2) such differences disappeared. In terms of QO2, the response to temperature compared well with that reported in a previous study (Dawirs 1983). Hence, we have no evidence of a clear antagonistic effect of temperature and salinity on the oxygen consumption rates of the megalopa in any of the studied populations. A possible explanation is that, at low temperatures, moderately low salinity (25 PSU) caused reduced feeding or assimilation rates, but larvae kept spending energy in sustaining essential functions and coping with osmotic stress. Alternatively, the ability to osmoregulate of the megalopae results in that oxygen consumption rate is little affected by the tested salinities, and variation in body mass or chemical composition was not detected. Osmoregulation is considered energetically costly which usually causes higher (but not lower) oxygen consumption rates (Dehnel 1960; Rivera-Ingraham et al 2016) in order to fuel the catabolism of osmotically active amino acids (Gilles 1973) and active salt exchange. Since at 25 PSU, megalopae of C. maenas from the Helgoland population osmoregulate, one would expect an effect of such salinity on oxygen consumption rates, as compared with 32.5 PSU where larvae do not need to regulate the haemolymph osmolality (Cieluch et al 2004). Since in addition, the osmoregulatory capacity of C. maenas megalopa from Helgoland increases with temperature (Torres et al. 2021a) one would expect that the effect of salinity on oxygen consumption rates would depend on temperature. Clearly, such effect should have occurred at least for the Helgoland population. Alternatively, the antagonistic effect of lower salinity and higher temperature on oxygen consumption rates may occur in the zoeal stages, which are practically osmoconformers (such as zoea II–IV). In crustaceans, the respiration sensitivity can change during ontogeny (Brown and Terwilliger 1999; Leiva et al. 2018). In zoeae I of C. maenas of the Helgoland population, Anger et al. (1998) found that low salinity depressed oxygen consumption rates. However, in their study, the effect of low salinity on such rates did not occur at the salinities tested by us (at T = 18 °C, VO2 = 0.040 μg O2 h−1 ind−1 at both 25 and 32 PSU in 2-day-old larvae: Table 1 in Anger et al. 1998). Instead, the drop in VO2 occurred at 20 PSU where we have recorded important mortality rates during the first and subsequent larval stages (Spitzner et al. 2019; Šargač et al. 2021). Perhaps oxygen consumption rates are depressed at salinities that are lower than those tested by us. Those salinities will be experienced in the Baltic Sea and might be experienced near the coast in the North Sea, in regions of freshwater influence. For instance, larvae produced from adult populations of the Wadden Sea, should experience low salinities near the mouth of estuaries such as, e.g. the Elbe and Wesser rivers (Bils et al. 2012, Fig. 2). Maps on larval spatial distribution and overlap with salinity fields are missing. However, biophysical models suggest an exchange of larvae between, e.g. Helgoland and Wadden Sea populations exposed to brackish waters from the Elbe and Wesser (Moksnes et al. 2014).

The differences between populations in how temperature modulates the effect of lower salinity on physiological rates, highlight the necessity to quantify among population responses to climate-driven environmental change. We found that survival and developmental time of larvae from the Cádiz population showed a pattern consistent with TMLS as found for Helgoland (Spitzner et al. 2019; Šargač et al. 2021) and the Irish Sea (survival to zoea II: Torres et al. 2020). At higher temperatures, osmoregulatory abilities might have been enhanced via various mechanisms operating on the molecular level (Pörtner et al. 2017) and (for the Helgoland population) those temperatures can amplify the capacity of C. maenas zoea I and megalopae to osmoregulate (Torres et al. 2021a). Most likely, the same capacity is present in the Cádiz population.

Another important difference was in the instantaneous growth rates, lower in Cádiz than in Helgoland and Kerteminde (Fig. S1). Such lower growth rates explain the reduced biomass of the megalopa because developmental times differed little among populations (Fig. S2). Moreover, freshly hatched Cádiz larvae had lower initial biomass than Helgoland, but slightly higher than Kerteminde (Fig. S3). In addition, in larvae from Cádiz, the mitigation effect of high temperature on growth (and on body mass) was weaker than the one found in the Helgoland population. This could be explained by the fact that larvae from the Cádiz population experience salinities that are higher than the 32.5 PSU used in our experiment. Moreover, due to habitat differences, some physiological process could be more sensitive to lower salinity (experimental salinities) than in, e.g. the Helgoland population. While larvae drifting from the Helgoland population into the German Bight are likely to be exposed to moderately low salinities (25–33 PSU), those of Cádiz are likely to experience salinities > 35 PSU in the Gulf of Cádiz (Criado-Aldeanueva et al. 2006). Furthermore, crabs from the local populations of Cádiz can be found in the coastal lagoons and ponds where salinity can reach values of 65 PSU (Drake and Arias 1995). Hence, the continuously high survival in all the treatments, despite the reduced biomass gain, shows high larval tolerance of the Cádiz population.

Effects of high temperature

An important difference between Cádiz and the other populations was the increased survival at the highest temperature, a pattern that was not present in Helgoland nor in Kerteminde. Thus, it appears that the populations differ in the temperature at which the metabolic demand exceeds the animal’s capacity to supply oxygen to tissues and cause a progressive decline in performance (Deutsch et al. 2015; Pörtner 2010). We were not expecting important differences in the response of larvae to high temperatures because, in Cádiz, C. maenas larvae are released earlier than in Helgoland and Kerteminde (January–April vs. April-June) and hence temperature differences may not be so contrasting (Cádiz: ~ 12–20 ºC vs Helgoland: ~ 9–18 ºC). Early release in Cádiz is consistent with data distribution of C. maenas larvae in the Cádiz Bay (November–May, with a peak in February: Drake et al. 1998) and in Portugal, where early zoea are present in the water column already in February (Queiroga et al. 1994) and arrival of megalopae into estuaries occurs already in April and May (Rey et al. 2016). However, the larval season in Cádiz includes April–May, when temperatures in the inlets can be well above 20 ºC (Drake et al. 1998) and those of the Gulf of Cádiz reach 20 ºC (Criado-Aldeanueva et al. 2006). In addition, environmental conditions affecting resource allocation to embryos may have driven the observed differences among populations in the larval response to salinity and temperature.

Differences between the Cádiz and Helgoland populations point to either local adaptation or some form of plasticity (Marshall et al. 2010; Fischer et al. 2011). Our data do not let us infer whether a single or both mechanisms drive the differences. Adult C. maenas exhibit seasonal changes in respiration rates (Klein Breteler 1975) and can become acclimated in warmer waters (Newell 1973). In other crustaceans, populations from different climatic regions can adjust metabolic rates to maintain function at the prevailing temperature (Saborowski et al. 2002). Plasticity, in the form of post-zygotic maternal effects, impact responses to lower salinity and increased temperature in C. maenas: lower salinity experienced during embryogenesis pre-empt the development of mitigation effects in the larval stages of a population of the Irish Sea (Torres et al. 2020), Helgoland and Kerteminde (Šargač et al. 2021). Perhaps higher salinities experienced during embryogenesis also pre-empt the development of mitigation responses, in the Cádiz population, although they exist at 20 PSU in terms of survival. We also expect that genetic differences between populations underpin the differences in the responses: in a genetic study covering the European coast, Roman and Palumbi (2004) pointed towards slight levels of population structuring for C. maenas, between the Atlantic and the North Sea, which may be associated with differentiation in genes driving both thermal and salinity tolerance. Further genetic studies are needed to determine the role of genetic diversity in driving the diversity of responses existing along the distribution range.

Portfolio of responses

From the Cádiz Bay to the Baltic Sea, C. maenas exhibit an important diversity in larval responses to high temperature and in the joint effect of salinity and temperature on physiological rates and how such rates are related as linked traits. Such diversity represents part of the portfolio of responses which is relevant for understanding response to climate-driven environmental conditions and the potential for invasion. In first instance, the diversity of responses affects the potential to buffer the meta-population from environmental variation. The general trend for marine species, under the effect of climate change, will be to shift their distribution upward or poleward related to sea warming (Chust et al. 2014). For the case of dispersive larvae, the portfolio effect (Schindler et al. 2015) can manifest in how individuals survive the larval phase, connect, and replenish populations as temperature increases. For instance, in the light of warming, southern populations may subsidise northern populations if the European coast is sufficiently connected by larval dispersal and if the direction of transport is poleward. Such connection may occur over the time scale of several generations as a stepping-stone process, as it is likely to have occurred with exotic species invading the North Sea from the English Channel (e.g. Giménez et al. 2020). However, whether such subsidy occurs, also depends on the mechanisms driving the phenotypic variation in the joint effects of temperature and salinity and on the local habitat characteristics. If southern populations were to be all-similar to that of Cádiz (i.e., more sensitive to lower salinity than that of Helgoland), the expansion may be limited if (1) the source of the variation is genetic (instead of adaptive plasticity) and (2) habitats of high latitude are consistently characterised by moderately low salinities. Overall, our work emphasises the need for incorporating multiple populations when quantifying the impact of multiple climate drivers on marine populations.

The second important point concerns the capacity of C. maenas as a global invader and as one of the top 100 species (Roman & Palumbi 2004; Lowe et al. 2000; Young and Elliott 2020). Invasion depends in first instance on transport pathways (e.g. on shipping routes) but also on the capacity of individuals to tolerate the conditions of new habitats (Tepolt and Somero 2014) for which our results are relevant. Our results suggest that different local populations within the home range may maximise the invasion potential in different world regions or locations along a latitudinal gradient. For example, two different invasions occurred in the Atlantic coast of North America: the poleward side of the invaded range is originated from an invasion from Northern Europe and the southern range originated by an invasion from South-Central Europe (Roman 2006; Darling et al. 2008; Jeffery et al. 2018). Invasions to different regions and the current structure appears to be driven by temperature tolerance (Compton et al. 2010; Jeffery et al. 2018), perhaps contributed by differences in larval thermal tolerance, in addition to that of the adults.