Ecology

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More

All of the statistical analyses were conducted using the R programming language version 4.0.5 within the RStudio IDE version 1.4.124,25. Data visualization and processing were performed with the ‘tidyverse’ collection, ‘foreach’ and ‘doParallel’ R packages26,27,28. Geographical Information System (GIS) operations on raster and vector files were conducted using the ‘sf’, ‘exactextractr’ and ‘raster’ R packages29,30,31.Data sources and pre-processingBiodiversity dataWe used the North American Breeding Bird Survey dataset as our source of biodiversity data due to its long temporal coverage and spatial extent14,32. The BBS is composed of bird species abundance records collected since 1966 from over 4,000 survey routes across the countries of Mexico, USA and Canada. For this study we focused solely on routes in the USA, due to their longer time dimension. Data collection follows public access roads that are 24.5 miles long (circa 39.2 km) using a point count protocol whereby routes are surveyed every half-mile (800 m) for a total of 50 stops. At each stop, observers stand for 3 min and record the species and the abundance of every bird seen or heard within 400 m of their location. The routes are surveyed by volunteers with experience in bird observation, and surveys are conducted from late April to July to capture the peak of the breeding season.We selected the years 2001 and 2016 as the two timepoints of our analysis. This 15-year timeframe corresponded to the longest possible timespan for which land cover data products were available at high spatial resolution18. Before analysis, we subset the BBS dataset by removing routes that had incomplete survey lengths (less than 50 point count stops, indicated by the RouteTypeDetailID field value being less than 2 in the extracted BBS dataset) or that were surveyed under adverse weather conditions such as high wind and rain (as indicated by the Run Protocol ID field being equal to 1), which could affect bird occurrence and detectability. Following this filtering process, the total number of BBS routes analysed was 960 (Extended Data Fig. 1).For higher precision when inferring the relationship between avian diversity and environmental variables, we subdivided each route into five segments of equal length, consisting of 10 count locations each. This approach was motivated by the need to more closely associate bird communities with the land cover composition in the area in which they are found. To minimize the spatial autocorrelation between adjacent segments and avoid overlaps in landscapes analysed, we filtered the data to keep only the first, third and fifth segment of each route. These segments therefore formed our sampling unit used in all analyses.We recognized that environmental conditions and stochastic trends in populations could introduce variability in biodiversity calculated from bird community data. We therefore extracted, for each segment and each species, the average population count across a 3-year period centred on our two timepoints (2000, 2001, 2002 and 2015, 2016, 2017)33. We then calculated the mean abundance of each species across these 3 years.The effect of observer experience34,35,36 was accounted for by sourcing the observer ID responsible for each route at each timepoint and including it as a random effect in the legacy model (see ‘Model development’ section). We also controlled for the time of day as it is plausible to expect visibility and avian species activity patterns to vary between early morning and later parts of the day. Time of day for each segment was calculated by averaging across the start and end time data entries associated with each route, and then including this as a covariate in both the legacy and equilibrium models (see ‘Model development’ section). However, we did not model detectability issues associated with traffic noise and disturbance for two reasons. First, all BBS surveys are conducted along public access roads with a vehicle, so the disturbance is expected to be similar across sites. Second, previous studies have found no clear evidence for noise being the main cause for reduced bird abundance near roads37.Following these procedures, our processed BBS dataset included entries of mean abundances of each species for a total of 2,880 segments, corresponding to segment 1, 3 and 5 of 960 routes (Extended Data Figs. 1 and 2). For each segment, at each timepoint we calculated different measures of alpha diversity following the Hill numbers framework38. We then selected to use the effective number of species at q = 1, calculated as the exponential of the Shannon–Wiener Index38. The effective number of species at q = 1 sits at the theoretical half-way point between the classic species richness measure that accounts only for the absolute number of species (q = 0) and the Berger-Parker dominance index (q = infinity), which instead only reflects the most common species. Thus, the effective number of species is a robust alternative to species richness, which does not take account of species rarity or detectability and can thus lead to biased biodiversity estimates16,17.Land cover and environmental dataLand cover data for the US for our focal years of 2001 and 2016 were sourced from the open-access NLCD CONUS products developed by the US Geological Survey (USGS)18,39. The NLCD products are high-resolution (30 m pixel dimensions) classified raster files covering the land area of the whole USA. This dataset provides us with the opportunity to look at finely gridded spatio-temporal changes in a landscape over a relatively long timeframe of 15 years, while utilizing data collected and analysed with the same methods (for example, land use classification algorithms).To reduce the number of potentially collinear explanatory variables included in our models, we aggregated the land cover variables provided by the NLCD dataset. We summarized these to five land cover categories: ‘urban’ (an aggregate of the Developed-Open Space (subclass 21), Developed-Low Intensity (22), Developed-Medium Intensity (23) and Developed-High Intensity classes); ‘forest’ (an aggregate of the Deciduous Forest (41), Evergreen Forest (42) and Mixed Forest (43) classes); ‘grassland’ (an aggregate of the Shrub (52), Grassland/Herbaceous (71) and Pasture/Hay (81) classes); ‘cropland’ (cultivated Crops (82) subclass) and ‘wetland’ (an aggregate of the Woody Wetland (90) and Herbaceous Wetland (95) classes). The Perennial Ice/Snow (12), Open Water (11) and Barren Land (31) classes were excluded from the analysis as they were very uncommon in our dataset. The distribution and total area of the land cover categories across the US are shown in Supplementary Figs. 1 and 2. Temperature data were sourced from the 30 arc-seconds gridded PRISM climate database19 and were extracted as the mean across May and June for each group of years from which bird abundances were taken.We first sampled the landscape surrounding each segment using a range of buffer shapes and sizes, and then selected the buffer type on the basis of the capacity of each buffer type to explain the response variable. The types of buffers that we explored were: a circular buffer around the centroid of the polygon defined by the vertices of each segment (4,000 m radius) and a series of three buffers around the segment line (500 m, 2,000 m and 4,000 m radius). The best fit was given by the smallest buffer size of 500 m, shown in Extended Data Fig. 2, which also coincides with the BBS protocol effective counting distance of 400 m and more closely reflects the size of bird territories14. Land cover variables were computed as a percentage of the total buffer area. Change in percentage points for each land cover type between the 2 years was computed by subtracting the values at the two timepoints. A change product is also provided by the USGS databases40, but it does not meet our needs because it considers land cover changes based on a ranking. Nonetheless, a comparison of urban land cover change between the timepoints showed a similar result (Supplementary Fig. 4). Land cover data were processed geospatially using the NAD 83 Conus Albers Coordinate Reference Systems projection, EPSG 5070.Model developmentTheoretical backgroundWe developed a statistical model that conceptualized extinction debts and colonization credits by combining the following two concepts: (1) the settled biodiversity of avian communities in a given landscape composition (that is, a system at equilibrium) and (2) the lagged response in the species diversity in a given landscape due to recent land cover changes (that is, a system moving to a new equilibrium). We reasoned that, given enough time, and with no further changes in land cover, the effective number of species at a given location would eventually equilibrate. The equilibrium distribution of the effective number of species emerges with the waning of the legacy effect of previous landscape compositions in encouraging or impeding the recruitment and survival of particular species. We did not model these ecological mechanisms directly, but instead expressed the equilibrium of the effective number of species, and the rate of approach to this equilibrium, as empirical functions of environmental covariates. It is important to keep in mind that during a finite time interval following environmental change, it is possible that our observations of effective number of species represent a system in a transitory state towards its new equilibrium. Yet, environmental changes may occur at rates that never allow the system to equilibrate. Although the equilibration processes are latent (that is, not amenable to direct observation), the combination of equilibrium and temporal legacy components into an integrated model, applied to a dataset with extensive environmental replication (due to spatial expansiveness), has allowed us to retrieve distributions for all relevant model parameters (see below).Model overviewThe observed effective number of species Rs,t at site s in year t for t = t1,t2 is modelled as a normally distributed variate with mean μs,t and standard deviation σ$$R_{s,t} approx mathrm{Normal}left( {mu _{s,t},sigma } right)$$
(1)
We assume that, under landscape change, the system is in a state of flux and that the data are from observations witnessing the transition between two (unattained) equilibria. The expected state of the system at any given point in time, μs,t, was formulated as a mixture of past and future equilibrium distributions (that is, a weighted average of the two distributions, where the weights are given by the complementary proportions ω and 1 − ω)$$mu _{s,t} = fleft( {x_{s,t_2};beta } right)omega left( {{Delta}x_{s,t_1,t_2};gamma } right) + fleft( {x_{s,t_1};beta } right)left( {1 – omega left( {{Delta}x_{s,t_1,t_2};gamma } right)} right)$$
(2)
Here, the function f describes the equilibrium distribution of the effective number of species as a function of the configuration of the local environment, captured in covariates xs,t. The weighting function ω depends on covariates ys,t derived from the difference in the local land cover between 2016 and 2001 (that is, it is a function of the land cover change that has taken place). The mixture weights ω and (1 − ω) determine the relative importance of the two equilibrium distributions (past or current). If ω = 1, the interpretation is that the new equilibrium distribution has been completely attained, and thus the current (2016) effective number of species is entirely explained by the current (2016) land cover. Conversely, if ω = 0, the current effective number of species is entirely explained by the past (2001) land cover. The vectors of parameters β and γ, presented in equation (2), are inferred from model fitting.We also augmented equation (2) with a function g of static covariates and random effects z that we expect to have an impact on the effective number of species. Thus, the model comprised equilibrium components, a temporal legacy component and static covariates:$$mu _{s,t} = fleft( {x_{s,t_2};beta } right)omega left( {{Delta}x_{s,t_1,t_2};gamma } right) + fleft( {x_{s,t_1};beta } right)left( {1 – omega left( {{Delta}x_{s,t_1,t_2};gamma } right)} right) + gleft( {z_s;alpha } right)$$
(3)
in which (fleft( {x_{s,t};beta } right)) are the equilibrium components for the two timepoints, (omega left( {Delta x_{s,t_1,t_2};gamma } right)) is the temporal legacy component, and (gleft( {z_s;alpha } right)) is the function that captures the static covariates and random effects, with α being the estimated static covariates parameter effects.Equilibrium componentsWe defined the equilibrium distribution of the effective number of species at a given timepoint as a function (fleft( {x_{s,t};beta } right)) of land cover. This function describes the expected effective number of species at location s, given sufficient time for the community to adapt to the given land cover composition. We now describe this function in more detail.The equilibrium component was formulated as a log-linear model comprising a total of I = 5 environmental covariates (the percentage cover of five landscape classes: urban, forest, grassland, wetland and cropland), using 2nd-order polynomial terms, captured by the coefficient j, to account for optima in effective number of species along each of the five environmental dimensions:$$fleft( {x_{s,t}} right) = {mathrm{exp}}left( {beta _0 + mathop {sum }limits_{i = 1}^{I = 5} mathop {sum }limits_{j = 1}^{J = 2} beta _{i,j}x_{i,s,t}^j} right)$$
(4)
In equation (4), the β parameters capture the effect of covariates on the equilibrium and are assumed to be the same for each environmental composition. A simplifying assumption necessary for the application of this model is that the effective number of species had equilibrated at the first timepoint. As data become available for more years in the future, the influence of this assumption on the model results will diminish and more accuracy will be achievable with multiple timepoints.To allow for conditionality in the effects of one land cover variable on the response of the effective number of species to another land cover variable, we extended this function with pairwise interaction terms k between all the linear terms for land cover variables and pairwise linear-quadratic terms, as follows:$$fleft( {x_{s,t}} right) = {mathrm{exp}}left( {beta _0 + mathop {sum }limits_{i = 1}^{I = 5} mathop {sum }limits_{j = 1}^{J = 2} beta _{0,i,j,}x_{i,s,t}^j + mathop {sum }limits_{i = 1}^{I = 4} mathop {sum }limits_{k = i + 1}^{K = 5} beta _{1,i,k}x_{i,s,t}x_{k,s,t}} right)$$
(5)
Temporal legacy componentThe main covariates, ({Delta}x_{i,s,t_1,t_2}), for the part of the model that captures the temporal legacy, (omega left( {{Delta}x_{s,t_1,t_2};gamma } right)), are derived from the change in land cover (({Delta}x_{i,s} = x_{i,s,t_2} – x_{i,s,t_1})) between the two timepoints$$x_{i,s,} = left{ {begin{array}{*{20}{l}} {x_{1,i,s} = left| {{Delta}x_{i,s}} right|,} hfill & {x_{2,i,s} = 0,} hfill & {{mathrm{if}},{Delta}x_{i,s} < 0} hfill \ {x_{1,i,s} = 0,} hfill & {x_{2,i,s} = {Delta}x_{i,s},} hfill & {{mathrm{otherwise}}} hfill end{array}} right.$$ (6) where ({Delta}x_{s,t,z}) is a vector of the ith environmental change variable (that is, urban, forest, grassland, wetland, cropland) at site s and for directionality z. The effect of these covariates on the mixture weights is given by:$$omega left( {{Delta}x_{s,t_1,t_2};gamma } right) = {mathrm{exp}}left( {mathop {sum }limits_{i = 1}^{I = 5} - gamma _{i,z}{Delta}x_{z,s,i}} right)$$ (7) This formulation weights the contribution that the environmental variables at the two timepoints have on the current effective number of species, as a function of the magnitude and directionality of change in each type of land cover covariate. The γ parameters, and subsequently the temporal legacy component, are allowed via the inclusion of the environmental change data ({Delta}x_{s,t,z}), to account for the distance between the land cover at the two timepoints, therefore quantifying how far the initial community would need to travel to reach equilibrium in 2016 as a function of the type, magnitude and directionality of change. It should be noted that our model, in equation (3), is only implicitly related to the speed with which the effective number of species reacts to environmental changes. Instead, it quantifies how much further it would still have to travel to reach the expected equilibrium associated with the current configuration of the landscape.Static covariatesAs described in model equation (3), we included a function of static covariates to which we can expect the effective number of species to respond without lags relating to the past landscape. We added a linear and quadratic fixed effect for temperature in 2016 to control for any difference in the effective number of species related to climatic characteristics and to allow for a parabolic relationship to be expressed (optima either at mean or extremes values). We also controlled for the heterogeneity of a landscape by including the effective number of land cover types, computed in the same way as the effective number of species, as a fixed effect40. A fixed effect for time of day, reflecting the time at which each segment was surveyed, was included to correct for differences in species detectability between early morning and later parts of the day41. An observer-level random effect was also added to control for variation between observers35,36 and partly account for between-route variation, given that we would expect observers who collect data from multiple routes to do so within a relatively small area. Spatial autocorrelation of the effective number of species was tested for all segments at once and by different radiuses for neighbour inclusion (500 m, 1,000 m, 5,000 m, 10,000 m, 100,000 m), using the Moran’s I statistic42. Spatial autocorrelation was not corrected for because Moran’s I was not significant at any spatial scale (P  > 0.05). Pseudo-replication between neighbouring segments was avoided by considering segments 1, 3 and 5, whose land cover buffers did not overlap (Extended Data Fig. 2).Model fittingThe model was fitted within a Bayesian framework using a Hamiltonian Markov chain Monte Carlo algorithm implemented in the STAN programming language43 version 4.3.0 and the ‘cmdstanr’ R package version 2.26.144.We ran 4 chains, sampling for 1,000 iterations with a burn-in period of 500 iterations each. These numbers of iterations were sufficient to achieve chain convergence. The STAN sampling was run on four parallel threads on a multi-core Intel i7 – 8750H processor with a maximum clock speed of 4.1 GHz.For the purposes of Bayesian inference, all slope parameters associated with the equilibrium component equation (5) and the static additive terms were assigned an unbiased prior (beta _{i,j} approx Nleft( {0,1} right),{mathrm{and}},z_s approx Nleft( {0,1} right),) where N is normal, with the aim of shrinking the parameter estimated towards 0 (that is, no covariate effect). A gamma distributed prior, with shape and rate 0.001, was assigned to the standard deviation of the random effect. For the following known and expected relationships, we also truncated the range of parameter values by bounding the upper or lower limits of the prior/posterior distributions. Intercept and standard deviation of the observer random effect were bounded below by 0. Linear effects for the environmental covariates and temperature were bounded below at 0, while their quadratic counterparts were bounded above at 0. Interaction terms were not limited. The temporal legacy component parameters were given a uniform (U) prior (gamma _isim Uleft( {0,1} right)), bounded between 0 and 1 to act as a weighting proportion between the present and the past. The upper bound on the gamma parameters to 1 does not bias us towards an increased contribution of the past land cover, but instead provides a more conservative approach.Model diagnostics were conducted by assessing chain convergence visually through trace plots, as well as statistically by employing the Gelman-Rubin test, which compares the estimated between-chain and within-chain variances45. Chain autocorrelation and the associated effective sample size were also monitored. In the case of low effective sample size, the chains were extended until the effective sample size exceeded a threshold value of 400. The marginal posterior distribution for each parameter was visualized via a density plot to check for multimodality.Model selection was conducted to inform choice of the size and shape of the land cover buffer around each sampled segment. We did so by comparing values of the Watanabe-Akaike Information Criterion leave-one-out (WAIC)-loo information criterion46 of four different models, each computed using land cover data calculated with two different buffer options of various sizes: a circular buffer around the centroid of the polygon defined by the vertices of each segment (4,000 m radius) and a series of buffers around the segment line (500 m, 2,000 m and 4,000 m radius). This approach was implemented through the ‘loo’ R package version 2.1, which provides an improvement on the original WAIC by including diagnostic measures around the point-wise log-likelihood value estimated around each sample draw47.Visualization of model predictionsA map of the USA (Fig. 1) was produced to represent the predicted extinction debts and colonization credits (that is, positive or negative distance in the effective number of species from the expected equilibria). The map was produced on a hexagonal grid at a spatial resolution of 10 km vertex-to-opposite-vertex, with each hexagon covering a total of 86 km2. Values of extinction debt and colonization credit were calculated by subtracting the predicted effective number of species produced by the model (equation 3) from the predicted effective number of species at equilibrium in 2016 (that is, when the legacy component equals 1). To correctly propagate and represent uncertainty in the extinction debts and colonization credits presented, this process was repeated 1,000 times for predictions originating from different draws from the posterior distribution. Uncertainty in the form of the geometric coefficient of variation, calculated as (root {2} of {{e^{left( {mathrm{log}left( {sigma + 1} right)^2} right)} – 1}}) where σ is the standard deviation, is mapped in Extended Data Fig. 4a. Extended Data Fig. 4 also includes a copy of Fig. 1 (Extended Data Fig. 4b) for reference, alongside upper (Extended Data Fig. 4c) and lower (Extended Data Fig. 4d) credible intervals.Over/underestimation values of biodiversity that could arise by neglecting debts and credits were computed as the difference between the effective numbers of species predicted by the equilibrium model and the legacy model, multiplied by 100 and then divided by the predicted effective number of species under the legacy model. This calculation results in a percentage measurement of the extent to which (in relative terms) the current effective number of species under- or overestimates the diversity that a given landscape can sustain at equilibrium.To further validate our predicted extinction debts and colonization credits, we compared the direction of the expected changes with the recorded difference in effective numbers of species between 2016 and 2019 (the latest year for which data are available). To do so, we sourced bird abundances from the North American BBS dataset14,32 for the year 2019 and conducted the same data processing as described above for the other two timepoints. We then conducted a Pearson correlation test to assess how well the observed change followed the model-predicted one. We are nevertheless aware that a 3-year timeframe is unlikely to be large enough for debts and credits to fully manifest.Plots were also generated to describe the behaviour of the mixture weight, ω (equation 7), which captures the contribution (weighting) of the landscape composition in determining the effective numbers of species at the two timepoints (Fig. 2 in the main text). Values of ω across the whole spectrum of plausible land cover change values (that is, from −100 to +100) were simulated by averaging over 10,000 draws from the posterior distribution of each γ parameter. Credible intervals were measured by taking the 95% range of the 10,000 draws.Explaining spatial variation in debts and creditsThe extinction debts and colonization credits predicted for the contiguous USA states were further modelled to identify which past land cover changes were the main drivers of the delayed biodiversity changes in USA bird communities. We considered the values of debts or credits associated with the 92,000 individual 86 km2 hexagons (Fig. 1) as a response variable. We then specified a Gaussian linear model including the magnitude of each land cover change as explanatory covariates. Positive and negative changes in each covariate were treated as separate linear components to differentiate their effects. The model was fitted to 1,000 sets of debts and credits, each originating from predictions based on independent draws from the posterior distribution. For each generalized linear model (GLM) fit, we then subsequently sampled each parameter distribution another 1,000 times and extracted the summarized parameter estimates from a total of 100,000 values. Model coefficients and their resulting uncertainty from the above process are presented in Fig. 4 and in more detail as part of Supplementary Table 3.Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Pollen-inferred landscape change and pre-industrial demographyRecently, data derived from tree rings or ice cores have been employed to approximate changes in human economic activity related to past epidemics, as well as to warfare and climatic variability46,47. However, none of these proxies is directly related to human demography or provides a basis to estimate variation in the Black Death’s mortality on a regional scale across Europe (to date only a single archaeological study using pottery as a proxy for demographic change on the national level, focusing on just a single country—England—has appeared48).In recent years, pollen data have been proven to be closely related to demographic variability. Most importantly, detailed comparisons of historical documentary data on population trends and landscape changes as revealed by pollen data have been carried out on a local scale and a close link between changes in European pollen data and changes in European local demography over the past millennium has been demonstrated on multiple occasions, that is, during the period and region of our concern here49,50. A strong link between long-term demographic trends as visible in regional settlement numbers and macro-changes in land cover (deforestation/afforestation) have also been confirmed for ancient Greece51. Additionally, a recent publication successfully employed pollen data to test the extent of the mortality associated with the sixteenth-century Spanish and Portuguese empires’ colonization of tropical regions in the Americas and Asia52. However, as of now there is no method to quantify past demographic trends in absolute numbers based on palaeoecological data. Consequently, we also focus in this paper on relative changes in historical societies’ populations and test the now common idea that the Black Death caused enormous mortality across Europe (with many scholars now arguing for a mortality exceeding 30% and upwards of 50% of the population within a few years) (see also Fig. 1). Using our BDP approach, we conclude this hypothesis is not maintainable. Our evidence for demography-related landscape changes (or lack thereof) negates it.Our main indicator is cereal pollen. In pre-industrial economies, rural labour availability (hence rural population levels) and the spatial scale of cereal cultivation were directly related. An increase in the extent and intensity of cereal cultivation—as reflected in pollen data—would have required not only a predilection and demand for cereals, but also greater availability of labour and thus population growth or significant immigration. The maintenance of existing agricultural activity, in turn, would have required relatively stable population levels53,54,55. The uniform ~50% mortality postulated for the Black Death across Europe should have resulted in a large and significant decline of cereal cultivation and parallel forest regrowth across Europe, as previously demonstrated for mid-fourteenth-century Sweden26 and singular sites in some regions of western Europe56. This result agrees with the fact that Black Death mortality could be high among people at productive age, as illustrated for England57,58. Moreover, even in the case of England, a comparatively commercialized and adaptive rural economy in mid-fourteenth-century Europe, the loss of 50% of the population led to a significant decline in the total area under cultivation (as documented by heterogeneous written sources)59. In Italy, another well-developed economy at that time, the expansion of large estates following the Black Death also did not compensate for the general loss of cereal productivity60. This effect, high mortality driving arable contraction, must have been yet more pronounced in more subsistence-oriented and less adaptive economies, with limited surplus production, such as in regions of the Iberian Peninsula, Germany, Sweden and particularly east-central Europe. Importantly, palaeoecological evidence for arable contraction may be indicative, to some extent, of not only rural population decline but also urban population decline in the region, as there is evidence in some areas, following the pandemic, of rural-to-urban migration, of country-dwellers repopulating urban centres10. Possibly less common was intraregional rural migration, as marginal lands were abandoned for better quality soils, which were more likely to remain under cultivation26,61.Therefore, cereal pollen remains our most potent pollen indicator related to demographic changes in pre-industrial European societies. Other pollen indicators, reflecting rewilding and reforestation (secondary ecological succession) of cereal fields abandoned as a result of significant mortality, or the transformation of cereal fields into pastures, which required less rural labour and thus also could have been a response to high plague mortality, play a secondary role in our analysis and provide further support for our conclusions.BDP data collectionExisting online palynological databases (the European Pollen Database (EPD)62 www.europeanpollendatabase.net, and the Czech Quaternary Palynological Database (PALYCZ)63, https://botany.natur.cuni.cz/palycz/), as well as personal contacts of the study authors and a systematic publication search were employed to identify palynological sites in Europe reaching the required chronological and resolution quality for the study of the last millennium. In order to enable statistical analysis, we included only sites clustered in well-defined historical-geographical regions, excluding isolated sites even if the quality of a site’s data was very good. Data of sufficient quality and amount from regions for which the Black Death is well-studied, notably central and northern England and the Low Countries, is not presently available; to the best of our knowledge, for each of these regions there currently is not more than a single isolated site56, which does not allow for the application of statistical approaches.In total, 261 pollen records with the average temporal resolution of 58 years and 14C-age control (or varve chronology), have been collected. The age–depth models of the sequences have been provided by authors in original publications, by the EPD or developed through the Clam package (version 2.3.4) of R software for the purpose of this study. The analytical protocol for pollen extraction and identification is reported in the original publications. The Pollen Sum includes all the terrestrial taxa with some exceptions based on the selection done in the original publications. The full list of sequences, exclusions from the Pollen Sum, age-depth models and full references are reported in Supplementary Data 1.The taxa list has been normalized by applying the EPD nomenclature. In this respect, the general name Cichorioideae includes Asteraceae subf. Cichorioideae of the EPD and PALYCZ nomenclatures, which primarily refers to the fenestrate pollen of the Cichorieae tribe64. Ericaceae groups Arbutus unedo, Calluna vulgaris, Vaccinium and different Erica pollen types, whereas deciduous Quercus comprehends both Q. robur and Q. cerris pollen types65. Rosaceae refers to both tree and herb species of the family. Finally, Rumex includes R. acetosa type, R. acetosella, R. crispus type, Rumex/Oxyria and Urtica groups U. dioica type and U. pilulifera.BDP summary pollen indicatorsIn order to connect changes visible in the pollen data to human demographic trajectories, we assembled four summary pollen indicators that describe specific landscapes related to human activity. They reflect different degrees of demographic pressure on the landscape (cereal cultivation, pastoral activities, which are less-labour intensive than cereal cultivation, abandonment and rewilding) as well as different durations of land abandonment that might have occurred post-Black Death. Our indicators account for the fact that Europe is a continent rich in natural heritage, with a wide range of landscapes and habitats and a remarkable wealth of flora and fauna, shaped by climate, geomorphology and human activity. In order to ensure uniform interpretation of the indicators, we relied on criteria that can be applied to all European landscapes regardless of their local specificity. Cereals and herding are directly related to human activities and are barely influenced by spatial differences. More complex is the succession of natural plants with their ecological behaviour and inter-species competition. For this reason, we relied on existing quantitative indicators of plant ecology.The Ellenberg L – light availability indicator66 provides a measure of sunlight availability in woodlands and consequently of tree-canopy thickness, reflecting the scale of the natural regeneration of woodland vegetation after cultivation or pasture activities61. Nonetheless, ecological studies have suggested that geographic and climatic variability between different European regions can influence the Ellenberg indicator system67,68,69,70,71. The original indicators were primarily designed for Central Europe58, but several studies developed Ellenberg indicators for other regions, reflecting the specific ecology of the selected taxa (British Isles;72 Czech Republic;73 Greece;74 Italy;75 Sweden76). Plants with L values between 5 and 8 are listed in the fast succession indicator, the ones with L values ranging from 1 to 4 are included in the slow succession indicator. The result is the following list:1) Cereals: only cultivated cereals have been included: Avena/Triticum type, Cerealia type, Hordeum type, Secale. 2) Herding includes pastoral indicators linked to the redistribution of human pressure: Artemisia, Cichorioideae, Plantago lanceolata type, Plantago major/media type, Polygonum aviculare type, Rumex, Trifolium type, Urtica, Vicia type. 3) Fast Succession comprises indicators of relatively recent reforestation of cultivated land after abandonment: Alnus, Betula, Corylus, Ericaceae, Fraxinus ornus, Juniperus, Picea, Pinus, Populus, deciduous Quercus, Rosaceae. 4) Slow Succession includes indicators of secondary succession established after several decades of abandonment: Abies, Carpinus betulus, Fagus, Fraxinus, Ostrya/Carpinus orientalis, Quercus ilex type.In order to validate the indicators overcoming the regional limits of Ellenberg values, a different subdivision has been provided following the Niinemets and Valladares shade tolerance scale for woody species of the Northern Hemisphere77. The subdivision of taxa in the Fast and Slow succession indicators remains the same with only three changes: Fraxinus ornus and Picea move from Fast to Slow succession and Fraxinus from Slow to Fast succession. Extended Data Figs. 1 and 2 show that the two groupings yield the same results, which confirms the reliability of our indicators. There is only one clear exception (Russia), with one more region where smaller-scale diversion occurs for only one indicator, Slow Succession (Norway). The different indicator behaviour results from the different attribution of Picea in our two sets of succession indicators: at high latitude, Picea characterizes the final stage of the ecological succession and hence its different attribution results in different summary indicator values in Russia for the two stages of ecological succession, fast and slow.Please note our summary indicators are not designed to reflect the entirety of the landscape and reconstruct all of its different components. Rather, they are a means of approximating changes in the landscape related to the types of human activities, and their intensity, as much as they relate to demographic changes in human populations using and inhabiting these landscapes.BDP analytical statistical and spatial methodsTo control for local specificity, pollen percentages of every taxon from each pollen site were standardized. From the taxa percentage in a given year the arithmetic mean calculated for the observations from the period 1250–1450 was subtracted and the result divided by the standard deviation for the 1250–1450 period. Standardized taxa results were assembled for each site into four BDP summary indicators. Since each indicator has different numbers of taxa, the sum of standardized taxa values calculated for a given year and site was divided by the number of taxa in the indicator. For the purposes of replication, this standardized pollen dataset, comprising the four indicators for each sample from each site, is available as Supplementary Data 2.This dataset has been analysed in two ways, statistically and spatially.For the statistical approach, standardized regional indices of landscape transformation were created for each region by calculating the average value for all sites within the region, for each of the subperiods analysed in the study (1250–1350 and 1351–1450; 1301–1350 and 1351–1400; 1325–1350 and 1351–1375). Differences between means for each subperiod were measured by the use of the bootstrapping based on 10,000 resamples. The 90% and 95% confidence intervals were estimated with the bias-corrected and accelerated method (BCa)78. These results are visualized in Fig. 5 for the comparison of the subperiods of 1250–1350 versus 1351–1450, and in Supplementary Figs. 4 and 6 for the comparison of the subperiods of 1300–1350 versus 1351–1400 and 1325–1350 versus 1351–1375, respectively.For the spatial approach, we employed the Bayesian model AverageR developed within the Pandora and IsoMemo initiatives (https://pandoraapp.earth/) to map the distribution of pollen indices across Europe. AverageR is a generalized additive model that has been described previously79. It relies on a thin-plate regression spline80 to predict new, unseen data using the following model:$$Y_{i} = {g}( {{{{mathrm{longitude}}}},{{{mathrm{latitude}}}}} ) + {varepsilon}_{i}$$Where Yi is the independent variable for site i; g(longitude, latitude) is the spline smoother; and εi ∼ N(0, σε) is the error term.The spline smoother can be written as X × β where X is a fixed design matrix and β is the parameter vector. Surface smoothing is controlled by employing a Bayesian smoothing parameter estimated from the data and trades-off bias against variance to make the optimal prediction75. This parameter β is assumed to follow a normal distribution: β ~ N(0, 1 /δ × λ × P), where P is a so-called penalty matrix of the thin plate regression spline, which penalizes second derivatives81. The δ parameter is by default set to 1 but this can be adjusted to suit smoothing needs for each application. In our study δ was set at 0.9 to match the preferred spatial scale of analysis for our dataset (approx. 250 to 500 km).AverageR was employed to generate smoothed surfaces for three sets of temporal bins (1250–1350 versus 1351–1450, as well as 1300–1350 versus 1351–1400 and 1325–1350 versus 1351–1375) and for the four BDP indicators (Supplementary Figs. 3, 5 and 7). For the same indicator the difference between the two temporal bins was plotted (Fig. 5; Supplementary Figs. 4 and 6).Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Erwin, D. H. Curr. Biol. 25, R930–R940 (2015).CAS 
Article 

Google Scholar 
Rabosky, D. L. & Hurlbert, A. H. Am. Nat. 185, 572–583 (2015).Article 

Google Scholar 
Hughes, M., Gerber, S. & Wills, M. A. Proc. Natl Acad. Sci. USA 110, 13875–13879 (2013).CAS 
Article 

Google Scholar 
Novack-Gottshall, P. M. et al. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-021-01656-0 (2022).Article 

Google Scholar 
Wright, D. F. Sci. Rep. 7, 13845 (2015).
Google Scholar 
Deline, B. in Elements of Paleontology hcxz (Cambridge Univ. Press, 2021).Ausich, W. I., Kammer, T. W., Rhenberg, E. C. & Wright, D. F. Palaeontology 58, 937–952 (2015).Article 

Google Scholar 
Zamora, S. & Rahman, I. A. Palaeontology 57, 1105–1119 (2015).Article 

Google Scholar 
Sheffield, S. L. & Sumrall, C. D. Palaeontology 62, 163–173 (2018).Article 

Google Scholar 
Deline, B. et al. Curr. Biol. 30, 1672–1679 (2020).CAS 
Article 

Google Scholar 
Novack-Gottshall, P. M. Paleobiology 33, 273–294 (2007).Article 

Google Scholar 
Anderson, P. S. Paleobiology. 35, 321–342 (2009).Article 

Google Scholar 
Cole, S. R. & Hopkins, M. J. Sci. Adv. 7, eabf4027 (2021).Article 

Google Scholar 
Hopkins, M. J. & St. John, K. Syst. Biol. 70, 1163–1180 (2021).Article 

Google Scholar 
Blomberg, S. P., Garland, T. Jr & Ives, A. R. Evolution 57, 717–745 (2003).Article 

Google Scholar 
Wright, A., Wagner, P. J., & Wright, D. F. in Elements of Paleontology hcx2 (Cambridge Univ. Press, 2021).Cole, S. R., Wright, D. F. & Ausich, W. I. Palaeogeogr. Palaeoclimatol. Palaeoecol. 521, 82–98 (2019).Article 

Google Scholar  More

Isolation and cultivationIn total, 39 mixotrophic flagellates were investigated (Table 1). Thirty-three isolates were enriched and isolated from euphotic zone samples at Station ALOHA (22° 45′ N, 158° 00′ W) in February and May 2019. To select for mixotrophic grazers of Prochlorococcus, whole seawater was amended with K medium [28] (1/20 final concentration), and live Procholorococcus (MIT9301) was added as prey (~5 × 106 cells mL−1 final concentration). Enriched seawater samples were incubated under ~70 µmol photons m−2 s−1 irradiance on a 12 h:12 h light:dark cycle and monitored by microscopy daily up to five days. Each day, samples were serially diluted to extinction (9 dilution steps, 12 replicates per dilution) in 96-well plates in nutrient-reduced K medium (1/20 concentration) with a constant background of Prochlorococcus cells. Wells at the highest dilution showing growth of putative grazers were subjected to 3–6 further rounds of dilution to extinction. Four additional mixotrophs were isolated in full K medium using water from earlier cruises, and two (dictyochophyte strains UHM3021 described in [24] and UHM3050) were enriched in K minus nitrogen medium (K-N) without Prochlorococcus enrichment. All isolates were rendered unialgal, but not axenic, and maintained at 24 °C in K-N medium (~0.2 µM N) amended with Prochlorococcus prey, under the same light conditions as above. Dense Prochlorococcus cells grown in Pro99 medium [29], were harvested and concentrated through gentle centrifugation at 2000 RCF for 5 minutes and resuspended in fresh K-N medium to minimize nutrient carryover. To ensure their long-term accessibility, the isolates used in this study are being transferred to the National Center for Marine Algae and Microbiota at the Bigelow Laboratory for Ocean Science, East Boothbay, ME, USA (ncma.bigelow.org).Table 1 Mixotrophic flagellates investigated in this study.Full size table18S rRNA gene sequencing and phylogenetic analysisCells were harvested by centrifuging 25–50 mL dense cultures at 3000 RCF for 10 min at 4 °C. Genomic DNA was extracted from the pellets using the ZymoBIOMICS DNA Kit (Zymo Research, Irvine, CA, USA). A near-full-length section of the eukaryotic small-subunit ribosomal RNA (18S rRNA) gene was amplified by PCR with the Roche Expand High Fidelity PCR System (Sigma-Aldrich, St. Louis, MO, USA) using either forward primer 5′-ACCTGGTTGATCCTGCCAG-3′ and reverse primer 5′-TGATCCTTCYGCAGGTTCAC-3′ [30], or Euk63F 5′-ACGCTTGTCTCAAAGATTA-3 and Euk1818R 5′-ACGGAAACCTTGTTACGA-3′ [31]. Amplicons were purified using spin columns (DNA Clean & Concentrator-25; Zymo Research, Irvine, CA, USA) and sequenced (Sanger) using the same PCR amplification primers and an additional reverse primer 1174R, 5′-CCCGTGTTGAGTCAAA-3′ [32], when necessary to connect two ends. For phylogenetic analyses, similar sequences were retrieved from the PR2 database [33] based on BLAST similarity, and two environmental homologs (GenBank Acc. FJ537342 and FJ537336) were retrieved from NCBI GenBank for the undescribed haptophyte taxon, which was not affiliated with any reference sequence from the PR2 database. Sequence alignments including 39 isolates, 29 reference and 2 outgroup taxa were created with MAFFT v7.450 using the G-INS-i algorithm [34] in Geneious R11.1.5 (http://www.geneious.com) [35]. Terminal sites that lacked data for any of the sequences were trimmed and any sites with greater than 25% gaps were removed from the alignment, which generated a total sequence length of 1617 bases. Phylogenetic analysis was performed using MrBayes v3.2.6 in Geneious R11.1.5 [36] with two runs of four chains for 1,000,000 generations, subsampling every 200 generations with burn-in length 100,000, under the GTR substitution model. The Bayesian majority consensus tree was further edited within iTOL v5 [37]. All 18 S rRNA gene sequences were deposited in GenBank with accession numbers MZ611704–MZ611740; MN615710–MN615711.Microscopic observationThe average diameter of flagellates in the exponential growth phase (n = 20 cells per strain) was measured by transmitted light microscopy using image analysis software (NIS-Elements AR, Nikon, Minato City, Tokyo, Japan) calibrated with a stage micrometer. Equivalent spherical diameter (ESD) and biovolumes were calculated assuming spherical cells. Chloroplasts were visualized by autofluorescence under epifluorescence microscopy. An average ESD of 0.64 µm was used for Prochlorococcus prey [24].Visual evidence of phagocytosis was obtained by adding fluorescent beads (0.5 μm YG Fluoresbrite Microspheres; Polysciences) to each culture. Samples post incubation (~2 h) were fixed with an equal volume of 4% ice-cold glutaraldehyde, and subsamples (20 µL) were mounted on a glass slide under a coverslip. Paired images captured using epifluorescence and transmitted light microscopy (Olympus BX51 with Leica DFC 7000 T color digital camera) were overlain to identify cells with ingested beads.Grazing experimentsLong-term grazing experiments were conducted for all 39 grazers, and 31 were used to quantify grazing rates based on rates of disappearance of Prochlorococcus cells, which persist but do not readily grow in K-N medium [24]. Rates were not calculated for eight isolates (seven Florenciella and one DictyX) because they were sampled at a lower frequency. Fifteen isolates representing all genera (or approximately genus-level clades) were examined in more detail by replicating grazing experiments two times (marked in bold in Supplementary Table S1), while the remaining sixteen isolates were tested once to survey within- and across- genus variation. Prior to the experiments, all grazer cultures were maintained/acclimated in the experimental medium (K-N with prey). Experiments were initialized by inoculating late-exponential-phase grazers into fresh K-N medium at a final concentration of ~103 flagellates mL−1, and adding live, unstained prey at a final concentration of 2–3 × 106 Prochlorococcus mL−1. Grazers were incubated for 3–8 days in total, depending on how fast prey were ingested. To minimize carryover of dissolved nutrients and prey growth in the grazing experiments, Prochlorococcus were grown to stationary phase in Pro99 medium, then pelleted (2000 RCF for 3–5 min) and resuspended in fresh K-N medium prior to addition to control and experimental cultures. Control cultures of grazer without added prey and prey without grazers were included during each grazing experiment to confirm that grazer growth and prey removal were attributable to grazing. Cell concentrations of prey and grazers were measured every 12–24 h by flow cytometry of glutaraldehyde-fixed samples at final concentration of 0.5% (Attune NxT; Thermo Fisher Scientific, Waltham, MA, USA, and CytoFLEX; Beckman Coulter Life Sciences, Indianapolis, IN, USA). Populations were distinguished based on light scatter and pigment autofluorescence and occasionally confirmed with DNA fluorescence (stained post-sampling with DNA stain SYBR Green I). Ambient bacteria concentrations, monitored using DNA stains, were ≤10% of the added Prochlorococcus at the start and during most periods of all grazing experiments.Grazing rates and biovolume conversion efficiencyWe calculated ingestion rates, I (prey grazer−1 h−1) for each grazer as:$$I = frac{{P_t – P_{t + 1}}}{{G_{{{{{rm{avg}}}}}}(T_{t + 1} – T_t)}}$$
(1)
where (P_t) and (P_{t + 1}) are the prey abundance at sampling interval t and t + 1 (cells mL−1), Gavg is the arithmetic mean grazer abundance (cells mL−1) over the time interval, and ((T_{t + 1} – T_t)) is the time (h) between two sampling intervals. Clearance rates (C, nL grazer−1 h−1) were calculated by dividing the ingestion rate by average prey concentration over the same interval, and specific clearance rate (body volume grazer−1 h−1) was calculated by dividing the clearance rate by cellular biovolume (µm3) of each grazer. Equation (1) uses a linear approximation of prey and grazer trajectories over the sampling interval, which was appropriate for our data where change could be relatively linear, concave-up, or concave-down. Other commonly used ingestion rate calculations assume exponential prey decline [38] and/or exponential grazer growth [39].Clearance rates over time in each experiment were assessed visually to obtain a representative series of rates that minimized potential influence of modest prey growth/decline observed in control cultures (Supplementary Fig. S1), as well as potential slowing of ingestion as the grazer neared carrying capacity or depleted prey to a low concentration. For each experiment a contiguous set of relatively constant rates were used to calculate a mean clearance rate. This assessment sometimes excluded the first 12–24 h, but not when removal rates were particularly fast. Intervals when the grazer neared carrying capacity were also often excluded, if grazing rates slowed down. To assess whether clearance rates increased as prey were depleted we plotted clearance rate and ingestion rate as a function of prey concentration (Supplementary Figs. S2 and S3). In general, there was no relationship between clearance rate and prey concentration, and ingestion increased linearly with prey concentration. These patterns imply that the prey concentrations in this experiment did not saturate the ingestion rates of these grazers. Under non-saturating prey concentrations the average clearance rate over an experiment should be a good estimate of the maximum clearance rate (Cmax). Consistent with this interpretation, functional responses fit to these experiments yielded Cmax estimates that were similar to the reported average clearance rates. Because these experiments were not designed with a sufficient range of prey density to esimate functional responses, we do not report the Cmax estimates, but we note here that our reported average clearance rates may be useful as approximate Cmax numbers in future work.Six grazers representing three classes (three dictyochophytes, two haptophytes, and one chrysophyte) were further investigated to determine functional grazing responses using a wide range of initial prey densities (105–107 cells mL−1). Functional responses were modeled using the Holling type II curve, (I = frac{{I_{{{{{rm{max}}}}}}P}}{{P + frac{{I_{{{{{rm{max}}}}}}}}{{C_{{{{{rm{max}}}}}}}}}}), where I is the ingestion rate over a sampling interval (Eq. 1), Imax is the maximum ingestion rate, Cmax is the maximum clearance rate and P is the arithmetic mean prey density between two sampling points. This curve was fit to ingestion rate data using maximum likelihood with R package bbmle [40].For 31 isolates we calculated the amount of grazer biovolume created per prey biovolume consumed, using data from the same grazing experiments used to calculate grazing rates. It was calculated based on the following formula:$$E = frac{{(F_{{{{{rm{f}}}}}} – F_{{{{{rm{i}}}}}})}}{{(P_{{{{{rm{i}}}}}} – P_{{{{{rm{f}}}}}})}}frac{{B_{{{{{rm{F}}}}}}}}{{B_{{{{{rm{p}}}}}}}}$$
(2)
where Ff and Fi are the final and initial flagellate concentrations, Pi and Pf are initial and final prey concentrations in each culture, and BF and BP are the cellular biovolume of prey and grazer. We refer to the quantity E as the biovolume conversion efficiency, and we use it as an indicator of physiological differences among diverse mixotrophs. Note that biovolume conversion efficiency can be greater than 1, if prey have greater nutrient:biovolume than the grazer.Quantitative PCRReal-time, quantitative PCR (qPCR) was performed to quantify the 18S rRNA gene abundances of representative mixotroph groups discriminated at approximately the genus level, including Florenciella, Rhizochromulina and another undescribed clade within the class Dictyochophyceae; Chrysochromulina and another undescribed clade within the division Haptophyta; clade H in the class Chrysophyceae; and Triparma eleuthera and Triparma mediterranea in the class Bolidophyceae. Primers (Supplementary Table S2) were designed to target a short region (95–176 bases) of the 18S rRNA gene and meet basic criteria (≤2 °C difference in melting temperature between members of a pair, %G + C content between 45 and 65%, ≤1 degenerate position per primer, no predicted primer dimers). Sequences considered targets for a given primer set had ≤1 mismatch across both primers, which included all or most known members within the corresponding targeted clade. Members in the nearest non-targeted clade had ≥3 mismatches distributed across both primers. Efficiency and specificity of the synthesized primers (IDT Inc., Coralville, IA, USA) was tested by ensuring there was specific amplification (qPCR followed by melting curve analysis and gel electrophoresis) when using DNA from cultures within the targeted group and no amplification when using DNA from cultures close to, but outside of the targeted group (Supplementary Table S3). Empirical observations of amplification success using control cultures were used to infer whether species known only by environmental sequences were likely to amplify with a given primer set (Supplementary Fig. S4).In situ gene abundances were quantified in water samples collected from Station ALOHA at 5, 25, 45, 75, 100, 125, 150, and 175 m, during HOT cruise numbers 259 (Jan), 262 (Apr), 264 (Jul), and 266 (Oct) of 2014. Seawater (ca. 2 L) was filtered through 0.02 μm pore-size, aluminum oxide filters (Whatman Anotop, Sigma-Aldrich, Saint Louis, MO, USA) and stored at −80 °C. Genomic DNA of both grazer cultures and environmental samples was extracted (MasterPure Complete DNA and RNA Purification Kit; Epicentre) as described elsewhere [41]. Four replicated PCR reactions (10 μL) were carried out for each sample except for Triparma (duplicates) and consisted of 5 μl of 2× PowerTrack SYBR Green Master Mix (Thermo Fisher Scientific, USA), 10 ng environmental DNA, 500 nM of each primer, and nuclease-free water. Reactions were run on an Eppendorf Mastercycler epgradient S realplex2 real-time PCR instrument. Each run contained fresh serial dilutions (1–6 log gene copies) of target-specific, 750-bp synthetic standards (gBlocks, IDT) prepared in triplicate. The cycling program included an initial denaturation step of 95 °C for 2 min, followed by 40 cycles of 95 °C for 5 s and 55 °C for 30 s. Specificity of amplification was checked with a melting curve run immediately after the PCR program and occasionally, by gel electrophoresis. Amplification efficiencies ranged from 95% to 106% for all the primers.To convert gene copies to cell numbers, 18S rRNA gene copy number per cell−1 was determined for representative isolates in the seven targeted genera/clades. Known quantities of cultured cells (106–107 cells) from each isolate with 2–8 replicates were pelleted at 4000 RCF for 15 min at 4 °C. DNA was extracted from the pelleted cells (MasterPure Complete DNA and RNA Purification Kit, Lucigen), quantified by fluorometry (Qubit, Invitrogen) and the extract volume adjusted to achieve a DNA concentration of 10 ng µL−1. The expected number of eukaryotic cells µL−1 of extract was calculated as the difference between the total cells in the sample prior to centrifugation and in the supernatant afterward (as determined by flow cytometry) divided by the final extract volume. Copy number of the 18S rRNA genes µL−1 of extract was determined by qPCR with the appropriate group-specific primers. The resulting value of gene copies µL−1 was divided by the equivalent number of eukaryotic cells µL−1 in the extract (assuming 100% extraction efficiency) to derive minimum estimates of gene copies cell−1. An average value for representatives within each genus/clade (1–5 isolates) was used to calculate in situ cell concentrations for the genus. These derived in situ abundances were compared to flow cytometric counts of total photosynthetic picoeukaryotes at Station ALOHA obtained from the Hawai’i Ocean Time-series Data Organization and Graphical System (https://hahana.soest.hawaii.edu/hot/hot-dogs/).Global distribution revealed through Tara Oceans 18S rRNA metabarcodesTo estimate the relative abundance of the OTUs closely related to our diverse isolates on a broader geographic scale, we searched the 18S rRNA-V9 sequence data from the 0.8–5 µm fraction of surface water sampled at 40 stations by the Tara Oceans project (http://taraoceans.sb-roscoff.fr/EukDiv/). Reads for ‘Tara lineages’ with highest similarity (E-value < 10−15) to each of our targeted clades (Supplementary Table S1) were expressed as a fraction of total reads excluding dinoflagellates but included all other Tara Oceans phytoplankton ‘taxogroups’: Bacillariophyta, Bolidophyceae, Chlorarachnea, Chlorophyceae, Chrysophyceae/Synurophyceae, Cryptophyta, Dictyochophyceae, Euglenida, Glaucocystophyta, Haptophyta, Mamiellophyceae, Other Archaeplastida, Other Chlorophyta, Pelagophyceae, Phaeophyceae, Pinguiophyceae, Prasino-Clade-7, Pyramimonadales, Raphidophyceae, Rhodophyta and Trebouxiophyceae. Dinoflagellates were excluded because of the difficulty in assigning phototrophic vs. heterotrophic status to all taxa, and because nearly all dinoflagellate reads were from a single, poorly annotated OTU that was also highly abundant in larger size fractions. More

Simpson, G. G. Tempo and Mode in Evolution (Columbia Univ. Press, 1944).Losos, J. B. Adaptive radiation, ecological opportunity, and evolutionary determinism. Am. Nat. 175, 623–639 (2010).PubMed 
PubMed Central 

Google Scholar 
Erwin, D. H. Novelty and innovation in the history of life. Curr. Biol. 25, 930–940 (2015).
Google Scholar 
Novack-Gottshall, P. M. General models of ecological diversification. I. Conceptual synthesis. Paleobiology 42, 185–208 (2016).
Google Scholar 
Marshall, C. R. Explaining the Cambrian “explosion” of animals. Annu. Rev. Earth Planet. Sci. 34, 355–384 (2006).CAS 

Google Scholar 
Hughes, M., Gerber, S. & Wills, M. A. Clades reach highest morphological disparity early in their evolution. Proc. Natl Acad. Sci. USA 110, 13875–13879 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Hopkins, M. J. & Smith, A. B. Dynamic evolutionary change in post-Paleozoic echinoids and the importance of scale when interpreting changes in rates of evolution. Proc. Natl Acad. Sci. USA 112, 3758–3763 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Wright, D. F. Phenotypic innovation and adaptive constraints in the evolutionary radiation of Palaeozoic crinoids. Sci. Rep. 7, 13745 (2017).PubMed 
PubMed Central 

Google Scholar 
Cantalapiedra, J. L., Prado, J. L., Hernández Fernández, M. & Alberdi, M. T. Decoupled ecomorphological evolution and diversification in Neogene-Quaternary horses. Science 355, 627–630 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Crouch, N. M. A. & Ricklefs, R. E. Speciation rate is independent of the rate of evolution of morphological size, shape, and absolute morphological specialization in a large clade of birds. Am. Naturalist 193, E78–E91 (2019).
Google Scholar 
Slater, G. J. & Friscia, A. R. Hierarchy in adaptive radiation: a case study using the Carnivora (Mammalia). Evolution 73, 524–539 (2019).PubMed 
PubMed Central 

Google Scholar 
Ronco, F. et al. Drivers and dynamics of a massive adaptive radiation in cichlid fishes. Nature https://doi.org/10.1038/s41586-020-2930-4 (2020).Cole, S. R. & Hopkins, M. J. Selectivity and the effect of mass extinctions on disparity and functional ecology. Sci. Adv. https://doi.org/10.1126/sciadv.abf4072 (2021).Erwin, D. H. et al. The Cambrian conundrum: early divergence and later ecological success in the early history of animals. Science 334, 1091–1097 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
Smith, A. B., Zamora, S. & Álvaro, J. J. The oldest echinoderm faunas from Gondwana show that echinoderm body plan diversification was rapid. Nat. Commun. 4, 1385 (2013).PubMed 
PubMed Central 

Google Scholar 
Deline, B. et al. Evolution and development at the origin of a phylum. Curr. Biol. 30, 1672–1679 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Pisani, D., Feuda, R., Peterson, K. J. & Smith, A. B. Resolving phylogenetic signal from noise when divergence is rapid: a new look at the old problem of echinoderm class relationships. Mol. Phylogenet. Evol. 62, 27–34 (2012).PubMed 
PubMed Central 

Google Scholar 
Smith, A. B. & Zamora, S. Cambrian spiral-plated echinoderms from Gondwana reveal the earliest pentaradial body plan. Proc. R. Soc. B 280, 20131197 (2013).PubMed 
PubMed Central 

Google Scholar 
Zamora, S. et al. in Early Palaeozoic Biogeography and Palaeogeography. Memoirs 38 (eds Harper, D. A. T. & Servais, T.) 157–171 (Geological Society, London, 2013).Lefebvre, B. et al. in Early Palaeozoic Biogeography and Palaeogeography. Memoirs 38 (eds Harper, D. A. T. & Servais, T.) 173–198 (Geological Society, London, 2013).Novack-Gottshall, P. M. Using a theoretical ecospace to quantify the ecological diversity of Paleozoic and modern marine biotas. Paleobiology 33, 273–294 (2007).
Google Scholar 
Wagner, P. J. On the probabilities of branch durations and stratigraphic gaps in phylogenies of fossil taxa when rates of diversification and sampling vary over time. Paleobiology 45, 30–55 (2019).
Google Scholar 
Sumrall, C. D. & Wray, G. A. Ontogeny in the fossil record: diversification of body plans and the evolution of “aberrant” symmetry in Paleozoic echinoderms. Paleobiology 33, 149–163 (2007).
Google Scholar 
Guensburg, T. E. & Sprinkle, J. in The Ecology of the Cambrian Radiation (eds Zhuravlev, A. I. U. & Riding, R.) 428–444 (Columbia Univ. Press, 2001).Zamora, S., Deline, B., Álvaro, J. J. & Rahman, I. A. The Cambrian Substrate Revolution and the early evolution of attachment in suspension-feeding echinoderms. Earth Sci. Rev. 171, 478–491 (2017).
Google Scholar 
Lloyd, G. T., Guillerme, T., Sherratt, E. & Wang, S. C. Claddis: measuring morphological diversity and evolutionary tempo. R package version 0.4.0 https://cran.r-project.org/package=Claddis (2020).Lloyd, G. T., Wang, S. C. & Brusatte, S. L. Identifying heterogeneity in rates of morphological evolution: discrete character change in the evolution of lungfish (Sarcopterygii; Dipnoi). Evolution 66, 330–348 (2012).PubMed 
PubMed Central 

Google Scholar 
Wagner, P. J. Early bursts of disparity and the reorganization of character integration. Proc. R. Soc. B 285, 20181604 (2018).PubMed 
PubMed Central 

Google Scholar 
Lewis, P. O. A likelihood approach to estimating phylogeny from discrete morphological character data. Syst. Biol. 50, 913–925 (2001).CAS 
PubMed 
PubMed Central 

Google Scholar 
Lloyd, G. T. Estimating morphological diversity and tempo with discrete character–taxon matrices: implementation, challenges, progress, and future directions. Biol. J. Linn. Soc. 118, 131–151 (2016).
Google Scholar 
Harmon, L. J. et al. Early bursts of body size and shape evolution are rare in comparative data. Evolution 64, 2385–2396 (2010).
Google Scholar 
Knope, M. L., Heim, N. A., Frishkoff, L. O. & Payne, J. L. Limited role of functional differentiation in early diversification of animals. Nat. Commun. 6, 6455–6461 (2015).CAS 

Google Scholar 
Cole, S. R., Wright, D. F. & Ausich, W. I. Phylogenetic community paleoecology of one of the earliest complex crinoid faunas (Brechin Lagerstätte, Ordovician). Palaeogeogr. Palaeoclimatol. Palaeoecol. 521, 82–98 (2019).
Google Scholar 
Revell, L. J. phytools: phylogenetic tools for comparative biology (and other things). R package version 0.7-47 https://cran.r-project.org/package=phytools (2020).Stayton, C. T. The definition, recognition, and interpretation of convergent evolution, and two new measures for quantifying and assessing the significance of convergence. Evolution 69, 2140–2153 (2015).PubMed 
PubMed Central 

Google Scholar 
Erwin, D. H. A conceptual framework of evolutionary novelty and innovation. Biol. Rev. 96, 1–15 (2021).
Google Scholar 
Lloyd, G. T. Journeys through discrete-character morphospace: synthesizing phylogeny, tempo, and disparity. Palaeontology 61, 637–645 (2018).
Google Scholar 
Anderson, P. S. L. & Friedman, M. Using cladistic characters to predict functional variety: experiments using early gnathostomes. J. Vertebr. Paleontol. 32, 1254–1270 (2012).
Google Scholar 
Deline, B. et al. Evolution of metazoan morphological disparity. Proc. Natl Acad. Sci. USA 15, E8909–E8918 (2018).
Google Scholar 
Matzke, N. J. & Irmis, R. B. Including autapomorphies is important for paleontological tip-dating with clocklike data, but not with non-clock data. PeerJ 6, e4553 (2018).PubMed 
PubMed Central 

Google Scholar 
Laing, A. M. et al. Giant taxon–character matrices: the future of morphological systematics. Cladistics 34, 333–335 (2018).PubMed 
PubMed Central 

Google Scholar 
Bush, A. M. & Bambach, R. K. Paleoecologic megatrends in marine Metazoa. Annu. Rev. Earth Planet. Sci. 39, 241–269 (2011).CAS 

Google Scholar 
Knope, M. L., Bush, A. M., Frishkoff, L. O., Heim, N. A. & Payne, J. L. Ecologically diverse clades dominate the oceans via extinction resistance. Science 367, 1035–1038 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Pimiento, C. et al. Functional diversity of marine megafauna in the Anthropocene. Sci. Adv. 6, eaay7650 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing v.4.0.2 (R Foundation for Statistical Computing, 2020).Zamora, S. & Rahman, I. A. Deciphering the early evolution of echinoderms with Cambrian fossils. Palaeontology 57, 1105–1119 (2014).
Google Scholar 
Cohen, K. M., Finney, S. C., Gibbard, P. L. & Fan, J.-X. The ICS International Chronostratigraphic Chart. Episodes 36, 199–204 (2013).
Google Scholar 
Bapst, D. W., Bullock, P. C., Melchin, M. J., Sheets, H. D. & Mitchell, C. E. Graptoloid diversity and disparity became decoupled during the Ordovician mass extinction. Proc. Natl Acad. Sci. USA 109, 3428–3433 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Bapst, D. W. A stochastic rate-calibrated method for time-scaling phylogenies of fossil taxa. Methods Ecol. Evol. 4, 724–733 (2013).
Google Scholar 
Bapst, D. W. Assessing the effect of time-scaling methods on phylogeny-based analyses in the fossil record. Paleobiology 40, 331–351 (2014).
Google Scholar 
Bapst, D. W. & Hopkins, M. J. Comparing cal3 and other a posteriori time-scaling approaches in a case study with the pterocephaliid trilobites. Paleobiology 43, 49–67 (2016).
Google Scholar 
Paradis, E. et al. ape: analyses of phylogenetics and evolution. R package version 5.4 https://cran.r-project.org/package=ape (2020).Yang, Z., Kumar, S. & Nei, M. A new method of inference of ancestral nucleotide and amino acid sequences. Genetics 141, 1641–1650 (1995).CAS 
PubMed 
PubMed Central 

Google Scholar 
Guillerme, T. et al. Disparities in the analysis of morphological disparity. Biol. Lett. 16, 20200199 (2020).PubMed 
PubMed Central 

Google Scholar 
Brusatte, S. L., Montanari, S., Yi, H.-Y. & Norell, M. A. Phylogenetic corrections for morphological disparity analysis: new methodology and case studies. Paleobiology 37, 1–22 (2011).
Google Scholar 
Guillerme, T. & Cooper, N. Time for a rethink: time sub-sampling methods in disparity-through-time analyses. Palaeontology 61, 481–493 (2018).
Google Scholar 
Wills, M. A. Crustacean disparity through the Phanerozoic: comparing morphological and stratigraphic data. Biol. J. Linn. Soc. 65, 455–500 (1998).
Google Scholar 
Lehmann, O. E. R., Ezcurra, M. D., Butler, R. J. & Lloyd, G. T. Biases with the Generalized Euclidean Distance measure in disparity analyses with high levels of missing data. Palaeontology 62, 837–849 (2019).
Google Scholar 
Hopkins, M. J. & John, K. S. A new family of dissimilarity metrics for discrete character matrices that include inapplicable characters and its importance for disparity studies. Proc. R. Soc. B 285, 20181784 (2018).PubMed 
PubMed Central 

Google Scholar 
Chessel, D., Dufour, A. B. & Thioulouse, J. ade4: analysis of ecological data–exploratory and Euclidean methods in environmental sciences. R package version 1.7-15 https://cran.r-project.org/package=ade4 (2020).Maire, E., Grenouillet, G., Brosse, S. & Villéger, S. How many dimensions are needed to accurately assess functional diversity? A pragmatic approach for assessing the quality of functional spaces. Glob. Ecol. Biogeogr. 24, 728–740 (2015).
Google Scholar 
Mouchet, M. A., Villéger, S., Mason, N. W. H. & Mouillot, D. Functional diversity measures: an overview of their redundancy and their ability to discriminate community assembly rules. Funct. Ecol. 24, 867–876 (2010).
Google Scholar 
Novack-Gottshall, P. M. General models of ecological diversification. II. Simulations and empirical applications. Paleobiology 42, 209–239 (2016).Villéger, S., Novack-Gottshall, P. M. & Mouillot, D. The multidimensionality of the niche reveals functional diversity changes in benthic marine biotas across geological time. Ecol. Lett. 14, 561–568 (2011).PubMed 
PubMed Central 

Google Scholar 
Mouillot, D. et al. The dimensionality and structure of species trait spaces. Ecol. Lett. https://doi.org/10.1111/ele.13778 (2021).Wills, M. A. in Fossils, Phylogeny, and Form: An Analytical Approach (eds Adrain, J. M. et al.) 55–143 (Kluwer Academic/Plenum, 2001).Foote, M. Discordance and concordance between morphological and taxonomic diversity. Paleobiology 19, 185–204 (1993).
Google Scholar 
Ciampaglio, C. N., Kemp, M. & McShea, D. W. Detecting changes in morphospace occupation patterns in the fossil record: characterization and analysis of measures of disparity. Paleobiology 27, 695–715 (2001).
Google Scholar 
Hopkins, M. J. & Gerber, S. in Evolutionary Developmental Biology: A Reference Guide (eds de la Rosa, L. N. & Müller, G.) 1–12 (Springer, 2017).Villéger, S., Mason, N. W. H. & Mouillot, D. New multidimensional functional diversity indices for a multifaceted framework in functional ecology. Ecology 89, 2290–2301 (2008).PubMed 
PubMed Central 

Google Scholar 
Anderson, M. J., Ellingsen, K. E. & McArdle, B. H. Multivariate dispersion as a measure of beta diversity. Ecol. Lett. 9, 683–693 (2006).PubMed 
PubMed Central 

Google Scholar 
Laliberté, E. & Legendre, P. A distance-based framework for measuring functional diversity from multiple traits. Ecology 91, 299–305 (2010).PubMed 
PubMed Central 

Google Scholar 
Novack-Gottshall, P. ecospace: simulating community assembly and ecological diversification using ecospace frameworks. R package version 1.0.1 https://cran.r-project.org/package=ecospace (2015).Laliberté, E. & Shipley, B. FD: measuring functional diversity from multiple traits, and other tools for functional ecology. R package version 1.0-12 https://cran.r-project.org/package=FD (2014).Dineen, A. A., Roopnarine, P. D. & Fraiser, M. L. Ecological continuity and transformation after the Permo-Triassic mass extinction in northeastern Panthalassa. Biol. Lett. 15, 20180902 (2019).PubMed 
PubMed Central 

Google Scholar 
Efron, B. Bootstrap methods: another look at the jackknife. Ann. Stat. 7, 1–26 (1979).
Google Scholar 
Kowalewski, M. & Novack-Gottshall, P. in Quantitative Methods in Paleobiology. Short Courses in Paleontology 16 (eds Alroy, J. & Hunt, G.) 19–54 (Paleontological Society and Paleontological Research Institute, 2010).Lingoes, J. C. Some boundary conditions for a monotone analysis of symmetric matrices. Psychometrika 36, 195–203 (1971).
Google Scholar 
Legendre, P. & Anderson, M. J. Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. Ecol. Monogr. 69, 1–24 (1999).
Google Scholar 
Burnham, K. P. & Anderson, D. R. Model Selection and Multi-Model Inference: A Practical Information-Theoretic Approach 2nd edn (Springer, 2002). More

Antonelli, A., Smith, R. J. & Simmonds, M. S. J. Unlocking the properties of plants and fungi for sustainable development. Nat. Plants 5, 1100–1102 (2019).PubMed 

Google Scholar 
Chen, I.-C., Hill, J. K., Ohlemuller, R., Roy, D. B. & Thomas, C. D. Rapid range shifts of species associated with high levels of climate warming. Science 333, 1024–1026 (2011).CAS 
PubMed 
ADS 

Google Scholar 
IPBES. Summary for Policymakers of the Global Assessment Report on Biodiversity and Ecosystem Services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. https://doi.org/10.5281/zenodo.3553579 (2019).Pecl, G. T. et al. Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. Science. 355, eaai9214 (2017).PubMed 

Google Scholar 
Warren, R. et al. Quantifying the benefit of early climate change mitigation in avoiding biodiversity loss. Nat. Clim. Chang. 3, 678–682 (2013).ADS 

Google Scholar 
Destro, G. F. G., Fernandes, V., Andrade, A. F. A., De Marco, P. & Terribile, L. C. Back home? Uncertainties for returning seized animals to the source-areas under climate change. Glob. Change Biol. 25, 3242–3253 (2019).ADS 

Google Scholar 
Travis, J. M. J. et al. Dispersal and species’ responses to climate change. Oikos 122, 1532–1540 (2013).
Google Scholar 
IPCC. Summary for policymakers. In Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (2014).IPCC. Summary for policymakers. In Climate Change and Land: An IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems (eds Shukla, P.R. & J. Skea, E. C.) (2019).IPCC. Special Report on 1.5 degrees: Summary for Policymakers. In Global Warming of 1.5°C. An IPCC Special Report on the Impacts of global Warming of 1.5°C Above Pre-industrial Levels and Related Global Greenhouse Gas Emission Pathways, in the Context of Strengthening the Global Response to the Threat of Climate Cha (2018).Ulloa Ulloa, C. et al. An integrated assessment of the vascular plant species of the Americas. Science 358, 1614–1617 (2017).CAS 
PubMed 
ADS 

Google Scholar 
Myers, N., Mittermeier, R. A., Mittermeier, C. G., da Fonseca, G. A. B. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).CAS 
PubMed 
ADS 

Google Scholar 
Coradin, L., Siminski, A. & Reis, A. Espécies Nativas da Flora Brasileira de Valor Econômico Atual e Potencial – Plantas para o futuro – Região Sul. (Ministério do Meio Ambiente, 2011).Nair, P. K. R. An introduction to agroforestry (Springer, 1993).
Google Scholar 
Sinclair, F. L. A general classification of agroforestry practice. Agrofor. Syst. 46, 161–180 (1999).
Google Scholar 
Somarriba, E. Revisiting the past: An essay on agroforestry definition. Agrofor. Syst. 19, 233–240 (1992).
Google Scholar 
Cerda, R. et al. Contribution of cocoa agroforestry systems to family income and domestic consumption: Looking toward intensification. Agrofor. Syst. 88, 957–981 (2014).
Google Scholar 
Montagnini, F. Integrating Landscapes: Agroforestry for Biodiversity Conservation and Food Sovereignty Vol. 12 (Springer, New York, 2017).
Google Scholar 
Siddique, I., Dionísio, A. C. & Simões-Ramos, G. A. Construindo Conhecimentos Sobre Agroflorestas em Rede. (UFSC, 2017).Jose, S. Agroforestry for conserving and enhancing biodiversity. Agrofor. Syst. 85, 1–8 (2012).
Google Scholar 
Sistla, S. A. et al. Agroforestry Practices Promote Biodiversity and Natural Resource Diversity in Atlantic Nicaragua. PLoS One 11, e0162529 (2016).PubMed 
PubMed Central 

Google Scholar 
Santos, P. Z. F., Crouzeilles, R. & Sansevero, J. B. B. Can agroforestry systems enhance biodiversity and ecosystem service provision in agricultural landscapes? A meta-analysis for the Brazilian Atlantic Forest. For. Ecol. Manage. 433, 140–145 (2019).
Google Scholar 
Reppin, S., Kuyah, S., de Neergaard, A., Oelofse, M. & Rosenstock, T. S. Contribution of agroforestry to climate change mitigation and livelihoods in Western Kenya. Agrofor. Syst. 94, 203–220 (2020).
Google Scholar 
Marconi, L. & Armengot, L. Complex agroforestry systems against biotic homogenization: The case of plants in the herbaceous stratum of cocoa production systems. Agric. Ecosyst. Environ. 287, 106664 (2020).CAS 

Google Scholar 
Somarriba, E. et al. Carbon stocks and cocoa yields in agroforestry systems of Central America. Agric. Ecosyst. Environ. 173, 46–57 (2013).
Google Scholar 
De Stefano, A. & Jacobson, M. G. Soil carbon sequestration in agroforestry systems: A meta-analysis. Agrofor. Syst. 92, 285–299 (2017).
Google Scholar 
Gomes, L. C. et al. Agroforestry systems can mitigate the impacts of climate change on coffee production: A spatially explicit assessment in Brazil. Agric. Ecosyst. Environ. 294, 106858 (2020).
Google Scholar 
Kofsky, J., Zhang, H. & Song, B.-H. The Untapped Genetic Reservoir: The Past, Current, and Future Applications of the Wild Soybean (Glycine soja). Front. Plant Sci. 9, 285–299 (2018).
Google Scholar 
Lorenzi, H. Arvores Brasileiras. (Plantarum, 2016).Zwiener, V. P. et al. Planning for conservation and restoration under climate and land use change in the Brazilian Atlantic Forest. Divers. Distrib. 23, 955–966 (2017).
Google Scholar 
Zechini, A. A. et al. Genetic conservation of Brazilian Pine (Araucaria angustifolia) through traditional land use. Econ. Bot. 72, 166–179 (2018).
Google Scholar 
Donazzolo, J., Stefenon, V. M., Guerra, M. P. & Nodari, R. O. On farm management of Acca sellowiana (Myrtaceae) as a strategy for conservation of species genetic diversity. Sci. Hortic. (Amsterdam) 259, 108826 (2020).CAS 

Google Scholar 
Favreto, R., Mello, R. S. P. & de Moura Baptista, L. R. Growth of Euterpe edulis Mart (Arecaceae) under forest and agroforestry in southern Brazil. Agrofor. Syst. https://doi.org/10.1007/s10457-010-9321-z (2010).Article 

Google Scholar 
Siminski, A., dos Santos, K. L. & Wendt, J. G. N. Rescuing agroforestry as strategy for agriculture in Southern Brazil. J. For. Res. 27, 739–746 (2016).
Google Scholar 
da Silva, L. C. R., Machado, S. A., Galvão, F. & Filho, A. F. Floristic evolution in an agroforestry system cultivation in Southern Brazil. An. Acad. Bras. Cienc. https://doi.org/10.1590/0001-3765201620150026 (2016).Article 
PubMed 

Google Scholar 
Gomes, V. H. F. et al. Species distribution modelling: Contrasting presence-only models with plot abundance data. Sci. Rep. 8, 1003 (2018).PubMed 
PubMed Central 
ADS 

Google Scholar 
Guisan, A. & Thuiller, W. Predicting species distribution: Offering more than simple habitat models. Ecol. Lett. 8, 993–1009 (2005).PubMed 

Google Scholar 
Raes, N. & Aguirre-Gutiérrez, J. A Modeling Framework to Estimate and Project Species Distributions in Space and Time Pontocaspian biodiversity RIse and DEmise View project Current and Future Biodiversity Patterns in Mainland Southeast Asia View project. (2018).Brooks, T. M. et al. Measuring terrestrial area of habitat (AOH) and its utility for the IUCN red list. Trends Ecol. Evol. 34, 977–986 (2019).PubMed 

Google Scholar 
Soberon, J. & Peterson, A. T. Interpretation of models of fundamental ecological niches and species’ distributional areas. Biodivers. Inform. 2, 1–10 (2005).
Google Scholar 
Levis, C. et al. Persistent effects of pre-Columbian plant domestication on Amazonian forest composition. Science 355, 925–931 (2017).CAS 
PubMed 
ADS 

Google Scholar 
Reis, M. S. et al. Domesticated landscapes in Araucaria Forests, Southern Brazil: A multispecies local conservation-by-use system. Front. Ecol. Evol. 6, 1–14 (2018).
Google Scholar 
IUCN Standards and Petitions Committee. Guidelines for Using the IUCN Red List Categories and Criteria. Version 14. Prep. by Stand. Petitions Comm. (2019).Gomes, V. H. F., Vieira, I. C. G., Salomão, R. P. & ter Steege, H. Amazonian tree species threatened by deforestation and climate change. Nat. Clim. Chang. 9, 547–553 (2019).ADS 

Google Scholar 
Guo, Y. et al. Prediction of the potential geographic distribution of the ectomycorrhizal mushroom Tricholoma matsutake under multiple climate change scenarios. Sci. Rep. 7, 46221 (2017).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Rodrigues, P., Silva, J., Eisenlohr, P. & Schaefer, C. Climate change effects on the geographic distribution of specialist tree species of the Brazilian tropical dry forests. Braz. J. Biol. 75, 679–684 (2015).CAS 
PubMed 

Google Scholar 
Wilson, O. J., Walters, R. J., Mayle, F. E., Lingner, D. V. & Vibrans, A. C. Cold spot microrefugia hold the key to survival for Brazil’s critically endangered araucaria tree. Glob. Chang. Biol. 25, 4339–4351 (2019).PubMed 
ADS 

Google Scholar 
Cámara-Leret, R. et al. Climate change threatens New Guinea’s biocultural heritage. Sci. Adv. 5, eaaz1455 (2019).PubMed 
PubMed Central 
ADS 

Google Scholar 
Esser, L. F., Saraiva, D. D. & Jarenkow, J. A. Future uncertainties for the distribution and conservation of Paubrasilia echinata under climate change. Acta Bot. Brasilica 33, 770–776 (2019).
Google Scholar 
Lima, V. P., Marchioro, C. A., Joner, F., ter Steege, H. & Siddique, I. Extinction threat to neglected Plinia edulis exacerbated by climate change, yet likely mitigated by conservation through sustainable use. Austral Ecol. 45, 376–383 (2020).
Google Scholar 
Santini, L., Benítez-López, A., Maiorano, L., Čengić, M. & Huijbregts, M. A. J. Assessing the reliability of species distribution projections in climate change research. Divers. Distrib. 27, 1–16 (2021).
Google Scholar 
Raes, N. et al. Historical distribution of Sundaland’s Dipterocarp rainforests at Quaternary glacial maxima. Proc. Natl. Acad. Sci. 111, 16790–16795 (2014).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Vaz, Ú. L. & Nabout, J. C. Using ecological niche models to predict the impact of global climate change on the geographical distribution and productivity of Euterpe oleracea Mart. (Arecaceae) in the Amazon. Acta Bot. Brasilica 30, 290–295 (2016).
Google Scholar 
Sánchez-Fernández, D. et al. Thermal niche estimators and the capability of poor dispersal species to cope with climate change. Sci. Rep. 6, 23381 (2016).PubMed 
PubMed Central 
ADS 

Google Scholar 
de Lima, R. A. F. et al. How much do we know about the endangered Atlantic Forest? Reviewing nearly 70 years of information on tree community surveys. Biodivers. Conserv. 24, 2135–2148 (2015).
Google Scholar 
Ribeiro, M. C. et al. The Brazilian Atlantic Forest: A Shrinking Biodiversity Hotspot (Springer, New York, 2011).
Google Scholar 
Díaz, S. et al. Assessing nature’s contributions to people. Science 359, 270–272 (2018).PubMed 
PubMed Central 
ADS 

Google Scholar 
Siddique, I. et al. Woody species richness drives synergistic recovery of socio-ecological multifunctionality along early tropical dry forest regeneration. For. Ecol. Manag. 482, 118848 (2021).
Google Scholar 
Harvey, C. A. et al. Climate-smart landscapes: Opportunities and challenges for integrating adaptation and mitigation in tropical agriculture. Conserv. Lett. 7, 77–90 (2014).
Google Scholar 
Schneidewind, U. et al. Carbon stocks, litterfall and pruning residues in monoculture and agroforestry cacao production systems. Exp. Agric. 55, 452–470 (2019).
Google Scholar 
Dinesh, D., Campbell, B. M., Bonilla-findji, O. & Richards, M. 10 Best Bet Innovations for Adaptation in Agriculture: A supplement to the UNFCCC NAP Technical Guidelines. Working paper 215 (2017).Lin, B. B., Perfecto, I. & Vandermeer, J. Synergies between agricultural intensification and climate change could create surprising vulnerabilities for crops. Bioscience 58, 847–854 (2008).
Google Scholar 
Perfecto, I., John Vandermeer & Angus Wright. 2019 Nature’s Matrix: Linking Agriculture, Biodiversity Conservation and Food Sovereignty. (Routledge, 2019).Hannah, L. et al. 30% land conservation and climate action reduces tropical extinction risk by more than 50%. Ecography 43, 943–953 (2020).
Google Scholar 
Zizka, A. et al. Biogeography and conservation status of the pineapple family (Bromeliaceae). Divers. Distrib. 26, 183–195 (2020).
Google Scholar 
Elias, G. A., Lima, J. M. T. & dos Santos, R. Threatened flora from the State of Santa Catarina, Brazil: Arecaceae. Hoehnea 46, e322018 (2019).
Google Scholar 
Pimm, S. L. et al. The biodiversity of species and their rates of extinction, distribution, and protection. Science 344, 1246752–1246752 (2014).CAS 
PubMed 

Google Scholar 
Brancalion, P. H. S. et al. What makes ecosystem restoration expensive? A systematic cost assessment of projects in Brazil. Biol. Conserv. 240, 108274 (2019).
Google Scholar 
Crouzeilles, R. et al. There is hope for achieving ambitious Atlantic Forest restoration commitments. Perspect. Ecol. Conserv. 17, 80–83 (2019).
Google Scholar 
Magnago, L. F. S. et al. Would protecting tropical forest fragments provide carbon and biodiversity cobenefits under REDD+?. Glob. Chang. Biol. 21, 3455–3468 (2015).PubMed 
ADS 

Google Scholar 
Rodrigues, A. C., Villa, P. M. & Neri, A. V. Fine-scale topography shape richness, community composition, stem and biomass hyperdominant species in Brazilian Atlantic forest. Ecol. Indic. 102, 208–217 (2019).
Google Scholar 
de Lima, R. A. F. et al. The erosion of biodiversity and biomass in the Atlantic Forest biodiversity hotspot. Nat. Commun. 11, 6347 (2020).PubMed 
PubMed Central 
ADS 

Google Scholar 
Loreau, M. Reconciling utilitarian and non-utilitarian approaches to biodiversity conservation. Ethics Sci. Environ. Polit. 14, 27–32 (2014).
Google Scholar 
Berkes, F. & Folke, C. Linking social and ecological resilience and sustainability. In Linking Social and Ecological Systems. Management Practices and Social Mechanisms for Building Resilience (Cambridge University Press, Cambridge, 2000).Fernandes, R. C. & Piovezana, L. The Kaingang perspectives on land and environmental rights in the south of Brazil. Ambient. Soc. 18, 111–128 (2015).
Google Scholar 
Machado Mello, A. J. & Peroni, N. Cultural landscapes of the Araucaria Forests in the northern plateau of Santa Catarina, Brazil. J. Ethnobiol. Ethnomed. 11, 51 (2015).PubMed 
PubMed Central 

Google Scholar 
Thuiller, W., Guéguen, M., Renaud, J., Karger, D. N. & Zimmermann, N. E. Uncertainty in ensembles of global biodiversity scenarios. Nat. Commun. 10, 1446 (2019).PubMed 
PubMed Central 
ADS 

Google Scholar 
Zurell, D. et al. A standard protocol for reporting species distribution models. Ecography 43, 1261–1277 (2020).
Google Scholar 
Warren, D. L., Matzke, N. J. & Iglesias, T. L. Evaluating presence-only species distribution models with discrimination accuracy is uninformative for many applications. J. Biogeogr. 47, 167–180 (2020).
Google Scholar 
Leroy, B. et al. Without quality presence-absence data, discrimination metrics such as TSS can be misleading measures of model performance. J. Biogeogr. 45, 1994–2002 (2018).
Google Scholar 
Araujo, M. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42–47 (2007).PubMed 

Google Scholar 
Raes, N. & ter Steege, H. A null-model for significance testing of presence-only species distribution models. Ecography 30, 727–736 (2007).
Google Scholar 
Newbold, T. Future effects of climate and land-use change on terrestrial vertebrate community diversity under different scenarios. Proc. R. Soc. B Biol. Sci. 285, 20180792 (2018).Loiselle, B. A. et al. Avoiding pitfalls of using species distribution models in conservation planning. Conserv. Biol. 17, 1591–1600 (2003).Bean, W. T., Stafford, R. & Brashares, J. S. The effects of small sample size and sample bias on threshold selection and accuracy assessment of species distribution models. Ecography 35, 250–258 (2012).
Google Scholar 
Meyer, A. L. S., Pie, M. R. & Passos, F. C. Assessing the exposure of lion tamarins (Leontopithecus spp.) to future climate change. Am. J. Primatol. 76, 551–562 (2014).PubMed 

Google Scholar 
Araújo, M. B. & Pearson, R. G. Equilibrium of species’ distributions with climate. Ecography 28, 693–695 (2005).
Google Scholar 
Guillera-Arroita, G. et al. Is my species distribution model fit for purpose? Matching data and models to applications. Glob. Ecol. Biogeogr. 24, 276–292 (2015).
Google Scholar 
Bascompte, J., García, M. B., Ortega, R., Rezende, E. L. & Pironon, S. Mutualistic interactions reshuffle the effects of climate change on plants across the tree of life. Sci. Adv. 5, eaav2539 (2019).PubMed 
PubMed Central 
ADS 

Google Scholar 
Hoffmann, A. A. & Sgrò, C. M. Climate change and evolutionary adaptation. Nature 470, 479–485 (2011).CAS 
PubMed 
ADS 

Google Scholar 
Thuiller, W. et al. Predicting global change impacts on plant species’ distributions: Future challenges. Perspect. Plant Ecol. Evol. Syst. 9, 137–152 (2008).
Google Scholar 
Mayle, F. E. Millennial-scale dynamics of southern Amazonian rain forests. Science 290, 2291–2294 (2000).CAS 
PubMed 
ADS 

Google Scholar 
Bullock, J. M. et al. Human-mediated dispersal and the rewiring of spatial networks. Trends Ecol. Evol. 33, 958–970 (2018).PubMed 

Google Scholar 
Ordonez, J. C. Constraints and opportunities for tree diversity management along the forest transition curve to achieve multifunctional agriculture. Curr. Opin. Environ. Sustain. 6, 54–60 (2014).Levis, C. et al. How People Domesticated Amazonian Forests. Front. Ecol. Evol. 5, 171 (2018).Mendes, P., Velazco, S. J. E., de Andrade, A. F. A. & De Marco, P. Dealing with overprediction in species distribution models: How adding distance constraints can improve model accuracy. Ecol. Modell. 431, 109180 (2020).
Google Scholar 
GBIF. GBIF Occurrence. https://www.gbif.org, https://doi.org/10.15468/dl.vjezvb (2019)Carvalho, G. flora: Tools for Interacting with the Brazilian Flora 2020. R package version 0.3.0. (2017).Raes, N. Partial versus full species distribution models. Nat. Conserv. 10, 127–138 (2012).
Google Scholar 
Zizka, A. et al. CoordinateCleaner: Standardized cleaning of occurrence records from biological collection databases. Methods Ecol. Evol. 10, 744–751 (2019).
Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna (2020).Oliveira, U. et al. The strong influence of collection bias on biodiversity knowledge shortfalls of Brazilian terrestrial biodiversity. Divers. Distrib. 22, 1232–1244 (2016).
Google Scholar 
Daru, B. H. et al. Widespread sampling biases in herbaria revealed from large-scale digitization. New Phytol. 217, 939–955 (2018).PubMed 

Google Scholar 
Aiello-Lammens, M. E., Boria, R. A., Radosavljevic, A., Vilela, B. & Anderson, R. P. spThin: An R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography 38, 541–545 (2015).
Google Scholar 
Proosdij, A. S. J., Sosef, M. S. M., Wieringa, J. J. & Raes, N. Minimum required number of specimen records to develop accurate species distribution models. Ecography 39, 542–552 (2016).
Google Scholar 
Beaumont, L. J. et al. Which species distribution models are more (or less) likely to project broad-scale, climate-induced shifts in species ranges?. Ecol. Modell. 342, 135–146 (2016).
Google Scholar 
Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).
Google Scholar 
Fourcade, Y., Besnard, A. G. & Secondi, J. Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics. Glob. Ecol. Biogeogr. 27, 245–256 (2018).
Google Scholar 
Austin, M. P. & Van Niel, K. P. Improving species distribution models for climate change studies: Variable selection and scale. J. Biogeogr. 38, 1–8 (2011).
Google Scholar 
Woodward, F. I. Climate and Plant Distribution. (Cambridge Univ. Press., 1987).IUCN. Plant Growth Forms Classification Scheme. Version: 1.0. https://www.iucnredlist.org/resources/classification-schemes (2020).Dormann, C. F. et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46 (2013).
Google Scholar 
Fremout, T. et al. Mapping tree species vulnerability to multiple threats as a guide to restoration and conservation of tropical dry forests. Glob. Chang. Biol. 26, 3552–3568 (2020).PubMed 
ADS 

Google Scholar 
Naimi, B. Package ‘ usdm ’. R Topics Document (2015).Syfert, M. M. et al. Using species distribution models to inform IUCN Red List assessments. Biol. Conserv. 177, 174–184 (2014).
Google Scholar 
Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Modell. 190, 231–259 (2006).
Google Scholar 
Elith, J., Kearney, M. & Phillips, S. The art of modelling range-shifting species. Methods Ecol. Evol. 1, 330–342 (2010).
Google Scholar 
Muñoz-Pajares, A. J. et al. Niche differences may explain the geographic distribution of cytotypes in Erysimum mediohispanicum. Plant Biol. 20, 139–147 (2018).PubMed 

Google Scholar 
Peng, L.-P. et al. Modelling environmentally suitable areas for the potential introduction and cultivation of the emerging oil crop Paeonia ostii in China. Sci. Rep. 9, 3213 (2019).PubMed 
PubMed Central 
ADS 

Google Scholar 
Merow, C., Smith, M. J. & Silander, J. A. A practical guide to MaxEnt for modeling species’ distributions: What it does, and why inputs and settings matter. Ecography 36, 1058–1069 (2013).
Google Scholar 
Boucher-Lalonde, V., Morin, A. & Currie, D. J. How are tree species distributed in climatic space? A simple and general pattern. Glob. Ecol. Biogeogr. 21, 1157–1166 (2012).
Google Scholar 
Elith, J., Ferrier, S., Huettmann, F. & Leathwick, J. The evaluation strip: A new and robust method for plotting predicted responses from species distribution models. Ecol. Modell. 186, 280–289 (2005).
Google Scholar 
Jiménez-Valverde, A. & Lobo, J. M. Threshold criteria for conversion of probability of species presence to either-or presence-absence. Acta Oecologica 31, 361–369 (2007).ADS 

Google Scholar 
Betts, J. et al. A framework for evaluating the impact of the IUCN Red List of threatened species. Conserv. Biol. 34, 632–643 (2020).PubMed 
PubMed Central 

Google Scholar 
ter Steege, H. et al. Estimating the global conservation status of more than 15,000 Amazonian tree species. Sci. Adv. 1, e1500936 (2015).PubMed 
PubMed Central 
ADS 

Google Scholar 
Dauby, G. et al. ConR : An R package to assist large-scale multispecies preliminary conservation assessments using distribution data. Ecol. Evol. 7, 11292–11303 (2017).PubMed 
PubMed Central 

Google Scholar  More

Settlement stage fishesSpecimen collectionData for settlement stage tropical larval fishes includes 1372 swimming speed measurements, collected from >75 unique taxa across 35 families of fishes, most of which are coral reef associated as adults. The data are collected from five locations, including: South Caicos Island (Turks and Caicos Islands, Caribbean – TCI), Green Island or Magnetic Island (exact location not recorded for individual samples, Great Barrier Reef, Australia GI/MI), Lizard Island (Great Barrier Reef, Australia – LI), Calabash Caye (Turneffe Islands Atoll, Belize – BLZ), and Moorea, Society Islands (MOR) (Table 1).Table 1 Sampling locations for settlement stage Ucrit data. Included are the Region, year of collection, location and associated name, as well as the total number of recorded measurements (count).Full size tableData from LI and TCI were obtained almost exclusively from specimens collected using light traps, placed 100–500 m meters off the leeward side of the Island, near either the School for Field Studies facilities (TCI) or the Lizard Island Research Station (LI). An additional eight specimens of newly hatched Acanthochromis polyacanthus, a pomacentrid species which does not have a pelagic phase, were captured with nets on the Lizard Island reefs. Data from GI/MI were obtained using a combination of light traps, beach seines, fence and dip nets.For data collected in Moorea (Mor), specimens arriving over night from the open ocean and attempting to settle on the reef were captured in nets placed on the reef crest. In Belize (BLZ) specimens were collected using a variety of techniques including crest nets, channel nets, light traps and night-light lift nets, although most individuals were collected using light traps and crest nets. Crest net locations were those reported in29. Unless stated otherwise in the “notes” field of the “ucrit_sett_dat” data table (see Online-only Table 1), all U-crit measurements on individuals captured by light traps or crest nets were made on the morning of capture, usually within 6 or 12 hours (please refer to the original publications for methodological details). In a few cases, some individuals were kept in the laboratory for up to 2 days to study changes in swimming speed associated with settlement (see24,30), and here the post-settlement status of the larvae was recorded in the field “stage” in the “fish_id_dat” data table (see Online-only Table 1). Some specimens of the pomacentrid, Abudefduf saxatilis were collected with hand nets from a fish attracting device deployed over a seagrass bed from a dock and are best considered as early post-settlement individuals, although the time since settlement is unknown. All specimens of the labrid, Clepticus parrae were collected with hand-nets from deep fore-reefs. Although they had settled to the fore-reef an unknown period of time before capture, these individuals had yet to undergo complete metamorphosis. Data from such post-settlement individuals should be used with caution, as it is known that swimming performance in some species decreases markedly upon settlement11,24.Most data were collected during the summer months (May through September for TCI and BLZ, November through February for LI), but in MOR, the data came from winter (August). In Belize, data were collected in 2003, 2004 and 2005, totalling 401 U-crit measurements. Data from TCI were collected in 2003, from 109 individuals. Data from LI represented just over half of the settlement stage swimming data (556 measurements), and were collected in 2001, 2002, 2003, 2004 and 2005. A total of 144 measurement were available from GI/MI and were collected in 1992. The 152 U-crit measurements from MOR were from 2010.Captured settlement stage larvae/pelagic juveniles were held in fresh seawater for a minimum of 1–2 hours to reduce stress prior to swimming trials, either with an aeration stone in 24 L buckets (BLZ), or in flow-through seawater aquarium facilities at the Lizard Island Research Station (LI) and the Department of Environment and Coastal Resources at South Caicos (TCI).U-crit protocolAll swimming experiments were conducted at ambient seawater temperatures, which ranged between 25 °C and 30 °C, depending on location and date.Settlement-stage individuals were swum using one of several swimming flumes, including a single-lane swimming chamber11,22, a six-lane swimming chamber12 (see Fig. 1), or a three-lane swimming chamber modified from the design of12. All swimming chambers were constructed from transparent Plexiglass (internal dimensions of swimming area: 185 mm × 50 mm × 50 mm). A removable lid, sealed with an O-ring was used to introduce fish to, and remove them from the chambers. One section of flow straighteners, 45-mm long, was placed just after the inflow in order to reduce turbulence within the chamber. Fish were retained within the swimming area by two 4.0 mm mesh metal retaining fences, which were covered with a finer mesh when required for very small larvae.Fig. 1Design of the swimming channel used for settlement stage larvae, pelagic juveniles and settled juveniles. Shown are Side view (a) & Top view (b, modified from12). This swimming channel can be operated at up to 50 cm s−1. Smaller designs with three channels, or only 1 channel that could obtain higher speeds were used for swimming faster individuals. For data collected by Leis and colleagues, higher speeds were obtained by blocking some lanes using Perspex.Full size imageFlow was generated using a 2.4 Kw swimming pool pump (although the size of the pump varied across the studies), or plumbed into a laboratory seawater system. The speed was set using a protractor mounted on a gate valve and calibrated using the procedures described under the technical validation section below. Faster speeds were also calibrated using an inline blue-white F-300 series flow meter. Flumes were plumbed using union valves so they could be dismantled and easily relocated and installed in field locations. Because the pumps used to run the flumes can heat the water temperature over time, they were plumbed with a minimum reservoir volume of 70 L, with a constant flow through of fresh seawater. A mercury thermometer located in the reservoir was used to ensure temperature remained ambient during the swimming trials. Example field deployments of the swimming channels at various locations can be seen in Fig. 2.Fig. 2Examples of the ‘fast’ swimming flume setup at various field locations. Shown are the Lizard Island Research Station (Great Barrier Reef, Australia, a), the Department of Environment and Coastal Resources aquarium facilities at South Caicos (Turks and Caicos Islands, b), and the dock at the University of Belize Institute of Marine Field Studies at Calabash Caye (Belize, c).Full size imageU-crit was measured by placing specimens in the swimming flume and incrementally increasing water speed until the individual could no longer maintain position for the full-time increment interval. The exact experimental protocols differed slightly among the studies. For fish measured at Lizard Island in November 2000–January 2001, November–December 2001 and South Caicos Island, and most fish at Calabash Caye, the experimental protocol followed7, with speed increments equivalent to approximately three total standard body lengths per second (bls−1) with a time interval of 2 min. This protocol was adopted because settlement stage larvae can vary substantially in size and subsequently their swimming capacity, as swimming speeds are strongly controlled by body size4. Aligning the speed increments with the approximate size category of fishes ensured that the overall duration of the U-crit experiment was relatively similar. For fish measured at Lizard Island during December 2003, speed increments used were 1.6 cm s−1 with a time interval of 5 min. At Lizard Island in 2005, specimens of Amblyglyphidon curacao were subjected to speed increments of 4.2 cm s−1 at intervals of 5 minutes. In Moorea in 2010, all individuals were subjected to speed increments of 6.1 cm s−1 at intervals of 2 minutes. For fish measured at Green Island and Magnetic Island speed increments used were 5 cm s−1 with a time interval of 5 min. At Calabash Caye an experiment was conducted to examine the impact of time increments on U-crit measurements, and the experimental protocol was recorded in this instance.U-crit swimming speed was calculated following17:$$U mbox{-} crit=U+left(t/ti,ast ,Uiright)$$
(1)
Where:U = penultimate speed (speed increment for which the fish swam for the entire duration of the set time interval). Ui = the velocity increment (varied by the specific study).t = the time swum in the final velocity incrementti = the set time interval for each velocity increment (varied by the specific study).While the speed increments used varied across studies in this collated dataset, previous studies have found no effect of varying the length of the time interval (ti) in terms of the resulting swimming speed between fish swum at two minute intervals and those swum at 15 minute intervals for six reef fish species10.Sample handling and morphological measurementsAfter each trial specimens were anaesthetised in chilled water or using clove oil (depending on the location and according to the relevant ethics approvals) and some were photographed while still fresh to maintain body flexibility and to avoid issues with shrinkage due to dehydration associated with preserved samples. Following photographing, the samples were preserved in either 70% ethanol, 95% ethanol, or 10% buffered formalin.From digital images the ImageTool (UTHSCSA 2002) software was used for image analysis. Measurements made from digital images (where available) are shown in Fig. 3, and included: total length (from the outer edge of the caudal fin to the tip of the upper jaw), caudal fin length (from the tip of the caudal fin to the caudal peduncle), body depth (the vertical height of the fish measured at the deepest region), body area (the area of the fish in lateral view excluding the fins but including the head and gut region), propulsive area (the area of the fish including the fins but excluding the head and gut region), muscle area (the area of the fish excluding the fins and the head and gut region), caudal fin depth, caudal peduncle depth and caudal fin area. All measurements were taken to the nearest 0.1 mm. Body width (at the widest region, usually the head) was also measured to the nearest 0.1 mm using vernier callipers. In some cases total lengths (TL) were measured pre-trial using callipers (BLZ, 2003 and 2004). Body length (BL, which is equivalent to SL for postflexion stages) was measured using an ocular micrometer on a dissecting microscope in some studies23.Fig. 3Morphological measurements of settlement stage fishes. Measurements include: total length (TL; outer edge of the caudal fin to the tip of the upper jaw), caudal fin length (CFL; tip of the caudal fin to the caudal peduncle), body depth (BD; height at the deepest region), body area (BA; area in lateral view excluding the fins), propulsive area (PA; area including the fins (naturally fully extended) but excluding the head and gut region where they are inflexible or lack overlaying muscle and cannot be used for propulsion), muscle area (MA; area excluding the fins and the head and gut region), caudal fin depth (CFD; widest section when fully extended), caudal peduncle depth (CPD; height at the narrowest point between the caudal fin and the fish’s body) and caudal fin area (CFA; area with the caudal fins naturally fully extended). Callipers were used to measure head width. Adapted from8.Full size imageLarval development dataset (Australia)Rearing protocolData gathered using a combination of the ‘fast’ and ‘slow’ swimming chambers (see below) on swimming abilities throughout development are available for six species, including two damselfish – Pomacentrus amboinensis and Pomacentrus mollucensis (Pomacentridae; Pomacentrinae); two cardinalfish – Ostorhinchus (Apogon) compressus and Sphaeramia nematoptera (Apogonidae); and two anemone fish – Amphiprion percula and Amphiprion melanopus (Pomacentridae; Amphiprioninae). Note that the pomacentrids have demersal eggs, whereas the apogonids orally brood their eggs.Australian specimens for assessing swimming speeds throughout larval development were obtained mostly from larvae reared at the James Cook University Aquarium facility, from adult broodstock collected from the northern Great Barrier Reef. Adult brood stock were kept in outside aquaria ranging in size from 1000 to 3000 L. The temperature of aquaria was kept between 27 and 29.5 °C, with larvae reared in the Autumn and Winter of 1998. Brood stock were fed a diet of chopped pilchards, prawns and Ascetes twice per day. Eggs were obtained from spawning broodstock before dark on the night of hatching and transferred to a rearing tank. Once hatched, larvae were reared and maintained in 200 L (120 × 60 × 30 cm) black painted glass aquaria that were illuminated by four “daylight” fluorescent tubes. The larvae were maintained in a 14:10 light/dark photo-period at 27.5–29 °C. Cultures of the algae Nannochloropsis sp. were used to green the water during the day. This kept light at the right intensity to prevent “bashing” behaviour (young larvae have a tendency to continually butt their heads against surfaces if the water is clear and the light intensity is too bright). Larvae were fed a diet of >52 micron sieved wild caught plankton, which was occasionally supplemented by rotifers and Artemia spp. when necessary. Larvae were fed twice per day to maintain prey densities of between 2–6 individuals per ml. Examples of ontogenetic series obtained through these rearing methods, and showing pre- and post- flexions stages are show in Fig. 4.Fig. 4Examples of larval developmental series obtained for larvae reared at the James Cook University Aquarium facility. Showin are Amphiprion melanopus (a) and Sphariamia nemaptopera (b).Full size imageU-crit protocolSwimming experiments for older [i.e., postflexion] larvae (see Fig. 4(a)ii–iv and Fig. 4(b)iv–vii) were carried out using the flumes described above for settlement stage fishes. However, these flumes were unsuitable for measurement of swimming capabilities of the delicate younger [i.e., preflexion] larvae (see Fig. 4(a)i and Fig. 4(b)ii,iii). Several characteristics had to be addressed in order to design equipment suitable for the measurement of the swimming capabilities of very young larvae. These included:

The apparatus needed to produce slow flow rates while maintaining laminar flow and minimal boundary layer effects. This is because newly hatched larvae are small enough to effectively utilise the boundary layer, which is broader for slower moving water.

The apparatus had to provide an environment suitable for very young larvae as the trauma of transferring larvae between containers can be fatal. Accordingly, stress associated with sudden changes in light intensity or water quality was minimised by “greening” the water with algae and the use of dark or clear surfaces to avoid “bashing” behaviour. In addition, the apparatus had to be set up within the immediate vicinity of the rearing tanks (or possibly in a rearing tank) to minimise the distance larvae had to be moved.

Two swimming channels were designed and used for younger larvae that were able to meet these requirements. These were designed to operate at “slow” and “medium” speeds. Both channels were able to produce laminar flow at much slower speeds. They consisted of a much wider swimming area so that most of the water flow occurred away from the sides, maximising the area of water not influenced by boundary layer effects. Both were able to be placed in a rearing aquarium of “greened” water. This prevented the larvae from exhibiting “bashing” behaviour, minimised the distance that larvae had to be transferred and meant that there was no change in water quality between the experimental apparatus and rearing tanks (Fig. 5).Fig. 5Design of swimming channels for younger larvae. Shown are side view (a) & Btop view (b). Dimensions are for the “slow” and “medium” flumes respectively.Full size imageFish were retained by a 0.3 mm mesh at the end of the swimming channel for the “slow” chamber and a 1 mm mesh for the “medium” chamber. The “slow” channel was powered by an Eheim 2,000 L per hour pump and the “medium” channel was powered by two such Eheim pumps. The speed for both channels was set using a protractor mounted on a gate valve as for the “fast” swimming chamber used for older larvae, and calibrated according to the description below under technical validation.Each clutch of eggs from each species was raised from hatching through to settlement and experiments were performed periodically throughout this larval period, with sampling days depending on the species. The first swimming trial was conducted on day 1, approximately 12 hours after hatching. Three clutches of each species were used for each swimming trial for the species Pomacentrus amboinensis, Sphaeramia nematoptera and Amphiprion melanopus to ensure that any clutch effects were considered4. While multiple broodstock were available for each species, no record was made at the time from which broodstock the replicate clutches were obtained. For other species only a single clutch was available. In some cases these data included light trap caught specimens to supplement the latest settlement stage. At each experimental age for each clutch 8–12 fish were used in the swimming trials.Larvae were subjected to incremental increases in flow rates equivalent to approximately 3 body lengths (BL) every two minutes until they could no longer maintain position, as for the experimental protocol described above for settlement stage fishes and U-crit calculated as per Eq. 1. Aligning the speed increments with the approximate size category of fishes ensured that the overall duration of the U-crit experiment was relatively similar throughout ontogeny.Sample handling and morphological measurementsFish that were swum, or siblings from the same batch at the same age, were anaesthetised in chilled water then fixed in 10% buffered formalin. After 12–48 hours, larvae were transferred to 70% alcohol and stored. Morphological measurements were carried out by capturing the image of each fish using a stereo dissecting microscope linked to a video recorder. These images were then saved as files on computer. As for settlement stage larvae, the image analysis program ImageTool was then used to measure lengths and areas for different regions of the fish.Measurements were made of total length (from the tip of the caudal fin to the tip of the upper jaw), body depth (the height of the fish measured at the deepest region), body area (the entire area of the fish excluding the fins) and total propulsive area (the area of the fish including the fins but excluding the head and gut region). The regions measured for both pre-flexion and post flexion larvae can be seen in Fig. 6.Fig. 6Measurements made on developmental series larvae. This includes post-flexion larvae which have developed a true caudal fin supported by a hypural plate and discrete soft rays) (a); and pre-flexion larvae that had no hypural plate or soft rays, but a continuous rayless fin-fold from anus to nape (b).Full size imageLarval development dataset (Taiwan and France)Data on development of swimming in larvae of ten species of pelagic-spawning tropical species of commercially important fishes reared by aquaculturists in Taiwan6,27,28 and two tropical species that brood their eggs that were reared for the aquarium trade in France are included26,28 (see Table 2). In addition, very limited, previously unpublished data on larvae of three species of pelagic spawning, commercial species reared in Taiwan are included. The emphasis in these studies was on postflexion-stage larvae, but for some species, swimming data on preflexion and flexion-stage larvae are included. The ‘standard’ six-lane swimming chamber was used for these measurements of U-crit. For larger larvae of some species, half of the lanes were closed off to achieve the faster speeds that these larvae can achieve. Despite this adjustment, some individuals were able to swim faster than the fastest speeds the swim chamber could achieve. In these cases, the speeds are reported in the database as greater than the maximum chamber speed.Table 2 Species whose U-crit swimming ontogeny were studied using reared larvae in Taiwan and France. Included are the location, family, species, the number of specimens assayed (N), the size range of the specimens, and the associated publication of the original data (where available).Full size tableIn Taiwan, larvae were obtained from commercial aquaculture farms (~22.4°N, 120.6°E) SE of Kaohsiung, southern Taiwan, in May 2004 and in May and June 2005. Rearing conditions varied with species, but most were reared in outdoor concrete or earth ponds. Exceptions were Epinephelus spp., which were reared in indoor concrete tanks, and Chanos chanos and some Eleutheronema tetradactylum, which were reared in outdoor concrete tanks under shade cloth. In most cases, the larvae were provided with a “natural” food source (phytoplankton and zooplankton that were resident in the pond). The aquaculturists did not maintain breeding stock, but obtained the pelagic eggs for rearing from elsewhere. The larvae obtained from the aquaculturists were placed in oxygenated plastic bags placed in insulated boxes and transported about 1 h by road to the National Museum of Marine Biology and Aquarium (NMMBA), Kenting, Taiwan (~22.1°N, 120.7°E). In the laboratory the larvae were acclimated in 40 l aquaria filled from the NMMBA seawater system. Each aquarium was fitted with an aerator and kept at ca. 25 °C. The larvae were fed twice daily with live, newly hatched brine shrimp (Artemia nauplii) and 50% of the total volume of water was exchanged with fresh seawater. The aquaria were cleaned daily by suctioning debris off the bottom. The species studied in Taiwan were all native to the western central Pacific, but the original brood stock may not have been obtained locally. The U-crit measurements were made in a shaded outdoor area where large tanks were located to hold adult fishes intended for either research purposes or for addition to the large public aquarium that forms part of the NMMBA campus. The extensive seawater system of NMMBA was used to supply seawater directly into the swimming chamber on a flow through basis. In some cases, this resulted in fluctuations in the calibration of the swimming chamber, which, as a result was calibrated more frequently than was normally the case. Water temperature in the chamber was recorded for each run, and all were within the range of temperatures in the nearby ocean, or in a few cases, the aquaculture ponds from which the larvae were obtained. The swimming chamber time increment interval was five minutes, with an increase in speed at each increment that varied with the flow from the seawater system and the number of lanes open in the swimming chamber, ranging from 1.6 to 5.3 cm s−1.The larvae studied at Lautan Production, a small company located in Meze, France (42.4°N, 3.6°E) in September 2010 were of two species reared for the aquarium trade. Gramma loreto is native to the western tropical Atlantic, and the brood stock came from Cuba. Pseudochromis fridmani is found only in the Red Sea, but the origin of the brood stock is otherwise unknown. Both species produce ‘egg balls’ that are laid in crevices or small caves and tended by an adult until hatching. The eggs typically hatch at night with little remaining yolk and with no fin-ray development, but with mouth open and eyes pigmented. Recently hatched larvae were removed from the spawning tank to rearing tanks with constant illumination and ‘green water’ at a temperature of 26 °C to 28 °C. Cohort date is for the morning when the larvae were removed from the spawning tank. For the first 5 days rotifers were supplied, and from 6 days after hatch (DAH), the larvae were fed with Artemia nauplii all by Lautan employees. Temperatures in the swimming chamber were similar to those in the rearing tank. Larvae of about 5 mm to settlement size (10–12 mm) were used to measure U-crit. The swimming chamber time increment interval was two minutes, with an increase in speed at each increment of 3.2 cm s−1.Larvae from Taiwan were either preserved in 75% ethanol or in some cases in Bouins Solution for future histology research. Larvae from Lautan were preserved in 75% ethanol. Measurements were made within 24–48 hours after preservation. Body length (BL) was measured on all larvae using a dissecting microscope ocular micrometer: this is Notochord Length (tip of snout to tip of notochord) for preflexion and flexion-stage larvae, and Standard Length for postflexion larvae (tip of snout to end of hypural plate). For some larvae from Taiwan additional measurements were made using Scion Image for Windows (Beta 4.02, Scion Corporation, Frederick, MD): Total Length (tip of snout to tip of posterior-most fin), Total Lateral Area (including fins) and Propulsive Area (Fig. 6), the last as defined by4 (see Fig. 6).Ethic declarationsData collated here are from a large array of studies collected across a range of institutions and locations, and to our knowledge in all cases complied with the required ethics procedures at the relevant institution at the time of data collection. Portions of this work were carried out under Australian Museum Animal Care and Ethics Approval 01/01 (JML) and James Cook Ethics Approvals A202, 402 (RF). In France, research was carried out under permits issued by CNRS to the USR 3278 CNRS/EPHE team to conduct research experiments in the field and laboratory at all locations (under the “Hygiène et Sécurité” section). In Moorea, the research was carried out under permits issued by le Délégué Régional à la recherche et à la technologie de la Polynesie française. More