Assessing warming impacts on marine fishes by integrating physiology‐guided distribution projections, life‐history changes and food web dynamics

Rising ocean temperatures pose a continuing threat to marine fish communities. As warming has far‐reaching impacts at multiple ecological levels, incorporating multimodal data is necessary for more accurately forecasting the responses of species and communities to the warming ocean. Range shifts, life‐history changes, and alterations of trophic dynamics are three important aspects of warming impacts, yet there has not been a formal integration of all three aspects in the same analysis. Here, we present a novel framework that integrates species distribution projections, life‐history changes, and food web dynamics to assess warming impacts on marine fish communities. We first introduce a simple yet effective way of incorporating thermal physiological data into the species distribution model without the need to empirically measure thermal performance curves. We then use the dynamic size‐spectrum model as the modeling backbone to incorporate data from species distributions and population‐level life history analyses. With this framework, we evaluate how individual species are affected under two warming scenarios (RCP4.5 and RCP8.5). We also simulate large‐scale top‐down and bottom‐up perturbations to examine community resilience under rising temperatures. We find that warming generally reduces species biomass and shifts species size spectra towards larger individuals, even though the maximum size tends to decrease under warming. However, the exact responses to rising temperatures differ among species and do not exhibit strong correlations with species size and the pace of life history. More severe warming also renders the focal community more vulnerable to top‐down perturbations, even though the community remains sufficiently resilient overall. The complex nature of species and community responses result from the fact that distribution range, life history, and food web dynamics change with warming in different directions that may not be intuitive to predict a priori. Importantly, we show that neglecting changes in species distribution or life history will lead to biased assessment of species and community responses. Our analyses highlight trophic dynamics, species biomass, and community resilience as three emergent properties that our framework can uniquely quantify. This integrative framework is readily applicable to other communities of interest and can be scaled up for multi‐regional or global analyses.


| INTRODUC TI ON
Anthropogenic carbon emissions have caused sea temperatures to rise globally, and this alarming trend is projected to persist to the end of the century and beyond (Bindoff et al., 2019). Marine ecosystems, especially fish communities, have been essential in providing resources and services and are under threat from the disrupted thermal environment. Assessing warming impacts on marine fish communities therefore has become a pressing issue at the intersection of ecology, conversation and fisheries. Changes in sea temperature have far-reaching consequences at all ecological levels (Hollowed et al., 2013). At the individual level, warming can alter physiological processes (e.g. metabolic rate) and whole-organism performance (e.g. locomotor speed and endurance; Schulte et al., 2013). At the population level, rising temperatures can cause significant shifts in life histories (van Gils et al., 2016;Wang et al., 2020), size structure and abundance (Blanchard et al., 2012), and spatial distribution (Gervais et al., 2021). Beyond individual species, warming can affect species interactions, most notably trophic dynamics (Reum et al., 2020), and ecosystem functions (Dossena et al., 2012). These effects can interact with one another and propagate to cumulatively impact community composition and resilience (Doney et al., 2012).
Since the effects of warming are all-encompassing, assessing how species and communities might respond to rising temperatures requires an approach that integrates information and processes from the individual to community levels.
Life-history changes, range shifts and alterations of trophic dynamics are three important axes for assessing the effects of warming on species and communities (Hollowed et al., 2013). To these ends, statistical life-history models, species distribution models and size-based food web models have been developed as useful analytical tools and are widely implemented for multiple taxonomic groups.
Over the past few decades, there have been key advances in more rigorously assessing warming effects on these three fronts. Recent studies have demonstrated the utility of using trait variation across spatial scales to forecast temporal shifts in life histories under future climates (e.g. Audzijonyte et al., 2020;Wang et al., 2020). At the same time, much progress has occurred to allow for incorporating data encompassing multiple ecological scales and processes for generating warming-related predictions and hypotheses testing (Blanchard et al., 2012Clarke et al., 2021;Peterson et al., 2019;Stuart-Smith et al., 2017;Talluto et al., 2016;Wabnitz et al., 2018). Despite significant progress on each front, there has not been a formal attempt to integrate information from species distribution models, statistical life-history analyses and food web models within the same analytical framework. As each of the three approaches captures a unique facet of warming impacts, evaluating the effects of rising temperatures on species and communities will be incomplete and biased if insights from all three approaches are not taken into account simultaneously.
Here, we present an integrative framework that uses dynamic size-spectrum models as the backbone to unite species distribution projections and population-level life-history data ( Figure 1).
Additionally, we introduce a simple yet effective approach for approximating thermal tolerance curves for species, which enhances the performance of distribution projections. This approximation does not require an empirically measured relationship between temperature and fitness, which is challenging, if not impractical, to collect for most marine fishes. Instead, the only physiological variable required is the maximum or minimum critical temperature (CT max or CT min , respectively), which has been estimated for a large number of species (e.g. Bennett et al., 2018;Comte & Olden, 2017). In its full capacity, our framework allows the integration of multimodal data from the individual, species and community levels, including thermal tolerance, species distribution projections, life-history, diet and distribution-informed trophic interactions.
Based on projected RCP4.5 and RCP8.5 scenarios, the objectives of this study are (a) to examine changes in species size spectra (the distribution of biomass density with body size), abundance, and mean body size, and (b) to assess changes in community size spectra and evaluate the community's ability to recover from largescale disturbances. Specifically, we expect that warming impacts would be more severe in species with faster life histories (i.e. smaller species) and inhabiting higher latitudes (temperate species). We also expect that warming would render communities more vulnerable to large disturbances, with more severe warming having more pronounced impacts. We test these predictions with the fish community inhabiting the continental shelf of eastern United States we show that neglecting changes in species distribution or life history will lead to biased assessment of species and community responses. Our analyses highlight trophic dynamics, species biomass, and community resilience as three emergent properties that our framework can uniquely quantify. This integrative framework is readily applicable to other communities of interest and can be scaled up for multi-regional or global analyses.

