Ecology

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

1.Balmford, A. & Gaston, K. J. Why biodiversity surveys are good value How the ‘terror crocodile’ grew so big. Nature 398, 204–205 (1999).ADS 
CAS 
Article 

Google Scholar 
2.Rondinini, C., Marco, M. D., Visconti, P., Butchart, S. H. M. & Boitani, L. Update or outdate: Long-term viability of the IUCN Red List Making the Red List financially. Conserv. Lett. 7, 126–130 (2014).Article 

Google Scholar 
3.Braithwaite, M. E. & Walker, K. J. 50 Years of Mapping the British and Irish Flora 1962–2012 (Botanical Society of the British Isles, 2012).
Google Scholar 
4.Barrett, G., Silcocks, A., Barry, S., Cunningham, R. & Poulter, R. The New Atlas of Australian Birds (Royal Australasian Ornithologists Union, 2003).
Google Scholar 
5.Biodiversity Center of Japan. Animal Distribution Atlas of Japan (Ministry of Environment Japan, 2010).
Google Scholar 
6.McGowan, K. & Corwin, K. The Second Atlas of Breeding Birds in New York State (Comstock Publishing Associates, 2008).
Google Scholar 
7.Rhoads, A. F. & Klein, W. M. J. The vascular Flora of Pennsylvenia: Annoted Checklist and Atlas (American Philosophical Society, 1993).
Google Scholar 
8.Baker, H. et al. Population estimates of birds in Great Britain and the United Kingdom British Birds. Br. Birds 99, 25–44 (2006).
Google Scholar 
9.Bonn, A., Rodrigues, A. S. L. & Gaston, K. J. Threatened and endemic species: Are they good indicators of patterns of biodiversity on a national scale?. Ecol. Lett. 5, 733–741 (2002).Article 

Google Scholar 
10.Albuquerque, F. & Beier, P. Rarity-weighted richness: A simple and reliable alternative to integer programming and heuristic algorithms for minimum set and maximum coverage problems in conservation planning. PLoS ONE 10, e0119905 (2015).Article 

Google Scholar 
11.Williams, P. et al. A comparison of richness hotspots, rarity hotspots, and complementary areas for conserving diversity of British birds. Conserv. Biol. 10, 155–174 (1996).Article 

Google Scholar 
12.Platts, P. J. et al. Conservation implications of omitting narrow-ranging taxa from species distribution models, now and in the future. Divers. Distrib. 20, 1307-1320. (2014).Article 

Google Scholar 
13.Kujala, H., Moilanen, A. & Gordon, A. Spatial characteristics of species distributions as drivers in conservation prioritization. Methods Ecol. Evol. 9, 1121–1132 (2018).Article 

Google Scholar 
14.Kukkala, A. S. & Moilanen, A. Core concepts of spatial prioritisation in systematic conservation planning. Biol. Rev. Camb. Philos. Soc. 88, 443–464 (2013).Article 

Google Scholar 
15.Lawler, J. J., White, D., Sifneos, J. C. & Master, L. L. Rare species and the use of indicator groups for conservation planning. Conserv. Biol. 17, 875–882 (2003).Article 

Google Scholar 
16.Kujala, H., Lahoz-Monfort, J. J., Elith, J. & Moilanen, A. Not all data are equal: Influence of data type and amount in spatial conservation prioritisation. Methods Ecol. Evol. 9, 2249–2261 (2018).Article 

Google Scholar 
17.Grantham, H. S. et al. Diminishing return on investment for biodiversity data in conservation planning. Conserv. Lett. 1, 190–198 (2008).Article 

Google Scholar 
18.Margules, C. R. & Pressey, R. L. Systematic conservation planning. Nature 405, 243–253 (2000).CAS 
Article 

Google Scholar 
19.Moilanen, A., Wilson, K. A. & Possingham, H. Spatial Conservation Prioritization: Quantitative Methods and Computational Tools (Oxford University Press, 2009).
Google Scholar 
20.Akasaka, M., Kadoya, T., Ishihama, F., Fujita, T. & Fuller, R. A. Smart protected area placement decelerates biodiversity loss: A representation-extinction feedback leads rare species to extinction. Conserv. Lett. 10, 539–546 (2017).Article 

Google Scholar 
21.Boakes, E., McGowan, P. & Fuller, R. Distorted views of biodiversity: spatial and temporal bias in species occurrence data. PLoS Biol. 8, e1000385 (2010).Article 

Google Scholar 
22.Williams, P. H., Margules, C. R. & Hilbert, D. W. Data requirements and data sources for biodiversity priority area selection. J. Biosci. 27, 327–338 (2002).CAS 
Article 

Google Scholar 
23.Da Fonseca, G. A. B. et al. … following Africa’s lead Community groups could. Nature 405, 393–394 (2000).Article 

Google Scholar 
24.Possingham, H. P., Grantham, H. & Rondinini, C. How can you conserve species that haven’t been found? Commentary. J. Biogeogr. 34, 758–759 (2007).Article 

Google Scholar 
25.Ohlemuller, R. et al. The coincidence of climatic and species rarity: High risk to small-range species from climate change. Biol. Lett. 4, 568–572 (2008).Article 

Google Scholar 
26.Kier, G. et al. A global assessment of endemism and species richness across island and mainland regions. Proc. Natl. Acad. Sci. USA. 106, 9322–9327 (2009).ADS 
CAS 
Article 

Google Scholar 
27.Raes, N., Roos, M. C., Slik, J. W. F., Van Loon, E. E. & Steege, H. T. Botanical richness and endemicity patterns of Borneo derived from species distribution models. Ecography (Cop.) 32, 180–192 (2009).Article 

Google Scholar 
28.Tulloch, A. I. T., Mustin, K., Possingham, H. P., Szabo, J. K. & Wilson, K. A. To boldly go where no volunteer has gone before: predicting volunteer activity to prioritize surveys at the landscape scale. Divers. Distrib. 19, 465–480 (2013).Article 

Google Scholar 
29.Rodewald, A. D., Strimas-Mackey, M., Schuster, R. & Arcese, P. Tradeoffs in the value of biodiversity feature and cost data in conservation prioritization. Sci. Rep. 9, 5921 (2019).Article 

Google Scholar 
30.Rabinowitz, D. Seven forms of rarity. In The biological aspects of rare plant conservation. (ed. Synge, H.) 205–217 (Wiley, Chichester, 1981).
Google Scholar 
31.Ministry of Environment. Red Data Book 2014: Plants I (Gyousei, 2015).
Google Scholar 
32.Ministry of Environment. Red Data Book 2014: Plants II (Gyousei, 2015).
Google Scholar 
33.Yahara, T. et al. Red list of Japanese vascular plants: Summary of methods and results. Proc. Japan. Soc. Plant Taxon. 13, 89–96 (1998).
Google Scholar 
34.Ball, I. R., Possingham, H. P. & Watts, M. Marxan and relatives: Software for spatial conservation prioritisation. In Spatial Conservation Prioritisation: Quantitative Methods and Computational Tools (eds Moilanen, A. et al.) 185–195 (Oxford University Press, 2009).
Google Scholar 
35.Yoshioka, A., Akasaka, M. & Kadoya, T. Spatial prioritization for biodiversity restoration: A simple framework referencing past species distributions. Restor. Ecol. 22, 185–195 (2014).Article 

