Philippa Kaur

1.IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (eds Pörtner, H.-O. et al.) (IPCC, 2019).2.Doney, S. C. et al. Climate change impacts on marine ecosystems. Ann. Rev. Mar. Sci. 4, 11–37 (2012).Article 

Google Scholar 
3.Bindoff, N. L. et al. in Special Report on the Ocean and Cryosphere in a Changing Climate (eds Pörtner, H.-O. et al.) Ch. 5 (IPCC, 2019).4.Griffith, G. P., Fulton, E. A., Gorton, R. & Richardson, A. J. Predicting interactions among fishing, ocean warming, and ocean acidification in a marine system with whole-ecosystem models. Conserv. Biol. 26, 1145–1152 (2012).Article 

Google Scholar 
5.Fu, C. et al. Risky business: the combined effects of fishing and changes in primary productivity on fish communities. Ecol. Modell. 368, 265–276 (2018).Article 

Google Scholar 
6.Tittensor, D. P. et al. Integrating climate adaptation and biodiversity conservation in the global ocean. Sci. Adv. https://doi.org/10.1126/sciadv.aay9969 (2019).7.IPBES: Summary for Policymakers. In Global Assessment Report on Biodiversity and Ecosystem Services (eds Díaz, S. et al.) (IPBES Secretariat, 2019).8.Boyce, D. G., Lotze, H. K., Tittensor, D. P., Carozza, D. A. & Worm, B. Future ocean biomass losses may widen socioeconomic equity gaps. Nat. Commun. 11, 2235 (2020).CAS 
Article 

Google Scholar 
9.Payne, M. R. et al. Uncertainties in projecting climate-change impacts in marine ecosystems. ICES J. Mar. Sci. 73, 1272–1282 (2016).Article 

Google Scholar 
10.Eyring, V. et al. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. Geosci. Model Dev. 9, 1937–1958 (2016).Article 

Google Scholar 
11.Tittensor, D. P. et al. A protocol for the intercomparison of marine fishery and ecosystem models: Fish-MIP v1.0. Geosci. Model Dev. 11, 1421–1442 (2018).Article 

Google Scholar 
12.Lotze, H. K. et al. Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change. Proc. Natl Acad. Sci. USA 116, 12907–12912 (2019).CAS 
Article 

Google Scholar 
13.Bryndum-Buchholz, A. et al. Twenty-first-century climate change impacts on marine animal biomass and ecosystem structure across ocean basins. Glob. Change Biol. 25, 459–472 (2019).Article 

Google Scholar 
14.Bryndum-Buchholz, A. et al. Differing marine animal biomass shifts under 21st century climate change between Canada’s three oceans. Facets 5, 105–122 (2020).Article 

Google Scholar 
15.Bryndum-Buchholz, A. et al. Climate-change impacts and fisheries management challenges in the North Atlantic Ocean. Mar. Ecol. Prog. Ser. 648, 1–17 (2020).Article 

Google Scholar 
16.Ruane, A. C. et al. The vulnerability, impacts, adaptation and climate services advisory board (VIACS AB v1.0) contribution to CMIP6. Geosci. Model Dev. 9, 3493–3515 (2016).Article 

Google Scholar 
17.Kwiatkowski, L. et al. Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections. Biogeosciences 17, 3439–3470 (2020).CAS 
Article 

Google Scholar 
18.Séférian, R. et al. Tracking improvement in simulated marine biogeochemistry between CMIP5 and CMIP6. Curr. Clim. Change Rep. 6, 95–119 (2020).Article 

Google Scholar 
19.Meehl, G. A. et al. Context for interpreting equilibrium climate sensitivity and transient climate response from the CMIP6 Earth system models. Sci. Adv. 6, eaba1981 (2020).Article 

Google Scholar 
20.Tebaldi, C. et al. Climate model projections from the Scenario Model Intercomparison Project (ScenarioMIP) of CMIP6. Earth Syst. Dyn. 12, 253–293 (2021).Article 

Google Scholar 
21.Heneghan, R. F. et al. Disentangling diverse responses to climate change among global marine ecosystem models. Prog. Oceanogr. 198, 102659 (2021).Article 

Google Scholar 
22.Zelinka, M. D. et al. Causes of higher climate sensitivity in CMIP6 models. Geophys. Res. Lett. 47, e2019GL085782 (2020).Article 

Google Scholar 
23.Kwiatkowski, L. et al. Emergent constraints on projections of declining primary production in the tropical oceans. Nat. Clim. Change 7, 355–358 (2017).CAS 
Article 

Google Scholar 
24.Cabré, A., Marinov, I. & Leung, S. Consistent global responses of marine ecosystems to future climate change across the IPCC AR5 Earth system models. Clim. Dyn. 45, 1253–1280 (2015).Article 

Google Scholar 
25.Laufkötter, C. et al. Drivers and uncertainties of future global marine primary production in marine ecosystem models. Biogeosciences 12, 6955–6984 (2015).Article 

Google Scholar 
26.Doney, S. C. Plankton in a warmer world. Nature 444, 695–696 (2006).CAS 
Article 

Google Scholar 
27.Rykaczewski, R. R. & Dunne, J. P. Enhanced nutrient supply to the California Current Ecosystem with global warming and increased stratification in an Earth system model. Geophys. Res. Lett. 37, L21606 (2010).Article 

Google Scholar 
28.Laufkötter, C., John, J. G., Stock, C. A. & Dunne, J. P. Temperature and oxygen dependence of the remineralization of organic matter. Glob. Biogeochem. Cycles 31, 1038–1050 (2017).Article 
CAS 

Google Scholar 
29.Coll, M. et al. Advancing global ecological modeling capabilities to simulate future trajectories of change in marine ecosystems. Front. Mar. Sci. 7, 741 (2020).Article 

Google Scholar 
30.Hawkins, E. & Sutton, R. The potential to narrow uncertainty in regional climate predictions. Bull. Am. Meteorol. Soc. 90, 1095–1107 (2009).Article 

Google Scholar 
31.Frölicher, T. L., Rodgers, K. B., Stock, C. A. & Cheung, W. W. L. Sources of uncertainties in 21st century projections of potential ocean ecosystem stressors. Glob. Biogeochem. Cycles 30, 1224–1243 (2016).Article 
CAS 

Google Scholar 
32.Gaines, S. D. et al. Improved fisheries management could offset many negative effects of climate change. Sci. Adv. 4, eaao1378 (2018).Article 

Google Scholar 
33.The State of World Fisheries and Aquaculture 2020 (FAO, 2020).34.Dahlke, F. T., Wohlrab, S., Butzin, M. & Pörtner, H.-O. Thermal bottlenecks in the life cycle define climate vulnerability of fish. Science 369, 65–70 (2020).CAS 
Article 

Google Scholar 
35.Stuart-Smith, R. D., Edgar, G. J. & Bates, A. E. Thermal limits to the geographic distributions of shallow-water marine species. Nat. Ecol. Evol. 1, 1846–1852 (2017).Article 

Google Scholar 
36.Carozza, D. A., Bianchi, D. & Galbraith, E. D. Metabolic impacts of climate change on marine ecosystems: implications for fish communities and fisheries. Glob. Ecol. Biogeogr. 28, 158–169 (2019).Article 

Google Scholar 
37.du Pontavice, H., Gascuel, D., Reygondeau, G., Stock, C. & Cheung, W. W. L. Climate-induced decrease in biomass flow in marine food webs may severely affect predators and ecosystem production. Glob. Change Biol. 27, 2608–2622 (2021).Article 

Google Scholar 
38.Piroddi, C. et al. Effects of nutrient management scenarios on marine food webs: a pan-European assessment in support of the marine strategy framework directive. Front. Mar. Sci. 8, 179 (2021).Article 

Google Scholar 
39.Maury, O. An overview of APECOSM, a spatialized mass balanced ‘Apex Predators ECOSystem Model’ to study physiologically structured tuna population dynamics in their ecosystem. Prog. Oceanogr. 84, 113–117 (2010).Article 

Google Scholar 
40.Maury, O. & Poggiale, J. C. From individuals to populations to communities: a dynamic energy budget model of marine ecosystem size-spectrum including life history diversity. J. Theor. Biol. 324, 52–71 (2013).Article 

Google Scholar 
41.Carozza, D. A., Bianchi, D. & Galbraith, E. D. The ecological module of BOATS-1.0: a bioenergetically-constrained model of marine upper trophic levels suitable for studies of fisheries and ocean biogeochemistry. Geosci. Model Dev. 9, 1545–1565 (2016).Article 

Google Scholar 
42.Carozza, D. A. et al. Formulation, general features and global calibration of a bioenergetically-constrained fishery model. PLoS ONE 12, e0169763 (2017).Article 
CAS 

Google Scholar 
43.Cheung, W. W. L. et al. Building confidence in projections of the responses of living marine resources to climate change. ICES J. Mar. Sci. https://doi.org/10.1093/icesjms/fsv250 (2016).Article 

Google Scholar 
44.Cheung, W. W. L., Dunne, J., Sarmiento, J. L. & Pauly, D. Integrating ecophysiology and plankton dynamics into projected maximum fisheries catch potential under climate change in the Northeast Atlantic. ICES J. Mar. Sci. 68, 1008–1018 (2011).Article 

Google Scholar 
45.Blanchard, J. L. et al. Potential consequences of climate change for primary production and fish production in large marine ecosystems. Phil. Trans. R. Soc. B 367, 2979–2989 (2012).Article 

Google Scholar 
46.Christensen, V. et al. The global ocean is an ecosystem: simulating marine life and fisheries. Glob. Ecol. Biogeogr. 24, 507–517 (2015).Article 

Google Scholar 
47.Gascuel, D., Guénette, S. & Pauly, D. The trophic-level-based ecosystem modelling approach: theoretical overview and practical uses. ICES J. Mar. Sci. 68, 1403–1416 (2011).Article 

Google Scholar 
48.Petrik, C. M., Stock, C. A., Andersen, K. H., van Denderen, P. D. & Watson, J. R. Bottom-up drivers of global patterns of demersal, forage, and pelagic fishes. Prog. Oceanogr. 176, 102124 (2019).Article 

Google Scholar 
49.Jennings, S. & Collingridge, K. Predicting consumer biomass, size-structure, production, catch potential, responses to fishing and associated uncertainties in the world’s marine ecosystems. PLoS ONE 10, e0133794 (2015).Article 
CAS 

Google Scholar 
50.Heneghan, R. F. et al. A functional size-spectrum model of the global marine ecosystem that resolves zooplankton composition. Ecol. Modell. 435, 109265 (2020).CAS 
Article 

Google Scholar 
51.Dunne, J. P. et al. GFDL’s ESM2 global coupled climate–carbon Earth system models. Part I: physical formulation and baseline simulation characteristics. J. Clim. 25, 6646–6665 (2012).Article 

Google Scholar 
52.Dunne, J. P. et al. Carbon Earth system models. Part II: carbon system formulation and baseline simulation characteristics. J. Clim. 26, 2247–2267 (2013).Article 

Google Scholar 
53.Dufresne, J.-L. et al. Climate change projections using the IPSL-CM5 Earth system model: from CMIP3 to CMIP5. Clim. Dyn. 40, 2123–2165 (2013).Article 

Google Scholar 
54.Dunne, J. P. et al. The GFDL Earth System Model Version 4.1 (GFDL-ESM 4.1): overall coupled model description and simulation characteristics. J. Adv. Model. Earth Syst. 12, e2019MS002015 (2020).
Google Scholar 
55.Krasting, J. P. et al. NOAA-GFDL GFDL-ESM4 Model Output Prepared for MIP6 CMIP Historical Version 20190726 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.859756.John, J. G. et al. NOAA-GFDL GFDL-ESM4 Model Output Prepared for CMIP6 ScenarioMIP ssp585 Version 20180701 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.870657.Boucher, O. et al. Presentation and evaluation of the IPSL-CM6A-LR climate model. J. Adv. Model. Earth Syst. 12, e2019MS002010 (2020).Article 

Google Scholar 
58.Boucher, O. et al. IPSL IPSL-CM6A-LR Model Output Prepared for CMIP6 CMIP Version 20180727 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.153459.Boucher, O. et al. IPSL IPSL-CM6A-LR Model Output Prepared for CMIP6 CMIP Historical Version 20180103 (Earth System Grid Federation, 2018); https://doi.org/10.22033/ESGF/CMIP6.5195 More

