Introduction

As the majority of marine animals are water-breathing ectotherms, temperature and dissolved oxygen concentration in sea water represent key environmental variables that influence their fitness and geographic distribution (Pauly 2010; Ern 2019). Indeed, temperature and oxygen may impact the distribution of marine organisms more so than the availability of food resources (Kramer 1987; Pauly 2021). In turn, if global oceans continue to warm and deoxygenate (e.g., Schmidtko et al. 2017; Breitburg et al. 2018; Bindoff et al. 2019; Laufkötter et al. 2020), water-breathing ectotherms may experience deleterious physiological effects associated with thermal stress and hypoxia (Pörtner 2001, 2002). In general, the response of marine water-breathing ectotherms to warming and deoxygenation consists of shifting their distribution to areas with more suitable environmental conditions (Deutsch et al. 2015; Poloczanska et al. 2016). Increasingly, it is also suggested that a response to warming and deoxygenation is the reduction of body size (Gardner et al. 2011).

Ocean-warming-induced changes in body size are predicted to impact water-breathing ectotherms globally (Cheung et al. 2013). Changes in body size can result from the complex interaction of temperature and dissolved oxygen as warming-induced increases in oxygen demands may not be met by an increase in available dissolved oxygen supply (Pörtner 2001; Pörtner and Knust 2007; Pauly 2010, 2021). While reductions in body size have been widely observed for aquatic ectotherms (Daufresne et al. 2009), the mechanisms behind this response remain debated (Lefevre et al. 2017; Pauly 2021).

The temperature-size rule (TSR) describes how ectothermic species grow faster and mature at a smaller size when reared at warmer temperatures (Atkinson 1994). The TSR refers to the plasticity of individuals throughout ontogeny, and thus can display differing affects to both size at maturity and maximum asymptotic size (Hoefnagel et al. 2018; Verberk et al. 2021). In wild populations, clines in body size across latitudinal gradients, including temperature and dissolved oxygen concentrations, can resemble the body-size responses observed at different laboratory rearing temperatures (Horne et al. 2015), referred to in the present study as a temperature-size response. While temperature-size responses are widely observed, several explanations aim to resolve the mechanisms behind the interaction of environment and growth (Deutsch et al. 2015; Pörtner et al. 2017; Clarke et al. 2021; Pauly 2021; Verberk et al. 2021).

A temperature-size response may depend on the species body size (Rubalcaba et al. 2020), sensitivity to changes in temperature or dissolved oxygen concentration (Forster et al. 2011; Hoefnagel and Verberk 2015), life history (Weber et al. 2015; Audzijonyte et al. 2020), geography (Deutsch et al. 2020; Clarke et al. 2021), environmental niche (Rypel 2014), feeding strategy and behaviour (van Rijn et al. 2017), and physiology (Atkinson et al. 2006). In order to elucidate further the mechanisms behind a temperature-size response, it is useful to compare the degree to which responses differ between species, rather than testing if a species displays a response at all (Verberk et al. 2021).

In the present study, we compare the strength of the temperature-size response of species of fishes, crustacean, and squid in New Zealand waters. Our study relies on generalised additive models (GAMs) fitted to trawl research data. The four fish species included in this study are hoki, Macruronus novaezelandiae Hector 1871, Australasian snapper, Pagrus auratus Forster 1801, southern blue whiting, Micromesistius australis Norman 1937, and orange roughy, Hoplostethus atlanticus Collett 1889, while the two invertebrate species are New Zealand arrow squid, Nototodarus sloanii Gray 1849, and New Zealand scampi, Metanephrops challengeri Balss 1914. Hereafter, the species are called hoki, snapper, blue whiting, orange roughy, arrow squid, and scampi. Species nomenclature follows the World Register of Marine Species (Horton et al. 2021).

