Yannick Geiger

Thierry Achard

Aline Maisse-François

Stéphane Bellemin- Laponnaz

Hyperpositive nonlinear effects in asymmetric catalysis

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Asymmetric amplification, which coins the obtention of highly enantioenriched chiral products from material of low enantiomeric purity, is a phenomenon which has fascinated scientists for several decades as it probably plays an important role in the emergence of biological homochirality. [START_REF] Carroll | A new definition of life[END_REF][START_REF] Frank | On spontaneous asymmetric synthesis[END_REF][START_REF] Noorduin | Complete chiral symmetry breaking of an amino acid derivative directed by circularly polarized light[END_REF][START_REF] Viedma | Chiral Symmetry Breaking During Crystallization: Complete Chiral Purity Induced by Nonlinear Autocatalysis and Recycling[END_REF][START_REF] Soai | Asymmetric autocatalysis and amplification of enantiomeric excess of a chiral molecule[END_REF][START_REF] Soai | Asymmetric autocatalysis of pyrimidyl alkanol and related compounds. Self-replication, amplification of chirality and implication for the origin of biological enantioenriched chirality[END_REF][START_REF] Donna | The origin of biological homochirality[END_REF][START_REF] Brewer | Chiral encoding may provide a simple solution to the origin of life[END_REF][START_REF] Bada | Origins of homochirality[END_REF][START_REF] Siegel | Chemistry: Shattered mirrors[END_REF] In enantioselective catalysis, such asymmetric amplification is known as positive non-linear effect [(+)-NLE] 11-14 since Kagan's pioneering work on non-linear relationships between product enantiomeric excess (e.e.) and catalyst e.e., which followed up on Wynberg and Feringa's study on e.e.-dependent chemical reactivity. [START_REF] Wynberg | Enantiomeric recognition and interactions[END_REF] A striking example is Noyori's DAIB-catalysed addition of dialkylzincs to benzaldehyde, where an only 15% e.e. ligand lead to a 95% e.e. product. [START_REF] Kitamura | Enantioselective addition of dialkylzincs to aldehydes promoted by chiral amino alcohols. Mechanism and nonlinear effect[END_REF] In the case of the intriguing Soai reaction, it combines a positive non-linear effect with a reaction that proceeds autocatalytically, i.e. faster as more product -which is also the catalyst -is formed ('far-from-equilibrium' process). [START_REF] Soai | Asymmetric autocatalysis and amplification of enantiomeric excess of a chiral molecule[END_REF][START_REF] Soai | Asymmetric autocatalysis of pyrimidyl alkanol and related compounds. Self-replication, amplification of chirality and implication for the origin of biological enantioenriched chirality[END_REF] Kagan and co-workers were also the first to propose mathematical models which rationalise and quantify the appearance of NLEs that are the consequence of the presence of self-organised dimeric, trimeric and even tetrameric catalysts (i.e. ML2, ML3 and ML4 models). [START_REF] Guillaneux | Nonlinear effects in asymmetric catalysis[END_REF] On the other hand, Noyori's group developed a model adapted to the DAIB ligand and the fluctional nature of dialkylzinc chemistry, in which a monomeric DAIB-alkylzinc adduct catalyses the reaction, while homochiral and heterochiral dimeric aggregates act as resting species. [START_REF] Noyori | Self and nonself recognition of chiral catalysts: The origin of nonlinear effects in the amino-alcohol catalyzed asymmetric addition of diorganozincs to aldehydes[END_REF][START_REF] Kitamura | Quantitative Analysis of the Chiral Amplification in the Amino Alcohol-Promoted Asymmetric Alkylation of Aldehydes with Dialkylzincs[END_REF][START_REF] Kitamura | Asymmetric addition of dialkylzincs to benzaldehyde derivatives catalyzed by chiral β-amino alcohols. Evidence for the monomeric alkylzinc aminoalkoxide as catalyst[END_REF] In both models, (+)-NLEs are achieved when the aggregates consist preferably of both catalyst enantiomers in equal parts (heterochiral aggregate) and when those aggregates are less (or not at all) catalytically active. Under those conditions, the active scalemic (i.e. nonenantiopure) catalyst may become as selective as the enantiopure one.

Interestingly, Kagan also showed that it should be also possible to find chiral systems where the product e.e.s "take values much larger for partially resolved ligand than enantiomerically pure ligand!". [START_REF] Guillaneux | Nonlinear effects in asymmetric catalysis[END_REF] They predicted this behaviour as a particular case of the ML3 system, where both homochiral (RRR) and heterochiral catalyst trimers (RRS) catalyse, the latter being the most selective one, and coined the term of hyperpositive NLEs. [START_REF] Guillaneux | Nonlinear effects in asymmetric catalysis[END_REF] However, to the best of our knowledge, a concrete example of a hyperpositive NLE has never been observed.

Here we present a case of a hyperpositive NLE using an ephedrine-based ligand in the catalytic addition dialkylzincs to aldehydes. Correlation between the product e.e.

(eeP) and ligand e.e. (eeL) for the asymmetric addition of dimethylzinc (or diethylzinc) to benzaldehyde using (1R,2S)-N-benzylephedrine [(-)-NBE] as chiral auxiliary ligand showed that an non-enantiopure ligand increased the product e.e. compared to the enantiopure ligand. A mechanistic study was carried out in which the influence of catalyst concentration and temperature were investigated and the reaction's kinetics analysed. Rationalisation points towards this system being different from the ML3system theorized by Kagan, consisting of competition between both discrete and dimeric homochiral catalytic species, and therefore also differing from Noyori's model for aminoalcohol-catalysed dialkylzinc additions.

Results

Hyperpositive non-linear effect Fig. 1 displays the hyperpositive correlation between the product e.e. (eeP) and ligand e.e. (eeL) for the asymmetric addition of dimethylzinc to benzaldehyde using (1R,2S)-N-benzylephedrine [(-)-NBE] [START_REF] Coote | Enantiospecific Stereodivergent Synthesis of trans-and cis-N(2),3-Dimethyl-4phenyl-1,2,3,4-tetrahydroisoquinolines[END_REF] as chiral auxiliary ligand. While the enantiopure ligand induces a very poor enantioselectivity of only 16%, the progressive reduction of the enantiomeric excess of the ligand leads to an approximately exponential growth of the product e.e., reaching a maximum of 53% e.e. with only 5% eeL. Thus, using an enantio-impure ligand of only 5% e.e. increased the product e.e. by 37% compared to the enantiopure ligand, which is astonishing -when seeking a more stereoselective catalyst, typical logic would suggest searching for alternative molecular structures of the ligand, but never to reduce the enantiopurity of the ligand. The use of ZnEt2 as the dialkylzinc reagent leads to a less pronounced but still significant eeP increase when decreasing eeL (cf. Supplementary Fig. 1). Kagan's ML3 model implies that the NLE curve can be hyperpositive only in the upper half of the ligand e.e. scale, which does not match with our observations shown in Fig. 1. [START_REF] Guillaneux | Nonlinear effects in asymmetric catalysis[END_REF] Therefore, the ML3 model is not applicable here and we have to search for other explanations for this peculiar phenomenon.

Mechanistic investigations

One striking difference between the use of enantiopure and scalemic ligand is that upon addition of the dialkyl zinc reagents, a white precipitate was observed in the latter case whereas enantiopure solutions stay clear.

Filtration and hydrolysis of the precipitates obtained from 50% e.e. ligand samples yielded the starting ligand in its racemic form, which is most likely due to the formation of a heterochiral dimer. Thus, it becomes clear that in this case precipitation of an inactive heterochiral aggregate drives the overall system to a positive (+)-NLE, while leaving an enantioenriched Zn complex in solution. Decreasing the e.e. of the ligand would thus generate more heterochiral -and inactive -dimer and therefore less catalytically active species (reservoir effect). [START_REF] Guillaneux | Nonlinear effects in asymmetric catalysis[END_REF] In other words, reducing the eeL induces a decrease of the catalyst concentration. We hypothesized that the observed hyperpositive NLE could be related to the active catalyst concentration. To this end, we evaluated the product e.e. while varying the concentration of the enantiopure ligand (Fig. 2a). Indeed, a decrease of the catalyst loading led to an increase of eeP. Interestingly, the shape of the curve is consistent with the one seen in Fig. 1, which made us compare the two datasets: assuming that the heterochiral aggregate fully precipitates and does not interfere in the catalytic process, a system with e. g. originally 20 mol% catalyst loading of a 50 % eeL-ligand should be identical to a 10 mol% catalyst loading of enantiopure catalyst (cf. Supplementary Fig. 2 for a more detailed graphical representation). Thus, we could correlate the results from the enantiopure catalysts with a simulated e.e. if a scalemic ligand was used. The resulting graph from Fig. 2a was then superimposed with the hyperpositive NLE curve (Fig. 2b), leading to an excellent match between both curves. The catalyst loading screening models accurately the eeP increase at lower eeL and proves that the observed hyperpositive NLE in Fig. 1 is due to a change in the concentration of the soluble part of the catalyst. The same catalyst concentration effect was observed when using ZnEt2 as reagent and a perfect correlation was obtained with the measured non-linear effect; dilution of the catalytic system while keeping the catalyst loading at 20 mol% also led to an eeP increase (cf.

Supplementary Fig. 3). The results point toward a possible participation of several active species with different enantioselectivities. In order to gain more insight into the system, we performed a temperature screening of the reaction, investigated the kinetic order in catalyst and studied the eeP evolution over time. Increasing the reaction temperature at constant catalyst concentration lead to a surprising eeP increase, both with the enantiopure and the 20% e.e. enantioenriched ligand (Fig. 3), which is an unusual inverse temperature dependence. This is obviously accompanied by a shift to the more enantioselective catalytic species, which overcompensates the usual eeP loss at high temperature. This temperature effect is less pronounced with ZnEt2 and correlates with its weaker hyperpositive NLE (cf. Supplementary Fig. 4). The kinetic order in catalyst c was determined by Visual Time-Normalised Analysis (VTNA) of the kinetic profiles obtained from reactions with different catalyst loadings, monitored by FTIR (the dataset can be found in Supplementary Fig. 5). [START_REF] Burés | A Simple Graphical Method to Determine the Order in Catalyst[END_REF][START_REF] Burés | Variable Time Normalization Analysis: General Graphical Elucidation of Reaction Orders from Concentration Profiles[END_REF][START_REF] Nielsen | Visual kinetic analysis[END_REF] Simultaneous analysis of all kinetic profiles did not give a single value of c, which indicates a non-constant c. Therefore, a segmented analysis was conducted and the results are shown in Fig. 4. The plot c vs. catalyst loading gives an unusual, Ushaped curve which starts at c = 1 at low catalyst loading, which goes down to 0.4 (at 10-15 mol% catalyst) but then increases again at higher catalyst loading. Exponential or near-exponential curves have been observed in case of catalyst aggregation, [%]