Global statistical analysisOur first attempt to identify plausible effects of meteorological covariates on COVID-19 spread applied a comparative regression analysis. To this end, we focused on the exponential onset of the disease, as it is the epidemic phase that allows for a better comparison between countries or regions, without the confounding effect of intervention policies. We first determined, for each of the spatial units (either countries or NUTS (nomenclature of territorial units for statistics) 2 regions), the day in which 20 or more cumulative cases were officially reported. We then fitted the first-order polynomial function f(t) = x0 + rt for the next 20 days of log-transformed data, where t represents time (in days) and ({{x}_0}) is the value at initial condition t = 0. The r parameter can be understood as the exponential growth rate, and is then used to estimate the basic reproduction number (R0) using the estimated serial interval T for COVID-19 of 4.7 days53, such that R0 = 1 + rT (ref. 54). (We note that we are interested here in the relationship between the reproductive number and not in the actual inference of R0.) Once R0 was obtained for all our spatial units, we filtered our meteorological data to match the same fitting period (with a 10-day negative delay to account for an incubation and reporting lapse) for every spatial unit. To compute a single average of the meteorological variables per regional unit, we computed a weighted average on the basis of the population contribution of each grid cell to the total population of the region. We did so to have an aggregated value that would better represent the impact of these factors on the population transmission of COVID-19, as the same variation in weather in a high-density urban area is more likely to contribute to a change in population-level transmission than that of an unpopulated rural area. We then averaged the daily values of temperature and AH for each country and computed univariate linear models for each of these variables as predictors of R0. Given the somewhat arbitrary criteria to select the dates to estimate the R0 in each country, a sensitivity analysis was run to test the robustness of the regressions to changes in the related parameters. We tested 70 different combinations of two parameters: the total number of days used for the fit (18–27) and the threshold of cumulative COVID-19 cases used to select the initial day of the fit (15–45). We also calculated the weather averages by shifting the selected dates accordingly. Then, a linear model for each of the estimates was fitted for both T and AH. A summary of the distribution of parameter estimates (the regression slope coefficients and the R2 of the models) is shown in Extended Data Fig. 3.Bivariate time-series analysis with scale-dependent correlationsTo examine associations between cases and climate factors in more detail, SDC was performed on the daily time series of both COVID-19 incidence and a given meteorological variable. SDC is an optimal method for identifying dynamical couplings in short and noisy time series20,21. In general, Spearman correlations between incidence and a meteorological time series assess whether there is a monotonic relation between the variables. SDC analysis was specifically developed to study transitory associations that are local in time at a specified temporal scale corresponding to the size of the time intervals considered (s). The two-way implementation (TW-SDC) is a bivariate method that computes non-parametric Spearman rank correlations between two time series, for different pairs of time intervals along these series. Different window sizes (s) can be used to examine increasingly finer temporal resolution. The results are sensitive to the value of this window size, s, with expected significant and highest correlation values at the scale of the transient coupling between variables. Correlation values decrease in magnitude as window size increases, and averages are computed over too long a time interval. Values can also decrease and become non-significant for small windows when correlations are spurious. Here, the method was applied for windows of different length (from s = 75 to 14 days) and, despite a weekly cycle showing up in some cases for small s, results removing this cycle were robust. We therefore did not remove this cycle.The results are typically displayed in a figure with the following subplots: (1) the two time series, to the left and top of the matrix of correlation values, respectively; (2) the matrix or grid of correlation values itself in the center, with significant correlations colored in blue when positive and in red when negative, with rows and columns corresponding to the temporal localization of the moving window along the time series on the left and top, respectively; (3) a time series at the bottom, below this grid, with the highest significant correlations for a given time (vertically, and therefore for the variable that acts as the driver, here the meteorological time series). To read the results, one starts at the diagonal and moves vertically down from it to identify a given lag for which significant correlations are found (the closest to the main diagonal). In some of the SDC figures, the time intervals with high local correlations are highlighted with boxes. These intervals alternate with other ones (left blank) for which no significant correlation is found. All colored areas correspond to significance levels of at least P  fs/fr, where fs is the sampling rate and fr the minimum frequency. Another strategy is that M be large enough that the M-lagged vector incorporates the temporal scale of the time series that is of interest. The larger the M, the more detailed the resulting decomposition of the signal. In particular, the most detailed decomposition is achieved when the embedding dimension is approximately equal to half of the total signal length. A compromise must be reached, however, as a large M implies increased computation, and too large a value may produce mixing of components. SSA is especially well suited for separating components corresponding to different frequencies in nonlinear systems. Here, we applied it to remove the weekly cycle.MSDC analysisMSDC provides a scan of the SDC analyses over a range of different scales (here, S from 5 to 100 days at 5-day intervals), by selecting the maximum correlation values (positive or negative) closer to the diagonal. The goal is to consider the evolution of transient correlations at all scales pooled together in a single analysis. The MSDC plot displays time on the x axis and scale (S) on the y axis, and positive and negative correlations either jointly or separately. The rationale behind MSDC is that correlations at very small scales can occur by chance because of coincident similar patterns, but that as one moves up to larger scales (by increasing S), the correlation patterns that are spurious tend to vanish, whereas those reflecting mechanistic links increase in strength. This increase in correlation values should occur up to the real scale of interaction, decreasing afterwards. By ‘real’, we mean here the temporal scale covering the extent of the interaction between the driver and the response process (in this case, the response of disease transmission to a given climate factor). Thus, continuity of the same sign correlations together with transitions to larger values are indicative of causal effects, whereas the rapid vanishing of small-scale significant correlations signals spurious ones.Process-based modelDescriptionThe dynamical model is a discrete stochastic model that incorporates seven different compartments: S, E, I, C, Q, R and D. The model structure is illustrated in Fig. 4. The transition probabilities of the stochastic model are based on the corresponding rates of the transitions between classes in the deterministic (mean-field) model (specified in Fig. 4b). These probabilities are defined as follows. P(e) = (1.0 − exp(−β dt)) is the probability of infection exposure of the susceptible class, where β = (1/N)(βII + βQQ) is the infection rate (of the deterministic model). P(i) = (1.0 − exp(−γ dt)) is the probability that an new exposed individual becomes infectious, where γ denotes the incubation rate. P(r) = (1.0 − exp(−Λ dt)) is the recovery probability, where λ0(1 − exp(λ1t)) is the (deterministic) recovery rate. P(p) = (1.0 − exp(−α dt)) is the protection probability, where α = α0exp(α1t). P(d) = (1.0 − exp(−K dt)) is the mortality probability, with K = k0exp(k1t). P(re) = (1.0 − exp(−τ dt)) is the release probability from confinement, where τ = τ0exp(τ1t). Finally, P(q) = (1.0 − exp(−δ  dt)) is the detection probability, where δ is the quarantine rate (for example, at which infected individuals are isolated from the rest of the population).In the model, both infected non-detected and infected detected individuals can infect susceptible ones. In the model incorporating temperature in the transmission rate, the respective values of βI and βQ are calculated as follows:$${beta }_{I}(t)={beta }_{I},T_{mathrm{inv}}(t);quad {beta }_{Q}(t)={beta }_{Q},T_{mathrm{inv}}(t)$$where (T_{mathrm{inv}}=fleft(frac{1-T(t)}{bar{T}}right)), with (bar{T}) corresponding to the overall mean of the temperature time series and f(·) to a Savitzky–Golay filter, used to smooth the temperature series with a window size of 50 data points and a polynomial order of 3. When the infection rate is constant, we simply omit the temperature term. For further comparison, in a third model, β is specified with a sinusoidal function of period equal to 12 months and an estimated phase.The number of individuals transitioning from compartment i to j at time t are determined by means of binomial distributions P(Xi,P(y)), where Xi corresponds to one of the compartments S, E, I, Q, R, D, C, and P(y) to the respective transition probability defined above. Thus,

e(t) = P(S(t), P(e)), new exposed individuals at time t

p(t) = P(S(t), P(p)), protected individuals at time t

i(t) = P(E(t), P(i)), new infected not detected individuals at time t

q(t) = P(I(t), P(q)), new infected and detected individuals at time t

r(t) = P(Q(t), P(r)), total recovered individuals at time t

d(t) = P(Q(t), P(d)), total dead individuals at time t

re(t) = P(C(t), P(re)), individuals released from confinement at time t

Then, the final dynamics are given by the following equations:$$S(t)=S(t-{rm{d}}t)-e(t)-p(t)+re(t)$$$$E(t)=E(t-{rm{d}}t)+e(t)-i(t)$$$$I(t)=I(t-{rm{d}}t)+i(t)-q(t)$$$$Q(t)=Q(t-{rm{d}}t)+q(t)-r(t)-d(t)$$$$R(t)=R(t-{rm{d}}t)+r(t)$$$$D(t)=D(t-{rm{d}}t)+d(t)$$$$C(t)=C(t-{rm{d}}t)+p(t)-re(t)$$CalibrationThe model was implemented using Python and calibrated by means of the least squares algorithm of the scipy library. The error function minimized with this algorithm was obtained from the normalized residuals on the basis of total cases (Q + R + D) and deaths (D).To search parameter space, we ran 100 calibrations starting from different initial choices of parameter combinations. The tolerance for termination in the change of the cost function was set to 1 × 10−10. Tolerance for termination by the norm of the gradient was also set to 1 × 10−10, and the tolerance for termination by the change of the independent variables was set to 1 × 10−10. The solver was the lsmr method (which is suitable for problems with sparse and large Jacobian matrices) with a differential step of 1 × 10−5. With this configuration, each fitting run usually converged after ~500 iterations.ValidationTo compare the model including an effect of T in the transmission rate to those without it, we calculated the chi-square, Akaike information criterion (AIC) and Bayesian information criterion (BIC) indices for the residuals obtained from the optimization process. The resulting values are shown in Supplementary Table 1.Our choice of T to modulate the infection rate (β) instead of AH underlies the fact that the temporal dynamics of both factors roughly follow the same shape, with the advantage that T shows less oscillatory behavior than AH. This fact adds stability to the model when the inverse relationship is used in the calculation of β (Supplementary Information). This selection is further reinforced by the results from the SDC analyses, which yielded larger correlations for temperature, even when penalizing for the larger autocorrelation structure.Our choice to modulate β using T instead of AH follows from the fact that the temporal dynamics of both climate variables present roughly the same shape, with the advantage that T exhibits weaker oscillations. This less fluctuating pattern provides stability to the model fitting when the inverse relationship is used in the calculation of β (Supplementary Information). Additionally, the transient correlations obtained with SDC yielded higher values for T than for AH (even when accounting for concurrent levels of autoregression in the two variables). More

1.Dairy Production and Products: Milk and Milk Products (FAO, 2013); http://www.fao.org/dairy-production-products/production/dairy-animals/cattle/en/2.Background: Corn and Other Feedgrains (USDA ERS, 2018); https://www.ers.usda.gov/topics/animal-products/dairy/background/3.National Agricultural Statistics Service (US Department of Agriculture); https://www.nass.usda.gov/index.php4.Capper, J. L., Cady, R. A. & Bauman, D. E. The environmental impact of dairy production: 1944 compared with 2007. J. Anim. Sci. 87, 2160–2167 (2009).CAS 
Article 

Google Scholar 
5.Niles, M. T. & Wiltshire, S. Tradeoffs in US dairy manure greenhouse gas emissions, productivity, climate, and manure management strategies. Environ. Res. Commun 1, 075003 (2019).Article 

Google Scholar 
6.Field, T. G. & Taylor, R. E. Scientific Farm Animal Production: An Introduction, Eleventh Edition (Pearson, 2018).7.Fuquay, J. W. Heat stress as it affects animal production. J. Anim. Sci. 52, 164–174 (1981).CAS 
Article 

Google Scholar 
8.St-Pierre, N. R., Cobanov, B. & Schnitkey, G. Economic losses from heat stress by US livestock industries. J. Dairy Sci. 86, E52–E77 (2003).Article 

Google Scholar 
9.Kadzere, C. T., Murphy, M. R., Silanikove, N. & Maltz, E. Heat stress in lactating dairy cows: a review. Livest. Prod. Sci. 77, 59–91 (2002).Article 

Google Scholar 
10.Bouraoui, R., Lahmar, M., Majdoub, A., Djemali, M. & Belyea, R. The relationship of temperature–humidity index with milk production of dairy cows in a Mediterranean climate. Anim. Res. 51, 479–491 (2002).Article 