The inclusion of two invertebrate species in the present study enabled the comparison of variable temperature-size responses between taxonomic groups, which has traditionally been based on studies of fishes (Pauly 2021). Further, as the bulk of research focuses on body size at maturation and warming, we investigated the influence of temperature and dissolved oxygen concentration on maximum body length, as the mechanisms between these responses may be different (Hoefnagel et al. 2018). Moreover, considering the complex interaction of temperature and dissolved oxygen within a temperature-size response, we assessed the relative importance of model covariates in order to parse between temperature, dissolved oxygen concentration and geographic location (representing other, latent environmental variables) in determining maximum body length of the study species. Doing so may reveal the relative influence of temperature and oxygen in eliciting a temperature-size response in water-breathing ectotherms. Overall, results from related research have found that a negative temperature-size response (e.g., smaller in warmer temperatures) was strongest in larger or more active species (Forster et al. 2012; van Rijn et al. 2017); therefore, we hypothesise that larger study species will exhibit a relatively stronger temperature-size response than smaller species.

Methods

Study area

Our study area was the exclusive economic zone (EEZ) of New Zealand, which spans latitudes 30 to 55°S (Fig. 1). The region is characterised by the influence of several fronts, which result in a broad temperature gradient, as summer sea surface temperatures reach over 21 °C in the north and 13 °C in the south, while winter sea surface temperatures can be as low as 5 °C in the north and 2 °C in the south (Garner 1969). From the north north-west, the Tasman Front and the Subtropical Front direct warm water across both islands and most of the EEZ, while colder and more nutrient-rich water is brought by the Subantarctic Front from the south (Leathwick et al. 2006a) (Fig. 1).

Fig. 1
figure 1

modified from Stephenson et al. (2018), and was created using QGIS (QGIS 2022)

Map of the study area, indicating notable areas mentioned in the text, plus bathymetry, cold (violet arrows), and warm (ochre arrows) ocean currents, fronts (brown and violet shade), and the Exclusive Economic Zone (EEZ) of New Zealand (red line). This figure has been

Trawl surveys

The length data employed in analyses were derived from bottom trawl research surveys carried out between 1979 and 2021 (NIWA 2014; SWPRON 2017). These surveys cover the continental shelf surrounding New Zealand, as well as the extensive submarine plateaus throughout the area (Leathwick et al. 2006a). In order to minimise variation between gear and sampling techniques, data were filtered to only include samples gathered by bottom trawl from trawling events that were conducted over a maximum of 3.5 nautical miles. Since the objective of the long-term bottom trawl research surveys is to assess changes in fished stocks, the geographical distribution of sampling effort is biased towards the most productive areas, including the Chatham Rise, the Campbell Plateau, and the Challenger Plateau (Figs. 1 and 2; Francis et al. 2002).

Fig. 2
figure 2

Location of all trawling events included in the trawl research survey database between 1979 and 2021 in New Zealand waters (NIWA 2014; SWPRON 2017). This figure was created using QGIS (QGIS 2022)

Study species

From all species present in the bottom trawl research survey database (n = 253), we selected the six most abundant fish and invertebrate species in terms of the number of individuals and frequency of encounter per sampling event. Hoki is a relatively long-lived, bentho-pelagic fish that is most abundant between 200–800 m depth (Papa et al. 2021). Snapper is a coastal fish species found around rocky reef ecosystems at depths down to 200 m, but is most commonly encountered between 10–50 m (Papa et al. 2021). Blue whiting is a mesopelagic fish found along the outer shelf and slope of New Zealand waters, between 130–800 m depth (Niklitschek et al. 2010). Orange roughy is a slow-growing, long-lived, bathypelagic fish that is found along continental slope margins, seamounts and ridges, between 700–1500 m depth (Papa et al. 2021). Arrow squid is a large pelagic squid that can be found as deep as 600 m depth, although it is primarily found shallower than 300 m (Jackson et al. 2000). Lastly, scampi burrows in soft sediments around the continental shelf and slope between 140–500 m depth (Smith 1999; Yaldwyn and Webber 2011). In the present dataset, the mean depth of species recorded ranged from 58 m (snapper) to 947 m (orange roughy) (Fig. 3a), mean temperatures from 4 °C (orange roughy) to 16 °C (snapper) (Fig. 3b), and mean dissolved oxygen concentrations ranged from 6 mg/L (orange roughy) to 8 mg/L (blue whiting) (Fig. 3c).

Fig. 3
figure 3

