Nature Climate Change | Volume 12 | November 2022 | 1037–1044 1044 Article https://doi.org/10.1038/s41558-022-01490-7 34. Sunday, J. M. et al. Thermal-safety margins and the necessity of thermoregulatory behavior across latitude and elevation. Proc. Natl Acad. Sci. USA 111, 5610–5615 (2014). 35. Pincebourde, S. & Casas, J. Narrow safety margin in the phyllosphere during thermal extremes. Proc. Natl Acad. Sci. USA 116, 5588–5596 (2019). 36. Somero, G. N. The physiology of climate change: how potentials for acclimatization and genetic adaptation will determine ‘winners’ and ‘losers’. J. Exp. Biol. 213, 912–920 (2010). 37. Johansson, F., Orizaola, G. & Nilsson-Örtman, V. Temperate insects with narrow seasonal activity periods can be as vulnerable to climate change as tropical insect species. Sci. Rep. 10, 8822 (2020). 38. Sinclair, B. J. et al. Can we predict ectotherm responses to climate change using thermal performance curves and body temperatures? Ecol. Lett. 19, 1372–1385 (2016). 39. Davis, A. J., Jenkinson, L. S., Lawton, J. H., Shorrocks, B. & Wood, S. Making mistakes when predicting shifts in species range in response to global warming. Nature 391, 783–786 (1998). 40. Suttle, K. B., Thomsen, M. A. & Power, M. E. Species interactions reverse grassland responses to changing climate. Science 315, 640–642 (2007). 41. Gouhier, T. C., Guichard, F. & Menge, B. A. Ecological processes can synchronize marine population dynamics over continental scales. Proc. Natl Acad. Sci. USA 107, 8281–8286 (2010). 42. Harley, C. D. G. Climate change, keystone predation, and biodiversity loss. Science 334, 1124–1127 (2011). Publisher’s note Springer Nature remains neutral with
regard to jurisdictional claims in published maps and
institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/. © This is a U.S. Government work and not under copyright protection in the US; foreign copyright protection may apply 2022, corrected publication 2023

Nature Climate Change Article https://doi.org/10.1038/s41558-022-01490-7 Methods CMIP6 simulations We obtained CMIP6 climate simulations for the historical forcing period (1850–2014) and future emissions scenario SSP5-8.5 (2015– 2100) via the CMIP6 data portal (https://esgf-node.llnl.gov/search/ cmip6/). Eight models from CMIP6 (AWI-CM-1-1-MR, BCC-CSM2-MR, CESM2, EC-Earth3, INM-CM5-0, MPI-ESM1-2-HR, MRI-ESM2-0 and NorESM2-MM) were selected on the basis of availability of daily air temperature at surface (‘tas’) at a 100 km nominal resolution at the time of download. While ‘tas’ at sub-daily frequencies is available for some models, daily data was selected to maximize the ensemble size. We resampled all datasets to a common 1° by 1° grid spanning −90° to 90° latitude and 0° to 360° longitude, and to a standard cal - endar without leap years. Spatial regions were defined on the basis of latitude as Northern Hemisphere Extra-tropics (90° S to 30° S), Tropics (30° S to 30° N) and Southern Hemisphere Extra-tropics (30° N to 90° N). Statistical analyses of climate data Quantile regression. Trends in the percentile values of global and regional temperature distributions were computed via quantile regression. Quantile regression can comprehensively model hetero - geneous conditional distributions, where the relationship between the quantiles of the dependent variable and the independent variable is different from the relationship between the means of the depend - ent variable and the independent variable. We applied quantile regression to analyse trends with respect to time at various percentile values (P 2.5, P10, P20, P30, P40, P50, P60, P70, P80, P90, P97.5). Analyses were performed using the R package quantreg, with significance level α = 0.1 and the default Barrodale and Roberts method to return confidence intervals for the estimated parameters. T o obtain the ensemble mean trends, we calculated the mean slope, upper bound and lower bound across the eight climate models at each geographi - cal location, then computed spatial averages for the full globe and three latitudinal regions. Variance. Trends in the magnitude of temporal variation of air tem - perature were examined at each geographical location using a mov - ing window approach. First, temperature was detrended by fitting a piecewise linear regression against time with Python package pwlf at each geographical location and extracting the residuals. Then, the temperature time series were divided into 10 yr windows starting in years 1855 through 2085 so as not to combine historical and future simulations (pre- and post-01-01-2015), and the variance of daily air temperature was calculated for each window. Windows were selected with no overlap to avoid statistical issues due to non-independence of estimates taken from partially overlapping time windows20. Scale-specific variability. Scale-specific variability was quantified using time-frequency decomposition. Specifically, at each geographi- cal location, wavelet analysis was conducted on multi-model mean temperature using the R package biwavelet43. Wavelet analysis resolves both the time and frequency domains of a signal (here a time series) via the wavelet transform. This is achieved via the convolution of a mother wavelet function and a time series across a set of windows τ and scales s.
We chose to use the Morlet wavelet, which represents a sine wave modulated by a Gaussian function:44 ψ0 (t) = π−1/4eiωt 0 e−r2 /2 where i is the imaginary unit, t represents non-dimensional time, and ω0 = 6 is the non-dimensional frequency 3. The continuous wavelet transform of a discrete time series x (t) with equal spacing δt and length T is defined as the convolution of x (t) with a normalized Morlet wavelet:44,45 Wx (s, τ) = √√√ √ δt s T−1 ∑ t=0 x (t) ψ0 ∗ ( (t − τ) δt s ) where * indicates the complex conjugate. By varying the wavelet scale s (that is, dilating and contracting the wavelet) and translating along localized time position τ , one can calculate the wavelet coefficients Wx (s, τ) across the different scales s and positions τ. These wavelet coef- ficients can be used to compute the bias-corrected local wavelet power, which describes how the contribution of each frequency or period in the time series varies over time:44,46,47 W2 x (s, τ) = 2s |Wx (s, τ)|2 where 2s is the bias correction factor 46. The scale s of the Morlet wavelet is related to the Fourier frequency f:47,48 1 f