Google Scholar 
11.West, J. W. Effects of heat-stress on production in dairy cattle. J. Dairy Sci. 86, 2131–2144 (2003).CAS 
Article 

Google Scholar 
12.Vitali, A. et al. Seasonal pattern of mortality and relationships between mortality and temperature–humidity index in dairy cows. J. Dairy Sci. 92, 3781–3790 (2009).CAS 
Article 

Google Scholar 
13.Pragna, P. et al. Heat stress and dairy cow: impact on both milk yield and composition. Int. J. Dairy Sci. 12, 1–11 (2017).CAS 
Article 

Google Scholar 
14.Hoffmann, I. Climate change and the characterization, breeding and conservation of animal genetic resources. Anim. Genet. 41, 32–46 (2010).Article 

Google Scholar 
15.Qi, L., Bravo-Ureta, B. E. & Cabrera, V. E. From cold to hot: a preliminary analysis of climatic effects on the productivity of Wisconsin dairy farms. AgEconSearch https://doi.org/10.22004/ag.econ.172411 (2014).16.Bohmanova, J., Misztal, I. & Cole, J. B. Temperature–humidity indices as indicators of milk production losses due to heat stress. J. Dairy Sci. 90, 1947–1956 (2007).CAS 
Article 

Google Scholar 
17.Field, C. B. et al. Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation (IPCC, 2021); https://www.ipcc.ch/report/managing-the-risks-of-extreme-events-and-disasters-to-advance-climate-change-adaptation/18.Mueller, N. D. et al. Cooling of US Midwest summer temperature extremes from cropland intensification. Nat. Clim. Chang. 6, 317–322 (2016).ADS 
MathSciNet 
Article 

Google Scholar 
19.Seneviratne, S. I., Donat, M. G., Mueller, B. & Alexander, L. V. No pause in the increase of hot temperature extremes. Nat. Clim. Chang. 4, 161–163 (2014).ADS 
Article 

Google Scholar 
20.Dairy 2014: Dairy Cattle Management Practices in the United States, 2014 (USDA, APHIS, NAHMS, 2016); https://www.aphis.usda.gov/animal_health/nahms/dairy/downloads/dairy14/Dairy14_dr_PartI_1.pdf21.Mondaca, M. R. & Cook, N. B. Modeled construction and operating costs of different ventilation systems for lactating dairy cows. J. Dairy Sci. 102, 896–908 (2019).CAS 
Article 

Google Scholar 
22.Ferreira, F. C., Gennari, R. S., Dahl, G. E. & De Vries, A. Economic feasibility of cooling dry cows across the United States. J. Dairy Sci. 99, 9931–9941 (2016).CAS 
Article 

Google Scholar 
23.Hayhoe, K. et al. Emissions pathways, climate change, and impacts on California. Proc. Natl Acad. Sci. USA 101, 12422–12427 (2004).ADS 
CAS 
Article 

Google Scholar 
24.Klinedinst, P. L., Wilhite, D. A., Hahn, L. G. & Hubbard, K. G. The potential effects of climate change on summer seasonal dairy cattle milk production and reproduction. Clim. Chang. 23, 21–36 (1993).ADS 
Article 

Google Scholar 
25.Mauger, G., Bauman, Y., Nennich, T. & Salathé, E. Impacts of climate change on milk production in the United States. Prof. Geogr. 67, 121–131 (2015).Article 

Google Scholar 
26.Key, N. & Sneeringer, S. Potential effects of climate change on the productivity of U.S. dairies. Am. J. Agric. Econ. 96, 1136–1156 (2014).Article 

Google Scholar 
27.Ortiz-Bobea, A., Knippenberg, E. & Chambers, R. G. Growing climatic sensitivity of U.S. agriculture linked to technological change and regional specialization. Sci. Adv. 4, eaat4343 (2018).ADS 
Article 

Google Scholar 
28.Butler, E. E., Mueller, N. D. & Huybers, P. Peculiarly pleasant weather for US maize. Proc. Natl Acad. Sci. USA 115, 11935–11940 (2018).CAS 
Article 

Google Scholar 
29.Lobell, D. B., Schlenker, W. & Costa-Roberts, J. Climate trends and global crop production since 1980. Science 333, 616–620 (2011).ADS 
CAS 
Article 

Google Scholar 
30.Tigchelaar, M., Battisti, D. S., Naylor, R. L. & Ray, D. K. Future warming increases probability of globally synchronized maize production shocks. Proc. Natl Acad. Sci. U. S. A. 115, 6644–6649 (2018).ADS 
Article 

Google Scholar 
31.PRISM Climate Data (Oregon State Univ., 2019); http://www.prism.oregonstate.edu/32.Daly, C. et al. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol. https://doi.org/10.1002/joc.1688 (2008).33.National Research Council. Nutrient Requirements of Dairy Cattle, Seventh Revised Edition (National Academies Press, 2001).34.Auldist, M. J., Walsh, B. J. & Thomson, N. A. Seasonal and lactational influences on bovine milk composition in New Zealand. J. Dairy Res. 65, 401–411 (1998).CAS 
Article 

Google Scholar 
35.Lobell, D. B. Climate change adaptation in crop production: beware of illusions. Glob. Food Sec. 3, 72–76 (2014).Article 

Google Scholar 
36.Mukherjee, D., Bravo-Ureta, B. E. & De Vries, A. Dairy productivity and climatic conditions: econometric evidence from South-eastern United States. Aust. J. Agric. Resour. Econ. 57, 123–140 (2013).Article 

Google Scholar 
37.Milk Cost of Production Estimates: Cost-of-Production Estimates-2016 Base (USDA ERS, 2021); https://www.ers.usda.gov/data-products/milk-cost-of-production-estimates/milk-cost-of-production-estimates/#Milk38.Liang, X. Z. et al. Determining climate effects on US total agricultural productivity. Proc. Natl Acad. Sci. USA 114, E2285–E2292 (2017).CAS 
Article 

Google Scholar 
39.Malikov, E., Miao, R. & Zhang, J. Distributional and temporal heterogeneity in the climate change effects on U.S. agriculture. J. Environ. Econ. Manage. 104, 102386 (2020).Article 

Google Scholar 
40.MacDonald, J. M., Law, J. & Mosheim, R. Consolidation in U.S. Dairy Farming Economic Research Report No. 274 (ERS, USDA, 2020); https://www.ers.usda.gov/publications/pub-details/?pubid=9890041.Hemme, T. & Otte, J. Pro-Poor Livestock Policy Initiative Status and Prospects for Smallholder Milk Production a Global Perspective (Food and Agriculture Organization of the United Nations, 2010).42.Osei-Amponsah, R. et al. Heat stress impacts on lactating cows grazing Australian summer pastures on an automatic robotic dairy. Animals 10, 869 (2020).Article 

Google Scholar 
43.Chang-Fung-Martel, J., Harrison, M. T., Rawnsley, R., Smith, A. P. & Meinke, H. The impact of extreme climatic events on pasture-based dairy systems: a review. Crop Pasture Sci 68, 1158 (2017).Article 

Google Scholar 
44.Livestock Hot Weather Stress. Operations Manual (NOAA, 1976); https://scirp.org/reference/referencespapers.aspx?referenceid=191321645.Pinheiro J., Bates D., Debroy S. S. D. Linear and nonlinear mixed effects models, R package nlme version 3.1-152 (2021).46.Conley, T. G. GMM estimation with cross sectional dependence. J. Econom. 92, 1–45 (1999).MathSciNet 
Article 

Google Scholar 
47.Borchers, H. W. pracma: practical numerical math functions, version 2.2.9.1–393 (2019).48.Colin Cameron, A., Gelbach, J. B. & Miller, D. L. Robust inference with multiway clustering. J. Bus. Econ. Stat. 29, 238–249 (2011).MathSciNet 
Article 

Google Scholar 
49.Zeileis, A., Köll, S. & Graham, N. Various versatile variances: an object-oriented implementation of clustered covariances in R. J. Stat. Softw. https://doi.org/10.18637/jss.v095.i01 (2020). More

1.Gartner, T., Mulligan, J., Schmidt, R. & Gunn, J. Natural Infrastructure (World Resources Institute, 2013).2.McDonald, R. I. et al. Water on an urban planet: urbanization and the reach of urban water infrastructure. Glob. Environ. Change 27, 96–105 (2014).Article 

Google Scholar 
3.Vorosmarty, C. J. et al. Global threats to human water security and river biodiversity. Nature 467, 555–561 (2010).CAS 
Article 

Google Scholar 
4.Grill, G. et al. Mapping the world’s free-flowing rivers. Nature 569, 215–221 (2019).CAS 
Article 

Google Scholar 
5.Tessler, Z. D. et al. Profiling risk and sustainability in coastal deltas of the world. Science 349, 638–643 (2015).CAS 
Article 

Google Scholar 
6.Palmer, M. A. Water resources: beyond infrastructure. Nature 467, 534–535 (2010).CAS 
Article 

Google Scholar 
7.Michalak, A. M. Study role of climate change in extreme threats to water quality. Nature 535, 349–350 (2016).CAS 
Article 

Google Scholar 
8.McDonald, R. I., Weber, K. F., Padowski, J., Boucher, T. & Shemie, D. Estimating watershed degradation over the last century and its impact on water-treatment costs for the world’s large cities. Proc. Natl Acad. Sci. USA 113, 9117–9122 (2016).CAS 
Article 

Google Scholar 
9.Romulo, C. L. et al. Global state and potential scope of investments in watershed services for large cities. Nat. Commun. 9, 4375 (2018).Article 
CAS 

Google Scholar 
10.Tellman, B. et al. Opportunities for natural infrastructure to improve urban water security in Latin America. PLoS ONE 13, e0209470 (2018).Article 

Google Scholar 
11.United Nations World Water Assessment Programme/UN-Water The United Nations World Water Development Report 2018: Nature-Based Solutions for Water (UNESCO, 2018).12.Palmer, M. A., Liu, J., Matthews, J. H., Mumba, M. & D’Odorico, P. Manage water in a green way. Science 349, 584–585 (2015).CAS 
Article 

Google Scholar 
13.Ziv, G., Baran, E., Nam, S., Rodríguez-Iturbe, I. & Levin, S. A. Trading-off fish biodiversity, food security, and hydropower in the Mekong River Basin. Proc. Natl Acad. Sci. USA 109, 5609–5614 (2012).CAS 
Article 

Google Scholar 
14.Harrison, I. J. et al. Protected areas and freshwater provisioning: a global assessment of freshwater provision, threats and management strategies to support human water security. Aquat. Conserv. Mar. Freshw. Ecosyst. 26, 103–120 (2016).Article 

Google Scholar 
15.The World Database on Protected Areas (IUCN and UNEP-WCMC, 2017); http://www.protectedplanet.net16.Huber-Stearns, H. R., Goldstein, J. H., Cheng, A. S. & Toombs, T. P. Institutional analysis of payments for watershed services in the western United States. Ecosyst. Serv. 16, 83–93 (2015).Article 

Google Scholar 
17.Moran, E. F., Lopez, M. C., Moore, N., Müller, N. & Hyndman, D. W. Sustainable hydropower in the 21st century. Proc. Natl Acad. Sci. USA 115, 11891–11898 (2018).CAS 
Article 

Google Scholar 
18.Zheng, H. et al. Benefits, costs, and livelihood implications of a regional payment for ecosystem service program. Proc. Natl Acad. Sci. USA 110, 16681–16686 (2013).CAS 
Article 

Google Scholar 
19.Adamowicz, W. et al. Assessing ecological infrastructure investments. Proc. Natl Acad. Sci. USA 116, 201802883 (2019).Article 
CAS 

Google Scholar 
20.McDonald R. I. Conservation for Cities: How to Plan & Build Natural Infrastructure (Island Press, 2015).21.Grill, G. et al. An index-based framework for assessing patterns and trends in river fragmentation and flow regulation by global dams at multiple scales. Environ. Res. Lett. 10, 015001 (2015).Article 

Google Scholar 
22.Poff, N. L. & Schmidt, J. C. How dams can go with the flow. Science 353, 1099–1100 (2016).CAS 
Article 

Google Scholar 
23.Liu, J. & Yang, W. Integrated assessments of payments for ecosystem services programs. Proc. Natl Acad. Sci. USA 110, 16297–16298 (2013).CAS 
Article 

Google Scholar 
24.Muller, M., Biswas, A., Martin-Hurtado, R. & Tortajada, C. Built infrastructure is essential. Science 349, 585–586 (2015).CAS 
Article 