K E Y W O R D S
Bayesian species distribution model, dynamic size-spectrum model, resilience, thermal performance curve, trophic interaction (25°-45°N, 65°-80°W). We choose this exclusive economic zone for the availability of long-term species occurrence data (~50 years) and its history of human exploitation and economical importance.
Nevertheless, our analytical framework can be readily applied to any marine communities of interest, provided that sufficient data are available.

| MATERIAL S AND ME THODS
We integrated physiology-guided species distribution models (hereafter physiology SDMs), life-history data and dynamic size-spectrum models to explore how species size spectra, community structure and community resilience might change in a warming climate.
Because of this overarching focus, we used temperature as the sole environmental factor in our analyses. Even though other factors, such as net primary production, are also crucial in mediating climate impacts on communities, temperature has a unique importance in that it directly influences the potential distributions and life histories of species (Figure 1). Our method therefore offers a data-driven approach to model temperature dependence of physiological and ecological processes within an ecosystem.

| Species occurrences and temperature data
We obtained species occurrence data in our study area from the Northeast Fisheries Science Center (NEFSC)'s bottom trawl surveys from 1948 to 2008, with a total of 872 recorded fish species (Families Actinopterygii and Elasmobranchii) (Depres & Benson, 2020). In all, 37 species with at least 1,000 occurrence records and sufficient population-level life-history data (see Section 2.4) were used for subsequent SDMs. These species were present in at least 36 out of 50 survey years. Collection depth in the raw survey data was bracketed by a minimum and a maximum depth, we therefore estimated the collection depth (in meters) of each occurrence record as the mean of the minimum and maximum depth associated with the recorded collection.
We then converted the geographical coordinates of each record into 1° × 1° cells for a program by Webb et al. (2020), which matched the cell coordinate, collection depth and collection time (month within year) to the global gridded temperature dataset from the Institute of Atmospheric Physics (Cheng et al., 2017). We discarded records without matched temperature due to missing data or spurious depth values. We chose this temperature dataset over others with finer geographical resolution for its depth-calibration of sea temperature.
To assess the effects of warming, we used the year 2013 as the baseline for comparison (referred to as 'current climate' hereafter).
We obtained the annual maximum of sea surface and sea bottom temperatures (maxSST and maxSBT, respectively) in 2013 from the NOAA World Ocean Atlas 2013, as well as those under the RCP4.5 and RCP8.5 warming scenarios during the years 2090-2100 from the Bio-ORACLE models (Assis et al., 2017). RCP4.5 represents a scenario in which the greenhouse gas concentrations stabilize, whereas RCP8.5 represents a scenario of increasing emissions over time. The geographical resolution of these datasets was 0.25° × 0.25°, resulting in a total of 677 cells within the study area. We used temperature maxima instead of time-averaged temperatures when evaluating F I G U R E 1 The analytical framework of this study, illustrating the integration of species distribution data and life-history data into dynamic size-spectrum model. Variables and model outputs in red are those directly affected by temperature. The arrows indicate how outputs from species distribution model and life-history model serve as input variables in the sizespectrum model species presence because we wanted to take a more conservative approach in testing the hypothesis regarding species distributions.
Using maximal temperatures that represent the most challenging condition for species, we would be more confident in trusting the predicted presences from the SDMs.