Boxplot summary of sampling and environmental information, including the interquartile range, minimum and maximum values (lines), plus outliers (dots, 1.5 times the interquartile range from quartiles), for hoki, snapper, blue whiting, orange roughy, arrow squid and scampi. Summary information includes a mean depth (m), b temperature at maximum depth (°C) from Bio-Oracle, c dissolved oxygen concentration at maximum depth (mg/L) from Bio-Oracle, and d maximum length (cm) records (NIWA 2014; SWPRON 2017; Assis et al. 2018). This figure was created using R Studio (Team R Development Core 2021)

Data and initial exploratory analyses

The total number of length observations for all six study species in the trawl research survey database was 22,661. We defined the maximum length for each species’ record as the 95th quantile value of length measurements (cm) taken at each sampling station (i.e., trawling event). All fish species are reported as standard length (SL), arrow squid as mantle length (ML), and scampi as carapace length (CL) (Fig. 3d). Hoki standard lengths were converted from total length (Kloser et al. 2011), snapper from fork length (Ferrell and Sumpton 1997), and blue whiting from fork length (Cohen et al. 1990) using reported length-length linear relationships. Body length was considered the most robust response variable to estimate maximum body size, rather than weight. Related to feeding, energy storage, or reproduction, weight in marine ectotherms is known to fluctuate seasonally (Craig et al. 2000; Mello and Rose 2005), as well, length–weight relationships vary between species based on body shape and other physiological factors (Jisr et al. 2018).

We identified the Gamma distribution as the best distribution model for maximum length data using diagnostics from the fitdistrplus package (Delignette-Muller and Dutang 2015) in R (Team R Development Core 2021). Other information gathered from the bottom trawl research survey database included latitude and longitude, expressed in Universal Transverse Mercator (UTM) coordinates (i.e., eastings and northings) in the models, as well as mean depth (m) of trawl. Data for mean temperature and dissolved oxygen concentration at maximum depth layer were extracted from the online database Bio-Oracle (Assis et al. 2018) and associated with the latitude and longitude of samples. These data layers are produced using monthly averages of climate data from 2000 to 2014 (Assis et al. 2018). Dissolved oxygen concentration data from Bio-Oracle were converted from mols/m−3 to mg/L in line with unit conversions provided by the International Council for the Exploration of the Sea (ICES n.d.). In situ temperature was recorded during trawling events, but long-term mean temperature and dissolved oxygen concentration derived from Bio-Oracle are likely to better reflect prevailing climatology, and were, therefore, considered to be more relevant for our analyses.

Before model fitting, we explored the distribution of the response variable (maximum length) and explanatory covariates to identify potential outliers and subsequently removed them (Zuur et al. 2010). Following this data trimming, the total number of observations for the six study species was reduced from 22,661 to 22,631 (Fig. 4). We then explored the range of maximum length (cm), depth (m), and environmental covariates of the study species, as well as the relationship between environmental covariates and location. Thereafter, we evaluated the degree of collinearity between temperature, dissolved oxygen concentration, depth, eastings, and northings for each species (Supplementary Fig. 1). A Pearson’s correlation coefficient greater than 0.7 in absolute value was considered indicative of collinearity (Leathwick et al. 2006b; Dormann et al. 2013). This collinearity analysis was warranted as GAMs (and other regression methods) are sensitive to correlated continuous covariate variables (Guisan et al. 2002; Dormann et al. 2013). For comparison purposes, depth was omitted from GAMs as it was found to be collinear with temperature and dissolved oxygen concentration for hoki and snapper respectively (Supplementary Fig. 1). This allowed the same full factorial model for all species to be undertaken, whilst preserving the model form that incorporates these key environmental variables central to our hypotheses. Moreover, we expect variation in depth to be contained within the effect of geographic location as a model covariate.

Fig. 4
figure 4

Observed distribution of the data used in the present study for the six study species, including hoki (n = 9420), snapper (n = 1673), blue whiting (n = 1308), orange roughy (n = 6763), arrow squid (n = 2912), and scampi (n = 555) (NIWA 2014; SWPRON 2017). Each record indicates a trawl sample station, where the 95th quantile length measurement was extracted. This figure was created using QGIS (QGIS 2022)