Google Scholar 
25.Veldkamp, T. I. E. et al. Water scarcity hotspots travel downstream due to human interventions in the 20th and 21st century. Nat. Commun. 8, 15697 (2017).CAS 
Article 

Google Scholar 
26.Cohen, S., Kettner, A. J. & Syvitski, J. P. M. Global suspended sediment and water discharge dynamics between 1960 and 2010: continental trends and intra-basin sensitivity. Glob. Planet. Change 115, 44–58 (2014).Article 

Google Scholar 
27.Dottori, F. et al. Development and evaluation of a framework for global flood hazard mapping. Adv. Water Resour. 94, 87–102 (2016).Article 

Google Scholar 
28.Byers L. et al. A Global Database of Power Plants (World Resources Institute, 2018); https://www.wri.org/publication/global-power-plant-database29.Liu, J. Integration across a metacoupled world. Ecol. Soc. 22, 29 (2017).Article 

Google Scholar 
30.Vercruysse, K., Grabowski, R. C. & Rickson, R. J. Suspended sediment transport dynamics in rivers: multi-scale drivers of temporal variation. Earth Sci. Rev. 166, 38–52 (2017).Article 

Google Scholar 
31.Wu, X.-X., Gu, Z.-J., Luo, H., Shi, X.-Z. & Yu, D.-S. Analyzing forest effects on runoff and sediment production using leaf area index. J. Mt. Sci. 11, 119–130 (2014).Article 

Google Scholar 
32.Wang, Y. et al. Annual runoff and evapotranspiration of forestlands and non-forestlands in selected basins of the Loess Plateau of China. Ecohydrology 4, 277–287 (2011).CAS 
Article 

Google Scholar 
33.Bilotta, G. S. & Brazier, R. E. Understanding the influence of suspended solids on water quality and aquatic biota. Water Res. 42, 2849–2861 (2008).CAS 
Article 

Google Scholar 
34.Stickler, C. M. et al. Dependence of hydropower energy generation on forests in the Amazon Basin at local and regional scales. Proc. Natl Acad. Sci. USA 110, 9601–9606 (2013).CAS 
Article 

Google Scholar 
35.Maltby, E. & Acreman, M. C. Ecosystem services of wetlands: pathfinder for a new paradigm. Hydrol. Sci. J. 56, 1341–1359 (2011).Article 

Google Scholar 
36.Shuster, W. D., Bonta, J., Thurston, H., Warnemuende, E. & Smith, D. R. Impacts of impervious surface on watershed hydrology: a review. Urban Water J. 2, 263–275 (2005).Article 

Google Scholar 
37.Borrelli, P. et al. Land use and climate change impacts on global soil erosion by water (2015–2070). Proc. Natl Acad. Sci. USA 117, 21994–22001 (2020).CAS 
Article 

Google Scholar 
38.Blöschl, G. et al. Changing climate both increases and decreases European river floods. Nature 573, 108–111 (2019).Article 
CAS 

Google Scholar 
39.Symes, W. S., Rao, M., Mascia, M. B. & Carrasco, L. R. Why do we lose protected areas? Factors influencing protected area downgrading, downsizing and degazettement in the tropics and subtropics. Glob. Change Biol. 22, 656–665 (2016).Article 

Google Scholar 
40.Kremen, C. & Merenlender, A. M. Landscapes that work for biodiversity and people. Science 362, eaau6020 (2018).Article 
CAS 

Google Scholar 
41.Dinerstein, E. et al. A global deal for nature: guiding principles, milestones, and targets. Sci. Adv. 5, eaaw2869 (2019).CAS 
Article 

Google Scholar 
42.Liu, J. et al. China’s environment on a metacoupled planet. Annu. Rev. Environ. Resour. 43, 1–34 (2018).CAS 
Article 

Google Scholar 
43.Viña, A., McConnell, W. J., Yang, H., Xu, Z. & Liu, J. Effects of conservation policy on China’s forest recovery. Sci. Adv. 2, e1500965 (2016).Article 

Google Scholar 
44.Chen, C. et al. China and India lead in greening of the world through land-use management. Nat. Sustain. 2, 122–129 (2019).Article 

Google Scholar 
45.Ouyang, Z. et al. Improvements in ecosystem services from investments in natural capital. Science 352, 1455–1459 (2016).CAS 
Article 

Google Scholar 
46.Vörösmarty, C. J. et al. Ecosystem-based water security and the Sustainable Development Goals (SDGs). Ecohydrol. Hydrobiol. 18, 317–333 (2018).Article 

Google Scholar 
47.Liu, J. et al. Nexus approaches to global sustainable development. Nat. Sustain. 1, 466–476 (2018).Article 

Google Scholar 
48.Flörke, M., Schneider, C. & McDonald, R. I. Water competition between cities and agriculture driven by climate change and urban growth. Nat. Sustain. 1, 51–58 (2018).Article 

Google Scholar 
49.McDonald, R. I. et al. Urban growth, climate change, and freshwater availability. Proc. Natl Acad. Sci. USA 108, 6312–6317 (2011).CAS 
Article 

Google Scholar 
50.Willner, S. N., Otto, C. & Levermann, A. Global economic response to river floods. Nat. Clim. Change 8, 594–598 (2018).Article 

Google Scholar 
51.Cattaneo, A., Nelson, A. & McMenomy, T. Global mapping of urban–rural catchment areas reveals unequal access to services. Proc. Natl Acad. Sci. USA 118, e2011990118 (2021).CAS 
Article 

Google Scholar 
52.Zarfl, C., Lumsdon, A. E., Berlekamp, J., Tydecks, L. & Tockner, K. A global boom in hydropower dam construction. Aquat. Sci. 77, 161–170 (2015).Article 

Google Scholar 
53.Schneider, A., Friedl, M. A. & Potere, D. A new map of global urban extent from MODIS satellite data. Environ. Res. Lett. 4, 044003 (2009).Article 

Google Scholar 
54.Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from spaceborne elevation data. EOS 89, 93–94 (2008).55.Yang, H. et al. A global assessment of the impact of individual protected areas on preventing forest loss. Sci. Total Environ. 777, 145995 (2021).CAS 
Article 

Google Scholar 
56.Smith, A. et al. New estimates of flood exposure in developing countries using high-resolution population data. Nat. Commun. 10, 1814 (2019).Article 
CAS 

Google Scholar 
57.Best, J. Anthropogenic stresses on the world’s big rivers. Nat. Geosci. 12, 7–21 (2019).CAS 
Article 

Google Scholar 
58.Hanasaki, N. et al. An integrated model for the assessment of global water resources—Part 1: model description and input meteorological forcing. Hydrol. Earth Syst. Sci. 12, 1007–1025 (2008).Article 

Google Scholar 
59.Bondeau, A. et al. Modelling the role of agriculture for the 20th century global terrestrial carbon balance. Glob. Change Biol. 13, 679–706 (2007).Article 

Google Scholar 
60.Pokhrel, Y. N. et al. Incorporation of groundwater pumping in a global Land Surface Model with the representation of human impacts. Water Resour. Res. 51, 78–96 (2015).Article 

Google Scholar 
61.Wada, Y., Wisser, D. & Bierkens, M. F. P. Global modeling of withdrawal, allocation and consumptive use of surface water and groundwater resources. Earth Syst. Dyn. 5, 15–40 (2014).Article 

Google Scholar 
62.Müller Schmied, H. et al. Variations of global and continental water balance components as impacted by climate forcing uncertainty and human water use. Hydrol. Earth Syst. Sci. 20, 2877–2898 (2016).Article 

Google Scholar 
63.Sheffield, J., Goteti, G. & Wood, E. F. Development of a 50-year high-resolution global dataset of meteorological forcings for land surface modeling. J. Clim. 19, 3088–3111 (2006).Article 

Google Scholar 
64.Dirmeyer, P. A. et al. GSWP-2: multimodel analysis and implications for our perception of the land surface. Bull. Am. Meteorol. Soc. 87, 1381–1398 (2006).Article 

Google Scholar 
65.Weedon, G. P. et al. The WFDEI meteorological forcing data set: WATCH forcing data methodology applied to ERA-Interim reanalysis data. Water Resour. Res. 50, 7505–7514 (2014).Article 

Google Scholar 
66.Bingham, H. C. et al. Sixty years of tracking conservation progress using the World Database on Protected Areas. Nat. Ecol. Evol. 3, 737–743 (2019).Article 

Google Scholar 
67.ArcGIS Desktop: Release 10.3.1 (Environmental Systems Research Institution, 2015).68.Domisch, S., Amatulli, G. & Jetz, W. Near-global freshwater-specific environmental variables for biodiversity analyses in 1 km resolution. Sci. Data 2, 150073 (2015).CAS 
Article 

Google Scholar 
69.Bennett, G. & Ruef, F. Alliances for Green Infrastructure: State of Watershed Investment 2016 (Forest Trends, 2016).70.R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2017).71.Wellman, B. & Frank, K. in Social Capital: Theory and Research (eds Lin, N. et al.) 233–273 (Routledge, 2001).72.Bates, D., Machler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar  More

1.Barnosky, A. D. et al. Has the Earth’s sixth mass extinction already arrived?. Nature 471, 51–57 (2011).CAS 
PubMed 
Article 
ADS 
PubMed Central 

Google Scholar 
2.Zaffos, A., Finnegan, S. & Peters, S. E. Plate tectonic regulation of global marine animal diversity. Proc. Natl. Acad. Sci. 114, 5653–5658 (2017).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
3.Claramunt, S. & Cracraft, J. A new time tree reveals Earth historys imprint on the evolution of modern birds. Sci. Adv. 1, e1501005 (2015).PubMed 
PubMed Central 
Article 
ADS 
CAS 

Google Scholar 
4.Leprieur, F., Descombes, P., Gaboriau, T., Cowman, P. F. & Parravicini, V. Plate tectonics drive tropical reef biodiversity dynamics. Nat. Commun. 7, 11461 (2016).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
5.Ficetola, G. F., Mazel, F. & Thuiller, W. Global determinants of zoogeographical boundaries. Nat. Ecol. Evol. 1, 0089 (2017).Article 

Google Scholar 
6.Mazel, F. et al. Global patterns of β-diversity along the phylogenetic time-scale: The role of climate and plate tectonics. Glob. Ecol. Biogeogr. 26, 1211–1221 (2017).Article 

Google Scholar 
7.Hofreiter, M. & Stewart, J. Ecological change, range fluctuations and population dynamics during the pleistocene. Curr. Biol. 19, R584–R594 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Pellissier, L. et al. Quaternary coral reef refugia preserved fish diversity. Science 344, 1016–1019 (2014).CAS 
PubMed 
Article 
ADS 
PubMed Central 

Google Scholar 
9.Jaramillo, C. et al. Effects of rapid global warming at the paleocene-eocene boundary on neotropical vegetation. Science 330, 957–961 (2010).CAS 
PubMed 
Article 
ADS 
PubMed Central 

Google Scholar 
10.Svenning, J.-C., Eiserhardt, W. L., Normand, S., Ordonez, A. & Sandel, B. The influence of paleoclimate on present-day patterns in biodiversity and ecosystems. Annu. Rev. Ecol. Evol. Syst. 46, 551–572 (2015).Article 

Google Scholar 
11.Steeman, M. E. et al. Radiation of extant cetaceans driven by restructuring of the oceans. Syst. Biol. 58, 573–585 (2009).PubMed 
PubMed Central 
Article 

Google Scholar 
12.Antonelli, A. & Sanmartín, I. Mass Extinction, gradual cooling, or rapid radiation? reconstructing the spatiotemporal evolution of the ancient angiosperm genus hedyosmum (Chloranthaceae) using empirical and simulated approaches. Syst. Biol. 60, 596–615 (2011).PubMed 
Article 
PubMed Central 

Google Scholar 
13.Nee, S., May, R. M. & Harvey, P. H. The reconstructed evolutionary process. Philos. Trans. R. Soc. Lond. B 344, 305–311 (1994).CAS 
Article 
ADS 