Google Scholar 
36.Naidoo, R. et al. Integrating economic costs into conservation planning. Trends Ecol. Evol. 21, 681–687 (2006).Article 

Google Scholar 
37.Japan Statistics. Population of Japan: Final Report of the 2005 Population Census. (Statistics Japan, 2010).38.R Core Team. R 4.0.0. (R Foundation for Statistical Computing, 2020).39.Zweig, M. H. & Campbell, G. Receiver-operating characteristic (ROC) plots: A fundamental evaluation tool in clinical medicine. Clin. Chem. 39, 561–577 (1993).CAS 
Article 

Google Scholar 
40.Hermoso, V., Ward, D. P. & Kennard, M. J. Using water residency time to enhance spatio-temporal connectivity for conservation planning in seasonally dynamic freshwater ecosystems. J. Appl. Ecol. 49, 1028–1035 (2012).Article 

Google Scholar 
41.Beyer, H. L., Dujardin, Y., Watts, M. E. & Possingham, H. P. Solving conservation planning problems with integer linear programming. Ecol. Modell. 328, 14–22 (2016).Article 

Google Scholar 
42.Schuster, R., Hanson, J. O., Strimas-Mackey, M. & Bennett, J. R. Integer linear programming outperforms simulated annealing for solving conservation planning problems. PeerJ 8, e9258. https://doi.org/10.7717/peerj.9258 (2020).Article 

Google Scholar  More