Model fitting and evaluation

We utilised GAMs to investigate the influence of temperature and dissolved oxygen concentration on the species’ maximum lengths. A tensor product smooth between eastings and northings (the effect of geographic position) was included to serve as a proxy for latent environmental covariates that were otherwise not included in the model (e.g., depth) and account for broad-scale spatial autocorrelation in the length data (Wood 2006; Grüss et al. 2016, 2018a, 2021).

To facilitate equivalent comparisons between species, we initially fitted the same “full” Gamma GAM to maximum length data for each species using the “gam” function from the mgcv package in the R environment (Wood 2006; Team R Development Core 2021) following the equation,

$$\mathrm{g}\left(L\right)=te\left(X,Y\right)+s\left(\mathrm{Temperature}\right)+s(\mathrm{Oxygen})$$
(1)

where L is maximum length of the species under consideration, g is the log-link function between maximum length and each term on the right side of the equation, te(X,Y) is the tensor product smooth between eastings and northings, and s is the thin plate regression spline with shrinkage (bs = “ts” specified in the smooth function of “gam” within the mgcv library) fitted to temperature and dissolved oxygen concentration (Grüss et al. 2018a). Each thin plate regression spline was limited to four degrees of freedom (k = 4) to prevent overfitting and to preserve the interpretability of GAM results (Mannocci et al. 2017; Grüss et al. 2018a; Weijerman et al. 2019).

After running all GAMs using the same form (“full model”; Eq. 1) for all species, collinear or non-significant covariates were removed, including oxygen for snapper and arrow squid (non-significant), and northings for scampi (collinear with oxygen); then, models were re-run. We refer to the re-run models as the “final models.” The adjusted final models were compared to full models via Akaike’s information criterion (AIC) via the “logLik” function of the mgcv package in R (Wood 2006; Team R Development Core 2021). When running the final models, an extra penalty was applied to temperature and dissolved oxygen concentration as their smoothing parameter approached zero, which allowed for the complete removal of temperature or dissolved oxygen concentration from the model when its smoothing was equal to zero (Marra and Wood 2011; Grüss et al. 2014). Ultimately, using the results from final GAMs, we predicted the maximum length of each study species in light of its observed temperature range.

In order to acknowledge the potential effects of fishing pressure on reducing the maximum lengths of our study species, all final models were re-run twice, using records split into two periods based on their collection dates. These periods include a “past” fishing period (1979–1998) and a “recent” fishing period (1999–2021). Here, we assume that the past fishing period, and fishing effort preceding this period (“early” commercial fishing), were more likely to remove the largest individuals from stocks. In turn, we assume that there will no longer be a strong effect of fishing reflected in the predicted maximum body length of species across their temperature range within the recent period relative to the previous period. These assumptions are based on the observation that fished stocks are able to recover in maximum length following a reduction in exploitation rate (Worm et al. 2009).

An evaluation procedure similar to the leave group out cross validation (LGOCV) procedure was employed to validate full and final GAMs (Kuhn and Johnson 2013; Grüss et al. 2016) using two performance metrics: (1) the adjusted coefficient of determination (adjusted R2), which is indicative of the proportion of the variance in maximum length explained by GAMs (Grüss et al. 2018b), and Spearman’s rank correlation coefficients (Spearman’s rho’s) between the predicted and observed maximum length values. Briefly, our evaluation procedure for each model consisted of (1) splitting randomly the dataset of interest into a training dataset, which receives 70% of the data, and a test dataset, which receives 30% of the data; (2) repeating this ten times so that 10 training datasets and 10 test datasets are generated; (3) fitting a model to each of the 10 training datasets; and (4) performing an evaluation of the 10 fitted models with the test datasets using performance metrics (adjusted R2 and Spearman’s rho’s in our case) (Grüss et al. 2016, 2020; Egerton et al. 2021). We considered that a given model passed the evaluation test if the median adjusted R2 across all 10 replicates were greater than 0.1 (Legendre and Legendre 1998; Grüss et al. 2016) and if the median Spearman’s rho’s across the 10 replicates was significantly different from zero (Grüss et al. 2014, 2021; Weijerman et al. 2019; Egerton et al. 2021).