Google Scholar 
14.Morlon, H., Parsons, T. L. & Plotkin, J. B. From the cover: Reconciling molecular phylogenies with the fossil record. Proc. Natl. Acad. Sci. 108, 16327–16332 (2011).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
15.Condamine, F. L., Rolland, J. & Morlon, H. Macroevolutionary perspectives to environmental change. Ecol. Lett. 16, 72–85 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
16.Condamine, F. L. et al. Deciphering the evolution of birdwing butterflies 150 years after Alfred Russel Wallace Deciphering the evolution of birdwing butterflies 150 years after. Sci. Rep. 5, 11860 (2015).PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
17.Lagomarsino, L. P., Condamine, F. L., Antonelli, A., Mulch, A. & Davis, C. C. The abiotic and biotic drivers of rapid diversification in Andean bellflowers (Campanulaceae). New Phytol. 210, 1430–1442 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
18.Rolland, J. & Condamine, F. L. The contribution of temperature and continental fragmentation to amphibian diversification. J. Biogeogr. 46, 1857–1873 (2019).Article 

Google Scholar 
19.Gaboriau, T. et al. Ecological constraints coupled with deep-time habitat dynamics predict the latitudinal diversity gradient in reef fishes. Proc. R. Soc. B Biol. Sci. 286, 20191506 (2019).Article 

Google Scholar 
20.Descombes, P. et al. Linking species diversification to palaeo-environmental changes: A process-based modelling approach. Glob. Ecol. Biogeogr. 00, 1–12 (2017).
Google Scholar 
21.Rangel, T. F. et al. Modeling the ecology and evolution of biodiversity: Biogeographical cradles, museums, and graves. Science 361, 5452 (2018).Article 
CAS 

Google Scholar 
22.Pontarp, M. et al. The latitudinal diversity gradient: Novel understanding through mechanistic eco-evolutionary models. Trends Ecol. Evol. 34, 211–223 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
23.Cowman, P. F. Historical factors that have shaped the evolution of tropical reef fishes: A review of phylogenies, biogeography, and remaining questions. Front. Genet. 5, 1–15 (2014).Article 

Google Scholar 
24.Bellwood, D. R. et al. Evolutionary history of the butterflyfishes (f: Chaetodontidae) and the rise of coral feeding fishes. J. Evol. Biol. 23, 335–349 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Cowman, P. F. & Bellwood, D. R. Coral reefs as drivers of cladogenesis: Expanding coral reefs, cryptic extinction events, and the development of biodiversity hotspots. J. Evol. Biol. 24, 2543–2562 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Sorenson, L., Santini, F., Carnevale, G. & Alfaro, M. E. A multi-locus timetree of surgeonfishes (Acanthuridae, Percomorpha), with revised family taxonomy. Mol. Phylogenet. Evol. 68, 150–160 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
27.Dornburg, A., Moore, J., Beaulieu, J. M., Eytan, R. I. & Near, T. J. The impact of shifts in marine biodiversity hotspots on patterns of range evolution: Evidence from the Holocentridae (squirrelfishes and soldierfishes). Evolution 69, 146–161 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
28.Cowman, P. F. & Bellwood, D. R. The historical biogeography of coral reef fishes: Global patterns of origination and dispersal. J. Biogeogr. 40, 209–224 (2013).Article 

Google Scholar 
29.Lohman, D. J. et al. Biogeography of the Indo-Australian archipelago. Annu. Rev. Ecol. Evol. Syst. 42, 205–226 (2011).Article 

Google Scholar 
30.Gaboriau, T., Leprieur, F., Mouillot, D. & Hubert, N. Influence of the geography of speciation on current patterns of coral reef fish biodiversity across the Indo-Pacific. Ecography 41, 1295–1306 (2017).Article 

Google Scholar 
31.McManus, J. W. Marine speciation, tectonics and sea- level changes in Southeast Asia. Proc. Fifth Int. Coral Reef 4, 133–138 (1985).
Google Scholar 
32.Potts, D. C. Sea-level fluctuations and speciation in Scleractinia. Proc. Fifth Int. Coral Reef 4, 51–62 (1985).
Google Scholar 
33.Hou, Z. & Li, S. Tethyan changes shaped aquatic diversification. Biol. Rev. https://doi.org/10.1111/brv.12376 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
34.Stadler, T. Mammalian phylogeny reveals recent diversification rate shifts. Proc. Natl. Acad. Sci. USA 108, 6187–6192 (2011).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
35.Bellwood, D. R. & Wainwright, P. C. The history and biogeography of Fishes on Coral Reefs. in Coral Reef Fishes, Dynamics and Diversity in a Complex Ecosystem, 5–32 (2002).36.Williams, S. T. & Duda, T. F. Did tectonic activity stimulate Oligo-Miocene speciation in the Indo-West Pacific?. Evolution 62, 1618–1634 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
37.Renema, W. et al. Hopping hotspots: Global shifts in marine biodiversity. Science 321, 654–657 (2008).CAS 
PubMed 
Article 
ADS 
PubMed Central 

Google Scholar 
38.Tea, Y.-K. et al. Phylogenomic analysis of concatenated ultraconserved elements reveals the recent evolutionary radiation of the fairy wrasses (teleostei: labridae: cirrhilabrus). Syst. Biol. 1, 1–12 (2021).
Google Scholar 
39.Hall, R. Southeast Asia’s changing palaeogeography. Blumea J. Plant Taxon. Plant Geogr. 54, 148–161 (2009).Article 

Google Scholar 
40.Keith, S. A., Baird, A. H., Hughes, T. P., Madin, J. S. & Connolly, S. R. Faunal breaks and species composition of Indo-Pacific corals: The role of plate tectonics, environment and habitat distribution. Proc. Biol. Sci. 280, 20130818 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
41.Cowman, P. F. & Bellwood, D. R. Vicariance across major marine biogeographic barriers: Temporal concordance and the relative intensity of hard versus soft barriers. Proc. Biol. Sci. 280, 20131541 (2013).PubMed 
PubMed Central 

Google Scholar 
42.Price, S. A., Claverie, T., Near, T. J. & Wainwright, P. C. Phylogenetic insights into the history and diversification of fishes on reefs. Coral Reefs 34, 997–1009 (2015).Article 
ADS 

Google Scholar 
43.Bellwood, D. R., Goatley, C. H. R. & Bellwood, O. The evolution of fishes and corals on reefs: Form, function and interdependence. Biol. Rev. 92, 878–901 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
44.Bowen, B. W., Rocha, L. A., Toonen, R. J. & Karl, S. A. The origins of tropical marine biodiversity. Trends Ecol. Evol. 28, 359–366 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
45.Price, S. A., Tavera, J. J., Near, T. J. & Wainwright, P. C. Elevated rates of morphological and functional diversification in reef-dwelling haemulid fishes. Evolution 67, 417–428 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
46.Kiessling, W., Simpson, C., Beck, B., Mewis, H. & Pandolfi, J. M. Equatorial decline of reef corals during the last Pleistocene interglacial. Proc. Natl. Acad. Sci. USA. 109, 21378–21383 (2012).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
47.Floeter, S. R. et al. Atlantic reef fish biogeography and evolution. J. Biogeogr. 35, 22–45 (2007).
Google Scholar 
48.Riginos, C., Buckley, Y. M., Blomberg, S. P. & Treml, E. A. Dispersal capacity predicts both population genetic structure and species richness in reef fishes. Am. Nat. 184, 52–64 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
49.Rocha, L. A. & Bowen, B. W. Speciation in coral-reef fishes. J. Fish Biol. 72, 1101–1121 (2008).Article 

Google Scholar 
50.Tedesco, P. A., Paradis, E., EvEque, C. L. & Hugueny, B. Explaining global-scale diversification patterns in actinopterygian fishes. J. Biogeogr. 44, 773–783 (2016).Article 

Google Scholar 
51.Rosenzweig, M. L. Species Diversity in Space and Time (Springer, 1995).Book 

Google Scholar 
52.Kisel, Y. & Barraclough, T. G. Speciation has a spatial scale that depends on levels of gene flow. Am. Nat. 175, 316–334 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
53.Fine, P. V. A. & Ree, R. H. Evidence for a time-integrated species-area effect on the latitudinal gradient in tree diversity. Am. Nat. 168, 796–804 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
54.Jetz, W. & Fine, P. V. A. Global gradients in vertebrate diversity predicted by historical area-productivity dynamics and contemporary environment. PLoS Biol. 10, e1001292 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
55.Konow, N., Price, S., Abom, R., Bellwood, D. & Wainwright, P. Decoupled diversification dynamics of feeding morphology following a major functional innovation in marine butterflyfishes. Proc. Biol. Sci. 284, 20170906 (2017).PubMed 
PubMed Central 

Google Scholar 
56.Clements, K. D., German, D. P., Piché, J., Tribollet, A. & Choat, J. H. Integrating ecological roles and trophic diversification on coral reefs: Multiple lines of evidence identify parrotfishes as microphages. Biol. J. Linn. Soc. 120, 729–751 (2017).
Google Scholar 
57.Lobato, F. L. et al. Diet and diversification in the evolution of coral reef fishes. PLoS ONE 9, e102094 (2014).PubMed 
PubMed Central 
Article 
ADS 
CAS 

Google Scholar 
58.Siqueira, A. C., Morais, R. A., Bellwood, D. R. & Cowman, P. F. Trophic innovations fuel reef fish diversification. Nat. Commun. 11, 1–11 (2020).Article 
CAS 

Google Scholar 
59.Louca, S. & Pennell, M. W. Extant timetrees are consistent with a myriad of diversification histories. Nature 580, 502–505 (2020).CAS 
PubMed 
Article 
ADS 
PubMed Central 

Google Scholar 
60.Morlon, H., Hartig, F. & Robin, S. Prior hypotheses or regularization allow inference of diversification histories from extant timetrees. bioRxiv (2020).61.McCord, C. L. & Westneat, M. W. Phylogenetic relationships and the evolution of BMP4 in triggerfishes and filefishes (Balistoidea). Mol. Phylogenet. Evol. 94, 397–409 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
62.Santini, F. & Carnevale, G. First multilocus and densely sampled timetree of trevallies, pompanos and allies (Carangoidei, Percomorpha) suggests a Cretaceous origin and Eocene radiation of a major clade of piscivores. Mol. Phylogenet. Evol. 83, 33–39 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
63.Santini, F., Carnevale, G. & Sorenson, L. First multi-locus timetree of seabreams and porgies (Percomorpha: Sparidae). Ital. J. Zool. 81, 55–71 (2014).Article 

Google Scholar 
64.Müller, R. D., Sdrolias, M., Gaina, C., Steinberger, B. & Heine, C. Long-term sea-level fluctuations driven by ocean basin dynamics. Science 319, 1357–1362 (2008).PubMed 
Article 
ADS 
CAS 
PubMed Central 

Google Scholar 
65.Heine, C., Yeo, L. G. & Müller, R. D. Evaluating global paleoshoreline models for the Cretaceous and Cenozoic. Aust. J. Earth Sci. 62, 275–287 (2015).CAS 

Google Scholar 
66.Kleypas, J. A. & Mcmanus, J. W. Environmental Limits to Coral Reef Development : Where Do We Draw the Line ?. Am. Zool. 39, 146–159 (1999).Article 

Google Scholar 
67.Bugayevskiy, L. M. & Snyder, J. P. Map Projections: A Reference Manual (Taylor & Francis, London, 1995).
Google Scholar 
68.Chang, J., Rabosky, D. L. & Alfaro, M. E. Estimating diversification rates on incompletely sampled phylogenies: Theoretical concerns and practical solutions. Syst. Biol. 69, 602–611 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
69.Morlon, H. et al. RPANDA: An R package for macroevolutionary analyses on phylogenetic trees. Methods Ecol. Evol. 7, 589–597 (2016).Article 

Google Scholar  More