| Estimating species' thermal performance curves
A thermal performance curve describes the relationship between temperature and an organism's relative fitness (Angilletta, 2006;Huey & Kingsolver, 1989;Huey & Stevenson, 1979). Although the curve should be left-skewed (Bennett & Huey, 1990;Martin & Huey, 2008), in practice its shape is often approximated using a symmetric Gaussian function, which has been shown to strike the best balance between goodness-of-fit to the data and model complexity, requiring only three parameters to define (Angilletta, 2006).
We approximated the thermal performance curve for a species with the following information: the maximum critical temperature (CT max ), the relative fitness at CT max and the temperature at which its fitness is maximal (optimal temperature, T opt ). We used thermal affinity as a proxy for its T opt (Audzijonyte et al., 2020;Stuart-Smith et al., 2015. For each species, we first examined the distribution of temperature associated with each observation. We then removed observations outside the 2.5 and 97.5 percentiles and used the midpoint of the truncated distribution to represent the species' thermal affinity. CT max data were available for 15 species from Bennett et al. (2018) and Comte and Olden (2017). Since there is no theoretical or empirical consensus regarding the value of relative fitness at CT max , we set the relative fitness at CT max to be 0.05 (5% of the fitness at T opt ).

| Physiology-guided species distribution models
We constructed SDMs using a Bayesian hierarchical framework recently developed by Talluto et al. (2016). This new method utilized Bayesian generalized linear models (GLMs) with environmental variables as predictors and presence-absence data as the binary response variable to produce probabilistic estimates of species presence (i.e. suitability) under different environmental conditions (the normal SDM). Since our data were presence-only, we created pseudoabsence data for a species as the presence of all other species (regarding both location and depth) that did not share the same cell-depth combination with any presence records of the focal species. We considered depth when defining pseudoabsences as it is an important dimension along which temperature and ecology vary for fishes.
For a thoroughly surveyed area, this approach was effective in inferring pseudoabsences while taking into account both niche differences and spatial accessibility. If the thermal performance curve was available for a species, the SDM could produce a physiology-guided model (the physiology SDM). This method represented a significant breakthrough in implementing correlative SDMs with species physiology, offering the potential for more accurately predicting species distributions (Gamliel et al., 2020). We constructed both the normal SDM and the physiology SDM for the 15 species with CT max data.
We also tested whether the physiology SDMs were sensitive to the use of thermal affinity as T opt by constructing a separate set of SDMs whose T opt were set to be 30% higher than thermal affinity. The SDM results were qualitatively similar with respect to the actual T opt values (Table S1), indicating that the SDMs were robust to potential difference between thermal affinity and T opt . For the other species, only the normal SDM was constructed. As is common for large-scale SDMs, we did not consider species interactions and anthropogenic interferences such as fishing at this point of the analyses. However, these factors played a role in producing the presence-absence patterns in the data, and it was difficult to isolate their effects at this stage of the analyses. We explicitly incorporated these important factors in the size-spectrum models instead.
When constructing the normal SDMs, we used temperature and its square term as predictors in the Bayesian GLM to model the unimodal, nonlinear relationship between temperature and habitat suitability (Gamliel et al., 2020;Talluto et al., 2016). For pelagic and benthic species, we used sea surface and sea bottom temperatures, respectively. For benthopelagic species, we used the average between sea surface and sea bottom temperatures. We used uninformative Gaussian priors as in Talluto et al. (2016) for the regression coefficients (μ = 0, = 10,000). To estimate the posterior distribution of these coefficients, we used an MCMC algorithm with an initial 1,000 tuning samples, a burn-in period of 4,000 samples, and an additional 4,000 samples, from which the result of Bayesian GLM was derived. We discerned the convergence of the MCMC using the Gelman index; values lower than 1.2 were considered convergence (Carpenter et al., 2017;Gelman & Rubin, 1992). For species with CT max data, we used relative fitness directly as a measure of suitability and performed Bayesian beta regression over the range of 0-40°C (the physiology sub-model), using vague Gaussian priors (μ = 0, = 1,000). Posterior distributions of the regression coefficients in the physiology sub-model were then used as priors for the same coefficients in Bayesian GLMs when constructing the physiology SDMs (Talluto et al., 2016). The weight of the physiology submodel was set to be equal to the normal SDMs, although the results were robust to variation in this parameter (Table S1). Other settings were consistent with the normal SDMs.
We used a 70% random sample from the species' presencepseudoabsence data for model training, and the remaining 30% of the data for testing model performance. For species whose presence-pseudoabsence ratio was below 0.2, we randomly chose a portion of the pseudoabsence data such that the presencepseudoabsence ratio was 0.2 in both the training and testing datasets. We performed cross-validation 20 times per species and used the Boyce index to assess model performance (Hirzel et al., 2006).
Boyce index has values ranging from −1 to 1, with higher values indicating better SDM performance. We also calculated k max values from each round of cross validation, which is a threshold for converting the continuous suitability values into binary presence-absence outcomes (Hirzel et al., 2006). Habitats with a lower suitability than k max would be considered unsuitable for the focal species (i.e. predicted absence). Overall, our models achieved high predictive accuracies (Table S2).
We used all presence-pseudoabsence data under 0.25° × 0.25° spatial resolution to better capture distributional characteristics of species under current and future climates. We recognized that the spatial resolution was finer here compared to the dataset used to construct SDMs. Since the temperature-suitability relationship would not vary with spatial resolution, we believed this would not influence presence predictions from the SDMs. For each species, we used the average of 20 k max values from the SDM (normal or physiology) with better performance as the threshold for translating continuous suitability values to presence or absence.

| Statistical life-history model
We focused on three life-history traits: asymptotic weight, weight at maturity and von Bertalanffy growth coefficient, as these are important diagnostic traits of a species' life history and key parameters for the size spectrum model. As length was more often used to represent size than weight for most species, we began by compiling population-level data on asymptotic length, length at maturity, von Bertalanffy coefficient and mean environmental temperature from FishBase and references therein (www. fishb ase.se). We supplemented data from FishBase with Google Scholar literature search, using a combination of species names, body size, size at maturation and von Bertalanffy as keywords.
For each population, we used the mean between males and females whenever data from both sexes were available. We estimated length at maturity as 2/3 asymptotic body length when the former was not directly available (Froese & Binohlan, 2000).
When there were fewer than three population-level observations, we broadened the search to include data from other species in the same genus (i.e. congeners), regardless of their range of distribution. Consequently, we excluded 24 species from further analyses due to insufficient data. We then converted weight from length using equations from (Barnes et al., 2008). For species whose length-weight equations were not available, we used equations at the genus or family levels. If even such equations were not available, we converted weight from length with a generic equation: weight = 0.02 × length 3 (Barnes et al., 2008). Overall, we obtained data for 37 species from 444 populations. In all, 17 of these 37 species included population-level data from congeners. For each species, we first calculated the mean temperature from all the populations. We then defined temperature anomaly as the difference in temperature between each population from the population-wide mean temperature . We used temperature anomaly of each population as the fixed-effect predictor in the statistical models.
To quantify temperature effect on life-history traits, we constructed separate linear models using the lm or lmer functions in the r package lme4 (Bates & Mächler, 2014), with maximum mass, mass at maturity or von Bertalanffy coefficient as the response variable.
We constructed four linear models with temperature anomaly as a fixed factor and species as a random factor. Four different models were constructed, representing different hypotheses regarding the slope and the intercept between life-history traits and temperature anomaly between species (Table S3). The Akaike information criterion (AIC) indicated that the model with species-specific slopes and intercepts was the best for all three life-history traits (Table S3). We obtained species-specific slopes and intercepts between the three life-history traits and temperature anomaly from the three lmer linear models, representing the magnitude of life-history changes with temperature for each species. To quantify life-history changes between current and future climates, we averaged the temperatures across cells in which a species was present under different climate scenarios (current, RCP4.5 and RCP8.5). We then inferred changes in life-history traits under warmer temperatures using the speciesspecific equations from the statistical models, given the difference between the temperature, averaged over suitable habitats, under current and future climates.

| Dynamic size-spectrum model
Size-spectrum models are a type of physiologically structured models in which vital ecological, physiological and life-history processes of individuals are described in relation to its size (Andersen & Beyer, 2006;Hartvig et al., 2011, see also figure 1 in Scott et al., 2014). These individual-based processes can be scaled to simulate the transference and circulation of biomass through the food web to generate steady-state size spectra, which are distributions of biomass density against size (i.e. weight). We used the mizer package in r to construct size-spectrum models . mizer can build size-spectrum models of different complexity, ranging from single-species models, to simple community models where individuals are only described by their size without species-specific information, to multispecies models where different sets of parameters can be set for each species in the community. For this study, we used the multispecies version of the mizer model. In all, 19 species that had reliable SDM results and sufficient spatial overlap with other species in the community were used in the dynamic size-spectrum model (Table S4). The standard model building and calibration procedures in mizer, as well as associated data collection, were described in more detail in the Supplementary Methods. We followed the same procedure when constructing models for the current and future climates.
However, several aspects of model input were different for future models, which we explain below.
To construct dynamic size-spectrum models under future climates, we input new life-history variables based on results from the statistical life-history models. We note here that we were comparing time-averaged conditions of past and future instead of a projected time series. Dietary preference parameters ( and ) did not vary systematically with temperature and therefore had the same values under all climate scenarios (Barnes et al., 2010). Spatial overlaps in the interaction matrices were re-calculated based on predicted dis- Here we followed Pimm (1984), who defined resilience as the rate at which a system returns to equilibrium following disturbance. To estimate community resilience under warming, we simulated both a topdown and a bottom-up perturbation as in Blanchard et al. (2010).
The top-down perturbation was as a 50% mortality on individuals larger than 10 g; the bottom-up perturbation was a 50% increase in the density of background resources. Each perturbation event persisted for 1.2 months. We ran the simulations for 300 years until the community reached equilibrium, as confirmed with augmented Dickey-Fuller tests (tseries package in r). We then evaluated the time needed for community biomass to reach steady-state (return time in months, T R ) following Neubert and Caswell (1997): where TB(t) was the percent change in the total biomass density relative to the steady-state value prior to perturbations (TB 0 , which was set to 1). The integration was carried out for the entire time frame of the analysis (t from 0 to t max ).

| Assessing the importance of warminginduced changes in species distribution and life history in predicting community responses
We first used a model selection approach to assess the contributions of species distribution and life history in gauging the responses of species biomass to warming. We constructed the full model with changes in asymptotic weight, von Bertalanffy coefficient and habitat availability predicted from the SDMs as predictors. We then eliminated the predictors one at a time and observed changes in the AIC scores. If eliminating a predictor resulted in a significantly lower AIC score, it would mean that the predictor was not essential in predicting changes in species biomass under warming.
To further demonstrate the importance of incorporating both changes in species distribution and life history in the communitylevel analyses, we constructed two alternative models under the RCP4.5 and RCP8.5 scenarios. In one alternative model (M no.inter ), we built and calibrated the model in the same way as described previously, except that we used the current interaction matrices in all future scenarios. In the other alternative model (M no.lifehis ), life-history traits remained unchanged between current and future climates. We compared the community size spectra between the original model and the two alternative models to quantify the importance of considering changes in trophic interactions and life history in affecting model outcomes. We also compared community resilience between the original and the two alternative models to evaluate the potential biases in estimating community responses to warming.

| Warming caused shifts in species distribution and altered trophic interactions
The Bayesian SDM performed well in predicting presences during cross-validation; the Boyce index and k max values were consistent among iterations, indicating the robustness of the SDM results (Table S2). The physiology SDM outperformed the normal SDM in 9 out of 15 species for which both SDMs were constructed. Notably, approximated thermal performance curves enabled the construction of SDMs for two species for which parameters of normal SDMs could not be estimated due to insufficient observations (Anchoa mitchilli and Cynoscion regalis). A previous study using the same Bayesian hierarchical models had shown that physiology SDMs built with empirically measured thermal tolerance curves had an advantage over the normal SDM in predicting species presence (Gamliel et al., 2020). Our results further demonstrated that even with much limited information (only critical temperature and thermal affinity), the approximated thermal performance curves can substantially enhance physiology SDM construction and performance.
Under warming, 13 temperate species shrunk in range due to lower-latitude habitats becoming unsuitable (Figures S1-S3; Table   S5). In contrast, three temperate and two subtropical species expanded in range; one subtropical species (Alosa sapidissima) was even predicted to be present throughout the study area in all climate scenarios ( Figures S1-S3). Shifts in range also affected spatial overlap (1) |TB (t)| 2 dt, from 0.06 in 2013 (SD = 0.04) to 0.13 under RCP4.5 (SD = 0.05) and 0.17 under RCP8.5 (SD = 0.05). We are aware that the correlation between patterns of spatial overlap and patterns of ecological interactions is often not straightforward in reality (Blanchet et al., 2020).

However, shifts in spatial overlap between species under different climate scenarios does change trophic interactions in our models
and has important consequences for model outcomes (see below).

| Warming accelerates life history, but the magnitude of trait responses varies between species
Across species, the asymptotic weight and weight at maturity decreased with temperature, while the von Bertalanffy growth coefficient increased as temperatures rose (dashed lines in Figure 2), indicating less time was required for species to reach their maximal size. This is consistent with the general finding that warming speeds up life history, making individuals growing faster yet to a smaller maximal size . At the same time, trait responses to temperature anomalies exhibited remarkable species-specific variation ( Figure 2). For the von Bertalanffy growth coefficient, part of the variation could be attributed to size (i.e. asymptotic weight), as growth increased more with temperature in larger species (multiple R 2 = 0.18, p = 0.07). On the other hand, the magnitude of temperature dependence for asymptotic weight did not significantly covary with size (multiple R 2 = 0.002, p = 0.86).

| Warming impacted species biomass and altered species size spectra
Warming reduced overall biomass of the community, with more severe warming producing more pronounced impacts. The response of biomass to warming was more subtle and varied at the species level ( Figure 3). In 10 species, warming substantially reduced biomass; in eight other species, biomass increased under warming. Notably, biomass of three species (Glyptocephalus cynoglossus, Melanogrammus aeglefinus and Rostroraja eglanteria) responded to two warming scenarios in different directions; it increased under RCP4.5 but decreased under RCP8.5 (Figure 3). Rising temperatures also substantially altered size spectra in each species (Figure 3). One common outcome was an increase in relative biomass density for larger individuals, causing an upward shift in species size spectra ( Figure   S4). As informed by life-history data, the response of maximal size to warming differed between species, producing species-specific patterns; for most species, the maximal size decreased under warming, resulting in more compressed size spectra. Maximal size for several species, on the other hand, did not change drastically under warming-the shape of their size spectra was therefore relatively unaltered.
Warming also impacted species abundance and mean body size. The majority of species experienced a reduction in abundance and an increase in mean body size as temperatures rose. Changes in spatial overlap, species size spectra and species abundance due to warming particularly impacted larger individuals through compounding effects on trophic interactions. As temperatures rose, larger individuals in most species fed at a lower rate compared to the maximum prey intake given their body sizes and prey abundance ( Figures S5-S7). The decrease in fish biomass across the community also meant that individuals that fed on other species had to rely more on the background resources to satisfy their dietary needs ( Figures S8-S10). Examining warming impacts along the axis of habitat type, the three pelagic species (orange circles in Figure 4) were impacted to a lesser degree by warming compared to other species of similar sizes; changes in their abundance and mean body size were F I G U R E 2 Changes in asymptotic weight (a), weight at maturity (b) and the von Bertalanffy growth coefficient (B) with temperature anomaly (deviation from the mean temperature of a species, 0 on the x-axis). Each solid line represents a trait-temperature anomaly regression. Grey lines are demersal species, orange lines are pelagic species and purple lines are benthopelagic species. The dotted line is the community-wide regression. Note the extensive variation between species. Life-history data were from FishBase (www.fishb ase.se) and literature search  Figure 4). Since warming impact was also minor in a subtropical demersal species (Rostroraja eglanteria; open black circle in Figure 4) and a temperate pelagic species (Alosa pseudoharengus; filled orange triangle in Figure 4), it was likely that pelagic and subtropical species in our focal community fared better under warming. Benthopelagic species (purple circles in Figure 4), on the contrary, were among those whose abundance and mean body size were affected the most by warming.

