Testing Bergmann’s rule in marine copepods

Macroecological relationships provide insights into rules that govern ecological systems. Bergmann’s rule posits that members of the same clade are larger at colder temperatures. Whether temperature drives this relationship is debated because several other potential drivers covary with temperature. We conducted a near-global comparative analysis on marine copepods (97 830 samples, 388 taxa) to test Bergmann’s rule, considering other potential drivers. Supporting Bergmann’s rule, we found temperature better predicted size than did latitude or oxygen, with body size decreasing by 43.9% across the temperature range (-1.7 to 30ºC). Body size also decreased by 26.9% across the range in food availability. Our results provide strong support for Bergman’s rule in copepods, but emphasises the importance of other drivers in modify-ing this pattern. As the world warms, smaller copepod species are likely to emerge as ‘winners’, potentially reducing rates of fisheries production and carbon sequestration.


Introduction
Although much of ecology has focused on seemingly large differences among organisms and ecosystems, especially in terms of their spatial and temporal variation, the discipline has been criticised for lacking general unifying principles (Lawton 1999, Allen and Hoekstra 2015, O'Connor et al. 2019. Macroecology seeks to discover ecological principles by averaging over finer-scale variation to reveal large-scale statistical relationships Maurer 1989, Gaston andBlackburn 2008). One such unifying principle is the importance of body size. Body size plays a central role in the physiology and ecology of organisms, governing processes such as respiration, metabolism, movement and trophic interactions (Peters 1986, Woodward et al. 2005, Yvon-Durocher et al. 2011. As a result, body size is increasingly used in ecosystem models to generalise traits across vast numbers of taxa, particularly in the marine environment (Blanchard et al. 2017). Better understanding of global patterns of body size, and their drivers, will provide stronger unifying principles in ecology and help support development of ecosystem models.
One of the earliest relationships identified in ecology was formulated by Bergmann (1848) in German, although translation into English (Mayr 1956, James 1970 subsequently contributed to confusion surrounding its definition. Bergmann's rule has been tested and verified for many taxa, including mammals Maurer 1989, Ashton et al. 2000), birds (James 1970, Ashton 2002, fish (Wilson 2009, Saunders andTarling 2018), reptiles Feldman 2003, Angilletta et al. 2004), amphibians (Olalla-Tárraga and Rodríguez 2007), phytoplankton (Sommer et al. 2017), nematodes (Van Voorhies 1996), insects (Chown and Gaston 2010, Osorio-Canadas et al. 2016, Tseng and Soleimani Pari 2019, crustaceans (Manyak-Davis et al. 2013, Garzke et al. 2015, Leinaas et al. 2016) and planktonic ciliates (Wang et al. 2020). Yet, the taxonomic level at which Bergmann's rule applies is commonly debated. Some theorists consider it to be intraspecific (Ashton et al. 2000, Olalla-Tárraga 2011, whilst others consider it interspecific (Blackburn et al. 1999, Hessen et al. 2013, leading to some confusion in the literature. Further, since the pattern was first described in endotherms, some question its applicability to ectotherms (Pincheira-Donoso et al. 2008, Watt et al. 2010. In this study, we consider Bergmann's rule to be defined as Bergmann (1848) himself defined it: species from the same taxonomic clade (here subclass) are generally smaller in warmer regions and larger in cooler regions (i.e. we consider only the interspecific version of Bergmann's rule, and don't consider size differences within species, Blackburn et al. 1999). It is clear that regardless of the precise definition used, much can be gained by investigating Bergmann's rule (Olalla-Tárraga 2011) because spatial patterns in body size at all taxonomic levels strongly influence the ecology of a system (Peters 1986).
Commonly, temperature is considered the primary driver of Bergmann's rule, although some have argued that other drivers might modify the anticipated patterns across taxonomic groups (James 1970, Millien et al. 2006, Yom-Tov and Geffen 2011. This is because latitude (and therefore temperature) is confounded with light availability, oxygen concentration (in aquatic systems), predation rate and food availability (Ho et al. 2010). Of these, light is unlikely to be a direct driver of Bergmann's rule because many species that follow the pattern do not depend directly on light. However, primary productivity depends on light, which could in turn influence the food available for many groups. It is also likely that drivers of Bergmann's rule differ between endotherms and ectotherms. For endotherms, it is generally accepted that temperature is the selective mechanism driving Bergmann's rule, with species from cooler regions conserving heat by being larger and consequently having lower body surface-area-to-volume ratios (Mayr 1956, Hessen et al. 2013. But this is not true for ectotherms (Olalla-Tárraga et al. 2006, Watt et al. 2010. Ectotherms might benefit at cooler temperatures from increased cell sizes (Van Voorhies 1996, Hessen et al. 2013, Leinaas et al. 2016, selective protection from mortality or increased fecundity, all of which scale with body size (Yampolsky andScheiner 1996, Vinarski 2014). Alternatively, the pattern might emerge as a result of a confounded driver such as oxygen concentration. For example, the 'oxygen (limitation) hypothesis' suggests that the size of marine ectotherms is limited by concentrations of dissolved O 2 (Chapelle andPeck 1999, Spicer andMorley 2019). Because gas solubility and water temperature are inversely correlated, this would predict larger sizes at cooler temperatures (Forster et al. 2012, Rollinson andRowe 2018).
Food availability is another potential driver of Bergmann's rule, where more food can result in faster growth rates (Lin et al. 2013) and larger body sizes (Vidal 1980, Huston and Wolverton 2011, Andriuzzi and Wall 2018. Conversely, larger body sizes might be favoured in areas where food is scarce, because animals must forage further to find food (Belovsky 1997, Brown et al. 2017, or size might provide some other selective advantage. Additionally, latitudinal variation in diet quality could further influence size (Berrigan andCharnov 1994, Ho et al. 2010).
Latitudinal variation in predation rate is also a plausible driver of Bergmann's rule within a taxonomic group (Wallerstein and Brusca 1982, Angilletta et al. 2004, Manyak-Davis et al. 2013) because predation rate tends to decline from the equator to the poles (Freestone et al. 2011). Predation rate could affect the size and growth of communities in several ways: through evolution towards species that mature at varied sizes (Kiørboe 2011, Manyak-Davis et al. 2013; through selective predator behaviour (Kiørboe 2011); through predation-related mortality prior to maximum size (Angilletta et al. 2004); and through selective advantages of allocating energy to predator defences (Kiørboe 2011).
We focus on Bergmann's rule in marine pelagic copepods, arguably the most abundant multicellular organism on Earth (Schminke 2007). Copepods are the primary link between phytoplankton and fish in aquatic systems, and they play a central role in fisheries production (Verity and Smetacek 1996). They are crustaceans that swim weakly and thus drift in currents. Marine copepods are an ideal group for testing Bergmann's rule in ectotherms because they are widespread over diverse environments, from the poles to the equator. There are a large number of copepod species, facilitating a more robust test of Bergmann's rule. Moreover, copepods are sensitive to many plausible drivers of Bergmann's rule, including temperature (Vidal 1980, Miller andWheeler 2012), food availability (Rutherford et al. 1999, Miller andWheeler 2012), predation (Kiørboe 2011, Miller andWheeler 2012), oxygen concentration (Rollinson and Rowe 2018) and latitude (Tseng and Soleimani Pari 2019).
A search for Bergmann's rule in the literature on copepods returned no previous studies; however, Brun et al. (2016) investigated what has been called the temperature-size rule (Atkinson 1994) in marine copepods. This rule describes the plastic phenotypic response to temperature within a species, with warmer temperatures leading to smaller individuals (Diamond andKingsolver 2009, Ghosh et al. 2013). As the study actually compared the mean sizes of species and not individual sizes within a species, it effectively tested Bergmann's rule. Based on a varied dataset collected using many different nets, Brun et al. (2016) found that body size declined weakly at warmer temperatures, but there was little or no effect of temperature between 10 and 30°C. However, the temperature effect was dwarfed by the effect of food, which surprisingly led to a decline in body size with increasing food availability. This contrasts with many other studies that have shown increasing body size with food concentration in copepods (Vidal 1980) and other ectotherms (Pafilis et al. 2009, Huston and Wolverton 2011, Andriuzzi and Wall 2018. In a recent study, Evans et al. (2020) found marine copepods in the North Atlantic conform to Bergmann's Rule. Although they found that temperature (2-27°C) had a more profound relationship with body size than described in Brun et al. (2016), they did not consider the effect of food availability. Other studies have found the intraspecific version of Bergmann's rule (often called temperature-size rule or James's rule) holds true for several marine copepods species (Garzke et al. 2015, Leinaas et al. 2016. Here, we test whether food is a more important driver of body size than temperature in marine copepods, whilst considering other drivers of the potential relationships with body size, such as oxygen levels, which have not been investigated previously for copepods. We test these relationships simultaneously because they are likely partially confounded with one another, which could modify their perceived relationships with body size when tested independently. Further, we account for natural differences in size based on diet. These tests are facilitated at a nearglobal scale by virtue of the continuous plankton recorder (CPR) survey dataset, the largest (~100 000 samples), most consistent (collected using the same device), global dataset on marine copepods (Richardson et al. 2006, Batten et al. 2019. We used a spatial comparative analysis to identify statistical relationships over environmental gradients across space (Brown 1995, Gaston andBlackburn 2008).

Sample collection
We chose the CPR dataset, assembled by the Global Alliance of CPR Surveys (GACS 2011, Batten et al. 2019, because it provides the largest, consistent, most spatially-extensive, species-resolved plankton dataset (Richardson et al. 2006). Data were sourced from four surveys: the North Atlantic CPR Survey, the Scientific Committee on Antarctic Research (SCAR) Southern Ocean CPR Survey (Hosie 2020), the North Pacific CPR Survey and the Australian CPR Survey (Fig. 1).
Although CPR surveys have better coverage in polar and temperate regions than in subtropical and tropical regions, samples have been collected in waters as warm as 30°C. All surveys use similar methods to collect and count copepods (Reid et al. 2003, Richardson et al. 2006 for more details). Specifically, all CPRs use the same mesh size (270 µm), the same mesh material (silk), the same size mouth opening (1.61 cm 2 ), towed at the same depths (5-10 m) and have similar designs (Hosie et al. 2003, Reid et al. 2003. The CPR is mainly towed behind ships of opportunity on their normal trading routes, but also behind research vessels. Each CPR tow is usually up to 450 nautical miles. The internal silk roll that captures the plankton is cut up into samples representing either 5 or 10 nautical miles, and microscopic counts of copepods are converted to number per m 3 .

Copepod data
Copepods are identified to species whenever possible, and shared training and staff exchanges amongst surveys ensure comparable data. Names of all copepod taxa were updated using the World Register of Marine Species (WoRMS) (<www.marinespecies.org/>). Data were available for 388 taxa and from 97 830 samples.
To test Bergmann's rule, we calculated the mean length of copepods in each sample and related this to environmental drivers. Maximum and minimum lengths of different copepod taxa were obtained from the Marine Planktonic Copepod Database from the Observatoire Océanologique de Banyuls-sur-Mer (<https://copepodes.obs-banyuls.fr/en/ index.php>; Razouls et al. 2020) and from Richardson et al. (2006). Because juveniles are more difficult to identify to species and reliably assign a size than adults, they were not included in the estimate of mean copepod size for each sample, although adults tend to be more common in CPR samples (Richardson et al. 2006). Each taxon was assigned a single size (the midpoint between their minimum and maximum lengths). We calculated the mean length of copepods in a CPR sample by multiplying abundance of adults of each taxon in the sample by their assigned length, and then dividing their sum by the total abundance of all adults within the sample.
Bergmann's rule could also potentially be influenced by spatial differences in the trophic structure of communities because body size is linked to the ecological role of a species (Woodward et al. 2005, Yvon-Durocher et al. 2011. Because carnivorous copepods are generally larger than omnivorous ones (Mauchline 1998, Supporting information), which could affect observed relationships, we distinguished obvious differences in diets among taxa. We used a combination of dietary studies from the literature (Richardson and Schoeman 2004) and morphological differences in copepod mouthparts (Huys and Boxshall 1991) to assign each taxon to one of two diet categories: carnivore; or herbivore/omnivore (hereafter called omnivore). To calculate the proportion of omnivores we divided the summed abundance of omnivore within each sample by the total abundance of copepods within samples.

Environmental data
We used sea surface temperature (SST) as an estimate of ocean temperature for the near-surface CPR samples, and chlorophyll-a concentration (Chl-a) as a proxy for food availability to copepods. Chl-a is correlated with copepod growth and fecundity in herbivorous and omnivorous species (Richardson and Verheye 1998, Hirst and Bunker 2003, Bunker and Hirst 2004, and more Chl-a leads to more of these grazers and thus more carnivorous zooplankton (Richardson and Schoeman 2004). As an index of food availability to copepods on each sample, we considered using the phytoplankton colour index, which is a visual assessment of the greenness of the silk in four levels (no colour, very pale green, pale green and green) (Richardson et al. 2006). This index is useful for analysis when averaged in time and space (Raitsos et al. 2005(Raitsos et al. , 2014; however Chl-a from satellite is more accurate for individual samples (Richardson et al. 2006). Further, temperature is rarely measured on CPR tows, so we used remotely-sensed Chl-a and SST data to ensure measurements were consistent among samples. We matched the location and time of collection of CPR samples with estimates of SST and Chl-a averaged over eight days prior to sampling to limit loss of data due to cloud cover and to represent feeding conditions over the recent past. For SST, we used daily Group for High Resolution SST data, a cloudfree global product based on satellite, buoy and ship data, and interpolated at 0.2° × 0.2° resolution (<https://doi. org/10.5067/GHCMC-4FM02>). For surface Chl-a from September 1997 until the end of 2016, we used daily data from the European Space Agency Ocean Colour Climate Change Initiative data (<https://doi.org/10.5285/9c334fbe6 d424a708cf3c4cf0c6a53f5>). Because this is a merged product from MERIS, Aqua-MODIS, SeaWiFS and VIIRS satellites, these data provide better coverage than products from a single satellite, and are more accurate globally (Mélin et al. 2017). For 2017 and 2018, we combined data from Aqua-MODIS (<https://doi.org/10.5067/AQUA/MODIS/L3M/ CHL/2018>) and VIIRS (<https://doi.org/10.5067/NPP/ VIIRS/L3M/CHL/2018>), the only available satellites. Chl-a data were retrieved at 4-km resolution and aggregated to 0.2° × 0.2° resolution to limit loss of data due to cloud cover and match the resolution of SST data. As there are no robust satellite Chl-a data before 1997, all analyses used CPR samples collected after this time.
To investigate whether oxygen is a driver of Bergmann's rule (Forster et al. 2012, Rollinson andRowe 2018), we used data from the World Ocean Atlas 2018 (WOA-18) (<https:// doi.org/10.17616/R3JZ3J>). We used these relatively coarser-resolution data (monthly climatology averaged at 1° × 1° resolution) because we could not find an observed global oxygen product at a finer resolution. We matched the resolution of the predictors for this investigation, by using the WOA-18 1° × 1° resolution products for SST and oxygen, and by aggregating latitude to 1° × 1° resolution.
Although it is plausible that predation pressure may alter patterns in Bergmann's rule, we were unable to account robustly for this effect in our study. This is because there are no reliable estimates of predation pressure on copepods by fish or invertebrate predators, which are poorly estimated in CPR samples (Richardson et al. 2006).

Statistical analyses
To test Bergmann's rule, we used a model-building approach. We fitted a weighted generalised linear mixed-effect model (GLMM, Bates et al. 2011), with fixed effects for all predictors (SST, latitude, oxygen, Chl-a and the proportion of omnivores) (Supporting information). We included the proportion of omnivores in the model to account for omnivorous copepods generally being smaller than carnivorous ones (Supporting information). For these models, we transformed the proportion of omnivores using an arcsine-square-root transformation, and Chl-a using a square-root transformation, because these predictors were extremely skewed and thereby not adequately distributed across their range. Upon visual inspection of model residuals, their magnitude generally increased at larger fitted values. We thus used a gamma error structure with a log-link function (there were no zero values for size). Because estimates of mean length of copepods in a sample are more precise when there are more specimens in a sample, we tried several weighting schemes to account for this (including weighting points by their square-root and fourth-root). However, we found that this choice did not affect our results in any meaningful way, so we opted for an unweighted approach.
To account for temporal and spatial correlation in the data, we used four random effects associated with the time and location of sampling. The first was a random intercept for Longhurst Province, which represents a global bioregionalisation of the ocean based primarily on differences in physical oceanography (Longhurst et al. 1995;<www. marineregions.org>). This accounts for natural, historical and sampling differences among marine regions, which minimises bias associated with seasonal differences in sampling regions as well as any tendency to oversample different times of the day within regions. The second was a random intercept associated with differences amongst CPR surveys, which accounts for the use of different vessels, large-scale regional effects and any minor methodological differences. The third was a random intercept of Tow within Survey, which adjusts for both temporal and spatial differences amongst tows. These intercepts further reduce biases due to time of sampling, and due to seasonal effects (especially because tows are of relatively short duration when compared with the seasons). The last was a random slope for days elapsed on Tow within Survey, which accounts for spatiotemporal autocorrelation amongst samples on tows, which arises due to differing weather and collection conditions throughout the period of the tow. This helped remove any general linear trends over the period of the tow.
Due to the large number of samples (~100 000) and thus high statistical power, we did not assess the significance of predictors. Rather, we selected the best model using the Bayesian information criterion (BIC), which is based on the goodness of fit (log-likelihood ratio) relative to model complexity (number of parameters). BIC is suitable to fitting heuristic models with large numbers of observations; it more harshly penalises overfitting than the more commonlyused Akaike information criterion (Schwarz 1978, Aho et al. 2014. From a preliminary analysis we found that SST, latitude and oxygen were strongly correlated (all r > 0.59, Supporting information) and could not be included together in the model. We thus first assessed their relative importance in separate models with the other predictors (Chl-a, and proportion of omnivores), and then retained the most important variable amongst SST, latitude and oxygen for inclusion in subsequent analyses.
We used the pseudo R-squared described in Nakagawa and Schielzeth (2013) to estimate the proportion of variation explained by the gamma GLMMs (Nakagawa et al. 2017). To interpret the ecological relevance of the drivers, we evaluated effect sizes using expected relationships of the mean length with predictors (Nakagawa andCuthill 2007, Sullivan andFeinn 2012). We also evaluated ecological relevance by converting copepod length estimates to body mass (using wetweight (mg) = 0.03493 × length (mm) 2.9878 , Pearre Jr. 1980) because the food that is available to higher-level predators such as fish is related to body mass rather than length.
To assess the effect of influential points within the model, we performed sensitivity tests by iteratively identifying and removing groups of outliers and high-leverage observations; these procedures had no impact on overall results. Results presented thus include all available data. By mapping residuals globally both with and without random effects (Supporting information), we confirmed that the use of random effects reduced spatial autocorrelation among samples; plotting residuals through time within Surveys and on Tows confirmed that our models reduced temporal autocorrelation. The code for the statistical analyses is publicly available on GitHub (<www.github.com/maxcampb/ Bergmann-Rule-Copepods>).

Results
We found that the model containing SST instead of latitude or oxygen concentration had much lower BIC, so SST was used in all subsequent models ( Table 1).
The final model for mean copepod length included SST, Chl-a (square-root transformed) and proportion of omnivores as predictors, and random effects for Longhurst Province, Survey and Tow within Survey. All predictors reduced BIC and increased the pseudo R 2 (Table 1). Together, fixed effects explained 9.2% of the variance in mean length of copepods (Table 1). Random effects explained a further 40.5%, where 0.1% was attributed to Longhurst Province (Fig. 2a), 2.3% was attributed to Survey (Fig. 2c) and 38.3% to the Tow within Survey effects (Fig. 2b).
The mean length of copepods declined with warmer SST (Fig. 2d), where the mean length of copepods decreased by a factor of 0.982 for each 1°C increase in SST, which is equivalent to a 0.85 mm decrease across the temperature range (−1.7 to 30°C) or a 43.9% decrease in mean length. This equates to an approximate linear decrease in copepod mass of ~2.7% per °C, equivalent to an 82.2% decrease in mass across the entire temperature range. Using a linear approximation, we found for each 1°C increase in temperature, the mean length of copepods was ~0.026 mm shorter.
The mean length of copepods also declined with increased Chl-a (Fig. 2e), where for each square increase in Chl-a (based on square-root transformation), the mean length of copepods decreases by a factor of 0.899. This is equivalent to a 0.45 mm decrease across the entire Chl-a range (0.02-9.51 mg m −3 ) or a 26.9% decrease in mean length. This equates to an 60.7% decrease in copepod mass across the Chl-a range. Using a linear approximation, we estimate for each square increase in Chl-a that copepods are ~0.153 mm shorter.
The proportion of omnivores was included in the model to account for the effect of trophic role on copepod length. The mean length of copepods decreased in samples that had a higher proportion of omnivorous copepods (Fig. 2f ). The mean length was 1.65 mm or 51.8% smaller in samples comprising only omnivores when compared with samples comprising only carnivores. This equates to a 88.7% decrease in copepod mass.
Investigating the relative importance of fixed effects showed that both SST and Chl-a were of similar importance in determining the mean length of copepods, as indicated by their BIC. Removing either of these substantially degraded BIC (Table 1), while removing the proportion of omnivores within samples degraded the BIC the most (Table 1).

Discussion
We believe that our study provides a robust test of Bergmann's rule because it is the largest study to date (based on 97 830 samples and 388 taxa within an order), is near-global in its distribution, and it applies a statistical approach that allows us to disentangle multiple correlated predictors (which is a substantial challenge in comparative analyses). We found evidence that marine copepods follow Bergmann's rule, where we see a strong decline in copepod body size with temperature, Table 1.
Step 1 compares three models with the WOA-18 data, to determine whether latitude, oxygen concentration or SST is the best predictor of size (n = 90 665).
Step 2 compares five competing models to explain Bergmann's rule in marine copepods including only best temperature-dependent predictor (SST, n = 97 830). BIC was the basis for model selection, where lower values indicate a more parsimonious model. Additionally, log-likelihood represents the goodness of fit (higher is better), and degrees of freedom (df) used represents the complexity. The R 2 value reported is the marginal estimate of the psuedo-R 2 . All models have the same representation of random effects (Longhurst Province, Survey and Tow within Survey intercept and slope).
Step 1 models -WOA-18 data as reported by Evans et al. (2020). We also found evidence that is consistent with the hypothesis that temperature is more important than latitude or concentration of dissolved oxygen. Furthermore, food availability was of similar importance to temperature, contrary to the findings by Brun et al. (2016), who found that body size was much more strongly related to food availability than to temperature. Nevertheless, our results corroborate findings by Brun et al. (2016) that copepod size decreases where more food is available. This contrasts with most previous work on copepods (Vidal 1980), nematodes (Andriuzzi and Wall 2018), lizards (Pafilis et al. 2009) and mammals (Huston and Wolverton 2011), which found that body size increases with food availability.

Explanations for Bergmann's rule
Although the mechanism underlying Bergmann's rule in copepods is unclear, it is likely that the relationship with temperature is more than a spurious correlation driven by differences in food availability or differences in trophic structure (McNab 1971, Belovsky 1997, Brown et al. 2017. Even after these drivers are accounted for, a strong negative relationship between copepod size and temperature remains. It is plausible that reduced predation rates in cooler environments may alter this relationship (Wallerstein and Brusca 1982, Angilletta et al. 2004, Manyak-Davis et al. 2013, and more work is needed to uncover the relationship between size and predation pressure. We also found that dissolved oxygen concentration (which decreases with increasing temperature) does not seem to adequately account for changes in size when compared with effects of temperature, which is consistent with the hypothesis that oxygen limitation is not responsible for Bergmann's rule. Instead, it is more likely that copepod size is regulated by temperature or some other unmeasured variable confounded with temperature. A potential mechanism is the negative correlation between growth efficiency and temperature, so that colder waters could produce larger copepods (Ikeda et al. 2001, Isla et al. 2008).
Our results suggest that when the negative relationship between taxon body size and temperature is adjusted for, the relationship between taxon body size and food availability is also negative. The direction of this relationship seems counterintuitive because typically more food leads to faster growth (Lin et al. 2013) and larger size (Vidal 1980, Berrigan andCharnov 1994). A potential explanation for the negative relationship between copepod length and increasing Chl-a is that it might result from diel vertical migration and satiation. Larger copepod species vertically migrate extensively, moving between near-surface and deeper layers, whereas smaller species are unable to migrate far and spend most of their time near the surface (Hays et al. 1994). For larger species that can balance the tradeoff between greater food availability at the surface and the reduced predation pressure at depth through diel vertical migration, individuals that are satiated spend less time in surface waters (Huntley and Brooks 1982, Atkinson et al. 1992, Huggett and Richardson 2000, Hays et al. 2001). Thus, under poor food conditions (here low Chl-a), larger species will spend longer in surface waters when Chl-a is low, increasing the mean size of the community observed in the CPR, which samples in the top 10 m. By contrast, in a rich food environment (high Chl-a), copepods are more rapidly satiated and larger copepods can return to deeper waters sooner, where they are not captured by the CPR. Thus, the negative relationship between copepod length and Chl-a could be the result of a combination of copepod behaviour and how samples were collected. This would be an artefact of the study, and this hypothesis requires further testing. In our dataset this explanation is probably more likely than copepods being larger in response to having to forage further (and deeper) to find food in areas where food is scarce, as suggested in other studies (Belovsky 1997, Brown et al. 2017). An alternative explanation is that copepods might grow larger in response to seasonality of their food supply (Brun et al. 2016). For example, copepods grow larger in systems with short seasonal pulses of food (e.g. Chla) by accumulating lipids for survival when food is limited (Kattner et al. 2007). Thus, being larger and having greater reserves could allow better survival during periods without food.
It is always difficult disentangling highly correlated environmental variables with considerable spatiotemporal autocorrelation present. Although we were able substantially reduce the amount of spatiotemporal autocorrelation using random effects, we acknowledge that there is still a significant amount of spatial autocorrelation remaining in the residuals (Supporting information). Strictly speaking, this could confound the effects in our analyses, which could in term lead to selecting sub-optimal candidate models using BIC and lead to us making poor inferences about effects. However, we believe that our model is an adequate representation of reality, because our results remained relatively consistent regardless of the choices we made within the modelling and because we have a large number of samples across a vast spatial extent, so strongly confounded effects are less likely to be present.
Nevertheless, we believe we were able to estimate the relative contributions of the drivers of Bergmann's rule because of our large number of samples, and through the use of random effects to adjust for spatial and temporal variations. Using BIC allowed us to select robust models that explain the data well without overfitting (Schwarz 1978, Aho et al. 2014, which is important when analysing very large datasets. We are more confident in the relationship at temperatures typical of temperate regions (5-20°C), and at lower Chl-a, where we have most data. Unfortunately, we have a paucity of data in tropical regions, although these are the largest areas in the ocean.

Implications for climate change
Although comparative studies like ours does not necessarily mean the analogous patterns will emerge with warming through time, space-for-time substitution in ecological modelling remains a robust approach for providing insights into future changes (Blois et al. 2013). The effect of temperature seems profound when considered across its range, with a 43.9% decrease in copepod length from −1.7 to 30°C. Bergmann's rule thus suggests that as oceans warm under climate change, the size (and mass) of copepods is likely to decline (Walther et al. 2002). Under a high-emissions scenario (RCP8.5), SST is likely to warm by ~2.7°C in 2090-2099 (compared to 1990-1999) based on the mean of the Coupled Model Intercomparison Project 5 (Bopp et al. 2013), or warm by ~0.6°C under a low-emissions scenario (RCP2.6) (Bopp et al. 2013). Based on our statistical model, the effect of warming of ~2.7°C under RCP8.5 could equate to a decrease in body mass of copepods globally of ~7.3%. For RCP2.6, the ~0.6°C warming could equate to a decrease in body mass of copepods globally of ~1.6%. These estimates would translate to a similar decline in copepod biomass globally assuming abundance remains unchanged.
However, most Earth System models also project a decline in primary production and Chl-a (Bopp et al. 2013, Stock et al. 2014, Lefort et al. 2015, Galbraith et al. 2017, Woodworth-Jefcoats et al. 2017. The relationship between copepod length and Chl-a is more counter-intuitive, lesswell understood and less substantiated by other studies than the relationship between copepod length and SST. More empirical evidence is required to confirm how copepod size varies with Chl-a. Regardless, decreases in primary production of between 2% and 16% by 2100 are predicted under RCP8.5 (Lefort et al. 2015). Using net primary production estimates from Bopp et al. (2013) and the conversion to Chl-a from Marañón et al. (2014), we find that under the RCP8.5 Chl-a is projected to decrease globally by ~0.086 mg m −3 , and could lead to an increase in body mass of copepods globally by ~1.2%. Under the RCP2.6, Chl-a is projected to decrease globally by ~0.020 mg m −3 , and could lead to an increase in body mass of copepods globally by ~0.3%. Thus given that the relationship between copepod length and Chl-a is more than an artefact of our study, the combined effects of increased temperature and decreased Chl-a are likely to decrease global copepod biomass by ~6.2% under the RCP8.5, or decrease by ~1.3% under a RCP2.6. Current Earth System Models also project a future decline in zooplankton biomass (Woodworth-Jefcoats et al. 2017) -and copepods dominate zooplankton biomass (Verity andSmetacek 1996, Sommer et al. 2001) -by ~7.9% globally (Stock et al. 2014). This decline in copepod size and mass could negatively impact global fisheries (Sheridan and Bickford 2011). No Earth System Models consider the effect of Bergmann's rule on copepod size.
There could be several other important ecosystem consequences of copepod size following Bergmann's rule as the climate warms. Because swimming ability and thus the amplitude of their vertical migration is related to their size (Hays et al. 1994, Ohman andRomagnan 2016), a decline in body size with warming implies less extensive vertical migration. Thus, reductions in body size could potentially weaken the biological pump that transfers carbon from surface layers to the deep ocean (Cavan et al. 2019). Further, copepods significantly contribute to carbon exports via their sinking faeces and moults following ecdysis -at rates mostly determined by their body size (Stamieszkin et al. 2015). Therefore, reduction in copepod body size with warming could have significant ramifications for deep-ocean systems (Levin andLe Bris 2015, Sweetman et al. 2017) and for feedbacks to the climate system (Portner et al. 2019).

Final thoughts
There are now many studies that support Bergmann's rule -across marine (Saunders andTarling 2018, Wang et al. 2020) and terrestrial systems (Arnett and Gotelli 2003, Ho et al. 2010) -and across ectothermic (Olalla-Tárraga and Rodríguez 2007, Wilson 2009) and endothermic taxa (Ashton et al. 2000, Brown et al. 2017. Despite this scientific support, there remains limited uptake of Bergmann's rule -and other macroecological ideas -in modelling studies. To project changes in biodiversity, ecosystems and fisheries under climate change, a host of modelling approaches are increasingly being coupled with Earth System Models (Everett et al. 2017), including nutrient-phytoplankton-zooplankton models (Stock et al. 2014), population models (Feng et al. 2018), size-spectrum models (Carozza et al. 2019), end-toend ecosystem models (Griffith et al. 2011(Griffith et al. , 2012 and statistical models (Grieve et al. 2017). There is considerable opportunity to include well-tested macroecological principles such as Bergmann's rule in future modelling efforts focused on climate change. Our analysis shows that these principles could substantially influence future projections.
This study also highlights the utility of using large global datasets for testing macroecological theory. Datasets such as the CPR, which have been collected consistently for decades have predominantly been used to understand ecosystem dynamics or describe global change (Edwards et al. 2010). There is great potential for comparative analyses with similar consistent, global datasets. Further, with the advent and increasing accessibility of powerful statistical techniques such as GLMMs -that make it possible to test multiple predictors whilst adjusting for spatial and temporal autocorrelation -there is increasing opportunity for providing robust and nuanced tests of macroecological relationships through spatial comparative analyses (Bolker et al. 2009). We recommend that future studies appropriately account for spatial and temporal autocorrelation, and consider simultaneously testing multiple potential predictors to avoid spurious and confounded relationships, as both were common in the past.
There is still much to learn about Bergmann's rule. Future research could be directed towards testing the rule across varied taxonomic levels, detailed investigations of regional differences, and testing nonlinear relationships between size and drivers of Bergmann's rule. These studies could consider the trophic roles of species when investigating these relationships, particularly for groups with diverse roles such as copepods, because they can explain a large amount of variation in body size and thereby may confound the relationships of interest. With directed research in this area, we could get closer to understanding and resolving the many complexities of Bergmann's rule that have been debated for decades.