Detailed population structure of L1–4 and a hierarchical sub-lineage naming systemWe assembled a high-quality data set of whole genomes, antibiotic resistance phenotypes, and geographic sites of isolation for 9584 clinical Mtb samples (“Methods” section and Supplementary Data 1). Of the total, 4939 (52%) were pan-susceptible, i.e., susceptible to at least isoniazid and rifampicin (and all other antibiotics when additional phenotypic data were available), and 4645 (48%) were resistant to one or more antibiotics (Supplementary Fig. 1a). Using the 62 SNS lineage barcode6, 738 isolates were classified as L1 (8%), 2193 as L2 (22%), 1104 as L3 (12%) and 5549 as L4 (58%, Supplementary Fig. 1b). Among the 4939 pan-susceptible isolates, we identified high-quality genome-wide SNSs (83,735 for L1, 56,736 for L2, 76,817 for L3, and 185,622 for L4) that we used in building maximum-likelihood phylogenies for each major lineage (L1–4, “Methods” section). We computed an index of genetic divergence (FST) between groups defined by each bifurcation in each phylogeny. Sub-lineages were defined as monophyletic groups that had high FST ( >0.33) and were also clearly separated from other groups in principal component analysis (PCA, see “Methods” section). We also defined internal groups to sub-lineages (see “Methods” section): an internal group is a monophyletic group genetically divergent (by FST and PCA) from its neighboring groups, but has one or more ancestral branches that show a low degree of divergence or low support (bootstrap values). Internal groups do not represent true sub-lineages in a hierarchical fashion, but defining them allows us to further characterize the Mtb population structure. We provide code to automate all the steps described above. Our approach is scalable and can be used on other organisms (see “Methods” section).To better classify Mtb isolates in the context of the global Mtb population structure, we developed a hierarchical sub-lineage naming scheme (Supplementary Data 2) where each subdivision in the classification corresponds to a split in the phylogenetic tree of each major Mtb lineage. Starting with the global Mtb lineage numbers (e.g., L1), we recursively introduced a subdivision (e.g., from 1.2 to 1.2.1 and 1.2.2) at each bifurcation of the phylogenetic tree whenever both subclades sufficiently diverged. Formally, we defined these splits using bootstrap criteria, and independent validations by FST and PCA (see “Methods” section). Internal groups were denoted with the letter “i” (e.g., 4.1.i1). This proposed system overcomes two major shortcomings of the existing schemas: same-level sub-lineages are never overlapping (unlike the system of Stucki et al.8 sub-lineage 4.10 includes sub-lineages 4.7–4.9), and the names reflect both phylogenetic relationships and genetic similarity (unlike semantic naming such as the “Asia ancestral” lineage in the system of Shitikov et al.7). Further, this naming system can be standardized to automate the process of lineage definition. These advantages come at the price of long sub-lineage names in the case of complex phylogenies (e.g., for L4, sub-lineage 4.10 gets the lineage designation 4.2.1.1.1.1.1.1). For compatibility with naming conventions already in use and to keep names as short as possible, we designed a second, shorthand, naming system that expands the Coll et al. lineage schema by adding new subdivisions and differentiating between sub-lineages and internal groups. For instance, sub-lineage 4.3.1 is designated as 4.3.i1, informing the user that this is an internal group of sub-lineage 4.3. To simplify the use of the hierarchical naming schema and the updated shorthand schema, we provide a table that compares them side by side along with naming systems currently in use (Supplementary Data 2).Using the sub-lineage definition rules and the sub-lineage naming scheme described above, we characterized six previously undescribed sub-lineages of L1 (Fig. 1 and Supplementary Fig. 2); five of which expand the current description of 1.2. We also detected an internal group of 91 isolates (1.1.3.i1) characterized by a long defining branch in the phylogeny (corresponding to 82 SNSs), a high FST (0.48), and geographically restricted to Malawi (85/91, 93% isolates, Fig. 1 and Supplementary Fig. 3). We estimated the date of the emergence of the MRCA of such a group (see “Methods” section) and we found it to be between 1497 and 1754. We found four previously undescribed sub-lineages of L3 (Fig. 2 and Supplementary Fig. 4), revising L3 into four main groups, whereas previously only two partitions of one sub-lineage were characterized (3.1). We found that the latter two partitions are in fact internal groups of the largest sub-lineage (3.1.1) in our revised classification.Fig. 1: Phylogenetic tree reconstruction of lineage 1 (binary tree).Gray circles define splits where the FST (fixation index) calculated using the descendants of the two children nodes is greater than 0.33. The sub-lineages are defined by colored areas (blue: sub-lineages already described in the literature; green: sub-lineages described here; purple: internal sub-lineages). Source data are provided as a Source Data file.Full size imageFig. 2: Phylogenetic tree reconstruction of lineage 3 (binary tree).Gray circles define splits where the FST (fixation index) calculated using the descendants of the two children nodes is greater than 0.33. The sub-lineages are defined by colored areas (green: sub-lineages described here; purple: internal sub-lineages). Source data are provided as a Source Data file.Full size imageL2 is divided into two groups: proto-Beijing and Beijing with the latter in turn partitioned into two groups: ancient- and modern-Beijing7. Each one of these groups is characterized by further subdivisions (three for the ancient-Beijing group and seven for the modern-Beijing group; see Supplementary Fig. 4). We found a new sub-lineage (2.2.1.2, Fig. 3, and Supplementary Fig. 5) within the previously characterized ancient-Beijing group. However, genetic diversity within the modern-Beijing group (2.2.1.1.1) was lower than in the other L2 sub-lineages and the tree topology and FST calculations did not support further hierarchical subdivisions. Although we did find three internal groups of modern-Beijing: two undescribed and one that corresponds to the Central Asia group7. For L4, our results support a complex population structure with 21 sub-lineages and 15 internal groups. In particular, we found 11 previously undescribed sub-lineages and 5 internal groups that expand our understanding of previously characterized sub-lineages (e.g., 4.2.2; 4.2 in the Coll et al. classification) or that were not characterized since these isolates were simply classified as L4 (e.g., 4.11, Fig. 4, and Supplementary Fig. 6) using the other barcodes.Fig. 3: Phylogenetic tree reconstruction of lineage 2 (binary tree).Gray circles define splits where the FST (fixation index) calculated using the descendants of the two children nodes is greater than 0.33. The sub-lineages are defined by colored areas (blue: sub-lineages already described in the literature; green: sub-lineages described here; purple: internal sub-lineages). Source data are provided as a Source Data file.Full size imageFig. 4: Phylogenetic tree reconstruction of lineage 4 (binary tree).Gray circles define splits where the FST (fixation index) calculated using the descendants of the two children nodes is greater than 0.33. The sub-lineages are defined by colored areas (blue: sub-lineages already described in the literature; green: sub-lineages described here; purple: internal sub-lineages). Source data are provided as a Source Data file.Full size imageA new barcode to define L1–4 Mtb sub-lineages and a software package to type Mtb strains from WGS dataWe defined a SNS barcode for distinguishing the obtained sub-lineages (Supplementary Data 3). We characterized new synonymous SNSs found in 100% of isolates from a given sub-lineage, but not in other isolates from the same major lineage, compiling 95 SNSs into an expanded barcode (Supplementary Data 3). We validated the barcode by using it to call sub-lineages in the hold-out set of 4645 resistant isolates and comparing the resulting sub-lineage designations with maximum-likelihood phylogenies inferred from the full SNS data (Supplementary Figs. 7–10). A sub-lineage was validated if it was found in the hold-out data and formed a monophyletic group in the phylogeny. Considering the “recent” sub-lineages, i.e., the most detailed level of classification in our system, we were able to validate eight out of nine L1 sub-lineages including five out of six of the new sub-lineages described here, with the exception of 1.1.1.2. We validated all four new L3 sub-lineages, all five L2 sub-lineages including the one previously undescribed, and 16 of the 21 L4 sub-lineages including two described here. The sub-lineages we could not confirm were not represented by any isolate in the validation phylogenies. We did not observe any paraphyletic sub-lineages in the revised classification system.We developed fast-lineage-caller, a software tool that classifies Mtb genomes using the SNS barcode proposed above. For a given genome, it returns the corresponding sub-lineage as output using our hierarchical naming system in addition to four other existing numerical/semantic naming systems, when applicable (see “Methods” section). The tool also informs the user on how many SNSs support a given lineage call and allows for filtering of low-quality variants. The tool is generalizable and can manage additional barcodes defined by the user to type the core genome of potentially any bacterial species.Geographic distribution of the Mtb sub-lineagesNext, we examined whether certain sub-lineages were geographically restricted, which would support the Mtb-human co-evolution hypothesis, or whether they constituted prevalent circulating sub-lineages in several different countries (i.e., geographically unrestricted)8. We used our SNS barcode to determine the sub-lineages of 17,432 isolates (see “Methods” section) sampled from 74 countries (Supplementary Fig. 11 and Supplementary Data 4, 5). We computed the Simpson diversity index (Sdi) as a measure of geographic diversity that controls for variable sub-lineage frequency (see “Methods” section) for each well-represented sub-lineage or internal group (n  > 20). We hypothesized that geographically unrestricted lineages would have a higher Sdi. We found Sdi to correlate highly (⍴ = 0.68; p-value = 5.7 × 10−7) with the number of continents from which a given sub-lineage was isolated (Supplementary Fig. 12). The Sdi ranged between a minimum of 0.05 and a maximum of 0.72, with a median value of 0.46 (Fig. 5). The known geographically restricted sub-lineages8 had an Sdi between 0.28 and 0.5 (Fig. 5 and Supplementary Table 1), while the known geographically unrestricted sub-lineages8,9 had an Sdi between 0.55 and 0.61 (Fig. 5 and Supplementary Table 2). We found 11 sub-lineages/internal groups with Sdi 0.61 (Supplementary Table 4), i.e., more extreme than previously reported geographically restricted or unrestricted sub-lineages, respectively.Fig. 5: Histogram of the Simpson diversity index calculated for sub-lineages of lineages 1–4.A data set of 17,432 isolates from 74 countries was used to perform this analysis. Yellow triangles designate the Simpson diversity index values of sub-lineages designated as geographically restricted by Stucki et al. Light gray circles designate the Simpson diversity index values of sub-lineages designated as geographically unrestricted by Stucki et al. Source data are provided as a Source Data file.Full size imageWhile the currently known geographically restricted sub-lineages are all in L4, we found evidence of geographic restriction for two sub-lineages/internal groups of L1. The first, the L1 internal group 1.1.3.i1, showed a very low Sdi (0.06) and was only found at high frequency among the circulating L1 isolates in Malawi (Fig. 6). This finding is also in agreement with the L1 phylogeny (Fig. 1) that shows a relatively long (82 SNS) branch defining this group. The second geographically restricted L1 sub-lineage is 1.1.1.1 (Sdi = 0.12) that was only found at high frequency among circulating L1 isolates in South-East Asia (Vietnam and Thailand, Fig. 7). To exclude the possibility that these two groups appeared geographically restricted as a result of oversampling transmission outbreaks, we calculated the distribution of the pairwise SNS distance for each of these two sub-lineages. We measured a median SNS distance of 204 and 401, respectively, refuting this kind of sampling error for these groups (typical pairwise SNS distance in outbreaks 0.67 and results on L4 transmissibility below.Differences in transmissibility between the Mtb global lineagesThe observation that some lineages/sub-lineages are more geographically widespread than others raises the question of whether this results from differences in marginal transmissibility across human populations. On a topological level, we observed L2 and L3 phylogenies to be qualitatively different from those of L1 and L4 (Figs. 1–4): displaying a star-like pattern with shorter internal branches and longer branches near the termini. We confirmed this quantitatively by generating a single phylogenetic tree for all 9584 L1–4 isolates and plotting cumulative branch lengths from root to tip for each main lineage (Supplementary Fig. 20). Star-like topologies have been postulated to associate with rapid or effective viral or bacterial transmission e.g., a “super-spreading” event in outbreak contexts25. To compare transmissibility between the four lineages, we compared the distributions of terminal branch lengths expecting a skew toward shorter terminal branch lengths supporting the idea of higher transmissibility. We found L4 to have the shortest median terminal branch length, followed in order by L2, L3, and L1 (medians: 6.2 × 10−5, 8.2 × 10−5, 10.2 × 10−5, 17.5 × 10−5, respectively; all pairwise two-sided Wilcoxon rank-sum tests significant p-value < 0.001; Fig. 9). Shorter internal node-to-tip distance is a second phylogenetic correlate of transmissibility; the distribution of this measure across the four lineages revealed a similar hierarchy to the terminal branch length distribution (Supplementary Fig. 21). We also computed the cumulative distribution of isolates separated by increasing total pairwise SNS distance (Supplementary Fig. 22). The proportion of L4 isolates separated by More

Review of C–R stability theoryTo set the context for how the R–N module will be used to understand the dynamics of nutrient-limited ecosystem models, we first briefly review stability results from modular food web theory. We do this by laying out a set of examples that serve to illustrate that in general, strong C–R interactions promote oscillatory dynamics while carefully placed weak C–R interactions dampen them5. We begin with the Rosenzweig–MacArthur C–R system as our base C–R module (Fig. 1a). It is biologically supported and produces a range of biologically plausible dynamics5, making it an appropriate system for this analysis. It exhibits three different dynamical phases over a gradient of interaction strengths (energetically defined sensu Nilsson et al. 2018) such that increasing the attack rate (({a}_{{CR}})) increases interaction strength15 (Fig. 2). We use the return time after a small perturbation (i.e., eigenvalues) to highlight the natural stability trade-off that occurs as interaction strength is changed, (i.e., the “checkmark” stability pattern)5,6. Equations and parameters can be found in Supplementary Results 1A. We draw your attention to three notable dynamical phases of the C–R module. At low interaction strengths the dominant eigenvalue (({lambda }_{{max }})) is negative and real and the C–R module follows a monotonic return to a stable equilibrium (Fig. 2a). During this phase ({lambda }_{{max }}) decreases from 0 (i.e., where ({a}_{{{CR}}}) allows the consumer to persist) to more negative values and thus stronger interactions tend to increase stability (Fig. 2d, i). At moderate interaction strengths, there is a sudden shift to eigenvalues with a non-zero complex part and population dynamics overshoot the equilibrium (Fig. 2b). Increases in interaction strength then further excite population dynamics and we observe less stable dynamics across this phase (Fig. 2d, ii). Last, the system reaches a Hopf bifurcation where the dominant eigenvalue becomes positive, yielding sustained cycles or oscillations (Fig. 2c, d, iii). As interaction strength increases across this phase, it is difficult to determine stability from the magnitude of a positive dominant eigenvalue; however, destabilization with increased interaction strength is readily observed in that the cycles become increasingly larger oscillations with a high coefficient of variation (CV)5. Note that while the Rosenzweig–MacArthur C–R system is shown here under a single set of parameters, analysis of the Jacobian shows the qualitative results to be general5. Moreover, the qualitative stability pattern remains for a type I and type III functional response5.Fig. 1: C–R and R–N base modules.a Rosenzweig–MacArthur C–R module modelled with Holling type II functional response and logistic resource growth, where (R) is resource biomass and (C) is consumer biomass. Parameters: (r) is the intrinsic growth rate of (R), (K) is the carrying capacity of (R), ({a}_{{mathrm {CR}}}) is the attack rate of (C) on (R), (e) is the assimilation rate of (C), ({R}_{0}) is the half-saturation density of (C), ({m}_{R}) and ({m}_{C}) are the mortality rates of (R) and (C), respectively. b R–N module modelled with a Monod nutrient uptake equation and external nutrient input, where (N) is a limiting-nutrient pool and (R) is the resource biomass. Parameters: ({I}_{N}) is external nutrient input to (N), ({a}_{{RN}}) is nutrient uptake rate by (R), (k) is the half-saturation density of (R), ({l}_{N}) and ({l}_{R}) are nutrient loss rates from (N) and (R), respectively.Full size imageFig. 2: C–R checkmark stability response.d Local stability (real and complex parts of the dominant eigenvalue; ({lambda }_{{max }})) as a function of interaction strength (({a}_{{{mathrm {CR}}}})) for the Rosenzweig–MacArthur C–R module. Time series reflect dynamics associated with region i, ii, and iii, respectively, following a perturbation that removes 50% of consumer biomass: a Stable equilibrium; monotonic dynamics. b Stable equilibrium; overshoot dynamics. c Unstable equilibrium; limit cycle. Boldness of arrows indicates the strength of interaction (({a}_{{CR}})).Full size imageWe now couple C–R modules into higher order food web modules to demonstrate how the addition of weak and/or strong interactions to a system can be used to predict dynamics at steady state (Fig. 3), constituting the “algebra” of C–R modules. Equations and parameters can be found in Supplementary Results 1B–D. We start with the three trophic level food chain (Fig. 3a), consisting of two coupled C–R modules (i.e., C1-R and P–C1). Theory has tended to find two weakly interacting C–R modules to generally produce locally stable equilibria16 (Fig. 3a). Increasing the strength of the C1–R interaction causes it to act like an oscillator (see Fig. 2c, above), and with enough increase this underlying oscillation is reflected in the limit cycles of the entire food chain (Fig. 3b). If the P–C1 interaction is strengthened as well, we end up with two coupled oscillators—the recipe for chaos17,18 (Fig. 3c). As such, coupled strong interactions are not surprisingly the recipe for complex and highly unstable dynamics.Fig. 3: Algebra of C–R modules.Time series showing the general dynamical outcomes for the food chain and diamond module at steady state with varied combinations of C–R interaction strengths. a Weak–weak interaction; point attractor. b Strong–weak interaction; limit cycle. c Strong–strong interaction; chaos. d Strong–strong, weak interaction; limit cycle. e Strong–strong, weak–weak interaction; point attractor.Full size imageFollowing McCann et al.19, we now add a weakly coupled consumer C2 to the food chain system of Fig. 3c. This weak consumer essentially draws energy away from the strong P–C1–R pathway and in doing so partially mutes the coupled oscillators, bringing the dynamics back to a more even limit cycle (Fig. 3d) and under certain conditions can drive equilibrium dynamics19. Last, the predator is weakly coupled to C2, creating a strong and weak pathway. The second weak interaction further draws energy away from the strong pathway, muting the oscillators entirely and bringing the system in this example to a point attractor (Fig. 3e). These examples show that well placed weak interactions (i.e., non-oscillatory phases, Fig. 2a, b) can be used to draw energy away from strong pathways and act as potent stabilizers of potentially oscillatory pathways. Note that weak interactions play a similarly stabilizing role in the omnivory module20 and further, weak interactions have been shown to stabilize large food web networks4,6 suggesting the principles derived from modular theory scale up to whole systems. Taken altogether, the oscillatory nature of strong C–R interactions generally promotes oscillatory dynamics in higher order systems, while the careful placement of weak C–R interactions—which are monotonic in nature—act to dampen oscillations. Although not discussed to our knowledge, we conjecture that if a subsystem exists such that strong interactions lead to monotonic dynamics (i.e., without oscillatory decay), strong interactions in this case would serve as a potent stabilizer. Below, we show the R–N module appears to be such a case.R–N module and stabilityTowards understanding how the R–N subsystem may interact in a higher order system, we first briefly consider the stability of the R–N module alone (akin to what we discussed for the C–R module above). The R–N module consists of a resource that takes up nutrients according to a Monod-like growth term, is open to flows from the external environment as a result of geochemical processes, and nutrients are lost to the external environment according to a linear term11 (Fig. 1b). Performing a local stability analysis about the interior equilibrium reveals the R–N module to be locally stable for all biologically feasible parameterizations, as determined by the signs of the trace and determinant of the Jacobian matrix (see Supplementary Results 2B). We now perform further numerical and analytical analyses to understand how stability is influenced by interaction strength.As the maximum rate of nutrient uptake (({a}_{{RN}})) is increased (i.e., R–N interaction strength), stability is generally increased (Fig. 4d), with the real part of the dominant eigenvalue (({lambda }_{{max }})) tending from 0 (i.e., where ({a}_{{RN}}) allows the resource to persist) towards an asymptote of ({-l}_{R}) (see Supplementary Results 2C). Numerical analysis reveals that the asymptote at ({-l}_{R}) can be approached from above or below depending on the relative leakiness of the R and N compartments (i.e., the rate at which nutrients are lost to the external environment from compartment R (({l}_{R})) and N (({l}_{N}))). For ({l}_{N} , > , {l}_{R}) (Fig. 4d), the R–N module only follows a monotonic return to equilibrium as interaction strength is increased, with increased interaction strength only tending to increased stability (i.e., reduce return time). For ({l}_{N} < {l}_{R}) (Fig. 4d), the R–N module follows a monotonic return to equilibrium for weak (Fig. 4a) and strong (Fig. 4c) interaction strength, but modest overshoot dynamics are observed for intermediate interaction strength (Fig. 4b). Stability tended to increase with interaction strength for weak to intermediate interaction strength (i.e., dominant eigenvalue becomes more negative), then slightly decrease as interaction strength became strong. A special case exists when ({l}_{R}={l}_{N}) (Fig. 4d), where stability increases with interaction strength until ({lambda }_{{max }}) becomes locked in at ({-l}_{R}), indicating stability does not change regardless of any further increase in interaction strength. Overall, the R–N interaction tends to generally stabilize in all cases (dominant eigenvalue goes from zero to a more negative saturating value with monotonic dynamics), although there are some intermediate cases that produce complex eigenvalues that suggest population dynamic overshoot potential (Fig. 4b). Note that we obtain qualitatively similar results when implicitly strengthening the R–N interaction by increasing nutrient loading (see Supplementary Results 2D and Supplementary Fig. 1). Now, given the above framework for coupled C–R modules—where weak C–R interactions with underlying monotonic dynamics dampen the oscillatory potential of strong C–R interactions—the underlying monotonic dynamics of the R–N module suggest that R–N interactions ought to be stabilizing when coupled to strong C–R interactions. Further, the underlying increase in stability (i.e., more rapid return to equilibrium) as R–N interaction strength is increased suggests the stabilizing potential of the R–N module ought to increase as the interaction becomes stronger.Fig. 4: R–N stability response to increasing interaction strength.Time series showing R density following a perturbation that lowered R density to 50% of equilibrium density for a low (({a}_{{RN}}=0.8)), b intermediate (({a}_{{RN}}=1)), and c high maximum rate of nutrient uptake (({a}_{{RN}}=2.8)). d Local stability (dominant eigenvalue; ({lambda }_{{max }})) of the R–N subsystem as ({a}_{{RN}}) is increased for ({l}_{N} , > , {l}_{R}), ({l}_{N}={l}_{R}), and ({l}_{N} < {l}_{R}), where ({l}_{R}) and ({l}_{N}) are the rate at which nutrients are lost to the external environment from compartment R and N, respectively. Solid lines are real parts and dashed lines are complex parts of ({lambda }_{{max }}).Full size imageTo look into this conjecture, we first coupled R–N to multiple configurations of strong and expectantly oscillatory C–R interactions and increased R–N interaction strength (({a}_{{RN}})). Following this, we added nutrient cycling and repeated the experiment to demonstrate that our results can be generalized to nutrient-limited ecosystem models. The full equations and parameter values for each model are listed in Supplementary Results 3A–D and 4A, B. We begin with the C–R–N system, where C–R and R–N are coupled through R (Fig. 5a). The initial increase in ({a}_{{RN}}) implicitly strengthens the C–R interaction and fuels the oscillatory potential of C–R and cycles emerge almost immediately after C is able to persist. As ({a}_{{RN}}) is increased further the cycles disappear and we obverse a steep stabilization phase, followed by a modest period of destabilization. Adding a weakly coupled predator gives a similar outcome, although the system continually stabilizes as ({a}_{{RN}}) is increased (Fig. 5b). If the P–C interaction is strengthened (i.e., both C1–R and P–C1 are strong, the recipe for chaos), R–N is unable to dampen oscillations even with a strong interaction strength, although a strong interaction gives tighter bound cycles than a weak interaction (Fig. 5c). We next add a weakly coupled consumer to the nutrient-limited food chain with strong P–C1 and C1–R interactions (Fig. 5d). As seen previously, this interaction draws energy out of the strong pathway, partially muting oscillatory potential. Thus, the ability for a strong R–N interaction to once again return the system to a stable equilibrium is not surprising. Finally, we add a detrital compartment to show that strong R–N interactions remain potent stabilizers in the context of nutrient cycling (Fig. 6b) when compared to a nutrient-limited food chain without nutrient cycling (Fig. 6a).Fig. 5: Nutrient-limited food chain stability.a–d Non-equilibrium dynamics (log10(C1,max/C1,min)) and equilibrium stability (real part of the dominant eigenvalue; ({lambda }_{{max }})) of the C–R–N, P–C–R–N with a single oscillator, P–C–R–N with coupled oscillators, and P–C1–C2–R–N modules, respectively, as ({a}_{{RN}}) is varied.Full size imageFig. 6: Nutrient-limited ecosystem module stability.a, b Non-equilibrium dynamics (log10(Cmax/Cmin)) and equilibrium stability (real part of the dominant eigenvalue; ({lambda }_{{max }})) of the C–R–N nutrient-limited food chain model and the C–R–N–D nutrient-limited ecosystem model, respectively, as ({a}_{{RN}}) is varied.Full size imageNote that we repeat our analysis of higher order modules by implicitly increasing R–N interaction strength through nutrient loading (see Supplementary Results 3E and 4C and Supplementary Figs. 2 and 3). In all cases, increased nutrient loading led to less stable dynamics, consistent with DeAngelis’ (1992) paradox of enrichment finding where increased nutrient loading lead to destabilizing autotroph–herbivore oscillations. More

1.Binney, H. et al. Vegetation of Eurasia from the last glacial maximum to present: key biogeographic patterns. Quat. Sci. Rev. 157, 80–97 (2017).ADS 
Article 

Google Scholar 
2.Clark, P. U. et al. The Last Glacial Maximum. Science 325, 710–714 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
3.Bigelow, N. H. Climate change and Arctic ecosystems: 1. Vegetation changes north of 55°N between the last glacial maximum, mid-Holocene, and present. J. Geophys. Res. 108, https://doi.org/10.1029/2002jd002558 (2003).4.Graham, R. W. et al. Timing and causes of mid-Holocene mammoth extinction on St. Paul Island, Alaska. Proc. Natl Acad. Sci. USA 113, 9310–9314 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Stuart, A. J. Late Quaternary megafaunal extinctions on the continents: a short review. Geol. J. 50, 338–363 (2015).Article 

Google Scholar 
6.Koch, P. L. & Barnosky, A. D. Late Quaternary extinctions: state of the debate. Ann. Rev. Ecol. Evol. Syst. 37, 215–250 (2006).Article 

Google Scholar 
7.Rabanus-Wallace, M. T. et al. Megafaunal isotopes reveal role of increased moisture on rangeland during late Pleistocene extinctions. Nat. Ecol. Evol. 1, 0125 (2017).Article 

Google Scholar 
8.Mann, D. H., Groves, P., Kunz, M. L., Reanier, R. E. & Gaglioti, B. V. Ice-age megafauna in Arctic Alaska: extinction, invasion, survival. Quat. Sci. Rev. 70, 91–108 (2013).ADS 
Article 

Google Scholar 
9.Capo, E. et al. Lake sedimentary DNA research on past terrestrial and aquatic biodiversity: overview and recommendations. Quaternary 4, https://doi.org/10.3390/quat4010006 (2021).10.Edwards, M. E. et al. Metabarcoding of modern soil DNA gives a highly local vegetation signal in Svalbard tundra. Holocene 28, 2006–2016 (2018).ADS 
Article 

Google Scholar 
11.Hughes, P. D., Gibbard, P. L. & Ehlers, J. Timing of glaciation during the last glacial cycle: evaluating the concept of a global ‘Last Glacial Maximum’ (LGM). Earth Sci. Rev. 125, 171–198 (2013).ADS 
Article 

Google Scholar 
12.Willerslev, E. et al. Fifty thousand years of Arctic vegetation and megafaunal diet. Nature 506, 47–51 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
13.Rasmussen, S. O. et al. A new Greenland ice core chronology for the last glacial termination. J. Geophys. Res. 111, https://doi.org/10.1029/2005jd006079 (2006).14.Mangerud, J. The discovery of the Younger Dryas, and comments on the current meaning and usage of the term. Boreas 50, 1–5 (2020).Article 

Google Scholar 
15.Bauska, T. K. et al. Carbon isotopes characterize rapid changes in atmospheric carbon dioxide during the last deglaciation. Proc. Natl Acad. Sci. USA 113, 3465–3470 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Wesser, S. D. & Armbruster, W. S. Species distribution controls across a forest‐steppe transition: a causal model and experimental test. Ecol. Monogr. 61, 323–342 (1991).Article 

Google Scholar 
17.Rijal, D. P. et al. Sedimentary ancient DNA shows terrestrial plant richness continuously increased over the Holocene in northern Fennoscandia. Sci. Adv. 7, eabf9557 (2021).18.Birks, H. H. Aquatic macrophyte vegetation development in Kråkenes Lake, western Norway, during the late-glacial and early-Holocene. J. Paleolimnol. 23, 7–19 (2000).ADS 
Article 

Google Scholar 
19.Guthrie, R. D. Origin and causes of the mammoth steppe: a story of cloud cover, woolly mammal tooth pits, buckles, and inside-out Beringia. Quat. Sci. Rev. 20, 549–574 (2001).ADS 
Article 

Google Scholar 
20.Mann, D. H., Peteet, D. M., Reanier, R. E. & Kunz, M. L. Responses of an Arctic landscape to Lateglacial and early Holocene climatic changes: the importance of moisture. Quat. Sci. Rev. 21, 997–1021 (2002).ADS 
Article 

Google Scholar 
21.Ritchie, M. in Competition and Coexistence (eds Sommer, U. & Worm, B.) 109–131 (Springer, 2002).22.Signor, P. W., Lipps, J. H., Silver, L. & Schultz, P. in Geological Implications of Impacts of Large Asteroids and Comets on the Earth vol. 190 (eds Silver, L. T. & Schultz, P. H.) 291–296 (1982).23.Haile, J. et al. Ancient DNA reveals late survival of mammoth and horse in interior Alaska. Proc. Natl Acad. Sci. USA 106, 22352–22357 (2009).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Librado, P. et al. Tracking the origins of Yakutian horses and the genetic basis for their fast adaptation to subarctic environments. Proc. Natl Acad. Sci. USA 112, E6889–E6897 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Nikolskiy, P. A., Sulerzhitsky, L. D. & Pitulko, V. V. Last straw versus Blitzkrieg overkill: climate-driven changes in the Arctic Siberian mammoth population and the Late Pleistocene extinction problem. Quat. Sci. Rev. 30, 2309–2328 (2011).ADS 
Article 

Google Scholar 
26.Pavlov, P., Svendsen, J. I. & Indrelid, S. Human presence in the European Arctic nearly 40,000 years ago. Nature 413, 64–67 (2001).ADS 
CAS 
PubMed 
Article 

Google Scholar 
27.Kuzmin, Y. V. & Keates, S. G. Siberia and neighboring regions in the Last Glacial Maximum: did people occupy northern Eurasia at that time? Archaeol. Anthropol. Sci. 10, 111–124 (2016).Article 

Google Scholar 
28.Stuart, A. J. & Lister, A. M. Extinction chronology of the woolly rhinoceros Coelodonta antiquitatis in the context of late Quaternary megafaunal extinctions in northern Eurasia. Quat. Sci. Rev. 51, 1–17 (2012).ADS 
Article 

Google Scholar 
29.Chang, D. et al. The evolutionary and phylogeographic history of woolly mammoths: a comprehensive mitogenomic analysis. Sci. Rep. 7, 44585 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
30.Vartanyan, S. L., Arslanov, K. A., Karhu, J. A., Possnert, G. & Sulerzhitsky, L. D. Collection of radiocarbon dates on the mammoths (Mammuthus primigenius) and other genera of Wrangel Island, northeast Siberia, Russia. Quat. Res. 70, 51–59 (2017).Article 
CAS 

Google Scholar 
31.Rogers, R. L. & Slatkin, M. Excess of genomic defects in a woolly mammoth on Wrangel island. PLoS Genet. 13, e1006601 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
32.Zimov, S. A., Zimov, N. S., Tikhonov, A. N. & Chapin, F. S. Mammoth steppe: a high-productivity phenomenon. Quat. Sci. Rev. 57, 26–45 (2012).ADS 
Article 

Google Scholar 
33.Yurtsev, B. A. The Pleistocene “Tundra-Steppe” and the productivity paradox: the landscape approach. Quat. Sci. Rev. 20, 165–174 (2001).ADS 
Article 

Google Scholar 
34.Rybczynski, N. et al. Mid-Pliocene warm-period deposits in the High Arctic yield insight into camel evolution. Nat. Commun. 4, 1550 (2013).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
35.Reimer, P. J. et al. The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0–55 cal kBP). Radiocarbon 62, 725–757 (2020).CAS 
Article 