| Ignoring changes in trophic interactions and life history leads to biased estimates of community structure and resilience
For both the RCP4.5 and RCP8.5 scenarios, the model containing all life-history and species-distribution predictors was either the single best model or among equally optimal models based on the AIC scores (Table S6). For this analysis, we excluded weight at maturity as a predictor from all linear models to avoid multicollinearity, as it was highly correlated with the asymptotic weight (R 2 > 0.86 in all climate scenarios). Indeed, changes in asymptotic weight, von Bertalanffy coefficient and habitat availability were all poor predictors of biomass changes by themselves, with the exception of changes in habitat availability under RCP8.5 ( Figure 5). Not considering changes in life history and spatial overlap between species also biased the predicted species size spectra, even though the severity of such biases varied between species. In most cases, the biases were between twofold to 10-fold, but they could be as high as two orders of magnitude in the most extreme cases (Figures S11 and   S12). The biases in the estimated community size spectra resulted in different return time of the community after large-scale disturbances, except for bottom-up perturbations under RCP85 (Table 1).
Overall, the deviations were greater for RCP4.5 than for RCP8.5, and greater for top-down than for bottom-up disturbances. These results highlighted the importance of accounting for shifts in trophic interactions and life-history traits simultaneously when assessing community responses to warming.

| DISCUSS ION
In this study, we integrated information on physiology-guided species distributions, life-history changes and food web dynamics to F I G U R E 4 Relative changes in abundance (a, b) and mean body size (c, d) for each species under RCP4.5 and RCP8.5, compared to 2013. Colouring of the circles are the same as in Figure 2.  assess the response of our focal community to warming. We demonstrated the merit of our integrative approach by showing that examining these aspects in isolation would produce partial, sometimes even misleading, assessments. Several key points emerged from our analyses.

