dataesr/scientific-paragraphs-categorization huggingface dataset:
According to Figure  2 , this domain was more severely affected by the heatwave, and is represented by Figure  9c. Firstly, the large error between all the experiments and the observations for the start of the growing season is noticeable (Figure  9a ). Secondly, from March to June LDAS_HRES analysis as well as LDAS_fc_d2 and LDAS_fc_d8 are only slightly correcting for this issue. This is a known issue as already mentioned by  [18] , the CO 2 -responsive version of ISBA is such that during the growing phase, enhanced photosynthesis corresponds to a CO 2 uptake, which results in vegetation growth from a prescribed LAI minimum threshold (1 m 2 m -2 for coniferous forest or 0. 3 m 2 m -2 for other vegetation types). These thresholds are likely too low and are currently being revisited using the CGLS LAI long term dataset.
Combining ERA5 atmospheric re-analysis, HRES analysis and its forecast within LDAS-Monde is highly relevant to foster research for land applications at various timescales from daily to annual. The use of HRES data to force LDAS-Monde is very promising and it can be complemented by ECMWF 51-member ensemble forecasts (~18 km spatial resolution). Use of the ECMWF ensemble in LDAS-Monde could help to capture uncertainties in the representation of LSVs. It would open the possibility to anticipate the impact of heatwaves at monthly temporal scales using a probabilistic method. Moreover, since one member of the ensemble is similar to HRES at a coarser spatial resolution, and the ensemble is available up to 15-days lead time (twice a day and up to 45 days twice a week) it can be used to test longer range forecast of LSVs than when using HRES.
, 2007; Rignot and Kanagaratnam, 2006; Luckman et al. , 2006]. In comparison, bed lubrication accelerates outlet glaciers by only 8 -10% during months of peak melting  [Rignot and Kanagaratnam, 2006; Joughin et al. , 2008]. [3] The ice sheet mass loss was estimated by  Rignot and Kanagaratnam [2006]  (hereinafter referred to as RK) in 1996, 2000 and 2005 from a comparison of ice discharge, D, and surface mass balance, SMB. Ice discharge was obtained from measurements of ice motion by satellite interferometric synthetic-aperture radar data (InSAR) and ice thickness by airborne radio echo sounding. Surface mass balance was derived from a positive degree day model (PDD) with melt retention and a snow accumulation grid, both representing the 1961 -1990 average. The spatial estimates of SMB were combined with the temporal variability in whole ice-sheet SMB from  Hanna et al.
Despite mounting evidence of species redistribution as climate warms, our knowledge of the coupling between species range shifts and isotherm shifts is limited. Compiling a global geo-database of 30,534 range shifts from 12,415 taxa, we show that only marine taxa closely track the shifting isotherms. In the oceans, the velocity of isotherm shifts interacts synergistically with anthropogenic disturbances and baseline temperatures such that isotherm tracking by marine life happens either in warm and undisturbed waters (e. g. Central Pacific Basin) or in colder waters where human activities are more pronounced (e. g. North Sea). On lands, increasing anthropogenic activities and temperatures negatively impact the capacity of terrestrial taxa to track isotherm shifts in latitude and elevation, respectively.
To fill this knowledge gap, we compiled range shifts for marine and terrestrial taxa from an exhaustive literature review building on and updating the most recent syntheses on climaterelated range shifts (see section Materials and Methods below). Based on this review, we compiled a geo-database (Table  S1  and Figs. S1-S2 in the Supplementary Materials) of 30,534 range shifts covering four kingdoms (Bacteria, Plantae, Fungi and Animalia), 20 phyla, 56 taxonomic classes and 12,415 taxa (Fig. S3 ). Focusing on latitudinal and elevational gradients, we first provided robust estimates of the velocity of range shift for the 20 most studied taxonomic classes (Fig. 1 ) by fitting models accounting for different methodological attributes among studies  (11)  (Fig. S4 ). Models were arranged in a full factorial design of geographical gradient (latitude vs.
Marine taxa (~80% being ectotherms in the database; Fig. S2C ) have moved towards the poles at a mean (±s. e. m. ) pace of 5. 92±0. 94 km. yr -1 , almost six times faster than terrestrial taxa (1. 11±0. 96) (one-way ANOVA: F = 12. 68; df factor = 1; df residuals = 45; P = 0. 0009). These rates far exceed the mean latitudinal range shift (0. 61±0. 24) reported by the seminal synthesis  (1). Yet, we found lower velocities than the values reported by subsequent syntheses focusing chiefly on terrestrial taxa (1. 76±0. 29) (2)  or exclusively on marine taxa (7. 20±1. 35) (3). Contrarily to a previous report (3) but according to expectations from thermal tolerance limits of marine ectotherms  (14) , our results suggest that both the leading (6. 02±1. 77) and trailing (6. 49±2. 13) edge of marine taxa are equally responsive to warming (one-way ANOVA: F = 0.
Note that the same pattern is observed for the three other studied life forms of the terrestrial realm, but the intercept is lower (see Table  S5  and Fig. S13 ). 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 (see. Fig. S11  for marine cryptogams). 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 (see. Fig. S11  for marine cryptogams). 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.
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 non-introduced taxon; and (iii) only if studies were based on assessments covering at least two historical censuses. Hence, we excluded studies reporting distribution changes from a single census (synchronous approach) or based on historical patterns of species mortality obtained from climatic reconstructions only, without real occurrence or abundance data to confirm model outputs. Data were only included in our database if the range shift estimates were reported at the taxa (i. e. subspecies, species or genus) level in the original research. When data on species range shifts were not directly available in the text or in tables in the main text or in the appendices of the publication, we first contacted the corresponding authors and requested the data.
In case of no positive response, we used the "WebPlotDigitizer" program (https://automeris. io/WebPlotDigitizer/) to extract range shift values from figures, when possible. This drastic selection procedure to ensure data quality led to a total of 30,534 range shift estimates at the taxa level extracted from 258 publications (Table  S1 ). When range shifts were reported for more than one geographically distinct survey area or more than one period delimited by two censuses (e. g. more than one resurvey of historical data), we considered them as independent case studies (N = 325). For each of these 325 case studies, 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 (Figs. S1-S2 ). 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 value by the study duration between two consecutive censuses (ending year -starting year + 1) to assess the rate of range shift (ShiftR), or the velocity of biological shifts, in kilometer per year along the latitudinal gradient and in meter per year along the elevational gradient.
In addition to the rate or velocity of range shift at the taxa level, we also retrieved information at the case study level, including: the size of the study area (Area); the starting year of the study (Start); the ending year of the study (End); as well as several other methodological attributes known to potentially affect the velocity of range shift  (11). More specifically, we considered: (i) the number of taxa in a study (Ntaxa) (continuous variable ranging from 1 to 4426; median = 21; mean = 122); (ii) the frequency of sampling (Sampling) (factor variable with four levels: "continuous"; "irregular"; or comparison of "two"; or "multiple" periods); (iii) whether range shift estimates were generated from "occurrence" or "abundance" data (PrAb) (factor variable with two levels); (iv) 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); (v) 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 procedure were carried out to calculate range shifts on a balanced dataset; "modeled" when range shifts were obtained from model outputs; and "resurveyed" when range shifts were calculated from paired designs such as permanent plots); and (vi) whether the "significance" of range shift estimates were assessed or "not" in the original study (Signif) (factor variable with two levels).
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 (23), 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"  (24). 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.
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 (Table  S1 ), 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/). 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)  (6, 7)  for each case study, following the approach used by Burrows et al. (7). 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 ) (6). 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.
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 the elevational gradient (°C. m -1 ) was then used as an estimate of the local adiabatic lapse rate.
Additional drivers of range shifts As baseline temperature conditions may affect the rate at which species are shifting their distributions  (25)  and the rate at which communities are reshuffled  (19) , 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 BaselineT 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 BaselineT throughout the entire study area by averaging values across all spatial grid cells overlapping with the polygons delineating the study area.
For marine systems, we used the Global Map of Human Impact on Marine Ecosystems  (12) , 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.
Description: assessing geographic and taxonomic biases To evaluate spatial biases in the reporting of range shift rates, 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  (26). 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 (see Fig. S1 ). 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.
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 (Table  S2 ). 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  (4, 11, 22)  would be inappropriate.
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, sometimes involving several observations for a given family or genus or per level of a methodological variable, LMM is the most appropriate modelling approach  (27). 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. 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.
We first selected the random effect structure after removing singularity fits, using the exact same procedure as for the models used to estimate the mean velocity of range shift per taxonomic class, and then the fixed effect structure, keeping the previously identified random structure constant. All continuous variables (LatVeloT, EleVeloT, BaselineT and HFI) were standardized to z-scores using the "gscale" function  (30)  from the "jtools" package in the R programming language  (23). This function standardizes each value of a given variable by subtracting it from the mean and dividing it by two times, instead of one time, the standard deviation of the focal variable. 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  (30).
Geographical bias in our knowledge of climate-induced range shifts For latitudinal range shifts on lands and in the oceans, our analysis confirms an important knowledge gap in the tropics  (4, 23). Indeed, there is a strong spatial unbalance in the database, such that the most developed regions of the northern hemisphere are better represented in terms of latitudinal range shifts: Western Europe; North America; the North Atlantic Ocean; and the North Sea (Fig. S1 ). By contrast, the southern hemisphere is underrepresented with Australia and the Tasmanian Sea being the most studied regions. We observed a similar pattern for elevational range shifts on lands. Again, mountain ranges from the most developed countries are better covered (i. e. Rockies and European Alps) than elsewhere, with tropical mountains poorly represented.
Hence, the geographical coverage of the data on species range shifts is far from global  (23). Finally, there is also a geographical bias among the different life forms (i. e. ectotherms, endotherms, phanerogams and cryptogams) covered in the database. Ectotherms are overrepresented for latitudinal range shifts in marine systems, whereas data on latitudinal range shifts for endotherms (mostly birds) dominate in North America, Finland and Australia (Fig. S2 ) where active groups of researchers are working on well-maintained databases of bird censuses.
23 μmol kg -1 (based on replicated samples, n = 21). Aragonite saturation and pCO2 were calculated using pH and A T measurements and the seacarb package  54  in R with dissociation constants K 1 and K 2  55  , Kf 56 and Ks  57. Seawater chemistry results are presented in Supplementary Tables  1  and 2. Shell-size analysis. Shell size was determined as the maximum shell length parallel to the hinge using brightfield microscopy and image analysis in ImageJ software. All statistical analyses were conducted in R (v. 3. 5. 3). As larval shell length data did not pass normality tests (Shapiro-Wilk test), shell-size was logtransformed to allow parametric statistical analysis. We tested the effect of day, treatment, and the interaction of the two using linear-mixed effects models, with day and treatment as fixed effects and replicate bucket as a random effect (lmer).
We visualized patterns of genetic variation throughout the experiment with PCA (prcomp function in R). Prior to PCA, the allele frequency matrix was centered and scaled using the scale function in R. Only larval samples that encompassed the full phenotypic distribution within a particular bucket were included in this analysis. In other words, the rows of the allele frequency matrix corresponding to larval samples that were selectively segregated based on shell size were removed, and PCA was run using the day 0 larval population and larval samples collected from each treatment on days 6, 26, and 43. A separate PCA was then implemented using allele frequency data from the day 0 larval population and all day 6 larval samples, which included discrete size groups from each treatment. This analysis thus explicitly examined a genomic signature of the individuals that were phenotypically distinct.