Google Scholar 
36.Pedersen, M. W. et al. Postglacial viability and colonization in North America’s ice-free corridor. Nature 537, 45–49 (2016).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
37.Slon, V. et al. Neandertal and Denisovan DNA from Pleistocene sediments. Science 356, 605–608 (2017).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Lorenz, M. G. & Wackernagel, W. Adsorption of DNA to sand and variable degradation rates of adsorbed DNA. Appl. Environ. Microb. 53, 2948–2952 (1987).ADS 
CAS 
Article 

Google Scholar 
39.Meyer, M. & Kircher, M. Illumina sequencing library preparation for highly multiplexed target capture and sequencing. Cold Spring Harb. Protoc. 2010, pdb.prot5448 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
40.Willerslev, E., Hansen, A. J. & Poinar, H. N. Isolation of nucleic acids and cultures from fossil ice and permafrost. Trends Ecol. Evol. 19, 141–147 (2004).PubMed 
Article 

Google Scholar 
41.Alsos, I. G. et al. The treasure vault can be opened: large-scale genome skimming works well using herbarium and silica gel dried material. Plants 9, https://doi.org/10.3390/plants9040432 (2020).42.Hill, M. O. Diversity and evenness: a unifying notation and its consequences. Ecology 54, 427–432 (1973).Article 