| Warming effects are largely dominated by species-specific responses
Species distribution models can inform warming impacts on the spatial matches between species and the environment, as well as the connectivity among species (Hollowed et al., 2013). Under warming, species often exhibit range shifts due to niche-tracking, either to avoid local extinction or to occupy newly available habitats (Román-Palacios & Wiens, 2020). Our SDM results conformed to the general observations that temperate species tended to experience more severe habitat shrinkage while subtropical species expanded in range (Pinsky et al., 2019). We note, however, that our SDM results did not quantify the full extent of species range shifts. A regionalscale study that covers the leading and trailing edges of species distribution would be required for this purpose. The SDM results also allowed us to consider the degree of warming for each species in a more realistic manner. Instead of assigning the same amount of temperature increase to all species (e.g. 2°C), we estimated the extent of warming for each species based on projected presence.
Because of this, the temperatures that some species experienced did not necessarily increase more under the more severe warming scenario. In Glytpocephalus cynoglossus, for example, mean temperature across occupied habitats was actually higher under RCP4.5 than in RCP8.5. This seemingly paradoxical result was due to the fact that under RCP4.5 there were more habitats whose temperatures were near the upper inhabitable threshold of G. cynoglossus. Under more severe warming (RCP8.5), these habitats became uninhabitable, resulting in a smaller number of occupied habitats with a lower mean temperature.
The responses of key life-history traits to warming captured changes of species phenology in the future. Due to this importance, several studies have examined warming-related changes in life-history traits, most commonly body size (Radchuk et al., 2019;Teplitsky & Millien, 2014;Zhang et al., 2015). However, there was often no consistent trend across taxa from these studies. More recent studies that focused on life-history responses to warming in marine systems also discovered that, even though warming often sped up life history, taxon-and lineage-specific responses to warming tended to be more dominant (Audzijonyte et al., 2020;Wang et al., 2020). Our analysis of warming effects on fish life history was in agreement with this general finding-the direction and magnitude of warming-induced life-history changes differed markedly among species (Figure 2).