Chemicals and reagentsEOs of basil (Ocimum basilicum L.), bergamot (Citrus bergamia Risso & Poit), camphor [Cinnamomum camphora (L.) J. Presl.], cinnamon (Cinnamomum zeylanicum Blume), citronella [Cymbopogon nardus (L.) Rendle], clove (Eugenia caryophyllus Wight), eucalyptus (Eucalyptus globulus Labill.), jasmine (Jasminum officinale L.), lavender (Lavandula angustifolia Mill.), lemon grass [Cymbopogan citratus (DC.) Stapf], mentha (Mentha piperita L.), rosemary (Rosmarinus officinalis L.), patchouli (Pogostemon patchouli Benth), and wild turmeric (Curcuma aromatica Salisb.) were procured from Talent Technologies (Talent Technologies, Kanpur, India). Acetylcholinesterase (AChE) activity assay kit, Anti-OBP2A antibody, ELISA kits, 1,1-diphenyl-2-picrylhydrazyl (DPPH), radioimmunoprecipitation (RIPA) buffer and phosphate buffer saline (PBS) were purchased from Sigma Aldrich (Sigma Aldrich Chemical Co., St. Luis, USA). TRPV1 antibody was purchased from Santa Cruz (Santa Cruz, California, USA). 1-chloro-2,4-dinitrobenzene (CDNB) was purchased from Cayman (Cayman Chemical Company, Michigan, USA). Human normal lung cell line (L-132) was obtained from the National Centre for Cell Sciences (NCCS), Pune, India. High performance liquid chromatography (HPLC) grade acetone was purchased from Merck (Merck Pvt. Ltd., Mumbai, India). All other chemicals used were of the highest analytical grade available.Test insects5–7 days old adult female Ae. albopictus mosquitoes were housed at the laboratory insectary, Division of Pharmaceutical Technology, Defence Research Laboratory, Tezpur, Assam, India. Mosquitoes were reared by maintaining temperature at 27 ± 2 °C, relative humidity: 75 ± 5% RH and 14L:10D h of light–dark alternative cycles in standard-sized wooden cages (75 cm × 60 cm × 60 cm) with a sleeve opening on one side as described previously63. 10% sucrose solution ad libitum were provided for nourishment. Before testing, the mosquitoes were starved for 24 h.Screening of EOsDose response study was performed to evaluate the best oils among the fourteen EOs. This study was approved (approval number: 032/2021TMCH, 28/08/2018) by the Institutional Human Ethical Committee (IHEC), of the Tezpur Medical College & Hospital (TMCH), Tezpur, Assam, India, and all experiments were performed in accordance with relevant guidelines and regulations. Five volunteers are chosen, not allergic to mosquito bite and all volunteers provided written informed consent. A volunteer’s thigh was marked according to the door opening hole of the K&D module as described by Klun and Debboun64. It is made of Plexiglas and the base of the rectangular cage (26 cm × 5 cm × 5 cm) has six holes, each with rectangular 3 × 4 cm holes that are opened and closed by a sliding door (Supplementary Fig. S8: Provide the photograph of K&D module). The flexor region of the forearms of a human volunteer was outlined with four rectangular (3 cm × 4 cm) test areas. A volume of 25 µL of each concentration of the EOs in soybean oil (40, 4 and 0.4 µg/cm2) and 25 µL of the soybean oil (diluent) as control was applied to the marked areas. After air drying for 5 min, a K&D module with matching cut outs in its floor was placed over the treated areas, containing five nulliparous 5–7 days old female mosquitoes in each hole. The doors of the cells were opened and the number of mosquitoes biting in each cell was recorded within a 2 min exposure, after which the doors were closed. After completion of each observation, mosquitoes were freed by opening cells of the K&D module in a sleeved screened cage. For each test, fresh sets of mosquitoes are used. Five replications for each test were carried out. The efficacy of EOs were determined by the percentage repellency against mosquitoes, using the formula or Eq. (2) described by WHO46.$$% ;{text{repellency}} = frac{C – T}{C} times 100$$
(2)
where, C is the number of mosquitoes landing, or biting at the control area; T is the number of mosquitoes landing or biting at the treated area.Fourier transform-infra red spectroscopy (FT-IR)Study of chemical compatibility for each formulation ingredients are necessary. All formulation ingredients possess specific value of vibrational frequency and have varied functional groups in their chemical structures. For compatibility study, each EOs, excipients to be used in cream formulation, and their physical mixture was placed one by one over the sample plate of the FT-IR instrument (Bruker, ALPHA, Billerica, MA, USA). The covering probe was placed over the sample and IR spectra was obtained over a wavelength of 2.5–25 μm at room temperature. Functional groups possessed by each individual ingredient should be identical in their physical mixture which confirms their compatibility37.Thermogravimetric analysis (TGA)The thermal behaviour of citronella oil, clove oil, lemon grass oil, their mixture and EO-MRC were evaluated using a thermal analyser (TG 209 F1 Libra®, NETZSCH-Gerätebau GmbH, 95100 Selb, Germany). Approximately about 10 mg sample weight was placed in the crucible each time. Nitrogen was used as a shielding gas. Heating program was fixed as 30–600 °C at a rate of 10 °C/min.Formulation development and optimizationFor optimization, a 17-run, 3-factor, 3-level Box-Behnken design (BBD) was utilized. A second order polynomial model was constructed by quadratic response surface methodology (RSM) using Design-Expert software (Version 6.0.8, Stat-Ease Inc., USA). Total seventeen formulations were obtained using EO concentrations as dependent variables against complete protection time (CPT) as independent variable or response variable. Analysis of variance (ANOVA) was performed using the same software to obtain the most effective formulation.Preparation of creamPhase inversion temperature method was applied for the preparation of EO-based mosquito repellent cream (EO-MRC). About 50 g cream sample was prepared in order to get enough for performing the various qualitative and quantitative assay. The oil phase (phase B) was prepared by dissolving the oil soluble excipients, except phase A (mosquito repellent active ingredients) under mild heating at 200 rpm in a hot magnetic plate stirrer (Magnetic Stirrer IKA RCT basic) and heated to 65 °C. The aqueous phase was prepared by mixing various aqueous soluble ingredients (phase C) under gentle heating and stirring. Temperature of the aqueous phase was raised to 65 °C. Phase A was gently added to the oil phase at a stirring speed of 200 rpm and 55 ± 2 °C. The mixture was then emulsified by adding phase C slowly and kept for 1 h at a stirring rate of 800 rpm and 60 ± 2 °C. The formulated EO-MRC was then kept for natural cooling.Efficacy assessmentCPT of the developed cream (EO-MRC) formulation was carried out by arm in cage bioassay. 1 mL EO-MRC was applied to ≈ 600 cm2 area of the forearm skin between the wrist and elbow and 1 mL of the 12% N, N-di ethyl benzamide (DEBA) based marketed cream (DBMC) was compared on the other arm. Two mosquito cages (size: 40 × 40 × 40 cm) each containing 200–250 non-blood-fed female Ae. Albopictus were used. One cage is designated for testing the EO-MRC and the other for the positive control (DBMC). During testing, hands were protected by surgical gloves for which the mosquitoes cannot bite while the volunteer avoids movement of the arm. EO-MRC and DBMC treated arms were exposed for 3 min at 30 min intervals to determine landing and/or probing activity. A single landing or probing of mosquito within a 3 min test interval concludes the test. CPT was calculated as the time (min) required for the first mosquito landing or probing after repellent application to the treated area. The median CPT and confidence intervals were estimated from the Kaplan–Meier Survival Function46.Efficacy was correlated with DEBA based marketed cream (DBMC). The inclusion of the specific commercial product DBMC is for comparison and does not constitute any recommendations.CharacterizationGas chromatography-mass spectroscopy (GC–MS)Qualitative studyDifferent chemical components in fourteen EOs and the selected blend were identified by a GC–MS system of Agilent Technologies (5301 Stevens Creek Blvd. Santa Clara, CA 95051, United States). Test sample concentration of 500 μg/mL was prepared in GC grade acetone. A sample volume of 1 μL was introduced into the injector held at 250 °C. Oven temperature of 40–300 °C was programmed at 20 °C/min. Helium was used as carrier gas at flow rate 1 mL/min. The injector and detector temperature were set at 250 °C and 230 °C (quad) and 150 °C (core) respectively37. Standard C7–C30 saturated alkanes were purchased from Sigma Aldrich Chemicals Co., St. Louis, USA. Retention indices (RI) of the identified components were determined for identification of the detected components.% Assay by GC–MS studyCalibration samples of eugenol and citronellol were prepared by dissolving an appropriate amount in GC grade acetone to get concentrations of 62.5 μg/mL, 125 μg/mL, 250 μg/mL and 500 μg/mL. Test samples of EO-MRC, clove oil and citronella oil were prepared by dissolving a required amount in acetone to quantify the EO components in the final formulation. A sample volume of 1 μL was introduced into the injector as described in ‘Qualitative study’ section.Physicochemical parametersPhysical parameters of the EO-MRC and placebo formulations were determined in order to establish aesthetic compliance and consumer acceptability. To determine the viscosity, a programmable viscometer was used (Model: DV2T, Ametek Brookfield, Middleboro, MA, USA); combined with software Rheo3000, version 1.2.2019.1 [R]. Sample volume was fixed at 30 g and viscosities were determined at 10 rpm for 40 s at room temperature using a T-Bar spindle (B-92) (Helipath spindle set, Brookfield Engineering Labs. Inc). Density was determined by using a pycnometer. pH of EO-MRC was checked by using digital pH meter (Labman Scientific instruments, Tamil Nadu, India).Spread ability of EO-MRC was determined as per the method reported earlier by Sabale65. In brief, 1 g of EO-MRC was placed on 1 cm2 pre-marked circular area on the glass slide (7.5 cm × 2.5 cm). EO-MRC was compressed using another glass slide placed from edge to centre of primary slide. 200 g of commercial weight was placed on the set up and allowed the gel to spread for the period of 1 min. The spread diameter was calculated with the aid of graph paper and spread ability was evaluated using formula expressed as Eq. (3):$$mathrm{Spread, ability}=mathrm{m}times frac{mathrm{l}}{mathrm{t}}$$
(3)
where, m is the commercial weight placed on the setup; l is the length of cream spread; and t is the time.Safety assessmentCytotoxicity by MTT assayThe reduction of tetrazolium salts is now widely accepted as a reliable way to examine cell proliferation. The yellow tetrazolium MTT (3-(4,5-dimethylthiazolyl-2)-2,5-diphenyltetrazolium bromide) is reduced by metabolically active cells, in part by the action of dehydrogenase enzymes, to generate reducing equivalents such as NADH and NADPH. With the help of spectrophotometric means, the resulting intracellular purple formazan can be quantified. The assay measures the cell proliferation rate and conversely, when metabolic events cause apoptosis or necrosis, the reduction in cell viability66.Cells cultured in T-25 flasks were trypsinized and aspirated into a 5 mL centrifuge tube. Cell pellet was obtained by centrifugation at 3000 rpm. The cell count was adjusted, using DMEM HG medium, such that 200 μL of suspension contained approximately 10,000 cells. To each well of the 96 well microtiter plate, 200 μL of the cell suspension was added and the plate was incubated at 37 ℃ and 5% CO2 atmosphere for 24 h. After 24 h, the spent medium was aspirated. 200 μL of different test concentrations viz. 62 µg/mL, 125 µg/mL, 250 µg/mL, 500 µg/mL, and 1000 µg/mL, of EO-MRC were added to the respective wells. The plate was then incubated at 37 °C and 5% CO2 atmosphere for 24 h. The plate was removed from the incubator and the drug containing media was aspirated. 200 μL of medium containing10% MTT reagent was then added to each well to get a final concentration of 0.5 mg/mL and the plate was incubated at 37 ℃ and 5% CO2 atmosphere for 3 h. Without disturbing the crystals formed in the wells, culture medium was completely removed. 100 μL of solubilisation solution (DMSO) was added to each well and the plate was then gently shake in a rocking shaker (ROCKYMAX™, Tarsons, Kolkata, India) to solubilize the formed formazan. The absorbance was measured at a wavelength of 570 nm and also at 630 nm using a microplate reader. The percentage growth inhibition was calculated and concentration of EO-MRC needed to inhibit cell growth by 50% (IC50) was generated from the dose–response curve for the cell line.Animals and ethics statementAll experimenting protocols using animal were performed according to the “Principles of Laboratory Animal care” (NIH publication 85–23, revised 1985) and approved by the Institutional Animal Ethical Committee (IAEC) of Defence Research Laboratory (DRL), Tezpur, Assam, India (approval no. CPCSEA/DRL/Protocol no. 3, 20/06/2018). All studies involving animals are reported in accordance with the ARRIVE guidelines for reporting experiments involving animals67. All efforts were made during the study period to minimize the suffering of animals and to reduce the number of animals used.5–8 weeks old, about 210–250 g of male healthy adult Wistar rats (Rattus norvegicus) and young and healthy New Zealand albino rabbits (Oryctolagus cuniculus) were obtained from the institutional animal housing facility and allowed to acclimatize for 7 days prior to the study. Standard food and purified water ad libitum were provided in clean and hygienic condition at 22–25 ℃, 40–70% RH with 12 h light–dark cycles.Acute dermal irritation studyAcute dermal irritation study was conducted on healthy New Zealand albino rabbits following the OECD test guidelines 40468. Approximately 24 h before the test, fur was removed from the dorsal area of the trunk. 0.5 g EO-MRC, was directly applied to the skin and after 4 h exposure period, residual EO-MRC was removed by using water without disturbing the integrity of the epidermis and examined for signs of erythema and oedema, at 60 min, and then at 24 h, 48 h and 72 h after EO-MRC removal. Dermal reactions are graded and recorded according to the grades in the Table 8. As per the method described by Banerjee et al.69; primary irritation index (PII) was calculated. Further, we have followed the Draize method of classification for PII scoring as non-irritant (if PII  More

1.Ceballos, G., Ehrlich, P. R. & Raven, P. H. Vertebrates on the brink as indicators of biological annihilation and the sixth mass extinction. Proc. Natl Acad. Sci. USA 117, 13596–13602 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
2.Committee, T. I. R. L. The IUCN Red List of Threatened Species – Strategic Plan 2017-2020. Tech. Rep., IUCN (2017).3.Witmer, G. W. Wildlife population monitoring: some practical considerations. Wild. Res. 32, 259–263 (2005).
Google Scholar 
4.McEvoy, J. F., Hall, G. P. & McDonald, P. G. Evaluation of unmanned aerial vehicle shape, flight path and camera type for waterfowl surveys: disturbance effects and species recognition. PeerJ 4, e1831 (2016).PubMed 
PubMed Central 

Google Scholar 
5.Burghardt, G. M. et al. Perspectives–minimizing observer bias in behavioral studies: a review and recommendations. Ethology 118, 511–517 (2012).
Google Scholar 
6.Giese, M. Effects of human activity on Adelie penguin Pygoscelis adeliae breeding success. Biol. Conserv. 75, 157–164 (1996).
Google Scholar 
7.Köndgen, S. et al. Pandemic human viruses cause decline of endangered great apes. Curr. Biol. 18, 260–264 (2008).PubMed 

Google Scholar 
8.Weissensteiner, M. H., Poelstra, J. W. & Wolf, J. B. W. Low-budget ready-to-fly unmanned aerial vehicles: an effective tool for evaluating the nesting status of canopy-breeding bird species. J. Avian Biol. 46, 425–430 (2015).
Google Scholar 
9.Sasse, D. B. Job-related mortality of wildlife workers in the united states, 1937–2000. Wildl. Soc. Bull. 31, 1015–1020 (2003).10.Kays, R., Crofoot, M. C., Jetz, W. & Wikelski, M. Terrestrial animal tracking as an eye on life and planet. Science 348, aaa2478 (2015).11.Altmann, J. Observational study of behavior: sampling methods. Behaviour 49, 227–266 (1974).CAS 
PubMed 

Google Scholar 
12.Hodgson, J. C. et al. Drones count wildlife more accurately and precisely than humans. Methods Ecol. Evolution 9, 1160–1167 (2018).
Google Scholar 
13.Betke, M. et al. Thermal imaging reveals significantly smaller Brazilian free-tailed bat colonies than previously estimated. J. Mammal. 89, 18–24 (2008).
Google Scholar 
14.Rollinson, C. R. et al. Working across space and time: nonstationarity in ecological research and application. Front. Ecol. Environ. 19, 66–72 (2021).
Google Scholar 
15.Junker, J. et al. A severe lack of evidence limits effective conservation of the world’s primates. BioScience 70, 794–803 (2020).PubMed 
PubMed Central 

Google Scholar 
16.Sherman, J., Ancrenaz, M. & Meijaard, E. Shifting apes: Conservation and welfare outcomes of Bornean orangutan rescue and release in Kalimantan, Indonesia. J. Nat. Conserv. 55, 125807 (2020).
Google Scholar 
17.O’Donoghue, P. & Rutz, C. Real-time anti-poaching tags could help prevent imminent species extinctions. J. Appl. Ecol. 53, 5–10 (2016).PubMed 

Google Scholar 
18.Lahoz-Monfort, J. J. & Magrath, M. J. L. A comprehensive overview of technologies for species and habitat monitoring and conservation. BioScience biab073. https://academic.oup.com/bioscience/advance-article/doi/10.1093/biosci/biab073/6322306 (2021).19.Gottschalk, T., Huettmann, F. & Ehlers, M. Thirty years of analysing and modelling avian habitat relationships using satellite imagery data: a review. Int. J. Remote Sens. 26, 2631–2656 (2005).
Google Scholar 
20.Steenweg, R. et al. Scaling-up camera traps: monitoring the planet’s biodiversity with networks of remote sensors. Front. Ecol. Environ. 15, 26–34 (2017).
Google Scholar 
21.Hausmann, A. et al. Social media data can be used to understand tourists’ preferences for nature-based experiences in protected areas. Conserv. Lett. 11, e12343 (2018).
Google Scholar 
22.Sugai, L. S. M., Silva, T. S. F., Ribeiro, J. W. & Llusia, D. Terrestrial passive acoustic monitoring: review and perspectives. BioScience 69, 15–25 (2018).
Google Scholar 
23.Wikelski, M. et al. Going wild: what a global small-animal tracking system could do for experimental biologists. J. Exp. Biol. 210, 181–186 (2007).PubMed 

Google Scholar 
24.Belyaev, M. Y. et al. Development of technology for monitoring animal migration on Earth using scientific equipment on the ISS RS. in 2020 27th Saint Petersburg International Conference on Integrated Navigation Systems (ICINS), 1–7 (IEEE, 2020).25.Harel, R., Loftus, J. C. & Crofoot, M. C. Locomotor compromises maintain group cohesion in baboon troops on the move. Proc. R. Soc. B 288, 20210839 (2021).PubMed 
PubMed Central 

Google Scholar 
26.Farley, S. S., Dawson, A., Goring, S. J. & Williams, J. W. Situating ecology as a big-data science: current advances, challenges, and solutions. BioScience 68, 563–576 (2018).
Google Scholar 
27.Lasky, M. et al. Candid critters: Challenges and solutions in a large-scale citizen science camera trap project. Citizen Science: Theory and Practice 6, https://doi.org/10.5334/cstp.343 (2021).28.Hastie, T., Tibshirani, R. & Friedman, J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction (Springer, 2001).29.Christin, S., Hervet, É. & Lecomte, N. Applications for deep learning in ecology. Methods Ecol. Evolution 10, 1632–1644 (2019).
Google Scholar 
30.Kwok, R. Ai empowers conservation biology. Nature 567, 133–135 (2019).ADS 
CAS 
PubMed 

Google Scholar 
31.Kwok, R. Deep learning powers a motion-tracking revolution. Nature 574, 137–139 (2019).ADS 
CAS 
PubMed 

Google Scholar 
32.LeCun, Y., Bengio, Y. & Hinton, G. Deep learning. Nature 521, 436–444 (2015).ADS 
CAS 
PubMed 

Google Scholar 
33.Pichler, M., Boreux, V., Klein, A.-M., Schleuning, M. & Hartig, F. Machine learning algorithms to infer trait-matching and predict species interactions in ecological networks. Methods Ecol. Evolution 11, 281–293 (2020).
Google Scholar 
34.Knudby, A., LeDrew, E. & Brenning, A. Predictive mapping of reef fish species richness, diversity and biomass in Zanzibar using IKONOS imagery and machine-learning techniques. Remote Sens. Environ. 114, 1230–1241 (2010).ADS 

Google Scholar 
35.Rey, N., Volpi, M., Joost, S. & Tuia, D. Detecting animals in African savanna with UAVs and the crowds. Remote Sens. Environ. 200, 341–351 (2017).ADS 

Google Scholar 
36.Beery, S., Morris, D. & Yang, S. Efficient pipeline for camera trap image review. in Proceedings of the Workshop Data Mining and AI for Conservation, Conference for Knowledge Discovery and Data Mining (2019).37.Kellenberger, B., Marcos, D. & Tuia, D. When a few clicks make all the difference: improving weakly-supervised wildlife detection in UAV images. in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (2019).38.Schofield, D. et al. Chimpanzee face recognition from videos in the wild using deep learning. Sci. Adv. 5, eaaw0736 (2019).ADS 
PubMed 
PubMed Central 

Google Scholar 
39.Ditria, E. M. et al. Automating the analysis of fish abundance using object detection: optimizing animal ecology with deep learning. Front. Mar. Sci. 7, 429 (2020).
Google Scholar 
40.Kellenberger, B., Veen, T., Folmer, E. & Tuia, D. 21 000 birds in 4.5 h: efficient large-scale seabird detection with machine learning. Remote Sens. Ecol. Conserv. 7, 445–460 (2021).
Google Scholar 
41.Ahumada, J. A. et al. Wildlife insights: a platform to maximize the potential of camera trap and other passive sensor wildlife data for the planet. Environ. Conserv. 47, 1–6 (2020).MathSciNet 

Google Scholar 
42.Eikelboom, J. A. J. et al. Improving the precision and accuracy of animal population estimates with aerial image object detection. Methods Ecol. Evolution 10, 1875–1887 (2019).
Google Scholar 
43.Weinstein, B. G. A computer vision for animal ecology. J. Anim. Ecol. 87, 533–545 (2018).PubMed 

Google Scholar 
44.Valletta, J. J., Torney, C., Kings, M., Thornton, A. & Madden, J. Applications of machine learning in animal behaviour studies. Anim. Behav. 124, 203–220 (2017).
Google Scholar 
45.Peters, D. P. C. et al. Harnessing the power of big data: infusing the scientific method with machine learning to transform ecology. Ecosphere 5, art67 (2014).
Google Scholar 
46.Yu, Q. et al. Study becomes insight: ecological learning from machine learning. Methods Ecol. Evol. 12, 2117–2128 (2021).47.Lucas, T. C. D. A translucent box: interpretable machine learning in ecology. Ecol. Monogr. 90, https://doi.org/10.1002/ecm.1422 (2020).48.Reichstein, M. et al. Deep learning and process understanding for data-driven Earth system science. Nature 566, 195–204 (2019).ADS 
CAS 
PubMed 

Google Scholar 
49.Camps-Valls, G., Tuia, D., Zhu, X. X. & Reichstein, M. Deep Learning for the Earth Sciences: A Comprehensive Approach to Remote Sensing, Climate Science and Geosciences (Wiley & Sons, 2021).50.Karpatne, A. et al. Theory-guided data science: A new paradigm for scientific discovery from data. IEEE Trans. Knowl. Data Eng. 29, 2318–2331 (2017).
Google Scholar 
51.Oliver, R. Y., Meyer, C., Ranipeta, A., Winner, K. & Jetz, W. Global and national trends, gaps, and opportunities in documenting and monitoring species distributions. PLoS Biol 19, e3001336 https://doi.org/10.1371/journal.pbio.3001336 (2021).52.Beery, S., Wu, G., Rathod, V., Votel, R. & Huang, J. Context R-CNN: long term temporal context for per-camera object detection. in 2020 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 13075–13085 (2020).53.Norouzzadeh, M. S. et al. Automatically identifying, counting, and describing wild animals in camera-trap images with deep learning. Proc. Natl Acad. Sci. USA 115, E5716–E5725 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
54.Schneider, S., Taylor, G. W., Linquist, S. & Kremer, S. C. Past, present and future approaches using computer vision for animal re-identification from camera trap data. Methods Ecol. Evolution 10, 461–470 (2019).
Google Scholar 
55.Beery, S., Van Horn, G. & Perona, P. Recognition in terra incognita. in 2018 European Conference on Computer Vision (ECCV), 456–473 (2018).56.Sugai, L. S. M., Silva, T. S. F., Ribeiro Jr, J. W. & Llusia, D. Terrestrial passive acoustic monitoring: review and perspectives. BioScience 69, 15–25 (2019).
Google Scholar 
57.Wrege, P. H., Rowland, E. D., Keen, S. & Shiu, Y. Acoustic monitoring for conservation in tropical forests: examples from forest elephants. Methods Ecol. Evolution 8, 1292–1301 (2017).
Google Scholar 
58.Desjonquères, C., Gifford, T. & Linke, S. Passive acoustic monitoring as a potential tool to survey animal and ecosystem processes in freshwater environments. Freshw. Biol. 65, 7–19 (2020).
Google Scholar 
59.Davis, G. E. et al. Long-term passive acoustic recordings track the changing distribution of North Atlantic right whales (eubalaena glacialis) from 2004 to 2014. Sci. Rep. 7, 1–12 (2017).
Google Scholar 
60.Wood, C. M. et al. Detecting small changes in populations at landscape scales: a bioacoustic site-occupancy framework. Ecol. Indic. 98, 492–507 (2019).
Google Scholar 
61.Kahl, S., Wood, C. M., Eibl, M. & Klinck, H. Birdnet: a deep learning solution for avian diversity monitoring. Ecol. Inform. 61, 101236 (2021).
Google Scholar 
62.Stowell, D., Wood, M. D., Pamuła, H., Stylianou, Y. & Glotin, H. Automatic acoustic detection of birds through deep learning: the first bird audio detection challenge. Methods Ecol. Evolution 10, 368–380 (2019).
Google Scholar 
63.Ford, J. K. B. in Encyclopedia of Marine Mammals 253–254 (Elsevier, 2018).64.Hughey, L. F., Hein, A. M., Strandburg-Peshkin, A. & Jensen, F. H. Challenges and solutions for studying collective animal behaviour in the wild. Philos. Trans. R. Soc. B: Biol. Sci. 373, 20170005 (2018).
Google Scholar 
65.Williams, H. J. et al. Optimizing the use of biologgers for movement ecology research. J. Anim. Ecol. 89, 186–206 (2020).PubMed 

Google Scholar 
66.Korpela, J. et al. Machine learning enables improved runtime and precision for bio-loggers on seabirds. Commun. Biol. 3, 1–9 (2020).
Google Scholar 
67.Yu, H. An evaluation of machine learning classifiers for next-generation, continuous-ethogram smart trackers. Mov. Ecol. 9, 14 (2021).
Google Scholar 
68.Browning, E. et al. Predicting animal behaviour using deep learning: GPS data alone accurately predict diving in seabirds. Methods Ecol. Evolution 9, 681–692 (2018).
Google Scholar 
69.Liu, Z. Y.-C. et al. Deep learning accurately predicts white shark locomotor activity from depth data. Anim. Biotelemetry 7, 1–13 (2019).
Google Scholar 
70.Wang, G. Machine learning for inferring animal behavior from location and movement data. Ecol. Inform. 49, 69–76 (2019).
Google Scholar 
71.Wijeyakulasuriya, D. A., Eisenhauer, E. W., Shaby, B. A. & Hanks, E. M. Machine learning for modeling animal movement. PLoS ONE 30, e0235750 (2020).72.Linchant, J., Lisein, J., Semeki, J., Lejeune, P. & Vermeulen, C. Are unmanned aircraft systems (UASs) the future of wildlife monitoring? A review of accomplishments and challenges. Mammal. Rev. 45, 239–252 (2015).
Google Scholar 
73.Hodgson, J. C., Baylis, S. M., Mott, R., Herrod, A. & Clarke, R. H. Precision wildlife monitoring using unmanned aerial vehicles. Sci. Rep. 6, 1–7 (2016).
Google Scholar 
74.Mathis, A. et al. DeepLabCut: markerless pose estimation of user-defined body parts with deep learning. Nat. Neurosci. 21, 1281–1289 (2018).CAS 
PubMed 

Google Scholar 
75.Graving, J. M. et al. DeepPoseKit, a software toolkit for fast and robust animal pose estimation using deep learning. Elife 8, e47994 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
76.Mathis, A., Schneider, S., Lauer, J. & Mathis, M. W. A primer on motion capture with deep learning: principles, pitfalls, and perspectives. Neuron 108, 44–65 (2020).CAS 
PubMed 

Google Scholar 
77.Kellenberger, B., Marcos, D. & Tuia, D. Detecting mammals in UAV images: best practices to address a substantially imbalanced dataset with deep learning. Remote Sens. Environ. 216, 139–153 (2018).ADS 

Google Scholar 
78.Kellenberger, B., Veen., T., Folmer, E. & Tuia, D. 21,000 birds in 4.5 hours: efficient large-scale seabird detection with machine learning. Remote Sens. Ecol. Conserv. https://doi.org/10.1002/rse2.200 (2021).79.Andrew, W., Greatwood, C. & Burghardt, T. Aerial animal biometrics: individual Friesian cattle recovery and visual identification via an autonomous UAV with onboard deep inference. in International Conference on Intelligent Robots and Systems (IROS) (2019).80.Schroeder, N. M., Panebianco, A., Gonzalez Musso, R. & Carmanchahi, P. An experimental approach to evaluate the potential of drones in terrestrial mammal research: a gregarious ungulate as a study model. R. Soc. open Sci. 7, 191482 (2020).ADS 
PubMed 
PubMed Central 

Google Scholar 
81.Bennitt, E., Bartlam-Brooks, H. L. A., Hubel, T. Y. & Wilson, A. M. Terrestrial mammalian wildlife responses to Unmanned Aerial Systems approaches. Sci. Rep. 9, 1–10 (2019).CAS 

Google Scholar 
82.Deneu, B., Servajean, M., Botella, C. & Joly, A. Evaluation of deep species distribution models using environment and co-occurrences. in International Conference of the Cross-Language Evaluation Forum for European Languages, 213–225 (Springer, 2019).83.Zhu, X. et al. Deep learning in remote sensing: A comprehensive review and list of resources. IEEE Geosci. Remote Sens. Mag. 5, 8–36 (2017).
Google Scholar 
84.Guirado, E., Tabik, S., Rivas, M. L., Alcaraz-Segura, D. & Herrera, F. Whale counting in satellite and aerial images with deep learning. Sci. Rep. 9, 1–12 (2019).CAS 

Google Scholar 
85.Duporge, I., Isupova, O., Reece, S., Macdonald, D. W. & Wang, T. Using very-high-resolution satellite imagery and deep learning to detect and count African elephants in heterogeneous landscapes. Remote Sens. Ecol. Conserv. https://doi.org/10.1002/rse2.195 (2020).86.Fretwell, P. T. & Trathan, P. N. Discovery of new colonies by Sentinel2 reveals good and bad news for emperor penguins. Remote Sens. Ecol. Conserv. https://doi.org/10.1002/rse2.176 (2020).87.Brodrick, P. G., Davies, A. B. & Asner, G. P. Uncovering ecological patterns with convolutional neural networks. Trends Ecol. Evolution 34, 734–745 (2019).
Google Scholar 
88.Audebert, N., Le Saux, B. & Lefèvre, S. Deep learning for classification of hyperspectral data: a comparative review. IEEE Geosci. Remote Sens. Mag. 7, 159–173 (2019).
Google Scholar 
89.McKinley, D. C. et al. Citizen science can improve conservation science, natural resource management, and environmental protection. Biol. Conserv. 208, 15–28 (2017).
Google Scholar 
90.Wäldchen, J. & Mäder, P. Machine learning for image based species identification. Methods Ecol. Evolution 9, 2216–2225 (2018).MATH 

Google Scholar 
91.Torney, C. J. et al. A comparison of deep learning and citizen science techniques for counting wildlife in aerial survey images. Methods Ecol. Evolution 10, 779–787 (2019).
Google Scholar 
92.Parham, J., Crall, J., Stewart, C., Berger-Wolf, T. & Rubenstein, D. I. Animal population censusing at scale with citizen science and photographic identification. in AAAI Spring Symposium-Technical Report (2017).93.Kühl, H. S. & Burghardt, T. Animal biometrics: quantifying and detecting phenotypic appearance. Trends Ecol. Evolution 28, 432–441 (2013).
Google Scholar 
94.Yu, X. et al. Automated identification of animal species in camera trap images. EURASIP J. Image Video Process. 2013, 1–10 (2013).ADS 

Google Scholar 
95.Mac Aodha, O. et al. Bat detective–deep learning tools for bat acoustic signal detection. PLoS Computat. Biol. 14, e1005995 (2018).
Google Scholar 
96.Schindler, F. & Steinhage, V. Identification of animals and recognition of their actions in wildlife videos using deep learning techniques. Ecol. Inform. 61, 101215 (2021).97.Avise, J. C. Molecular Markers, Natural History and Evolution (Springer Science & Business Media, 2012).98.Vidal, M., Wolf, N., Rosenberg, B., Harris, B. P. & Mathis, A. Perspectives on Individual Animal Identification from Biology and Computer Vision. Integr. Comp. Biol. 61, 900–916 https://doi.org/10.1093/icb/icab107 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
99.Berger-Wolf, T. Y. et al. Wildbook: crowdsourcing, computer vision, and data science for conservation. Preprint at https://arxiv.org/abs/1710.08880 (2017).100.Parham, J. et al. An animal detection pipeline for identification. in IEEE Winter Conference on Applications of Computer Vision (WACV), 1075–1083 (IEEE, 2018).101.Weideman, H. et al. Extracting identifying contours for African elephants and humpback whales using a learned appearance model. in Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV) (2020).102.Brust, C.-A. et al. Towards automated visual monitoring of individual gorillas in the wild. in 2017 IEEE International Conference on Computer Vision Workshops (ICCVW), 2820–2830 (2017).103.Li, S., Li, J., Tang, H., Qian, R. & Lin, W. ATRW: a benchmark for Amur tiger re-identification in the wild. in 2020 ACM International Conference on Multimedia, 2590–2598 (2020).104.Bendale, A. & Boult, T. E. Towards open set deep networks. in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 1563–1572 (2016).105.Mathis, M. W. & Mathis, A. Deep learning tools for the measurement of animal behavior in neuroscience. Curr. Opin. Neurobiol. 60, 1–11 (2020).CAS 
PubMed 

Google Scholar 
106.Sanakoyeu, A., Khalidov, V., McCarthy, M. S., Vedaldi, A. & Neverova, N. Transferring dense pose to proximal animal classes. in 2020 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 5233–5242 (2020).107.Zuffi, S., Kanazawa, A., Jacobs, D. W. & Black, M. J. 3D menagerie: modeling the 3D shape and pose of animals. in 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 6365–6373 (2017).108.Biggs, B., Roddick, T., Fitzgibbon, A. & Cipolla, R. Creatures great and smal: recovering the shape and motion of animals from video. in 2018 Asian Conference on Computer Vision (ACCV), 3–19 (Springer, 2018).109.Biggs, B., Boyne, O., Charles, J., Fitzgibbon, A. & Cipolla, R. Who left the dogs out? 3D animal reconstruction with expectation maximization in the loop. in 2020 European Conference on Computer Vision (ECCV), 195–211 (Springer, 2020).110.Zuffi, S., Kanazawa, A., Berger-Wolf, T. & Black, M. J. Three-D safari: learning to estimate zebra pose, shape, and texture from images” in the wild”. in 2019 IEEE International Conference on Computer Vision (ICCV), 5359–5368 (2019).111.Wang, Y., Kolotouros, N., Daniilidis, K. & Badger, M. Birds of a feather: capturing avian shape models from images. in 2021 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 14739–14749 (2021).112.Haalck, L., Mangan, M., Webb, B. & Risse, B. Towards image-based animal tracking in natural environments using a freely moving camera. J. Neurosci. methods 330, 108455 (2020).PubMed 

Google Scholar 
113.Pettorelli, N. et al. Satellite remote sensing for applied ecologists: opportunities and challenges. J. Appl. Ecol. 51, 839–848 (2014).
Google Scholar 
114.Davies, A. B., Tambling, C. J., Kerley, G. I. H. & Asner, G. P. Effects of vegetation structure on the location of lion kill sites in African thicket. PLoS ONE 11, e0149098 (2016).PubMed 
PubMed Central 

Google Scholar 
115.Froidevaux, J. S. P., Zellweger, F., Bollmann, K., Jones, G. & Obrist, M. K. From field surveys to LiDAR: shining a light on how bats respond to forest structure. Remote Sens. Environ. 175, 242–250 (2016).ADS 

Google Scholar 
116.Risse, B., Mangan, M., Stürzl, W. & Webb, B. Software to convert terrestrial LiDAR scans of natural environments into photorealistic meshes. Environ. Model. Softw. 99, 88–100 (2018).
Google Scholar 
117.Haalck, L. & Risse, B. Embedded dense camera trajectories in multi-video image mosaics by geodesic interpolation-based reintegration. in 2021 IEEE Winter Conference on Applications of Computer Vision (WACV), 1849–1858 (2021).118.Schonberger, J. L. & Frahm, J.-M. Structure-from-motion revisited. in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 4104–4113 (2016).119.Mur-Artal, R. & Tardós, J. D. ORB-SLAM2: an open-source SLAM system for monocular, stereo, and RGB-D cameras. IEEE Trans. Robot. 33, 1255–1262 (2017).
Google Scholar 
120.Kuppala, K., Banda, S. & Barige, T. R. An overview of deep learning methods for image registration with focus on feature-based approaches. Int. J. Image Data Fusion 11, 113–135 (2020).ADS 

Google Scholar 
121.Lisein, J., Linchant, J., Lejeune, P., Bouché, P. & Vermeulen, C. Aerial surveys using an unmanned aerial system (UAS): comparison of different methods for estimating the surface area of sampling strips. Tropical Conserv. Sci. 6, 506–520 (2013).
Google Scholar 
122.Wu, C. Critical configurations for radial distortion self-calibration. in 2014 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 25–32 (2014).123.Ferrer, J., Elibol, A., Delaunoy, O., Gracias, N. & Garcia, R. Large-area photo-mosaics using global alignment and navigation data. in Mts/IEEE Oceans Conference, 1–9 (2007).124.Guisan, A. & Zimmermann, N. E. Predictive habitat distribution models in ecology. Ecol. Model. 135, 147–186 (2000).
Google Scholar 
125.Lehmann, A., Overton, J. M. & Austin, M. P. Regression models for spatial prediction: their role for biodiversity and conservation. Biodivers. Conserv. 11, 2085–2092 (2002).126.Breiman, L. Random forests. Mach. Learn. 45, 5–32 (2001).MATH 

Google Scholar 
127.Parravicini, V. et al. Global patterns and predictors of tropical reef fish species richness. Ecography 36, 1254–1262 (2013).
Google Scholar 
128.Smoliński, S. & Radtke, K. Spatial prediction of demersal fish diversity in the baltic sea: comparison of machine learning and regression-based techniques. ICES J. Mar. Sci. 74, 102–111 (2017).
Google Scholar 
129.Čandek, K., Čandek, U. P. & Kuntner, M. Machine learning approaches identify male body size as the most accurate predictor of species richness. BMC Biol. 18, 1–16 (2020).
Google Scholar 
130.Baltensperger, A. P. & Huettmann, F. Predictive spatial niche and biodiversity hotspot models for small mammal communities in Alaska: applying machine-learning to conservation planning. Landscape Ecol. 30, 681–697 (2015).131.Faisal, A., Dondelinger, F., Husmeier, D. & Beale, C. M. Inferring species interaction networks from species abundance data: a comparative evaluation of various statistical and machine learning methods. Ecol. Inform. 5, 451–464 (2010).
Google Scholar 
132.Van Horn, G. et al. The inaturalist species classification and detection dataset. in 2018 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 8769–8778 (2018).133.Copas, K. et al. Training machines to improve species identification using GBIF-mediated datasets. in AGU Fall Meeting Abstracts, Vol. 2019, IN53C–0758 (2019).134.Lennox, R. J. et al. A novel framework to protect animal data in a world of ecosurveillance. BioScience 70, 468–476 (2020).
Google Scholar 
135.Strubell, E., Ganesh, A. & McCallum, A. Energy and policy considerations for deep learning in NLP. in Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, 3645–3650 (2019).136.Samek, W., Montavon, G., Vedaldi, A., Hansen, L. K. & Müller, K.-R. Explainable AI: Interpreting, Explaining and Visualizing Deep Learning, Vol. 11700 (Springer Nature, 2019).137.Swanson, A. et al. Snapshot Serengeti, high-frequency annotated camera trap images of 40 mammalian species in an African savanna. Sci. data 2, 1–14 (2015).
Google Scholar 
138.de Lutio, R. et al. Digital taxonomist: identifying plant species in community scientists’ photographs. ISPRS J. Photogramm. Remote Sens. 182, 112–121 (2021).139.Mac Aodha, O., Cole, E. & Perona, P. Presence-only geographical priors for fine-grained image classification. in Proceedings of the IEEE/CVF International Conference on Computer Vision, 9596–9606 (2019).140.Gurumurthy, S. et al. Exploiting Data and Human Knowledge for Predicting Wildlife Poaching. in Proceedings of the 1st ACM SIGCAS Conference on Computing and Sustainable Societies, 1–8, https://doi.org/10.1145/3209811.3209879 (ACM, 2018).141.Datta, S., Anderson, D., Branson, K., Perona, P. & Leifer, A. Computational neuroethology: a call to action. Neuron 104, 11–24 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
142.Joska, D. et al. AcinoSet: a 3D pose estimation dataset and baseline models for Cheetahs in the wild. 2021 IEEE International Conference on Robotics and Automation (ICRA) Preprint at https://arxiv.org/abs/2103.13282 (IEEE, Xi’an, China, 2021).143.Chen, Q. & Koltun, V. Photographic image synthesis with cascaded refinement networks. in 2019 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 1511–1520 (2017).144.Lee, J., Hwangbo, J., Wellhausen, L., Koltun, V. & Hutter, M. Learning quadrupedal locomotion over challenging terrain. Sci. Robot. 5, eabc5986 (2020).145.Botella, C., Joly, A., Bonnet, P., Munoz, F. & Monestiez, P. Jointly estimating spatial sampling effort and habitat suitability for multiple species from opportunistic presence-only data. Methods Ecol. Evolution 12, 933–945 (2021).
Google Scholar 
146.Beery, S., Cole, E., Parker, J., Perona, P. & Winner, K. Species distribution modeling for machine learning practitioners: a review. in Proceedings of the 4th ACM SIGCAS Conference on Computing and Sustainable Societies (2021).147.Arzoumanian, Z., Holmberg, J. & Norman, B. An astronomical pattern-matching algorithm for computer-aided identification of whale sharks Rhincodon typus. J. Appl. Ecol. 42, 999–1011 (2005).
Google Scholar 
148.de Knegt, H. J., Eikelboom, J. A. J., van Langevelde, F., Spruyt, W. F. & Prins, H. H. T. Timely poacher detection and localization using sentinel animal movement. Sci. Rep. 11, 1–11 (2021).
Google Scholar 
149.Walter, T. & Couzin, I. D. TRex, a fast multi-animal tracking system with markerless identification, and 2D estimation of posture and visual fields. eLife 10, e64000 (2021).PubMed 
PubMed Central 

Google Scholar 
150.Kellenberger, B., Tuia, D. & Morris, D. AIDE: accelerating image-based ecological surveys with interactive machine learning. Methods Ecol. Evolution 11, 1716–1727 (2020).
Google Scholar 
151.Settles, B. Active learning. Synth. lectures Artif. Intell. Mach. Learn. 6, 1–114 (2012).MathSciNet 
MATH 

Google Scholar 
152.Ofli, F. et al. Combining human computing and machine learning to make sense of big (aerial) data for disaster response. Big Data 4, 47–59 (2016).PubMed 

Google Scholar 
153.Simpson, R., Page, K. R. & De Roure, D. Zooniverse: observing the world’s largest citizen science platform. in Proceedings of the 23rd International Conference on World Wide Web 1049–1054 (2014).154.Pocock, M. J. O., Roy, H. E., Preston, C. D. & Roy, D. B. The biological records centre: a pioneer of citizen science. Biol. J. Linn. Soc. 115, 475–493 (2015).
Google Scholar  More