In order to estimate the relative importance of temperature and dissolved oxygen concentration in our final GAMs, we followed the same approach as Grüss et al. (2016). For each species, this approach compares the predictions of the final model with the predictions of models after random permutation of the values of a given predictor (temperature, dissolved oxygen concentration, eastings, or northings) within the complete dataset fed into the GAM (“random GAMs”) (Thuiller et al. 2012; Grüss et al. 2016). Thereafter, the index of relative importance of the predictors (temperature, dissolved oxygen concentration, eastings, or northings) in explaining maximum length is obtained as one minus the Pearson’s correlation coefficient between the predictions of the final model and the predictions from the random models (Grüss et al. 2016; Dove et al. 2020; Bolser et al. 2020).

Results

The length data employed in our models were found to be representative of those observed for the study species, including adult size classes, for all species (Table 1). The temperature data employed in models were found to represent a portion of the reported temperature ranges of species (Supplementary Table 1). Across the study area, there was a north to south latitudinal gradient of temperature and dissolved oxygen concentrations from approximately 19 °C and 7.5 mg/L to 7.5 °C and 9.5 mg/L. In deeper water layers, these gradients were still evident, although with higher variability (18 °C and 4.9 mg/L to 1.5 °C and 9 mg/L; Supplementary Fig. 2).

Table 1 Observed and reported maximum length (cm) for hoki (standard length, SL), snapper (SL), blue whiting (SL), orange roughy (SL), arrow squid (mantle length, ML), and scampi (carapace length, CL). Reported maximum lengths: 1. Cohen et al. 1990 2. Randall et al. 1998, 3. Thomsen 1998, 4. Roper et al. 1984, 5. Cryer et al. 2005. Note: for comparison purposes, the reported maximum lengths for hoki, snapper, and blue whiting were converted to SL using published length-length linear relationships (see methods)

Temperature was found to have a significant effect on the maximum length of all six study species, while dissolved oxygen was found to be significant for three species: hoki, blue whiting, and orange roughy (Table 2). Following removal of collinear or non-significant covariates, the final GAMs for snapper, arrow squid and scampi differed from their original full models, while the final models for hoki, blue whiting, and orange roughy were identical to the original full models. Overall, the adjusted final models for snapper, arrow squid, and scampi provided little to no improvements to the performance of GAMs, as was confirmed by AIC comparison (Table 3), but the final model for scampi passed the evaluation tests while its original full model did not. Further, dissolved oxygen concentration was significant for New Zealand scampi following removal of northings as they were collinear with oxygen (Table 3). Thus, dissolved oxygen concentration was significant for four species within final models. As such, all final GAMs passed the evaluation tests, with the median adjusted R2 values ranging between 0.104 (orange roughy) and 0.44 (hoki), while median Spearman’s rho values ranged from 0.29 (arrow squid) to 0.55 (hoki) (Supplementary Fig. 3).

Table 2 Results for the full generalised additive models (GAMs) of all study species. Models express maximum length as a function of temperature, dissolved oxygen concentration, and an interaction term (tensor product smooth) between eastings and northings (Eq. 1)
Table 3 Results for the final generalised additive models (GAMs) for snapper, arrow squid, and scampi; final GAMs result from the removal of collinear or non-significant covariates from the full GAMs (given by Eq. 1). Also listed is the Akaike’s information criterion (AIC) for final and full GAMs. Note that the final GAM for scampi included the smooth effect of eastings instead of an interaction term between eastings and northings, as northings were found to be collinear with dissolved oxygen concentration and were, therefore, removed from the scampi GAM

The species differed in their thermal ranges and temperature-size responses, but all showed decreasing maximum length with temperature (Fig. 5). Temperature was generally a more important covariate than dissolved oxygen concentration in explaining maximum length in the final models, while eastings and northings were the most important covariates for all species except for hoki and scampi (Fig. 6). Finally, the predicted maximum lengths across the past fishing period (1979–1998) and those from the recent period (1999–2021) generally displayed the same relationship between maximum body length and temperature, with hoki the only species predicted to be smaller across the warmer extent of its distribution within the recent time period (Supplementary Fig. 4).