Temperature [°C] where either only the monomer or only the aggregate catalyses, but these do not contain minimum or maximum values. 21 In addition, the evolution of eeP over the course of the reaction was investigated. Fig. 5 displays eeP vs. time plots at an early stage of the reaction and with different catalyst loadings. For a given catalyst loading, the eeP values are not constant but increase as the reaction proceeds, a result that is observed in all cases. These observations indicate, in addition to the fact that the eeP is directly dependent on the catalyst loading, that the ratio of the different catalytic active species changes over time. Changing the catalyst loading seemingly does not affect significantly the way in 0,0

Catalyst order c

Catalyst loading [mol%] which this ratio evolves (at least in the early stage of the reaction) but it affects the initial distribution of the different catalysts. (standard reaction conditions as in Fig. 1). Some plots have been omitted for clarity; the full dataset can be found in Supplementary Fig. 6. Each point is the mean of three different experiments; the vertical bars depict standard deviations.

Discussion

We propose the reaction scheme shown in Fig. 6 that accounts for our observations.

Our results indicate the presence of two different enantiodivergent pathways: a monomeric (-)-NBE-ZnMe-catalysed mechanism and a possible dimeric (-)-NBE-ZnMe-catalysed mechanism that operate at a steady state. Assuming that the eeP produced by the monomer catalyst (eem) is higher than with the dimer (eed) and that both catalysts are in a dynamic equilibrium, then decreasing the catalyst concentration leads to an eeP increase as the equilibrium is shifted to the more enantioselective monomer. Decreasing eeL has the same effect as it causes the lowering of the effective catalyst in solution (because of the precipitating heterochiral aggregate) and therefore gives rise to the observed hyperpositive NLE. This is evidenced by the following points: 1) Chiral alkylzinc aminoalkoxides, issued from dialkylzinc and a chiral aminoalcohol, are known to generate homochiral as well as heterochiral dimers. [START_REF] Kitamura | Enantioselective addition of dialkylzincs to aldehydes promoted by chiral amino alcohols. Mechanism and nonlinear effect[END_REF] Since the latter are insoluble in our case, only monomeric and homochiral dimeric complexes (both consisting only of the major enantiomer) stay in solution. This was supported by [START_REF] Carroll | A new definition of life[END_REF] H DOSY NMR experiments [START_REF] Neufeld | Accurate molecular weight determination of small molecules via DOSY-NMR by using external calibration curves with normalized diffusion coefficients[END_REF][START_REF] Bachmann | Solution structures of alkali metal cyclopentadienides in THF estimated by ECC-DOSY NMR-spectroscopy (incl. software)[END_REF][START_REF] Bachmann | New External Calibration Curves (ECCs) for the Estimation of Molecular Weights in Various Common NMR Solvents[END_REF][START_REF] Kreyenschmidt | Molecular Weight Estimation of Molecules Incorporating Heavier Elements from van-der-Waals Corrected ECC-DOSY[END_REF] of enantiopure (-)-NBE reacted with ZnMe2 in toluene-d8, which indicated that the complexes have a higher molecular weight (MW) than expected from monomeric (-)-NBE-ZnMe over a broad temperature range (cf. Supplementary Fig. 7).

2) The monomer/homochiral dimer equilibrium is necessarily affected by the overall catalyst concentration as well as by the reaction temperature. The higher the catalyst concentration and the lower the temperature, the more it aggregates. This is seen in the upwards slope in Fig. 1a and 3, considering the homochiral dimer less enantioselective than the monomeric catalyst.

3) The impact of the catalyst loading on the product e.e. is already visible at the initial stage of the reaction, where other potential aggregates (e. g. catalyst/product-aggregates) play a minor role and therefore are not likely to be the origin of the eeP variation. 4) The intriguing catalyst order / catalyst loading dependence is also consistent with a model in which both monomer and homochiral dimer catalyse. If we consider both monomeric and dimeric catalysts to have an inherent first order kinetic behaviourwhich is the case for monomeric DAIB-ZnMe [START_REF] Kitamura | Enantioselective addition of dialkylzincs to aldehydes promoted by chiral amino alcohols. Mechanism and nonlinear effect[END_REF][START_REF] Kitamura | Quantitative Analysis of the Chiral Amplification in the Amino Alcohol-Promoted Asymmetric Alkylation of Aldehydes with Dialkylzincs[END_REF] then we would expect first order behaviour at very low and at very high catalyst loading, when either only the monomeric or only the dimeric catalyst is present. In the intermediate regime, the observed order is expected to change because (i) the overall kinetic constant of the system varies as a function of the dimeric catalyst's kinetic constant (k2) and the monomer's kinetic constant (k1) ratio, and (ii) because of the aggregation-induced decreasing amount of catalyst molecules. Both effects lead to a non-linear variation of the reaction rate and, therefore, to c ≠ 1.

In the light of Noyori's model for aminoalcohol-catalysed addition of dialkylzinc to aldehyde, these findings are notable because catalysis by aggregates was always excluded (except in certain conditions where the model didn't fit experimental results)

and this supposition has only been rarely questioned. [START_REF] Chaloner | Enantioselective reaction of diethylzinc with arenecarbaldehydes in the presence of ephedrine derivatives[END_REF][START_REF] Zhang | Dramatic temperature effect in asymmetric catalysis in the enantioselective addition of diethylzinc to aldehydes[END_REF][START_REF] Kang | The interaction of chiral amino thiols with organozinc reagents and aldehydes: a mechanism of amino thiol-catalyzed addition of organozinc reagents to aldehydes[END_REF] In conclusion, we have discovered an example of a hyperpositive NLE in the (-)-NBEcatalysed enantioselective addition of dialkylzincs to benzaldehyde. Catalyst loading screenings, temperature screenings, conversion vs time and eeP vs time plots indicate that the catalytic system consists of two enantiodivergent pathways, in which both monomeric and homochiral dimeric (-)-NBE-ZnMe complexes catalyse the reaction, with different enantioselectivities. Further work is in progress to gain a deeper understanding of the mechanism of the catalyst system. The dimeric catalyst RR is depicted as being in equilibrium with higher-order aggregates like trimers or tetramers, which have been observed with ligands derived from ephedrine and other aminoalcohols [START_REF] Coates | Reactions between some organozinc compounds and 2-dimethylaminoethanol, acetoxime, phenylisocyanate, and benzophenone. Some observations on the methylzinc methoxide tetramer[END_REF][START_REF] Boersma | Investigations on organozinc compounds X. Coordination chemistry of organozinc compounds RZnX: organozinc-oxygen derivatives of potentially bidentate ligands[END_REF][START_REF] Van Der Schaaf | Organozinc complexes with monoanionic chelating phenolates or 2-pyridylmethanolates. Molecular structure of [Zn(CH2SiMe3){OCH2(2-Py)}]4[END_REF][START_REF] Fitzpatrick | Thiol and disulfide derivatives of ephedra alkaloids 2: A mechanistic study of their effect on the addition of diethyl zinc to benzaldehyde[END_REF] and whose presence (and possible catalytic activity) we cannot fully exclude so far.

Fig. 1 |

 1 Fig. 1 | Hyperpositive non-linear effect. (1R,2S)-N-benzylephedrine [(-)-NBE]-catalysed

Fig. 2 |

 2 Fig. 2 | Impact of the catalyst loading in the enantiomeric excess of the product. a,

Fig. 3 |

 3 Fig. 3 | Temperature-dependence study. Addition of ZnMe2 to benzaldehyde in presence of (-)-NBE ligand as function of the temperature. : enantiopure (-)-NBE, : (-)-NBE with

Fig. 4 |

 4 Fig. 4 | Catalyst order as function of the catalyst loading. (-)-NBE-catalysed addition of

Fig. 5 |

 5 Fig. 5 | Evolution of the product e.e. as function of time and at different catalyst

Fig. 6 |

 6 Fig. 6 | Proposed enantiodivergent model for the (-)-NBE-catalysed addition of

Acknowledgments

The authors gratefully acknowledge the CNRS and the Ministère de l'Enseignement Supérieur et de la Recherche (MESR) for a PhD grant to Y.G. This work is a part of the Strasbourg IDEX program and we thank the NIE Labex for support (ANR-11-LABX-0058_NIE). We thank Prof Lutz H. Gade (University of Heidelberg) and Prof Joseph Moran (University of Strasbourg) for helpful comments on this work.

Data availability

Complete experimental procedures, compound characterisation data and kinetic rate profiles are available in the Supplementary Information.

Competing interests

The authors declare no competing interests.


Jonathan Lenoir
email: jonathan.lenoir@u-picardie.fr

Romain Bertrand

Lise Comte

Luana Bourgeaud

Tarek Hattab

Jérôme Murienne

Gaël Grenouillet

Species better track climate warming in the oceans than on land

of the paper: "Compiling a global geo-database of >30,000 range shifts, the authors show that marine species closely track shifting isotherms whereas terrestrial species lag behind, likely due to wider thermal safety margins and movement constraints imposed by human activities."