| Warming responses of trophic dynamics, species biomass and community resilience are three emergent properties that may not be intuitively predicted
Numerous studies have examined warming effects on communitylevel trophic dynamics, yet there was no universal trend between the degree of impact and trophic position. Kirby and Beaugrand (2009) found that warming enhanced bottom-up control in the North Sea, thereby more severely impacting top predators. By contrast, Kratina et al. (2012) found that warming intensified top-down control in freshwater communities. As in life-history traits, the exact consequences of warming on trophic dynamics appeared to be taxon specific and system specific, and it is not uncommon that warming led to an overall neutral effect at the community level (Kratina et al., 2012;Nelson et al., 2020). The interactions between bottom-up and topdown controls in complex food webs under warming have thus been considered an emergent property that may not be intuitively predicted a priori (Lynam et al., 2017). In our analyses, warming can alter trophic interactions through two mechanisms. First, changes in species size distributions affected prey suitability to predators.
Second, range shifts due to warming changed the encounter probabilities among species. Considering these two mechanisms, warming might be more impactful to higher level consumers feeding on other fish. Indeed, dietary shifts were more pronounced in higher level consumers in our focal community (e.g. Gadus morhua, Lophius

americanus, Merluccius bilinearis, Squalus acanthias, Urophycis chuss
and Zoarces americanus), which consumed less fish as temperatures rose ( Figures S8-S10). However, our model only accounted for a decrease in background resource level under warming without explicitly modelling changes in size distribution, taxonomic composition and spatial distribution of the underlying planktons. Therefore, the impacts on species lower in the food chain could be underestimated in our analyses.
However, how such reduction might relate to characteristics of species themselves, such as body size and habitat type, is not as frequently tested. We predicted that small, temperate species would be impacted more by warming. Contrary to our prediction, there was no clear size dependency of warming impacts, as far as total biomass and abundance were concerned. Smaller species with faster life histories did not suffer more reduction in biomass or abundance compared to larger species (Figures 3 and 4). Overall, the magnitude of warming effects varied more among species than it did along the size axis, indicating that species-level variation played a more important role. The response of a species' biomass and abundance to warming therefore could not be readily predicted from its life history alone, nor could it be gauged simply from the extent of habitat shrinkage.
The response of biomass to warming is therefore also an emergent property from our analytical framework. This emergent property is crucial in assessing community-level responses to warming and in evaluating the community's resilience to perturbations and its ecosystem functions. The highly species-specific responses to warming have led to another emergent property at the community level-the resilience to larger perturbations. Our analyses have shown that community resilience depended on type of perturbation, and that more severe warming did not necessarily render the community more vulnerable.