Fig. 5
figure 5

Maximum lengths predicted by final GAMs (solid line) with 95% confidence interval bands, as a function of temperature (°C) for hoki, snapper, blue whiting, orange roughy, arrow squid, and scampi. The dashed lines of matching colour represent the reported preferred temperature range for each species, including mean preferred temperature (point), as reported by FishBase (for fishes) (Froese & Pauly 2021; www.fishbase.org) and SeaLifeBase (for invertebrates) (Palomares & Pauly 2021; www.sealifebase.org). This figure was created using R Studio (Team R Development Core 2021)

Fig. 6
figure 6

Estimated relative importance of eastings, northings, dissolved oxygen concentration, and temperature (with 95% confidence intervals), in explaining maximum length predicted by the generalised additive models of a hoki, b snapper, c blue whiting, d orange roughy, e arrow squid, and f scampi. See the main text for a description of the calculation of indices of relative importance. This figure was created using R Studio (Team R Development Core 2021)

Discussion

We found that the maximum body length of all six study species displayed an inverse relationship with temperature. The largest species, hoki, displayed the strongest temperature-size response, with maximum length predicted to decrease by approximately 72% (from 85.5 cm to 24.1 cm SL) across the sampled temperature range (Fig. 5). The predicted decrease in maximum body length was similar for the remaining species, ranging between 30% for arrow squid and 18% for blue whiting (Fig. 5). Overall, the strength of the temperature-size response for the study species was likely underestimated, as our data covered only part of the study species’ known temperature and dissolved oxygen concentration ranges (Supplementary Table 1; Fig. 5). The two species with the widest observed temperature ranges, hoki and arrow squid, displayed the strongest reductions in maximum body length. Conversely, blue whiting maintained the narrowest observed temperature range, resulting in the smallest predicted reduction in maximum body length (Fig. 5).

Our results support our initial hypothesis that the largest species display the strongest temperature-size response. Despite a similar temperature range extent, hoki’s relative reduction in maximum body length was two times stronger than that of arrow squid. Experimental results from Forster et al. (2012) found that the strength of a temperature-size response in water-breathing ectotherms increased with their maximum body lengths. Further, when investigating the temperature-size responses of 74 fish species across the Mediterranean Sea, van Rijn et al. (2017) concluded that active large species were those that displayed the largest reductions in maximum body size with temperature. While species’ traits were not considered in the models developed in the present study, hoki is by far the largest of our study species. Species’ traits, life histories and geographic variables may influence the strength of temperature-size responses, and, thus, should be considered in future models.

According to the gill-oxygen limitation theory (GOLT), marine ectotherms reach an asymptotic body size due to the limitations of gill surface area to body weight ratio, when the two-dimensional respiratory surface (gill) area is unable to provide sufficient oxygen required for the growth and function of three-dimensional bodies (Pauly 2010, 2021; Pauly and Cheung 2018). Therefore, one mechanism under the GOLT would see species reduce their maximum body length (i.e., limit growth) with increasing temperatures, in order to increase the gill surface area to body weight ratio and, thus, service temperature-dependent oxygen demands (Pauly and Cheung 2018). Another mechanism, the oxygen- and capacity-limited thermal tolerance (OCLTT) concept, states that species have upper thermal limits, whereby an animal’s capacity to adequately supply oxygen to maintain performance diminishes as it approaches this limit (Pörtner 2001, 2002). According to the OCLTT, an animal’s ability to supply oxygen does not scale sufficiently with body size relative to oxygen demand, whilst increasing temperatures and oxygen demand will compound oxygen limitation (Pörtner et al. 2017; Audzijonyte et al. 2019).

While we were unable to mechanistically investigate the variable temperature-size responses between species, our results suggest that large-bodied species’ temperature- and size-dependent growth efficiencies were constrained under increasing temperatures (Pörtner 2001; Kozłowski et al. 2004; Pauly 2010; Verberk et al. 2011; Hoefnagel et al. 2018; Rubalcaba et al. 2020), as the assumed increases in oxygen demand (with temperature) were not met despite dissolved oxygen concentration being above known critical levels for marine organisms (Cook and Herbert 2012; Brennan et al. 2016; Shi et al. 2021).