4πs ω0 +√2+ω2 0 . When ω0 = 6, the scale s is approximately equal to the reciprocal of the Fourier frequency f, so period p ≈ s. The local wavelet power spectrum can then be visualized via heat maps and contour plots45,47. From the resulting local wavelet power spectrum heat map with time on the x axis, period (scale) on the y axis and power on the z axis, scale-averaged wavelet power was computed at annual (between 3 d and 2 yr) and multi-annual (between 2 yr and 30 yr) periodicities. This was achieved by taking the weighted sum of the local wavelet power across all scales for each time location τ:44,47 W−2 x (τ) = δjδt Cδ J ∑ j=0 ||Wx (sj, τ)| | 2 sj where Cδ = 0.776 for the Morlet wavelet, δJ represents the spacing between successive scales and δt represents the spacing between suc- cessive time locations 44. Scale-averaged power was then regressed against time using Generalized Least Squares (GLS) regression for the period 1850–2100 at each geographic location. T o determine the robustness of results to the choice of period for scale averaging, we also performed analysis of trends separately at interannual (between 2 yr and 7 yr) and multi-annual (between 7 yr and 30 yr) scales and found qualitatively similar results. Temporal autocorrelation. The temporal autocorrelation of air temperature was quantified by calculating the spectral exponent at each geographical location20. As described above, temperature was detrended by fitting a piecewise linear regression at each geographical location and extracting the residuals. The detrended temperature was divided into 10 yr windows starting in years 1855 through 2085. Fourier transforms of each time series were computed via fast Fourier trans- form using the Python package NumPy. Periodograms were prepared with frequency on the x axis and power spectral density on the y axis. The spectral exponent, β, was calculated as the slope of the regression line relating log-transformed power to log-transformed frequency. β expresses the relative contributions of frequencies to the power spectrum. In the case of equal contribution from all frequencies, β = 0. Greater contribution from low frequencies than from high frequencies results in a more negative value of β and indicates greater temporal autocorrelation in the time domain. Analysis of decadal trends. For each climate model, GLS regression was used to detect statistically significant trends (P < 0.05) in variance and temporal autocorrelation with respect to time in the presence of potentially autocorrelated residuals. T o measure inter-model agree- ment, we calculated the multi-model mean trend as the mean of trends calculated for each of the 8 models at each geographic location, then assessed the proportion of models that agreed with the sign of the multi-model mean trend. Inter-model agreement was considered as statistically significant at the α = 0.1 level on the basis of a binomial test. ANCOVA was used to quantify the relationship between temporal Nature Climate Change Article https://doi.org/10.1038/s41558-022-01490-7 autocorrelation and time while accounting for potential differences between land and sea environments. Statistically significant main effects and interactions were reported for P < 0.05. Modelling temperature impacts on ecology Thermal tolerance data. We obtained experimentally derived ther- mal tolerance parameters for a set of terrestrial ectotherms (n  = 38) published by Deutsch et al. 5 and used them to predict physiological response to CMIP6-simulated temperature. Deutsch et al.5 gathered data from 31 thermal performance studies published between 1974 and 2003 based on a collection of insects from 35 different locations. For each species, experimental intrinsic growth rates at multiple tem- peratures were used to fit a TPC, yielding least-squares estimates of key parameters such as critical thermal maximum (CTmax), optimum temperature (Topt), and sigma (σ ). We used a numerical scheme to reconstruct the curves whereby the rise in performance up to Topt was modelled as Gaussian and the decline beyond Topt was quadratic5,49: P (T) = ⎧⎪ ⎨⎪ ⎩ exp [− ( T−Topt 2σ ) 2 ] for T ≤ Topt 1 − ( T−Topt Topt −CTmax ) 2 for T > Topt . (1) This allowed negative growth rates to arise at high temperatures, but growth rates were bound at zero at low temperatures. Negative performance values indicate that mortality surpasses reproduction rates. Because P(T) is capped at 1 under this numerical scheme, P(T) represents the relative fitness of each species based on its normalized maximum growth rate. However, scaling this relative or normalized maximum growth rate by two orders of magnitude (that is, by a factor of 0.1 or 10.0) had limited quantitative and no qualitative impact on our results (Extended Data Fig. 4). Overall, increasing the growth rate scaling factor had no impact on population stability but promoted extinction risk. Isolation of temperature aspects. T o isolate projected changes in mean temperature and variability, we transformed the future (2050– 2100) time series using z-score normalization. Using this approach, we modified projected time series to match the historical (1950–2000) mean and/or standard deviation. Working in 10 yr moving windows between 2050 and 2100, each series x i with mean m 1 and standard deviation s1 was transformed to series y i with mean m2 and standard deviation s2: yi = m2 + (xi − m1) s2 s1 . (2) According to the scenario, m 2 and s2 were alternatively defined as (1) high emissions scenario mean and standard deviation (‘Mean, variance and autocorrelation’), (2) high emissions scenario mean and historical standard deviation (‘Mean and autocorrelation’), (3) histori- cal mean and high emissions scenario standard deviation (‘Variance and autocorrelation’) and (4) historical mean and standard deviation (‘ Autocorrelation’). High emissions scenario statistics refer to the properties of future series x i and confer no change to that aspect of the time series. Population dynamical modelling. T o model the effects of tempera- ture change on the stability and extinction probability of global ecto- therm populations, temperature dependence was integrated in the growth rate term of a population dynamical model50. While more com- plex synecological models can capture a range of community-level effects including competition and predation, we chose to model first order autecological dynamics to produce foundational insights about the role of temperature fluctuations on single-species population dynamics. Specifically, we used the r − α logistic growth model to simulate temperature-dependent growth and negative density-dependence: dN dt = N (rt − αN) (3) with population size N, time t, temperature-dependent growth rate rt, and self-regulation in the form of intraspecific competition α. This r − α logistic model is easily interconvertible with the classical r − K formula- tion (r/α = K), but has the advantage of handling negative values of r without issues51. This approach is sensitive to the effects of tempera- tures at and above the critical thermal maximum, which can yield negative growth rates that are important for determining population dynamics as well as long-term fitness. We extracted time series of daily temperature at the source loca- tions for each species from the ensemble of eight climate simulations. Daily intrinsic growth rates were computed from temperature using equation (1), incorporated into the r − α logistic growth model depicted in equation (3), and the model was then numerically solved using the explicit Runge-Kutta method of order 5(4) implemented in the Python SciPy package to obtain daily population densities. Rather than delin- eating active periods, which may shift under climate change, we con- sidered the full year to account for potential changes in fitness due to shifts in activity. The sensitivity of the results to strong (α = 1) and weak (α = 0.1) self-regulation was examined and found to be extremely limited (Extended Data Fig. 5). We also assessed the sensitivity of our results to absolute rather than relative or normalized growth rates by scaling rt by a factor of 0.1 or 10 in our simulations. Scaling rt by two orders of mag- nitude in this manner had very little quantitative and no qualitive impact on our results. This suggests that the effects of temperature fluctuations on changes in the spatiotemporal distribution of population abundance, stability and extinction were not contingent upon the use of relative fitness (that is, normalized growth rate) versus absolute fitness (that is, growth rate scaled by a factor of 0.1 or 10). These sensitivity analyses also showed that our results are robust to temperature-mediated changes in the maximum instantaneous growth rate52,53. Analysis of population changes. T o quantify temperature-driven changes in ecological stability and extinction probability, we compared population sizes and dynamics between a historical period (1950– 2000) and a future period (2050–2100). Here we defined latitudinal regions according to traditional delineations in ecology: Northern Hemisphere Extra-tropics, 60° S to 23° S; Tropics, 23° S to 23° N; and Southern Hemisphere Extra-tropics, 23° N to 60° N. Population abundance was computed as the mean population
size (N) for a time period. Population stability was computed as the inverse of the coefficient of variation, or mean population divided
by population standard deviation. Percent changes in population size and stability were computed for each of the climate models as (future − historical) /historical × 100% and plotted without outliers in
Fig. 4. Statistically significant changes in population abundance and stability between the historical and future periods were identified via the Mann–Whitney U-test, with the eight models as replicates. Extinction probability was quantified as the proportion of ensem- ble simulations for which the population declined to zero during a 50 yr simulation. Changes in extinction probability were calculated as the difference between future and historical extinction probabilities. Statistically significant changes in extinction probability were identi- fied on a regional basis via the Mann–Whitney U-test. Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article. Nature Climate Change Article https://doi.org/10.1038/s41558-022-01490-7 Data availability The CMIP6 simulation data used in this paper are available via the data portal https://esgf-node.llnl.gov/search/cmip6/. The ecology data are available for download at https://doi.org/10.1073/pnas.0709472105. Code availability The code can be accessed on GitHub at https://github.com/KateDuffy/ climate-change-ecology54. References 43. Gouhier, T. C., Grinsted, A. & Simko, V. R package biwavelet: Conduct Univariate and Bivariate Wavelet Analyses (Version 0.20.17). 44. Torrence, C. & Compo, G. P. A practical guide to wavelet analysis. Bull. Am. Meteorol. Soc. 79, 61–78 (1998). 45. Grinsted, A., Moore, J. C. & Jevrejeva, S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Process. Geophys. 11, 561–566 (2004). 46. Liu, Y., Liang, X. S. & Weisberg, R. H. Rectification of the bias
in the wavelet power spectrum. J. Atmos. Ocean. Technol. 24, 2093–2102 (2007). 47. Cazelles, B. et al. Wavelet analysis of ecological time series. Oecologia 156, 287–304 (2008). 48. Maraun, D. & Kurths, J. Cross wavelet analysis: significance testing and pitfalls. Nonlinear Process. Geophys. 11, 505–514 (2004). 49. Huey, R. B. & Stevenson, R. D. Integrating thermal physiology and ecology of ectotherms: a discussion of approaches. Am. Zool. 19, 357–366 (1979). 50. Mallet, J. The struggle for existence: how the notion of carrying capacity, K, obscures the links between demography, Darwinian evolution, and speciation. Evol. Ecol. Res. 14, 627–665 (2012). 51. Vasseur, D. A. Theoretical Ecology: Concepts and Applications (eds. McCann, K. S. & Gellner, G.) 243–262 (Oxford Univ. Press, 2020). 52. Estay, S. A., Clavijo-Baquet, S., Lima, M. & Bozinovic, F. Beyond average: an experimental test of temperature variability on the population dynamics of Tribolium confusum. Popul. Ecol. 53, 53–58 (2011). 53. Bozinovic, F. et al. The mean and variance of environmental temperature interact to determine physiological tolerance and fitness. Physiol. Biochem. Zool. 84, 543–552 (2011). 54. Duffy, K. climate-change-ecology (GitHub, 2022); https://github. com/KateDuffy/climate-change-ecology Acknowledgements This work was primarily supported by the National Science Foundation (NSF) grant CCF-1442728 while K.D. was a PhD
student at the SDS Lab in Northeastern University. Additional support was provided for K.D. and A.R.G. by NSF SES-1735505
and for T.C.G. by NSF OCE-2048894. We acknowledge the background support from a previous NSF Expeditions in Computing grant (award no. 1029711) and an ongoing DOD Strategic Environmental Research and Development Program funding (no. RC20-1183). K.D. and A.R.G. acknowledge support from the NASA Ames Research Center. Author contributions K.D., T.C.G. and A.R.G. conceived, designed and refined the project. K.D. performed the data analysis and modelling. K.D., T.C.G. and A.R.G. interpreted the results. K.D. wrote the manuscript with contributions from T.C.G. and A.R.G. Competing interests The authors declare no competing interests. Additional information Extended data is available for this paper at
https://doi.org/10.1038/s41558-022-01490-7. Supplementary information The online version
contains supplementary material available at
https://doi.org/10.1038/s41558-022-01490-7. Correspondence and requests for materials should be addressed to Kate Duffy. Peer review information Nature Climate Change thanks David Vasseur and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Reprints and permissions information is available at
www.nature.com/reprints.