| Utility of the framework in conservation and management
The size-spectrum models have proven a valuable tool for conservation and management in fisheries science (Canales et al., 2020;Datta & Blanchard, 2016;Szuwalski et al., 2017;Wo et al., 2020;Zhang et al., 2016Zhang et al., , 2018. Our framework can further enhance the potential of size-spectrum models for more applied inquiries. In our models,

| Model limitations and future directions
Having demonstrated the utility of our integrative analyses, here we list several limitations to our framework and point out directions for further research. First, the SDMs could not differentiate between range shifts due to thermal niche tracking and those due to changes in underlying resources. The impact of changes in local resources was potentially taken into account simultaneously by the SDMs and the size-spectrum model. In this regard, the SDMs and the size-spectrum models may not be truly independent from each other. In our analysis, we used only temperature when constructing SDMs, which should have partially alleviated the issue of redundant information. We caution future studies that use our framework to be thoughtful about including information on local resource levels (e.g. primary production) when constructing SDMs for the analyses.
Second, out statistical life-history models assumed that temperature was the only factor affecting life-history trajectories of species. This is obviously an over-simplification. Anthropogenic factors, most notably fishing-related selection, have been shown to drive lifehistory evolution in marine species (Kuparinen et al., 2016). However, our framework can readily accommodate these factors by including them as predictors in the life-history analyses to better predict lifehistory changes as both temperature and fishing pressure change in the future. In addition, performing statistical analyses on life-history traits required sufficient population-level data. When such data are not available, we recommend two alternative approaches. The first one is to estimate key life-history traits for which data are lacking using equations developed by Thorson et al. (2017). Another alternative is to use temperature-dependent scaling in metabolic rate to infer changes in life-history traits. The package therMizer, which im- Third, changes in suitable habitats from the SDMs could potentially constrain the output from size-spectrum models, as we used this information to produce a first-pass, ball-park biomass predictions under future climates for mizer. We did not detect this issue for biomass predictions under RCP4.5, as changes in suitable habitats was a poor predictor of species biomass. This constraint could be more prominent under more extreme warming scenarios (RCP8.5), when species shifted more in range. Our biomass predictions can be We advocate our analytical framework as an innovative tool for evaluating warming responses in other communities. When applied to a multi-regional or global dataset, our framework can test large-scale hypotheses regarding warming impacts on species. For example, how would warming effects vary among fishes inhabiting different places in the water column? How might warming differentially affect fishes occupying different ecological niches and places in the phylogenetic tree? Questions like these represent exciting and fruitful avenues for future research.

ACK N OWLED G EM ENT
We would like to thank Julia Blanchard and the size-spectrum read-

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.