We found that temperature was a more important covariate than oxygen for five of the six study species, with blue whiting the exception (Fig. 6). Across the study area, dissolved oxygen levels were at or near saturation (Supplementary Fig. 2). Therefore, limitations of available dissolved oxygen were unlikely to drive a temperature-size response in species. This conclusion adheres to the proposed mechanisms including the GOLT, and the OCLTT, in that relative reductions in maximum size (i.e., limited growth) are likely required in order to meet the oxygen demands brought about by temperature-induced increases in metabolic activity of water-breathing ectotherms (Pörtner and Peck 2010; Pauly 2021). Furthermore, the inclusion of two invertebrate species in the present study allows for further generalisation of temperature-size responses in marine ectotherms, which has been primarily developed with regard to fishes (Pauly 2010, 2021). As crustaceans and cephalopods generally have a higher tolerance of hypoxia (Shi et al. 2021), future research should continue to compare temperature-size responses between taxonomic groups, and the interactive effect of temperature and dissolved oxygen concentration on the metabolic scope of water-breathing ectotherms. As such, it may be that species that are less tolerant of low oxygen levels may experience a stronger temperature-size response with warming, even if dissolved oxygen levels remain at or near saturation.

Unsurprisingly, the most important variable in explaining maximum body length in our models was generally geographic position (eastings and northings), as this represents latent ecological and environmental variables (Wood 2006; Grüss et al. 2016, 2018a, 2021; Bolser et al. 2020; Egerton et al. 2021). These variables could relate to variation in depth, the availability of suitable habitat, ecological competition, prey availability, size selective predation or mortality, and/or fishery effects (Connolly and Roughgarden 1999; Tu et al. 2018; Penn et al. 2018; Deutsch et al. 2020; Kulatska et al. 2021). These factors associated with a species’ distribution may influence the temperature-size response. By including the interaction term between eastings and northings in our models, we were able to compare the importance of temperature and dissolved oxygen concentration in explaining maximum length relative to that of other, unmeasured variables.

It is important to note that all of the species included in the present study are targeted by fisheries in New Zealand waters. The present study uses only survey data, and we did not include fishing pressure in our models. Overall, the maximum body lengths of our study species may have been truncated, as fisheries preferentially target the largest individuals in a population, such that our ability to display the extent of a temperature-size relationship may be confounded by fishing pressure (Cheung et al. 2013; Tu et al. 2018). Yet, the maximum body length ranges for all of our study species approximate adult sizes (Table 1; Fig. 3d). Further, by comparing more recent survey data (1999–2021) to past survey data (1979–1998; Supplementary Fig. 4), we assume that the influence of fishing on removing large individuals from the population no longer influences the maximum length of species across our study area (i.e., beyond previous, higher intensity fishing era). Given that we found few differences in our predicted trends between time periods, we conclude that our results are not severely confounded by fishing effort in our study, but most likely due to the wide environmental and temporal gradient of the study.

In conclusion, our results support the hypothesis that larger-bodied aquatic ectotherms experience the strongest temperature-size responses. They also support expectations from the GOLT and OCLTT theories that increase in oxygen demand may be size- and temperature-dependent. Climate change has already resulted in measurable impacts on marine organisms and ecosystems. This includes the poleward shift in distribution of thousands of species as they track their thermal affinity (Chaudhary et al. 2021), to the increasing number and spatial coverage of marine areas experiencing episodic or chronic hypoxia throughout the last century (Limburg et al. 2020). The observed and projected changes to ocean temperature and dissolved oxygen concentration will likely tax the metabolic demands of water-breathing ectotherms across depths and latitudes (Oschlies 2021), but to varying degrees between species, size-classes or functional groups (Audzijonyte et al. 2020; Rubalcaba et al. 2020). Thus, our results extend the findings of previous studies proposing that larger-bodied species may experience the most significant temperature-size responses with future ocean warming.