There is mounting evidence of species redistribution as climate warms. Yet, our knowledge of the coupling between species range shifts and isotherm shifts remains limited. Here we introduce BioShifts, a global geo-database of 30,534 range shifts. Despite a spatial imbalance towards the most developed regions of the northern hemisphere and a taxonomic bias towards the most charismatic animals and plants of the planet, data shows that marine species are better at tracking isotherm shifts, and move towards the pole six times faster than terrestrial species. More specifically, we find that marine species closely track shifting isotherms in warm and relatively undisturbed waters (e.g. Central Pacific Basin) or in cold waters subject to high human pressures (e.g. North Sea). On land, human activities impede the capacity of terrestrial species to track isotherm shifts in latitude, with some species shifting in the opposite direction to isotherms. Along elevational gradients, species follow the direction of isotherm shifts but at a pace that is much slower than expected, especially in areas with warm climates. Our results suggest that terrestrial species are lagging behind shifting isotherms more than marine species, which is likely related to the interplay between wider thermal safety margin of terrestrial vs. marine species and more constrained physical environment for dispersal in terrestrial vs. marine habitats.

The redistribution of life on Earth in response to climate change [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Chen | Rapid range shifts of species associated with high levels of climate warming[END_REF][START_REF] Poloczanska | Global imprint of climate change on marine life[END_REF][START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF] is now considered a global change driver on its own with far-reaching implications for ecosystem and human health [START_REF] Pecl | Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being[END_REF] .

Managing the consequences of climate-driven species redistributions requires a better understanding of the capacity of species of various taxonomic groups and from different habitats to shift their distribution ranges and track shifting isotherms. As climate warms, isotherms are shifting poleward and upslope to cooler latitudes and elevations in most parts of the world, generating spatially-structured patterns in the velocity of isotherm shifts [START_REF] Loarie | The velocity of climate change[END_REF][START_REF] Burrows | The pace of shifting climate in marine and terrestrial ecosystems[END_REF] . Marine organisms seem to closely track this complex mosaic of climate velocities [START_REF] Pinsky | Marine taxa track local climate velocities[END_REF] . However, the pattern is less clear for terrestrial organisms [START_REF] Chen | Rapid range shifts of species associated with high levels of climate warming[END_REF] . Evidence suggests that biotic responses on land are lagging behind the velocity of climate change, particularly for long-lived species and poor-dispersers [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF][START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF] .

To date, a comprehensive analysis of the coupling between the velocity of species range shifts and the velocity of isotherm shifts across biological systems (i.e. terrestrial vs. marine) and life forms (e.g. ectotherms vs. endotherms) is still lacking [START_REF] Pinsky | Greater vulnerability to warming of marine versus terrestrial ectotherms[END_REF][START_REF] Pinsky | Climate-driven shifts in marine species ranges: scaling from organisms to communities[END_REF] .

To fill this knowledge gap, we compiled data on range shifts for marine and terrestrial species in both the northern and southern hemispheres from an exhaustive literature review building on and updating the most recent syntheses on climate-related range shifts [START_REF] Chen | Rapid range shifts of species associated with high levels of climate warming[END_REF][START_REF] Poloczanska | Global imprint of climate change on marine life[END_REF][START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF] (BioShifts [START_REF] Comte | BioShifts: a global geo-database of climate-induced species redistribution over land and sea[END_REF] ; see Data Availability for download) (Supplementary Fig. 1). The BioShifts geodatabase encompasses 30,534 range shifts documented from 258 peer-reviewed studies [START_REF] Comte | BioShifts: a global geo-database of climate-induced species redistribution over land and sea[END_REF] , spanning a total of: 12,415 harmonized species names; four kingdoms (Bacteria, Plantae, Fungi and Animalia); 20 phyla; and 56 classes. It also contains several methodological attributes (e.g. study area, study period, sampling effort, data quality), that can be used to account for methodological variations in meta-analyses or quantitative reviews such as ours [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] . Based on this geo-database, we first carefully assessed the geographical and taxonomic biases [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF][START_REF] Feeley | Most 'global' reviews of species' responses to climate change are not truly global[END_REF] impeding our knowledge of climate change effects on species redistribution (Extended Data 1-3). We then provided robust estimates of the velocity of latitudinal and elevational range shifts for the 20 most studied taxonomic classes (Fig. 1), with the aim to compare our estimates with former estimates from the existing literature [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Chen | Rapid range shifts of species associated with high levels of climate warming[END_REF][START_REF] Poloczanska | Global imprint of climate change on marine life[END_REF] while accounting for potential methodological biases [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] .

To do so, we fitted several linear mixed-effects models (LMMs) with methodological attributes treated as random effects [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] , all arranged in a full factorial design of geographical gradient (latitude vs. elevation) [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Chen | Rapid range shifts of species associated with high levels of climate warming[END_REF] × biological system (marine vs. terrestrial) [START_REF] Feeley | Most 'global' reviews of species' responses to climate change are not truly global[END_REF] × hemisphere (North vs. South) [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF] × positional parameter (centroid vs. margins) [START_REF] Poloczanska | Global imprint of climate change on marine life[END_REF][START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF] (Supplementary Table 1). We expected (1) faster velocities of latitudinal range shifts in the oceans than on land 3 , (2) faster velocities of range shifts in the northern hemisphere than in the southern hemisphere due to the interhemispheric asymmetry in the rate of climate warming over the twentieth century 16 and (3) different velocities of shifts across the distribution range. For each taxonomic class, we thus estimated the velocity of range shift separately for the centroid and the margins of the range, making the distinction between the trailing and leading edge to test for different types of range shifts [START_REF] Poloczanska | Global imprint of climate change on marine life[END_REF][START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF][START_REF] Rumpf | Elevational rear edges shifted at least as much as leading edges over the last century[END_REF][START_REF] Freeman | Expanding, shifting and shrinking: The impact of global warming on species' elevational distributions[END_REF] : trailing-edge contractions; leading-edge expansions; lockstep shifts involving both trailing-edge contractions and leading-edge expansions together with a displacement at the centroid of the range; and lean range shifts involving stable margins with the centroid shifting within the existing range.

Then, we assessed the coupling between the velocity of isotherm shifts and the velocity of range shifts at the species level, along the latitudinal gradient in marine systems and separately for the latitudinal and elevational gradients in terrestrial systems (Extended Data 4-6). As before, we controlled for varying methodologies [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] , using LMMs, and tested for two-way interaction terms between the velocity of isotherm shifts and: (i) life-form categories (ectotherms, endotherms, phanerogams, cryptogams); (ii) baseline temperatures (the historical mean annual temperature regime prior to the baseline survey); and (iii) the standardized human footprint index (a scaled variable summarizing the direct and indirect human pressures on both terrestrial and marine environments: human population density; buildings; roads; agricultural lands; pollution; commercial shipping; industrial fishing; ocean acidification; etc.) [START_REF] Halpern | A Global Map of Human Impact on Marine Ecosystems[END_REF][START_REF] Venter | Global terrestrial human footprint maps for 1993 and 2009[END_REF] . Among animals, we distinguished ectotherms from endotherms to test the hypothesis that ectotherms better track shifting isotherms than endotherms. Theory and evidence suggest that ectotherms are more sensitive to temperature fluctuations than endotherms [START_REF] Paaijmans | Temperature variation makes ectotherms more sensitive to climate change[END_REF] and conform more closely to their physiological limits of thermal tolerance, especially so for marine ectotherms [START_REF] Sunday | Thermal tolerance and the global redistribution of animals[END_REF] , and are thus more likely to closely track shifting isotherms [START_REF] Pinsky | Marine taxa track local climate velocities[END_REF] . For chlorophyllous organisms, we distinguished phanerogams or seed-bearing plants (angiosperms and gymnosperms) from cryptogams reproducing by spores (ferns, mosses, lichens and algae) to test the hypothesis that greater dispersal abilities in cryptogams allows to better track shifting isotherms [START_REF] Angert | Do species' traits predict recent shifts at expanding range edges?[END_REF] . As historical temperature regime may affect the rate at which species are shifting their distributions [START_REF] Guo | Land-use change interacts with climate to determine elevational species redistribution[END_REF][START_REF] Bertrand | Ecological constraints increase the climatic debt in forests[END_REF] , we tested whether the coupling between species range shifts and isotherm shifts varied with baseline temperatures. Similarly, we accounted for potential effects of anthropogenic disturbances on the magnitude of the coupling. We expected that land-use intensity in terrestrial systems may impede species range shifts through habitat fragmentation [START_REF] Warren | Rapid responses of British butterflies to opposing forces of climate and habitat change[END_REF] , while exploitation of marine resources (e.g. demersal or bottom trawling) may accelerate the relocation of exploited stocks through massive population die-offs and crashes followed by local extinctions [START_REF] Engelhard | Climate change and fishing: a century of shifting distribution in North Sea cod[END_REF] .

Results and discussion

Geographical, taxonomic and methodological biases matter. We found a strong spatial imbalance in the data towards the most developed regions of the northern hemisphere (Extended Data 1-2) and a clear taxonomic bias towards the most charismatic animals (Aves, Actinopterygii, Amphibia, Mammalia) and plants (Magnoliopsida, Liliopsida) (Extended Data 3). This supports former claims that global meta-analyses on species range shifts are not truly global [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF][START_REF] Feeley | Most 'global' reviews of species' responses to climate change are not truly global[END_REF] and that most species remain understudied, while others attract most of the public, scientific and government attention [START_REF] Troudet | Taxonomic bias in biodiversity data and societal preferences[END_REF] . In addition to these geographic and taxonomic biases, differences in methodological attributes among studies play a key role in the observed variation in the velocity of range shifts among the 12,415 species included in BioShifts [START_REF] Comte | BioShifts: a global geo-database of climate-induced species redistribution over land and sea[END_REF] (Fig. 2). Most of this variation was explained by methodological attributes, which contributed from 6 to 82% (mean = 36%, median = 35%) of the total variation (Supplementary Table 1). By contrast, differences among taxonomic classes and positions at range margins (trailing edge vs. leading edge) contributed only 0 to 50% (mean = 10%, median = 7%) of the total variation. These findings confirm the importance of accounting for varying methodologies in meta-analyses [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] . However, contrary to former meta-analyses arguing against the use of studies reporting range shifts for a single or a handful number of species because such studies are more likely to select responsive species 1-3 , we found no relationship between sample size and the velocity of range shifts (Supplementary Fig. 2). Hence, we recommend future meta-analyses on range shifts to consider all available information instead of constraining the analyses to the set of studies that focused on multiple species (e.g. more than 3 species).

Marine organisms move faster than terrestrial organisms. Once methodological attributes were accounted for, the variation in the estimated velocity of range shifts among taxonomic classes and positions at range margins was still fairly large, ranging from 3.20 m.yr -1 downslope for freshwater fishes (trailing edge, northern hemisphere) to 12.39 m.yr -1 upslope for amphibians (trailing edge, southern hemisphere) and from 6.52 km.yr -1 equatorward for reptiles (trailing edge, northern hemisphere) to 18.54 km.yr -1 poleward for insects (centroid of the range, northern hemisphere) (Fig. 3). Marine species (~80% being ectotherms in the database, Extended Data 2) have moved towards the poles at a mean (±s.e.m.) pace of 5.92±0.94 km.yr -1 (one-sample Student's t-test: t = 6.26; df residuals = 23; P = 2.20 × 10 -6 ), almost six times faster than terrestrial species (one-way ANOVA: F = 12.68; df factor = 1; df residuals = 45; P = 8.88 × 10 - 4 ). This mean velocity far exceeds the one reported by the first synthesis (0.61±0.24 km.yr -1 ) 1 but is very similar in magnitude and direction to the mean velocity reported by a more recent synthesis focusing exclusively on marine species (7.20±1.35 km.yr -1 ) 3 . Importantly, the mean velocity of latitudinal range shifts we found for terrestrial systems (1.11±0.96 km.yr -1 ) was nonsignificantly different from zero (one-sample Student's t-test: t = 1.15; df residuals = 22; P = 0.25). This contradicts a former synthesis from 2011 reporting a mean positive velocity of latitudinal range shifts across a wide range of taxonomic groups (1.76±0.29 km.yr -1 ) 2 . Although the authors of this synthesis chiefly focused on terrestrial taxonomic groups, they also included data on range shifts from several marine taxonomic groups (molluscs and algae) in their analyses, which could explain the discrepancy. Along the elevation gradient, we found that terrestrial species have shifted upslope at a mean pace of 1.78±0.41 m.yr -1 (one-sample Student's t-test: t = 4.33; df residuals = 36; P = 1.13 × 10 -3 ), slightly faster than what was previously reported (1.22±0.18 m.yr -1 ) 2 .

Although we found a tendency towards faster latitudinal range shifts in the northern hemisphere (4.24±0.70 km.yr -1 ; one-sample Student's t-test: t = 4.66; df residuals = 36; P = 4.28 × 10 -5 ) as opposed to the southern hemisphere (1.07±0.34 km.yr -1 ; one-sample Student's t-test: t = 1.43; df residuals = 9; P = 0.19), the difference was not significant (one-way ANOVA: F = 3.08; df factor = 1; df residuals = 45; P = 0.09). More data on species range shifts are thus clearly needed in the southern hemisphere (see geographical biases in Extended Data 1-2) to be able to confirm or infirm our hypothesis of faster range shifts in the northern hemisphere related to the faster rates of climate warming over there compared with the southern hemisphere [START_REF] Friedman | Interhemispheric Temperature Asymmetry over the Twentieth Century and in Future Projections[END_REF] .

Regarding the position within the range, terrestrial taxa seem to have relatively stable latitudinal distributions, showing no clear signal of range shift at the trailing edge (-0.17±1.61 km.yr -1 ; one-sample Student's t-test: t = -0.10; df residuals = 5; P = 0.92), the centroid of the range (2.41±2.45 km.yr -1 ; one-sample Student's t-test: t = 0.98; df residuals = 7; P = 0.36) or the leading edge (0.81±0.65 km.yr -1 ; one-sample Student's t-test: t = 1.24; df residuals = 8; P = 0.25). By contrast, marine species seem to be very sensitive to warming, showing trailing-edge contractions (6.49±2.13 km.yr -1 ; one-sample Student's t-test: t = 3.04; df residuals = 7; P = 0.02), leading-edge expansions (6.02±1.77 km.yr -1 ; one-sample Student's t-test: t = 3.40; df residuals = 8; P = 9.32 × 10 -3 ) and poleward shifts at the centroid of the range (5.13±0.41 km.yr -1 ; one-sample Student's t-test: t = 12.54; df residuals = 6; P = 1.57 × 10 -5 ). Our results also indicate that the leading and trailing edge of marine species are equally sensitive to warming (one-way ANOVA: F = 0.03; df factor = 1; df residuals = 15; P = 0.87), which is consistent with expectations from thermal tolerance limits of marine ectotherms [START_REF] Sunday | Thermal tolerance and the global redistribution of animals[END_REF] (though in contrast to a previous report [START_REF] Poloczanska | Global imprint of climate change on marine life[END_REF] ). In turn, this suggests that marine species are moving in lockstep 4 towards the poles. Similarly, along elevational gradients, the trailing and leading edge of terrestrial species have moved towards the summits at a comparable mean pace of 2.34±0.67 m.yr -1 and 2.15±0.60 m.yr -1 , respectively (one-way ANOVA: F = 0.03; df factor = 1; df residuals = 23; P = 0.87). This indicates that terrestrial species are moving in lockstep towards mountain summits, which is very consistent with two recent syntheses concluding on symmetric boundary shifts in mountains [START_REF] Rumpf | Elevational rear edges shifted at least as much as leading edges over the last century[END_REF][START_REF] Freeman | Expanding, shifting and shrinking: The impact of global warming on species' elevational distributions[END_REF] .

Note, however, that the mean upslope shift was significant at the leading edge (one-sample Student's t-test: t = 6.19; df residuals = 12; P = 4.65 × 10 -5 ), but only marginally significant at both the trailing edge (one-sample Student's t-test: t = 2.07; df residuals = 11; P = 0.06) and the centroid of the range (one-sample Student's t-test: t = 2.13; df residuals = 11; P = 0.06).

Marine species are better at tracking isotherm shifts.

Assessing the degree of coupling between species range shifts and isotherm shifts (Extended Data 4-6), we found that marine species better track isotherm shifts in latitude than terrestrial species (Figs. 4-5). For marine systems, our best model explained 33% of the total variation in the velocity of species range shifts (Fig. 4a). Only 4% of the total variation was related to fixed effects, namely the velocity of isotherm shifts (VIS), standardized human footprint index (HFI), baseline temperatures (BT), life forms (LF) and synergistic effects between VIS and HFI or BT (Fig. 4b), whereas 29% was explained by random effects or methodological attributes (Supplementary Table 2). Again, this strongly supports the idea that varying methodologies in estimates of climate-driven biological responses can contribute to most of the explained variation and need to be explicitly considered in quantitative reviews [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] . Noteworthy, we found that faster climate velocities combined with higher human pressures in the oceans (e.g. commercial shipping, industrial fishing, ocean acidification) or warmer sea surface temperatures during the baseline survey increases the velocity of species range shifts along the latitudinal gradient for both marine ectotherms and cryptograms (Figs. 5cd, Extended Data 7, Supplementary Table 2). More specifically, we found that marine species closely track shifting isotherms either in initially warm and undisturbed waters (e.g. Central Pacific Basin) [START_REF] Halpern | A Global Map of Human Impact on Marine Ecosystems[END_REF] or in initially cold waters where human activities are more pronounced (e.g. Norwegian Sea, North Sea and English Channel) (Fig. 6c, Extended Data 8). This pattern is unlikely to result from a collinearity issue between the velocity of isotherm shifts and the standardized human footprint index (R 2 = 0.05) (Supplementary Fig. 3). Instead, it may stem from the combination of two processes. First, marine species are living closer to their upper thermal limits in the tropics, where sea surface temperatures are the highest, thus increasing the likelihood of local extirpations at their trailing edges as climate warms [START_REF] Pinsky | Greater vulnerability to warming of marine versus terrestrial ectotherms[END_REF] . Second, lower constraints on dispersal and colonization in the oceans (as opposed to terrestrial habitats) 3 may help species to rapidly shift their distribution towards the newly available habitats. By contrast, at high latitudes where the thermal safety margin of marine species is larger [START_REF] Pinsky | Greater vulnerability to warming of marine versus terrestrial ectotherms[END_REF] , climate warming alone is unlikely to explain isotherm tracking. Instead, anthropogenic activities (e.g. fishing pressure and pollution in the North Sea) may render populations more sensitive to climate change by reducing abundance and density, truncating the age distribution and leading to the depletion of fish stock at the trailing edge of their range [START_REF] Engelhard | Climate change and fishing: a century of shifting distribution in North Sea cod[END_REF] . In parallel, successful management actions at higher latitudes, such as along the Norwegian's coastlines in the Norwegian Sea and the Barents Sea, combined with climate warming, may increase population sizes of commercial fishes at the leading edge of their range [START_REF] Kjesbu | Synergies between climate and management for Atlantic cod fisheries at high latitudes[END_REF] , thus promoting successful colonization.

Unlike in the oceans, the degree of coupling between the velocity of species range shifts and the velocity of isotherm shifts is comparatively poor on land (Fig. 4). Again, this inconsistency with a former synthesis 2 from 2011 may stem from the fact that we here analyzed the coupling separately for the marine and terrestrial systems, a distinction that was not possible in 2011 due to a lack of data at that time, for marine systems. Our best model explained 47% of the total variation in the velocity of species range shifts along the latitudinal gradient on land (Supplementary Table 2), of which the largest proportion was explained by varying methodologies among studies (Fig. 4a). Among fixed effects, we found that range shifts were best explained by differences between life forms and a negative interaction term between the velocity of isotherm shifts and the standardized human footprint index (explaining 8% of the total variation; Fig. 4b). Such antagonistic effect between climatic and human-related drivers suggests that habitat loss and fragmentation -associated with high population densities and other human activities such as agricultural practices -in the lowlands, combined with limited species' dispersal abilities -relative to the speed at which isotherms are shifting along the latitude gradient [START_REF] Schloss | Dispersal will limit ability of mammals to track climate change in the Western Hemisphere[END_REF] likely impede the capacity of terrestrial taxa to track shifting isotherms. Again, this pattern is unlikely to result from a collinearity issue between the velocity of isotherm shifts and the standardized human footprint index (R 2 = 0.09) (Supplementary Fig. 3). Interestingly, we found that when exposed to a high degree of anthropogenic disturbances (HFI > 0.3), terrestrial species tend to shift in the opposite direction to isotherms (i.e. HFI conditions for which both velocities show opposite signs in Fig. 5b and Extended Data 9), most likely due to local extinction processes at the leading edge. For instance, a previous study showed that during 1970-1999, habitat loss and degradation led to a decline in the distribution sizes of three-quarters of butterfly species that approach their northern climatic range margins in Britain, outweighing the climate-induced species range shifts that were expected from climate warming [START_REF] Warren | Rapid responses of British butterflies to opposing forces of climate and habitat change[END_REF] . This is consistent with the general idea that land-use and climate change may act as opposing forces on species distribution changes. In addition, air conducts heat 25 times less effectively than water [START_REF] Pinsky | Climate-driven shifts in marine species ranges: scaling from organisms to communities[END_REF] , which makes terrestrial species, in general, less sensitive than marine species to temperature fluctuations and thus less likely to move as a direct response to climate warming [START_REF] Pinsky | Greater vulnerability to warming of marine versus terrestrial ectotherms[END_REF] . The availability of thermal microrefugia (e.g. shaded environments) on land may also allow species to more easily regulate their body temperature (e.g. microhabitats may allow terrestrial ectotherms to increase their thermal safety margin by 3°C on average as compared with marine ectotherms [START_REF] Pinsky | Greater vulnerability to warming of marine versus terrestrial ectotherms[END_REF] ). Hence, we confirm that isotherm tracking is very unlikely for terrestrial taxa living in the lowlands [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF][START_REF] Bertrand | Ecological constraints increase the climatic debt in forests[END_REF] .

Along elevational gradients, the best model explained 11% of the total variation in the velocity of range shifts (Fig. 4a; Supplementary Table 2) and showed that the velocity of isotherm shifts interacts with both baseline temperatures and life forms (explaining 2% of the total variation; Figs. 5a-6a, Extended Data 9-10). Contrary to the latitudinal gradient, the standardized human footprint index was not selected as a meaningful explanatory variable in the best model. Noteworthy, we found a better coupling between the velocity of species range shifts and the velocity of isotherms shifts for ectotherms in cold environments (i.e. close to mountaintops). The geographic isolation and habitat area constraints specific to mountaintops (e.g. sky islands) may exacerbate local extinction events through reduced population sizes as climate warms and habitat area shrinks, thus paying off part of the climatic debt for ectotherms living close to mountaintops. For instance, mass extinction events associated with climate warming and pathogen outbreaks have already been reported for several amphibian species endemic to mountainous regions [START_REF] Pounds | Widespread amphibian extinctions from epidemic disease driven by global warming[END_REF] . For endotherms, phanerogams and cryptogams, the slope of the relationship between the velocity of species range shifts and the velocity of isotherm shifts in mountainous systems is negative, especially under warm climates (Extended Data 9). However, we found that velocity values for range shifts along elevational gradients are always positive, except under very warm baseline temperature conditions (BT > 20°C). This indicates that endotherms, phanerogams and cryptogams are in general shifting their elevational ranges upslope to track shifting isotherms but are consistently lagging behind climate change. Isotherms may be shifting upslope at a pace that is simply too fast for species with limited dispersal abilities and long life spans, such as trees, to keep pace. Additionally, in the tropics, the higher importance of biotic interactions 32 may further impede the rate of range shifts over what is expected from climate change alone [START_REF] Hillerislambers | How will biotic interactions influence climate change-induced range shifts?[END_REF] . But most importantly, the global climatic grids that are currently available, and that we used here, may still be too coarse in spatial resolution (1 km 2 at best) to allow a reliable assessment of the true velocity of isotherm shifts experienced along mountain slopes. Indeed, the topoclimatic and microclimatic heterogeneity that is available across few metres in mountain systems, something that is not accounted for here, may provide a strong spatial buffer against climate warming [START_REF] Lenoir | Local temperatures inferred from plant communities suggest strong spatial buffering of climate warming across Northern Europe[END_REF] , allowing species to shift at relatively small spatial distances and seemingly "stay" [START_REF] Graae | Stay or go -how topographic complexity influences alpine plant population and community responses to climate change[END_REF] relative to the velocity of isotherm shifts that is measured at a coarser spatial resolution. Hence, the slow velocities of species range shifts that we observed in mountainous areas could also be the result of local compensation effects involving short distance escapes and species persistence within microrefugia.

General implications. To conclude, the coupling between species range shifts and isotherm shifts is not uniform across biological systems, confirming the lags observed in the biotic responses of terrestrial organisms to climate change [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF][START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a continental scale[END_REF] . Noteworthy, we demonstrate complex interactions between the velocity of climate warming, the degree of human pressures on the environment, historical temperature regimes and species characteristics. We suggest that commercial fishing may speed up the displacement of marine species distribution through resource depletion and population crashes at the trailing edge, whereas low constraints on dispersal in the oceans may allow marine species living close to their upper thermal limits to better track climate warming at the leading edge. On land, habitat loss and fragmentation due to land-use changes may impede the ability of terrestrial species to track shifting isotherms. These complex interactions need to be accounted for to improve scenarios of biodiversity redistribution and its consequences on human well-being [START_REF] Pecl | Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being[END_REF] under future climate change. The fact that marine species better track climate warming than terrestrial species also suggests that biodiversity redistribution will have more immediate and far-reaching consequences in the oceans than on land. For instance, community reshuffling [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF][START_REF] Bertrand | Ecological constraints increase the climatic debt in forests[END_REF] and the "tropicalisation" of temperate ecosystems [START_REF] Vergés | The tropicalization of temperate marine ecosystems: climate-mediated changes in herbivory and community phase shifts[END_REF][START_REF] Vergés | Tropicalisation of temperate reefs: Implications for ecosystem functions and management actions[END_REF] is likely to happen much faster in marine than in terrestrial systems, with more direct and rapid consequences on ecosystem health and functioning (e.g. increasing fish herbivory in kelp forests [START_REF] Vergés | Long-term empirical evidence of ocean warming leading to tropicalization of fish communities, increased herbivory, and loss of kelp[END_REF] ) as well as on the valuable services (e.g. fishery) and disservices (e.g. coastal erosion) oceans can provide.

However, it is important to bear in mind that our findings, as well as former syntheses on the topic, are still dependent on data availability and thus suffer from severe taxonomic and geographic biases. Despite a broad taxonomic coverage of the tree of life (Extended Data 3), species range shifts recorded in BioShifts 13 cover only 0.6% of the described biodiversity on Earth (N = 2,094,892 taxa). Besides, it is noteworthy that species range shifts in the southern hemisphere and in tropical regions in general are underrepresented. These limitations may affect our perception of species redistribution, and by consequence challenge global biodiversity conservation efforts [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF][START_REF] Feeley | Most 'global' reviews of species' responses to climate change are not truly global[END_REF] . It is thus more important now than ever to continue to study and document range shifts in areas and for taxonomic groups that have been so far somewhat neglected. Our database on species range shifts provides solid foundations to build a truly global monitoring of species redistribution. We thus call for future research perspectives linking our database on species range shifts with existing but scattered databases on species traits [START_REF] Kattge | TRY -a global database of plant traits[END_REF][START_REF] Kissling | Establishing macroecological trait datasets: digitalization, extrapolation, and validation of diet preferences in terrestrial mammals worldwide[END_REF][START_REF] Meiri | Traits of lizards of the world: Variation around a successful evolutionary design[END_REF][START_REF] Oliveira | a global database for amphibian ecological traits[END_REF][START_REF] Frimpong | Fish Traits: A Database of Ecological and Life-history Traits of Freshwater Fishes of the United States[END_REF][START_REF] Wilman | EltonTraits 1.0: Species-level foraging attributes of the world's birds and mammals[END_REF] to 308 improve our ability to anticipate biodiversity redistribution under climate change.

309

Methods

Literature search. We reviewed the scientific and peer-reviewed literature reporting climatedriven range shifts under contemporary climate change. By contemporary climate change, we here mean the period stretching from the beginning of the 19th century and onwards. As a general approach, we started from the reference lists of the most recent meta-analyses and syntheses on the topic 2-4 that we completed by regularly searching the scientific literature published between 2014 and 2018, following the same protocol as in Lenoir & Svenning 4 (Supplementary Fig. 1). Because of the clear focus on latitudinal and elevational range shifts in the scientific literature and the lack of information on the other geographical dimensions [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF] , we excluded several reports focusing exclusively on bathymetric or longitudinal range shifts. Broad inclusion criteria comprised studies: (i) focusing on relatively recent (since 1850s) distribution changes; (ii) based on occurrence or abundance data of at least one species; and (iii) only if studies were based on assessments covering at least two historical censuses with a minimum of 10 years between censuses. Hence, we excluded studies reporting distribution range changes from a single census (synchronous approach comparing data from different ontogenetic life stages of the same species like seedlings vs. adult trees) or based on historical patterns of species mortality obtained from climatic reconstructions only, without real occurrence or abundance data from at least two different time periods to confirm model outputs. We also excluded studies focusing exclusively on distributional range changes of invasive alien species. This selection procedure led to a total of 258 published and peer-reviewed studies for which we could extract data on species range shifts [START_REF] Pinsky | Marine taxa track local climate velocities[END_REF][START_REF] Engelhard | Climate change and fishing: a century of shifting distribution in North Sea cod[END_REF] .

We used Google Sheets to store the raw data on species range shifts in a dynamic and common file that we shared among co-authors, while always keeping a regular copy of the database saved on several computers to ensure backups. Once studies were clearly identified and stored as ".pdf" files in a common folder in Google Docs, each co-author picked studies, one by one, and entered data manually in the database. Some of the ".pdf" files were carefully annotated to help us quickly identify and recover any meaningful information in the main text or display items (e.g. tables or figures). When data on species range shifts were not directly available in the main text, in tables or in the appendices of the publication, we first contacted the corresponding authors and requested the data. In case of no positive response from the original authors and when data on range shifts could be extracted from published figures, we used the "WebPlotDigitizer" program (https://automeris.io/WebPlotDigitizer/). When range shifts were reported for more than one geographically distinct survey area or between more than two censuses (e.g. more than one resurvey of historical data), we considered them as independent case studies (N = 325). This data entry procedure led to a total of 30,534 range shift estimates at the species level (see Data Availability to access the database: BioShifts [START_REF] Comte | BioShifts: a global geo-database of climate-induced species redistribution over land and sea[END_REF] ).

Range shift estimates, as reported by the original authors, were coded as positive values if poleward in latitude, or upward in elevation, and negative otherwise (equatorward and downward). When the authors reported horizontal range shifts with both the magnitude and direction (i.e. azimuth) values, we used trigonometric relationships to transform these values into latitudinal range shifts for consistency with the main bulk of data available in the scientific literature. Next, we divided each range shift estimate by the study duration between two consecutive censuses (ending year -starting year + 1) to assess the rate or velocity of range shift (ShiftR), in kilometer per year along the latitudinal gradient and in meter per year along elevational gradients. In addition to the velocity of range shift at the species level, we also retrieved information at the case study level (N = 325), including methodological attributes known to potentially affect the velocity of range shift [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] : the starting year of the study (Start); the ending year of the study (End); the size of the study area (Area); the number of taxa in a study (Ntaxa) (continuous variable ranging from 1 to 4426; median = 21; mean = 122); the frequency of sampling (Sampling) (factor variable with four levels: "continuous"; "irregular"; or comparison of "two"; or "multiple" periods); whether range shift estimates were generated from "occurrence" or "abundance" data (PrAb) (factor variable with two levels); the spatial resolution of the raw data used to estimate range shifts (Grain) (factor variable with three levels: "fine" for data based on GPS coordinates with a spatial resolution lower than 10 km; "coarse" for data based on range maps or atlas grids with a spatial resolution greater than 100 km; and "medium" for intermediate situations); the quality of the approach used to estimate range shifts (Quality) (factor variable with four levels: "low" when no data cleaning procedures were performed before computing range shifts; "balanced" when data cleaning or resampling procedures were carried out to calculate range shifts on a balanced dataset; "modeled" when range shifts were obtained by computing the difference in the position of a given range parameter estimated from species distribution models (SDMs) independently calibrated during at least two different time periods (note that in this case SDM outputs represent the realized and not the potential species distribution for a given time period); and "resurveyed" when range shifts were calculated from paired designs such as permanent plots); and whether the "significance" of range shift estimates were assessed or "not" in the original study (Signif) (factor variable with two levels). To improve the balance in the number of observations among levels of a given factor variable, we merged some levels with poor data coverage together for the Sampling and Quality variables. For instance, the levels "continuous" and "irregular" were merged together with the level "multiple" such that Sampling was used in our analyses as a factor variable with two levels: "two" vs.

"multiple". Regarding the Quality variable, we merged the level "resurveyed" together with the level "balanced" such that Quality was used in our analyses as a factor variable with three levels: "low"; "balanced"; and "modeled". Still at the case study level (N = 235), we digitized the study region in Google Earth and used the resulting polygons to retrieve spatial information such as the total area covered by the study. If no clear maps delineating the study area was reported in the original study (e.g. map displaying the study region), we used national geographic boundaries or any meaningful spatial information from the text to delineate the study area. All spatial polygons were used to produce a geo-database (Extended Data 1-2).

Taxonomic harmonization. Before undertaking any taxonomic harmonization procedure, the last version of our database, dated April 2018, contained 13,570 entries of taxa at any taxonomic rank up to the genus level (i.e. subspecies, species and genus). Using the R programming language 301 , we assembled an R script in order to retrieve, for each taxonomic entry, the most recent accepted name and its associated classification. After a visual inspection for obvious syntax correction, three steps of taxonomic verification were performed. First, names were searched in the National Center for Biotechnology Information (NCBI) taxonomy database using the function "classification" from the R package "taxize" 302 . Then, in the same way, any remaining taxonomic entity not found in NCBI was verified with the Integrated Taxonomic Information System (ITIS) database. The full taxonomic classification was also retrieved during these two steps. Third, the last remaining taxonomic entities not found in NCBI and ITIS were checked using the Global Biodiversity Information Facility (GBIF) database, using the function "name_backbone" in the R package "rgbif". If we found a match, the corrected taxonomic entity was again checked in NCBI and ITIS by undergoing the previously mentioned procedure once again to retrieve a reliable taxonomic classification. Finally, only names at the species and the genus level were kept for the analyses (subspecies being aggregated at the species level).

Following this taxonomic harmonization procedure, the final number of taxa names in the database was reduced to 12,415. Because marine and terrestrial taxa shift at different rates and directions to potentially track the complex mosaic of local climate velocities 8 , we calculated the observed local velocity of temperature change (i.e. the spatial shift of isotherms over time) [START_REF] Loarie | The velocity of climate change[END_REF][START_REF] Burrows | The pace of shifting climate in marine and terrestrial ecosystems[END_REF] for each case study, following the approach used by Burrows et al. [START_REF] Burrows | The pace of shifting climate in marine and terrestrial ecosystems[END_REF] . We divided the temporal change in annual mean temperature observed over the studied period (°C.yr -1 ) by the corresponding spatial gradient (°C.km -1 or °C.m -1 ) as a measure of the velocity of temperature change (km.yr -1 or m.yr -1 ) [START_REF] Loarie | The velocity of climate change[END_REF] . The temporal gradient was calculated using time-series data from the CRU covering the period 1901-2016 at a spatial resolution of 0.5° (about 55 km at the equator) and from the Met office Hadley Centre observations datasets covering the period 1870-2018 at a spatial resolution of 1° (about 111 km at the equator) for terrestrial and marine studies, respectively. To do so, we regressed annual mean temperature (°C) values for all years throughout the study period as well as the two preceding years against time (yr) using linear regressions. When the starting year was prior to 1901 or 1870 for terrestrial and marine systems, respectively, we started the time series in 1901 or 1870 depending on the climate series. The slope parameter (°C.yr -1 ) of this model was then used as an estimate of the temporal gradient. For the sake of comparison with the rate of range shift usually calculated along the latitudinal and elevational gradients, we calculated the spatial gradient of annual mean temperature along the latitudinal (km.yr -1 ) and along elevational (m.yr -1 ) gradients, separately. This allowed us to assess both the latitudinal and elevational velocity of temperature change (LatVeloT and EleVeloT). To assess the latitudinal spatial gradient of annual mean temperature across land and sea, we used spatial grids from WorldClim and BIO-ORACLE, respectively, at a spatial resolution of 5 arc-minute (about 9.2 km at the equator). The WorldClim grid of annual mean temperature was downloaded at the finest spatial resolution, which is 30 arc-second (about 1 km at the equator), but aggregated at 5 arc-minute to be consistent with the spatial resolution of sea surface temperatures. Latitudinal spatial gradients were calculated as in Burrows et al. [START_REF] Burrows | The pace of shifting climate in marine and terrestrial ecosystems[END_REF] based on a 3 × 3 neighborhood sub-grid with the centre cell being the focal cell and its eight neighboring cells used to calculate the difference in temperatures for each northern and southern (resp. southern and northern in the southern hemisphere) pairs divided by the distance between them. Average differences (°C.km -1 ) for the focal centre cell were calculated, excluding any missing values (usually along coastlines), using weightings of 1 and 2 for cells diagonal and adjacent, respectively, to the focal centre cell. For the elevation gradient, we used the temperature data from the WorldClim grid of annual mean temperature at the finest spatial resolution (30 arc-second which is about 1 km at the equator) and calculated the spatial gradient across each case study using a linear model relating annual mean temperature (the response variable) to both elevation and latitude (the explanatory variables). We used latitude as a covariate in this model to account for the latitudinal variation in temperature observed within studies covering large spatial extents, i.e. elevation values close to the equator are not directly comparable, in terms of temperature, to elevation values close to the poles. The coefficient parameter along elevational gradients (°C.m -1 ) was then used as an estimate of the local adiabatic lapse rate. For the study areas that were larger in extent than the spatial resolution of the temperature grids, we computed the mean values of LatVeloT or EleVeloT throughout the entire study area by averaging values across all spatial grid cells overlapping with the polygons delineating the study area.

Additional drivers of range shifts.

As baseline temperature conditions may affect the rate at which species are shifting their distributions [START_REF] Guo | Land-use change interacts with climate to determine elevational species redistribution[END_REF] , we extracted annual mean temperature values during the year of the initial census (Start) as well as the two preceding years and calculated the mean (hereafter BT in °C). For terrestrial and marine systems, we used time-series data from CRU and the Met office Hadley Centre observations datasets, respectively. When the initial census of a given publication source was prior to 1901 or 1870 for terrestrial and marine systems, respectively, we used the oldest years available from the time series to compute baseline temperature conditions. Similar to climate velocity variables, when the study areas were larger in extent than the spatial resolution of the temperature grids, we computed the mean values of BT throughout the entire study area by averaging values across all spatial grid cells overlapping with the polygons delineating the study area.

As anthropogenic disturbances such as land-use intensity or industrial fishing may act as confounding factors on the velocity of range shift [START_REF] Guo | Land-use change interacts with climate to determine elevational species redistribution[END_REF] , we retrieved information on anthropogenic impacts for both the terrestrial and marine environment. For terrestrial systems, we downloaded the Global terrestrial Human Footprint maps for the year 2009 [START_REF] Venter | Global terrestrial human footprint maps for 1993 and 2009[END_REF] . These maps, at a spatial resolution of 30 arc-second (about 1 km at the equator), provide remotely-sensed and bottom-up survey information on eight variables measuring the direct and indirect human pressures on the environment: (1) the extent of built environments; (2) human population density; (3) electric infrastructure (night-light time); (4) crop lands; (5) pasture lands; (6) roads; (7) railways; and (8) navigable waterways acting like roads for people to access natural resources. All eight pressure variables were scaled by the original authors based on their degree of influence on the terrestrial environment. For instance, human population density and night-time lights were scaled between 0 and 10 while roads were scaled between 0 and 8. Scores for each of the eight individual threats were then summed and weighted by the original authors to make a composite map of the global human footprint index ranging from 0 to 50. For marine systems, we used the Global Map of Human Impact on Marine Ecosystems [START_REF] Halpern | A Global Map of Human Impact on Marine Ecosystems[END_REF] , also available at 30 arc-second resolution (about 1 km at the equator). This gridded dataset provides a cumulative impact score ranging from 0.01 to 90.1 for the minimum and maximum value, respectively. It was developed on the basis of expert judgment, to estimate ecosystem-specific impacts with respect to 17 anthropogenic drivers of ecological change (e.g. commercial shipping, demersal and pelagic fishing, ocean acidification, pollution). To allow comparison between terrestrial and marine systems, we rescaled both indices between 0 and 1 (standardized human footprint index or standardized HFI) and computed the mean per study area. The original authors have extensively validated HFI values against satellite imagery, yielding high confidence they represent conditions of human pressure on the environment [START_REF] Venter | Global terrestrial human footprint maps for 1993 and 2009[END_REF] .

Description: assessing geographic and taxonomic biases.

To evaluate spatial biases in the reporting of species range shift, we built 2° × 2° gridded maps, on top of which we overlaid the digitized polygons associated with the observations gathered in the database for both the terrestrial and marine realm, and separately for latitudinal and elevational range shifts. For each 2° × 2° grid cell, we also computed the relative proportion of ectotherms vs. endotherms for animals and phanerogams vs. cryptogams for plants and plant-like life forms (e.g. lichens and algae). We distinguished ectotherm from endotherm life-forms due to their contrasting sensitivity to temperature fluctuations in the environment, with ectotherms being unable to directly regulate their body temperatures as opposed to endotherms. The distinction between phanerogams and cryptogam life-forms allowed to contrast between two different reproduction strategies among chlorophyllous organisms: the evolved seed-bearing plants (angiosperms and gymnosperms) vs.

the other plant-like life forms reproducing by spores (ferns, mosses, lichens and algae). We then generated cartograms using the diffusion-based method for producing density-equalizing maps 303 . The number of range shift rates per 2° × 2° grid cell (i.e. sample size) was used to distort the map: the bigger the grid cell, the larger the sample size (Extended Data 1). We additionally estimated the phylogenetic coverage of the range shift database with respect to the whole tree of life described in the Open Tree of Life (https://tree.opentreeoflife.org) collapsed at the level of taxonomic classes and the total number of species recorded in the Catalogue of Life (http://catalogueoflife.org/).

Detection: estimating the velocity of range shifts per taxonomic class. Data coverage in our

database is very unbalanced between: the marine vs. terrestrial realm; the northern vs. southern hemisphere; and the margins vs. centroid of the species range (Supplementary Table 1). Besides, data on species range shifts do not even exist for some taxonomic classes in some of the combination of realm × hemisphere × position in the species range. For instance, dicots (Magnoliopsida) are exclusively terrestrial organisms while cartilaginous fishes (Chondrichthyes) almost exclusively live in marine habitats except for a few sharks and rays living in freshwater habitats during all or part of their lives. Thus, a single model to estimate the velocity of range shifts per taxonomic class while accounting for methodological biases [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF][START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF][START_REF] Feeley | Most 'global' reviews of species' responses to climate change are not truly global[END_REF] would be inappropriate. Hence, we divided latitudinal range shifts (N = 16,952) into a full factorial design 304 with eight experimental units based on all possible combinations of levels across three factor variables: biological system (marine vs. terrestrial); hemisphere (north vs. south); and range position (centroid vs. margins). We did the same for elevational range shifts (N = 13,582) except that there were only four possible experimental units (i.e. terrestrial systems only). To ensure robust fit, we further focused on taxonomic classes with more than 30 observations per experimental unit (N = 20 taxonomic classes fulfilling this sample size criterion) (Fig. 1, Supplementary Table 1), which reduced our sample size to 16,399 and 13,341 observations for latitudinal and elevational range shifts, respectively. Among the 12 possible combinations, only one combination (latitude × margins × terrestrial × south) could not be fulfilled due to a lack of range shift data (N = 8). This resulted in a total of 11 sub-models (i.e. factorial models) (Supplementary Table 1).

For each of the 11 factorial models for which the data were available, we built a linear mixed-effects model (LMM) relating the velocity of species range shift (ShiftR) for a given taxon (i.e. the response variable) against taxonomic Class, a factor variable with as many levels as the number of taxonomic classes within the focal experimental unit (e.g. Amphibia vs. Aves for latitudinal range shifts at the centroid of the distribution in terrestrial systems of the southern hemisphere) (Supplementary Table 1). Note that if a given factorial model only had data for one unique taxonomic class (e.g. Actinopterygii for latitudinal range shifts at the centroid of the distribution in marine systems of the southern hemisphere) (Supplementary Table 1), then the variable Class was not included in the fixed effects of the LMM. For the five LMMs focusing on the rate of range shift at the margins of the distribution, we added an extra factor variable (Margin) with two levels ("leading" vs. "trailing" edge) in the fixed effects, to provide robust estimates of the rate of range shift at both the leading and trailing edges. Given the complex structure of the database, involving repeated observations per taxonomic units (e.g. family, genus) or methodological levels, LMM is the most appropriate modelling approach 304 . This allowed to provide estimates of the velocity of range shifts per taxonomic class that are representative across all levels of a given methodological variable rather than providing estimates for each level separately, while accounting for taxonomic non-independence. More specifically, we included Genus as a random intercept term nested (or not: in case of singularity fit) within Family to account for potential taxonomic autocorrelation in the residuals of the models. In addition, because the different methodological approaches used in the scientific studies may also contribute to a non-negligible fraction of the variation in range shifts [START_REF] Brown | Ecological and methodological drivers of species' distribution and phenology responses to climate change[END_REF] , we considered several methodological variables as random intercept terms in the LMMs (Area, Start, Ntaxa, Sampling, PrAb, Grain, Quality and Signif). To be included in the random part of our LMMs, the continuous variables Area, Start and Ntaxa were first transformed into factor variables with four levels each, using the respective quantiles as cutting points. Then, for each factorial model separately, we selected only the set of uncorrelated variables with at least two levels having more than four observations. We used the "lmer" function from the "lme4" package 305 in the R programming language 301 . We used a model selection procedure where the best random effect structure was identified by testing all the combinations of random factors and selecting the one with the lowest Akaike information criterion with small-sample correction (AICc). To compare AICc values among candidate models, we set the restricted maximum likelihood argument to "FALSE" in the "lmer" function (i.e. REML = FALSE for maximum likelihood) 305 . To ensure robust estimations, all the singular fits were removed from the list of candidate models prior to model selection. In case of singular fits across all candidate models, we used case study (Source) as the unique random intercept term. If the random intercept term Source also led to a singular fit, then we used a linear regression model (LM). For each of the LMMs (or LMs in case of singular fits for all the candidate models) focusing on the velocity of range shift at the margins of the distribution, we also included an interaction between Margin and Class that we tested against a model without the interaction term in the fixed effects based on the AICc value. When the absolute difference in AICc value between these two candidate models was greater than two, we selected the model with the lowest AICc value. Otherwise, in case of equivalent competing models, we selected the one with the interaction effect between Margin and Class considering that it allows flexible range shift estimations at the trailing and leading edge. Once the best LMM was selected for each factorial model (Supplementary Table 1), we set REML to TRUE 305 to extract coefficient estimates among the different levels of the factor variables Class and Margin. To test whether the estimated rate of range shift for a given taxonomic class and at a given position within the range was significantly different from zero, we reran each of the 11 selected best models using a bootstrap approach (N = 5,000 iterations). From these 5,000 estimates, we computed the mean and median velocity of range shift as well as the standard deviation and 95% confidence interval per taxonomic class. Finally, to assess the overall goodness-of-fit of the different factorial models, as well as to compare the relative importance of biological versus methodological effects on the rate of range shift, we computed the marginal (i.e. variance explained by the fixed effects) and conditional (i.e. variance explained by both the fixed and random effects) R 2 values 306 of each bootstrap iteration and for each factorial model using the "r.squaredGLMM" function from the "MuMIn" package in the R programming language 301 . Attribution: coupling between species' range shifts and isotherms' shifts. We assessed the coupling between the velocity of species range shifts and the velocity of isotherm shifts (VIS) using LMMs built separately for the latitudinal and elevational gradients. We specified the velocity of species latitudinal (km.yr -1 ) or elevational (m.yr -1 ) range shifts as the response variable and either the latitudinal or elevational VIS (LatVeloT / EleVeloT; continuous variables) as the main explanatory variable. To account for potential interacting effects on the relationship between the velocity of range shifts and VIS, we added several covariates in our models: baseline temperature (BT; a continuous variable); standardized human footprint index (standardized HFI; a continuous variable representing human pressures on the environment bounded between 0 and 1); and life forms (LF) (a factor variable with 4 levels: ectotherm, endotherm, cryptogam, phanerogam). As temperature regimes and human pressures on the environment are not directly comparable between lands and oceans, we further modeled the coupling between the velocity of species latitudinal range shifts and VIS in latitude (LatVeloT) separately for the marine and terrestrial realm. We tested for all two-way interaction terms between each covariate (BT, HFI and LF) and VIS (either LatVeloT or EleVeloT). We also tested for a unimodal relationship between the estimated rates of range shifts and baseline temperature conditions (BT) using a second-order polynomial term. The variables Position within the range (a factor variable with 3 levels: trailing edge, centroid and leading edge) and Hemisphere (a factor variable with 2 levels: North vs. South) were not incorporated as covariates in the models as both variables had no effect to explain the variation in the rates of latitudinal and elevational range shifts per taxonomic class (Supplementary Table 2).

Similar to the LMMs developed at the taxonomic class, the aforementioned explanatory variables were used as fixed effects in LMMs, whereas the methodological attributes (Area, Start, Ntaxa, Sampling, PrAb, Grain, Quality and Signif) were used as random intercept terms.

Starting from the beyond optimal model (full model with all fixed effects) 305 separately for the velocity of latitudinal range shifts in marine and terrestrial systems as well as for the velocity of elevational range shifts, we tested all model combinations and selected the best model based on the lowest AICc value, setting REML to FALSE (i.e. maximum likelihood) for model selection and then to TRUE to estimate the coefficients once the best model was selected 305 . We first selected the random effect structure after removing singular fits, using the exact same procedure as for the models used to estimate the mean velocity of range shift per taxonomic class. We then selected the fixed effect structure, keeping the previously identified random structure constant.

All continuous variables (LatVeloT, EleVeloT, BT and HFI) were standardized to z-scores using the "gscale" function 307 from the "jtools" package in the R programming language 301 . This function standardizes each value of a given variable by subtracting it from the mean and dividing it by two times the standard deviation of the focal variable (instead of one time as more commonly done). This rescaling formula is recommended over the traditional formula of dividing by one time the standard deviation because it allows direct comparisons of model coefficients with untransformed binary predictors 307 . For the sake of consistency, we focused on the set of species belonging to the taxonomic classes with more than 30 observations, resulting in 16,521 (1,403 marine vs. 15,118 terrestrial) and 13,459 observations for latitudinal and elevational range shifts, respectively. The 95% confidence intervals around each of the estimated coefficients were calculated using bootstraps (N = 5,000 iterations), similar to the models used to estimate the mean velocity of range shift per taxonomic class.

Finally, to illustrate the capacity of species to track the shifting isotherms, we mapped the predicted slopes for each combination of the predictors identified in the best models, separately for latitudinal (marine or terrestrial) and elevational range shifts. A slope value of one between the velocity of species range shifts and the velocity of isotherm shifts indicates a perfect coupling with species closely tracking the shifting isotherms. To do so, we built a 2° × 2° gridded map, on top of which we overlaid the digitized polygons associated with each observation used in the previous models. We then generated cartograms using the diffusion-based method for producing density-equalizing maps 303 . As before (see section entitled "Detection: assessing geographic and taxonomic biases"), the number of range shift rates per grid cell (i.e. sample size) was used to distort the map: the bigger the grid cell, the larger the sample size. Finally, we tested whether the slope estimated for each 2° × 2° grid cell (i.e. according to the grid-specific baseline temperature and the standardized human footprint index) significantly differed from a value of one (indicating a perfect coupling), based on 5,000 bootstrap iterations.

several grants from the Agence Nationale de la Recherche (TULIP ANR-10-LABX-41, CEBA ANR-10-LABX-25-01).

Competing Interests Statements

Authors declare no competing interests. Reprints and permissions information is available at www.nature.com/reprints or to methodological attributes (random effect terms in the models) for each of the 10 factorial models for which we had data (Supplementary Table 1). Each factorial model represents a combination of positional parameter (Cen: centroid; Mrg: margins) × spatial gradient (L: latitude; E: elevation) × biological systems (M: marine; T: terrestrial) × hemisphere (N: north; S: south). Note that the "Margin" factor variable with two levels (leading edge vs. trailing edge) was only tested in model combinations focusing on margins (Mrg). Error bars represent the distribution of 5,000 bootstrap iterations. suggest that species are not tracking the shifting isotherms. Note that negative slopes do not necessarily mean that species are shifting in the opposite direction to isotherm shifts (see Fig. 5).

The number of range shift estimates (i.e. sample size) in each grid cell was used to distort the map: the bigger the grid cell, the larger the sample size. Grid cells with a black and bold border display areas where species are closely tracking the shifting isotherms (i.e. slope parameter not significantly different from 1 based on 5,000 bootstrap iterations).

  Climate velocity. Using the spatial information obtained from the digitized polygons as well as the temporal information (Start and End years) available from each of the 258 publication sources, we retrieved basic temperature information to calculate the velocity of temperature change throughout the study period. Terrestrial climate data were obtained from WorldClim v. 1.4 (http://www.worldclim.org/) and the Climate Research Unit (CRU) TS v. 3.23 (https://crudata.uea.ac.uk/cru/data/hrg/) while marine climate data were obtained from BIO-ORACLE (http://www.bio-oracle.org/) and the Met office Hadley Centre observations datasets (https://www.metoffice.gov.uk/hadobs/hadisst/).

Fig. 1│Taxonomic

  Fig. 1│Taxonomic coverage. Number of taxa (in parenthesis), in log scale, per taxonomic class: from the least (top) to the most (bottom) studied taxonomic class. Only taxonomic classes with more than 30 observations per factorial model are displayed.

Fig. 2│Sources of

  Fig. 2│Sources of variation in species range shifts. Proportion of explained variation either related to the taxonomic class and position at the range margin (fixed-effect terms in the models)

Fig. 3│Mean velocity

  Fig. 3│Mean velocity of species range shifts per taxonomic class. Estimated velocity of range shift per taxonomic class (i.e. effect size) in km.yr -1 and m.yr -1 for (a) latitudinal and (b) elevational range shifts, respectively, after accounting for methodological variation. Outputs are displayed for all possible combinations of positional parameter (TE: trailing edge; CE: centroid vs. LE: leading edge) × hemisphere (N: north; S: south) × biological systems (M: marine; T: terrestrial). Violin plots represent the distribution of 5,000 bootstrap iterations. Stars show significant deviations from zero shift (*: P < 0.05; **: P < 0.01; ***: P < 0.001).

Fig. 4│Degree of

  Fig. 4│Degree of coupling between species range shifts and isotherm shifts. Models outputs in terms of (a) proportion of explained variation and (b) effect size related to the velocity of isotherm shifts (VIS), baseline temperatures (BT + BT 2 ), standardized human footprint index

Fig. 5│Main determinants

  Fig. 5│Main determinants of the velocity of species range shifts. Results are displayed along (a) elevational and (b, c, d) latitudinal gradients for both the (a, b) terrestrial and (c, d) marine realms. Panel a shows the interaction effect between baseline temperatures and the velocity of isotherm shifts in elevation for ectotherms. Panel b shows the interaction effect between the standardized human footprint index and the velocity of isotherm shifts in latitude for terrestrial ectotherms. Panel c shows the interaction effect between baseline temperatures and the velocity of isotherm shifts in latitude for marine ectotherms while setting the standardized human footprint index to its median value in the database. Panel d shows the interaction effect between the standardized human footprint index and the velocity of isotherm shifts in latitude for marine ectotherms while setting baseline temperatures to the median value in the database. The two white lines and the white hatching represent the range of conditions for which marine ectotherms closely track the shifting isotherms in latitude (i.e. slope parameter not significantly different from 1 based on 5,000 bootstrap iterations). Note that negative slopes do not necessarily indicate species range shifts in the opposite direction to isotherm shifts, unless the signs of the two estimates (for a given combination of baseline temperatures and standardized human footprint index) are opposite.

Fig. 6│Maps of

  Fig. 6│Maps of the degree of coupling between the velocity of species range shifts and the velocity of isotherm shifts. Cartograms show the predicted slope coefficient between the velocity of species range shifts and the velocity of isotherm shifts per 2° × 2° grid cell along (a)

  shift (km/year)

Acknowledgments

We acknowledge authors who kindly sent us their data on species range shifts estimates. In particular, we are thankful to Kristin Kleisner and Christopher Hassall who kindly provided data on behalf of the NOAA Northeast Fisheries Science Center, The Nature Conservancy, the British Arachnological Society and the Spider Recording Scheme. Finally, we wish to acknowledge

Data Availability

The data supporting the findings of this study are available in the BioShifts geo-database in the Figshare Digital Repository 13 available at https://figshare.com/s/ebd19485a00757ababb0

Code Availability R scripts used in the analyses are available at https://figshare.com/s/ebd19485a00757ababb0

Author Contribution

JL, LC, JM and GG initiated and conceived the project idea. LC and JL built the general structure of the database. GG, LC, RB, TH and JL reviewed the scientific literature and filled the database throughout the project duration. GG ensured data curation. LB and LC carried out the taxonomic harmonization of the database with help from JM. TH linked the taxonomic backbone of the database to the Open Tree of Life (https://tree.opentreeoflife.org) and to Catalogue of Life (http://catalogueoflife.org/) to produce a visualization of the phylogenetic coverage of the database. GG, LC, JL and RB prepared the set of methodological variables included as covariates in the subsequent analyses. RB and JL analyzed the data with help from LC, LB and GG. TH, RB and JL produced all the figures. JL wrote the manuscript with contribution from all co-authors. JL and RB contributed equally to this work.