Google Scholar 
43.Koleff, P., Gaston, K. J. & Lennon, J. J. Measuring beta diversity for presence-absence data. J. Anim. Ecol. 72, 367–382 (2003).Article 

Google Scholar 
44.Dixon, P. VEGAN, a package of R functions for community ecology. J. Veg. Sci. 14, 927–930 (2003).Article 

Google Scholar 
45.Grootes, P. M. & Stuiver, M. Oxygen 18/16 variability in Greenland snow and ice with 10−3- to 105-year time resolution. J. Geophys. Res. Oceans 102, 26455–26470 (1997).ADS 
CAS 
Article 

Google Scholar 
46.Andersen, K. K. et al. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. Nature 431, 147–151 (2004).ADS 
CAS 
PubMed 
Article 

Google Scholar 
47.Stuiver, M. & Grootes, P. M. GISP2 oxygen isotope ratios. Quat. Res. 53, 277–284 (2017).Article 
CAS 

Google Scholar 
48.Johnsen, S. J. et al. The δ18O record along the Greenland Ice Core Project deep ice core and the problem of possible Eemian climatic instability. J. Geophys. Res. Oceans 102, 26397–26410 (1997).ADS 
CAS 
Article 

Google Scholar 
49.Fuhrer, K., Neftel, A., Anklin, M. & Maggi, V. Continuous measurements of hydrogen peroxide, formaldehyde, calcium and ammonium concentrations along the new grip ice core from summit, Central Greenland. Atmos. Environ. A 27, 1873–1880 (1993).ADS 
Article 

Google Scholar 
50.Mayewski, P. A. et al. Major features and forcing of high-latitude northern hemisphere atmospheric circulation using a 110,000-year-long glaciochemical series. J. Geophys. Res. Oceans 102, 26345–26366 (1997).ADS 
CAS 
Article 

Google Scholar 
51.Alley, R. B. et al. Abrupt increase in Greenland snow accumulation at the end of the Younger Dryas event. Nature 362, 527–529 (1993).ADS 
Article 

Google Scholar 
52.Holden, P. B. et al. PALEO-PGEM v1.0: a statistical emulator of Pliocene–Pleistocene climate. Geosci. Model Dev. 12, 5137–5155 (2019).ADS 
CAS 
Article 

Google Scholar 
53.Karger, D. N. et al. Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4, 170122 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
54.Martindale, A. et al. Canadian Archaeological Radiocarbon Database (CARD 2.1) (Laboratory of Archaeology at the University of British Columbia, and the Canadian Museum of History, accessed 6 February 2020).55.Vermeersch, P. M. Radiocarbon Palaeolithic Europe database: a regularly updated dataset of the radiometric data regarding the Palaeolithic of Europe, Siberia included. Data Brief 31, 105793 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Rue, H., Martino, S. & Chopin, N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J. R. Stat. Soc. B 71, 319–392 (2009).MathSciNet 
MATH 
Article 

Google Scholar 
57.Lindgren, F. & Rue, H. Bayesian spatial modelling with R-INLA. J. Stat. Softw. 63, 1–25 (2015).Article 

Google Scholar 
58.Martiniano, R., De Sanctis, B., Hallast, P. & Durbin, R. Placing ancient DNA sequences into reference phylogenies. Preprint at https://doi.org/10.1101/2020.12.19.423614 (2020).59.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).60.Gu, Z., Eils, R. & Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32, 2847–2849 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Wang, Y. et al. Supporting Data for: Late Quaternary Dynamics of Arctic Biota from Ancient Environmental Metagenomics https://dataverse.no/privateurl.xhtml?token=86979109-5605-43b5-b3fb-f470d85b114c (2021).62.Theodoridis, S. et al. Climate and genetic diversity change in mammals during the Late Quaternary. Preprint at https://doi.org/10.1101/2021.03.05.433883 (2021). More