Ecology

Study areaAll simulations were run at a 100 × 100 m resolution for the entire European Alps, which cover ~200,000 km². Elevations reach 4,810 m above sea level at the highest peak (Mont Blanc, elevational data were obtained from ref. 44). Mean annual temperature ranges from approximately −13 up to 16 °C and annual precipitation sums reach up to ~3,600 mm (climatic conditions were obtained from WorldClim45).Species dataTrue presences/absences were derived from complete species lists of 14,040 localized plots covering subalpine and alpine non-forest vegetation of the Alps, compiled from published46 and unpublished data sources. For more information see the supplementary information in ref. 21.Data on demographic rates as well as dispersal parameters were taken from ref. 21, Supplementary Table 1. Detailed values are provided in Supplementary Table 1.Environmental variablesCurrent climate dataMaps of current climatic conditions were generated on the basis of mean, minimum and maximum monthly temperature obtained from WorldClim45 and monthly precipitation sums derived from ref. 47 at a resolutions of 0.5 arcmin and 5 km, respectively. Temperature and precipitation data were downscaled to 100 m as described in ref. 21 and using ordinary kriging with elevation as covariable. As the reference periods of these two datasets do not match (temperature values represent averages for 1950–2000, while precipitation covers 1970–2005) temperature values were adapted using the E-OBS climate grids available online (www.ecad.eu/download/ensembles/download.php). We used these spatially refined temperature and precipitation grids to derive maps of mean annual temperature and mean annual precipitation sum. We chose only two climatic variables to keep models simple and, therefore, interpretation of results more straightforward. In fact, the climatic drivers of population performance and distribution can be more complex48 and vary among species, life history stages and vital rates49. However, as correlations between different descriptors of temperature (such as coldest month or warmest month temperature, Pearson correlation of 0.84) as well as between precipitation variables are high in the European Alps, we chose mean annual temperature and mean annual precipitation sum as they give the most basic description of how climatic conditions change over geographical and elevational gradients.Future climate dataMonthly time series of mean temperature as well as precipitation sums predicted for the twenty-first century were extracted from the Cordex data portal (http://cordex.org). We chose two IPCC5 scenarios from the RCP families representing mild (RCP 2.6) and severe (RCP 8.5) climate change to consider the uncertainty in the future climate predictions. Both scenarios were generated by Météo-France/Centre National de Recherches Météorologiques using the CNRM-ALADIN53 model, fed by output from the global circulation model CNRM-CM5 (ref. 50). The RCP 2.6 scenario assumes that radiative forcing reaches nearly 3 W m−2 (equal to 490 ppm CO2 equivalent) mid-century and will decrease to 2.6 W m−2 by 2100. In the RCP 8.5 scenario, radiative forcing continues to rise throughout the twenty-first century and reaches >8.5 W m−2 (equal to 1,370 ppm CO2 equivalent) in 210024.These time series were statistically downscaled (delta method) by (1) calculating differences (deltas) between monthly temperature and precipitation values hindcasted for the current climatic conditions (mean 1970–2005) and forecasted future values at the original spatial resolution of 11′; (2) spatially interpolating these differences to a resolution of 100 × 100 m using cubic splines and (3) adding them to the downscaled current climate data separately for each climatic variable21,36. Subsequently, we calculated running means (averaged over 35 years) for every tenth year (2030, 2040 through to 2080) for each climatic variable and finally derived, on the basis of the monthly data, mean annual temperature and mean annual precipitation sums for these decadal time steps. The application of 35-yr running means ensures that we use the same time interval for parameterization and prediction. Furthermore, for long-lived species such as alpine plants, running means over long time intervals appear most appropriate to characterize climatic niches33.Soil dataIn addition to the climatic data we used a map of the percentage of calcareous substrate within a cell (5′ longitude × 3′ latitude dissolved to 100 m resolution; further referred to as soil) as described in the supplementary information of ref. 21.Environmental suitability modellingWe parameterized logistic regression models (LRMs) with a logit link function using species presence/absence data as response and the three environmental (two bioclimatic and one soil) variables as predictors. All predictor variables entered the model as second-order polynomials in agreement with the standard unimodal niche concept.From the models, we also derived a threshold value to use for translating predicted probability of occurrence (as a measure of site suitability) into predicted presence or absence of each species at a site (called occurrence threshold, OT, henceforth). The threshold was defined such that it optimized the true skills statistic (TSS), a measure of predictive accuracy derived from comparing observed and predicted presence–absence maps51.Genetic model and niche partitioningSpecies-specific suitability curves for the annual mean temperature gradient derived from the LRMs were partitioned into suitability curves of ecologically different haplotypes of a species as defined by allelic variation in up to three diploid loci (Extended Data Fig. 6). The number of alleles was varied (n = 1, 2, 3, 5 and 10 alleles) as was the proportion of the entire species’ (temperature) niche covered by each haplotype. Models with more than one locus were run with diallelic loci, as otherwise computational efforts would have increased excessively (for each haplotype the number of seeds, juveniles and adults has to be stored and all seeds have to be distributed separately). Each combination of haplotype number and allelic niche size used in a particular simulation is further referred to as setting. Species-specific suitability curves along the other two dimensions (precipitation and soil) remained unpartitioned to ease interpretation of results. The implications of relaxing this assumption by modelling niche partitioning along different environmental gradients are hard to predict. Outcomes would probably depend on the direction and strength of individual specialization along these gradients, whether they are positively or negatively correlated, as well as on how both temperature and precipitation patterns will be affected by climate change. As a consequence, the patterns we found could be re-enforced, for example when the replacement of cold-adapted haplotypes within the species elevational range is further delayed, for example, because those haplotypes adapted to warmer conditions can cope less well with higher precipitation at higher elevations. Vice versa, maladaptation to the warming temperatures might be attenuated, for example, if temperature increase is paralleled by precipitation decrease and emerging drought stress. If, in this case, haplotypes from lower elevations can better cope with both higher temperatures and less water availability than those of median elevations, they may replace the latter faster at these median elevations and hence shorten the phase of maladaptation.Allelic effects were assumed to define the temperature optimum additively. Hence, the heterozygotes’ optimum is always exactly between the optima of the two corresponding homozygotes, corresponding to a codominant genetic model. Finally, all haplotypes corresponding to one setting were assumed to have constant (temperature) niche size, defined as a proportion (ω = 50%, 75% and 100%, for one haplotype only 100%) of the entire species’ (temperature) niche. The temperature niche was computed as the difference between the upper and lower temperature values at which the LRM-derived suitability curve predicted a suitability equal to OT (with precipitation and soil set to the respective optima of the species, also derived from the LRMs). To derive the same geographic distribution under current climatic conditions for each setting, the union of the niches of all haplotypes of one set has to approximate the niche of the single-species model (Extended Data Fig. 6). Note, however, that, the aspired equality of niches is impossible, as the niches resulting from a logistic regression with quadratic terms are always elliptic in shape. Therefore, the optima of the two extreme homozygotes (that is, those carrying haplotypes adapted to the coldest or warmest margins of the entire species’ niche) are fixed at: niche limits ± 1/2 × ω × niche size and all other optima are distributed between them at equal distances (Extended Data Fig. 6 gives a schematic illustration). As a consequence, the performance of the extreme haplotypes, both coldest and warmest, is modelled to be somewhat higher towards the cold and warm margins of the temperature niche, respectively, compared with the performance modelled for the species without intraspecific niche partitioning (compare the overlap of the black with the red and blue curves in Extended Data Fig. 6a). However, as haplotype number did not affect modelled range loss (compare Table 1), this marginal mismatch does not apparently impact our results. Homozygotes were ordered from the cold-adapted A1A1 up to the warm-adapted AnAn.Finally, the suitability curves of the genotypes were assumed to have the same value at their optimum as the species-specific suitability curve at this point (Extended Data Fig. 6).Artificial landscapesArtificial landscapes were defined as a raster of 50 × 112 cells (of 100 × 100 m). These rasters were homogeneous with respect to precipitation and soil carbon conditions which were set to the values optimal for each species according to the LRMs. With respect to temperature, by contrast, we implemented a gradient across the raster running from the minimum –9.1 °C to the maximum +3.8 °C temperature for which the LRM predicts values >OT across all six species. Buffering by 1 °C at both limits was done to avoid truncating simulation results. Further 4 °C at the lower limit were added to consider simulated temperature increase (below). The final temperature range represented by the raster ran from −14.1 to +4.8 °C. Temperature increased linearly along this axis by an increment of 0.171 °C per cell, derived from a 20° slope and a temperature decrease of 0.5 °C per 100 m of elevational change. Along the 50-cell breadth of the landscape, temperature values were kept constant. On the basis of these grids, we implemented a moderate and a severe climate change scenario, characterized by temperature increases of 2 and 4 °C, respectively, over 80 yr. Therefore, temperature of each raster cell increased annually by 0.025 and 0.05 °C, respectively.Simulations of spatiotemporal range dynamicsCATS21 is a spatially and temporarily explicit model operating on a two-dimensional grid (of 100 m mesh size in this case). CATS combines simulations of local species’ demography with species’ distribution models by scaling demographic rates relative to occurrence probabilities (suitabilities) predicted for the respective site or grid cell by the latter. Dispersal among grid cells is modelled as a combination of wind, exozoochoric and endozoochoric (that is, animal dispersal via attachment to the fur or ingestion and defecation, respectively) dispersal of plant seeds. Time proceeds in annual steps. The annually changing occurrence probabilities at each site affect the demography of local populations and hence, eventually, the number of seeds that are produced in each grid cell in the respective year. As a consequence, local populations grow or decrease, become extinct or establish anew and hence the species’ distribution across the whole study area changes as a function of the changing climate.Demographic modellingClimate-dependence of local demography was modelled by linking demographic rates (seed persistence, germination, survival, flowering frequency, seed yield and clonal reproduction) and carrying capacity to occurrence probabilities predicted by LRMs by means of sigmoidal functions. Furthermore, all rates were fixed at OT at a value ensuring stable population sizes; for more information see refs. 21,33. Demographic rates were confined by zero and a species-specific maximum value (Supplementary Table 1), which was assumed to be the same for all genotypes of a species. As a corollary, the demographic rates of all genotypes of a species are the same under optimal environmental conditions but their actual values for a particular site in a particular year differ due to different temperature optima of genotypes. In addition, germination, survival and clonal reproduction were modelled as density-dependent processes to account for intraspecific competition between genotypes. In our application, for all density-dependent functions, the species abundance is defined as the sum of all adult individuals in a given cell, regardless of their genotypes. Density dependence is commonly achieved by multiplying rates with (C – N)/C, where N is the species abundance and C is the (site- and genotype-specific) carrying capacity. This correction for density dependence causes the functions to drop to zero when N approximates C. To avoid the subsequent collapse of population sizes, we defined density-dependent rates as (C – N)/C × rate() + N/C × rate(OT), which ensures stable population sizes at densely populated sites occupied by only one genotype. To account for uncertainty in parameters of demographic rates, we assigned each species two value sets representing the upper and lower end of a plausible range of values on the basis of information derived from databases and own measurements (Supplementary Table 1).The simulations allowed for cross-pollination between genotypes. We used the relative amount of flowers (genotype-specific flowering frequency as defined by the sigmoid curve for the given suitability in the given raster cell for the given year × number of adults of that genotype in the population of that cell) to derive an estimate of the haplotype frequencies in the total pollen produced by the population within a grid cell. For the multiallelic case we allowed for recombination between loci with a recombination rate of 0.1%. These frequencies were set equal to the probability that particular haplotypes are transmitted to each year’s seed yield by pollination. Spatial pollen dispersal was accounted for in the following way: in each cell, 5% of the pollen involved in producing the annual seed yield, was assumed to stem from outside the respective raster cell. The proportions of different haplotypes in this 5% were derived from the overall pollen frequencies in all cells within a 700 m radius around the target cell (following estimates in ref. 52). Subsequently, produced seeds of each genotype were divided into resulting genotypes regarding the adult’s haplotype composition and the haplotype frequencies in the cells’ entire pollen load.Dispersal modellingFor wind dispersal of plant species we parameterized the analytical WALD kernel53 on the basis of measured seed traits and wind speed data from a meteorological station in the Central Alps of Austria. Exozoochorous and endozoochorous plant kernels were parameterized on the basis of correlated random walk simulations for the most frequent mammal dispersal vector in the study area, the chamois (Rupicapra rupicapra L.). For more details, see ref. 33. To account for uncertainties in species-specific dispersal rates, the proportion of seeds dispersed by the more far-reaching zoochorous kernels was assumed either as high (1–5%) or low (0.1–0.5%), setting upper and lower boundaries of a credible range of the dispersal ability of species.Simulation set up and simulation initializationTo test for the effects of climate change on genetic diversity in 2080, we ran CATS over the period 2000 to 2080 for each of the six study species across the entire Alps under a full factorial combination of (1) three niche sizes (50%, 75% and 100%); (2) six numbers of haplotypes (equal to two, three, five and ten alleles for one locus and four and eight for the diallelic two- and three-locus models, respectively); (3) three climate scenarios (current climate, RCP 2.6 and RCP 8.5); and (4) two sets of demographic and dispersal parameters. As a ‘control’ we also ran simulations for all climate scenarios and the two demographic and dispersal parameter sets for a setting with one genotype filling the whole niche of the species. To account for stochastic elements in CATS four replications were run for each combination of ‘treatments’.For simulations in artificial landscapes we used, instead of RCP 2.6 and RCP 8.5, ‘artificial’ climatic scenarios with an assumed warming of 2 and 4 °C, respectively, and no change in precipitation.All simulation runs were started with homozygotic individuals only. As initial distribution, for each simulation run each cell predicted to be environmentally suitable to the species (that is, occurrence probability of species >OT)—and within the real distribution range of the species28 (not relevant for simulations in artificial landscapes, of course)—was assumed to be occupied by an equal number of adults of each (homozygotic) genotype, with the total sum equal to the carrying capacity of the site. To accommodate this arbitrary within-cell genetic mixture of homozygotes (and missing heterozygotes) to actual local conditions we started simulations of range dynamics with a burn-in phase of 200 yr, run under constant current climatic conditions. To have a smooth transition from the burn-in phase under current climate (corresponding to the climate of the years 1970–2005; see current climate data) to future climate projections (starting with 2030) and to derive annual climate series, climate data were linearly interpolated between these two time intervals.Statistical analysisWe evaluated the contribution of climate scenario, haplotype number and haplotype niche size to overall species’ range change as well as the change in the frequency of the warm-adapted haplotype by means of linear models. In these models, log(range size in 2080/range size in 2000) and log(percentage of warm-adapted haplotype in 2080/percentage of warm-adapted haplotype in 2000), averaged over the four replicates and the two demographic and dispersal parameter sets, were the response variables. For the analysis of the change of the warm-adapted haplotype simulation settings with 100% niche size were ignored, as in this setting all haplotypes have the same temperature optimum (that is, neutral genetic variation). Approximate normality of residuals was confirmed by visual inspection.As indicator of the ‘topographic opportunity’ remaining to the species in the real world we calculated the area colonizable at elevations higher than those already occupied at the end of the simulation period. We therefore drew a buffer of 1 km around each cell predicted to be occupied in 2080 and then summed the area within these buffers at a higher elevation than the focal, occupied cell. Overlapping buffer areas were only counted once. This calculation was done for simulations conducted with one full-niche genotype per species only.Sensitivity analysisWe interpret the simulated relative decrease of warm-adapted haplotypes mainly as an effect of (1) the unrestricted expansion of cold-adapted haplotypes at the leading edge and (2) the resistance of the locally predominating haplotype that becomes increasingly maladapted with progressive climate warming, to recruitment by better-adapted haplotypes from below that are either rare or not present at all initially. However, the results achieved, and our conclusions, may be sensitive to assumptions about particular parameter values. Parameters that control the longevity of adult plants, and indirectly the rate of recruitment of new individuals, as well as those controlling gene flow via pollen (instead of seeds) may be particularly influential in this respect. We additionally ran simulations on artificial landscapes under alternative values of these parameters. In particular, we set the maximum age of plants to 10 yr instead of 100 yr and raised the proportion of locally produced pollen assumed to be transported up to 700 m to 10%. Both of these values are thus probably set to rather extreme levels: a maximum age of 10 yr is much shorter than the 30–50 yr assumed to be the standard age of (non-clonal) alpine plants31; and a cross-pollination rate between cells of 10% is high given that among the most important alpine pollinators only bumblebees are assumed to transport pollen >100 m regularly54,55. We add that we ran these additional simulations only in combination with the demographic species parameters set to high values because a short life span combined with low-level demographic parameters did not allow for stable populations of some species, even under current climatic conditions.For individual species, adapting plant age and cross-pollination rate between cells (Extended Data Fig. 7), did change the magnitude of loss of the warm-adapted haplotype. Nevertheless, for all of them the warm-adapted haplotype still became rarer when climate warmed and this effect increased with the level of warming. We are confident that our conclusions are qualitatively insensitive to variation of these parameters within a realistic range.Finally, in the simulations where we assumed a multilocus structure of the temperature niche, the recombination rate may also affect simulation results because it determines the rate by which new haplotypes can emerge. We also tested sensitivity of our simulations to doubling the recombination rate to 0.2%. Again, we found that a higher recombination rate had little qualitative effect on the results. Quantitatively, it resulted in a slightly pronounced relative decrease of the warmth-adapted haplotype in most species (Extended Data Fig. 8).Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Concentrations and sources of organic compoundsMid- and long-chain n-alkanesThroughout most sections, the biomarker composition of the Hala Hu record is dominated by long-chain n-alkanes (Fig. 2a). In most samples, concentrations were highest for nC31 (ca 5–35 µg/g d.w.) and decreasing with chain-length (Supplementary Fig. S3–1). Generally, the concentrations were low in the glacial period, and increased between ca. 10 and 8 ka cal BP, while gradually decreasing after ca. 5 ka cal BP. These alkanes can be attributed to vascular higher plants from the lake catchment27, which is characterized by alpine meadows. Mid-chain nC23 and nC25-alkanes are frequently attributed to aquatic macrophytes28. However, n-alkane patterns of aquatic and terrestrial plants overlap. Moreover, we consider the contribution from macrophyte to the sedimentary n-alkane pool at the coring location to be minor, because of the specific n-alkane pattern of the samples, the overall low concentrations of mid-chain n-alkanes, and the deep water depth of the core site. An exception is a section in the core dated to ca. 15–14 ka cal BP, that has concentrations from ca. 20 up to 50 µg/g d.w. for individual mid-chain n-alkanes (Fig. 2a). This can be explained by enhanced contribution from submerged aquatic plants29 when the lake level was much lower than present9,23. Alternatively, a contribution from other mid-chain producers, such as Sphagnum species30,31, is possible during phases of lake regression and the potential formation of peatlands around the lake shore.Fig. 2: Summary of measured proxy parameters in cores H7 and H11 (overlap 7.4−9.1 kyrs BP).a δD values and concentrations of mid- and long-chain n-alkanes. Arrows (ASM) indicate the maximum strength of the summer monsoon over Asia. b Concentrations of aquatic biomarkers (C20-HBI, alkenones, n-alkenes) and of microbial derived PMI; alkenone index (Uk´3738) and C37-alkenone δD values. c Concentrations of firemarkers (ΣM, L, G). d Titanium, sulfur, strontium, and calcium contents from XRF-scanning22,23. e Lake level reconstructed from ostracod assemblages9. Dark dashed interval 8.4−8.0 kyrs BP indicates mass flow layer. Light and medium grey shaded areas mark episodes of late glacial and mid-Holocene regime shifts within the aquatic ecosystem.Full size imageAlgal biomarkersAside from n-alkanes, a range of compounds of mostly aquatic origin can be identified in the aliphatic and ketone fraction of the sediment extracts. Unsaturated mid-chain alkenes nC21:1, nC23:1, nC25:1, and nC27:1 were abundant in traces in large parts of the core, but exhibit very high concentrations ( >100 µg/g d.w.) in the glacial period in a single sample at 17.5 ka cal BP and between 15 and 14 ka cal BP. Relatively high concentrations up to 25 µg/d.w. were also observed in the mid-Holocene sequence from ca. 9 to 5 ka cal BP (Fig. 2b). The origin of those compounds is not fully resolved and so far n-alkenes have not been detected in samples from aquatic and terrestrial vascular plants on the TP. However, algae such as eustigmatophytes and chlorophytes have been suggested as possible precursors in the open freshwater Lake Lugu, southeastern TP32, and in Lake Challa, eastern Africa33,34. Therefore, we consider phytoplankton as the most likely source of those compounds in Hala Hu.It is notable that different n-alkene distribution patterns are visible throughout the core, with a predominance of nC27:1 and nC25:1 in the mid-Holocene section and nC23:1 and nC25:1 in other core sections (Supplementary Fig. S3–1). This indicates either different source organisms for at least nC23:1 compared to nC27:1 (supported by more depleted δD-values, see below and Supplementary Fig. S3–3) or a change towards the synthesis of longer chain-lengths by the same source organisms, triggered by changing environmental conditions.The C20-highly branched isoprenoid (HBI) compound is another observed phytoplankton biomarker, which is widely absent in the glacial sequences but shows traces throughout the Holocene with peak abundances (up to 30 µg/g d.w) in the mid-Holocene (8−5 ka cal BP; Fig. 2b). As often been found in cyanobacterial and algal mats (e.g.,35,36,37), it has recently been assigned as a trophic indicator derived from diatoms in lake systems38.Other observed biomarkers of algal origin are alkenones, which are primarily produced by haptophytes, even though a variety of species are possible precursors39,40 The summed concentrations of C37-, C38-, and C39-compounds were low (300 ng/g d.w.; ca. 7.8–6.3 ka cal BP) (Fig. 2b). Here, the di- and tri-unsaturated homologues of all chain lengths appear to elute as peaks undisturbed by contamination, while the tetra-unsaturated alkenones show co-elution with a fatty acid ethyl ester (FAEE) at least for the C37-compound.Pentamethylicosane (PMI) is a compound that was detected in most samples in minor concentrations (400 ng/g d.w.) during an episode in the glacial period (ca 16.6−14 ka cal BP) and during the late Holocene (4.5 ka cal BP to present). This compound has been assigned to microorganisms, i.e bacteria and archaea, often related to the methane cycle41,42,43,44.Fire markersAnother analyzed biomarker group, anhydrous sugars levoglucosan (L), galactosan (G), and mannosan (M) are generated by combustion and pyrolysis of cellulose and hemicellulose, thus, are often referenced as “pyromarkers” or “firemarkers”45,46,47. They occur in low concentrations ( More

1.Science Plan and Implementation Strategy IGBP Report No. 53/IHDP Report No. 19 (Global Land Project, 2005).2.Millennium Ecosystem Assessment—Ecosystems and Human Well-Being: Desertification Synthesis Encyclopedia of the Anthropocene vols 1–5 (MEA, 2017).3.Huang, J., Yu, H., Guan, X., Wang, G. & Guo, R. Accelerated dryland expansion under climate change. Nat. Clim. Change 6, 166–171 (2015).Article 

Google Scholar 
4.Pimm, S. L. The complexity and stability of ecosystems. Nature 307, 321–326 (1984).Article 

Google Scholar 
5.Tilman, D. & Downing, J. A. Biodiversity and stability in grasslands. Nature 367, 363–365 (1994).Article 

Google Scholar 
6.Belnap, J. Surface disturbances: their role in acceleration desertification. Environ. Monit. Assess. 37, 38–57 (1995).Article 

Google Scholar 
7.Zhao, Y., Jia, R. L. & Wang, J. Towards stopping land degradation in drylands: water-saving techniques for cultivating biocrusts in situ. Land Degrad. Dev. 30, 2336–2346 (2019).Article 

Google Scholar 
8.Rodriguez-Caballero, E. et al. Dryland photoautotrophic soil surface communities endangered by global change. Nat. Geosci. 11, 185–189 (2018).CAS 
Article 

Google Scholar 
9.Coe, K. K. & Sparks, J. P. Physiology-based prognostic modeling of the influence of changes in precipitation on a keystone dryland plant species. Oecologia 176, 933–942 (2014).Article 

Google Scholar 
10.Ferrenberg, S., Tucker, C. L. & Reed, S. C. Biological soil crusts: diminutive communities of potential global importance. Front. Ecol. Environ. 15, 160–167 (2017).Article 

Google Scholar 
11.Belnap, J. & Gillette, D. A. Soil surface disturbance: impacts on potential wind erodibility of sand desert soils in SE Utah, USA. Land Degrad. Dev. 8, 355–362 (1997).Article 

Google Scholar 
12.Rutherford, W. A. et al. Albedo feedbacks to future climate via climate change impacts on dryland biocrusts. Sci. Rep. 7, 44188 (2017).13.Duniway, M. C. et al. Wind erosion and dust from US drylands: a review of causes, consequences, and solutions in a changing world. Ecosphere 10, e02650 (2019).14.Ferrenberg, S., Faist, A. M., Howell, A. & Reed, S. C. Biocrusts enhance soil fertility and Bromus tectorum growth, and interact with warming to influence germination. Plant Soil 429, 77–90 (2018).CAS 
Article 

Google Scholar 
15.Eldridge, D. J. et al. The pervasive and multifaceted influence of biocrusts on water in the world’s drylands. Glob. Change Biol. 26, 6003–6014 (2020).16.Ferrenberg, S., Reed, S. C. & Belnap, J. Climate change and physical disturbance cause similar community shifts in biological soil crusts. Proc. Natl Acad. Sci. USA 112, 12116–12121 (2015).CAS 
Article 

Google Scholar 
17.Reed, S. C. et al. Changes to dryland rainfall result in rapid moss mortality and altered soil fertility. Nat. Clim. Change 2, 752–755 (2012).CAS 
Article 

Google Scholar 
18.Concostrina-Zubiri, L. et al. Biological soil crusts across disturbance-recovery scenarios: effect of grazing regime on community dynamics. Ecol. Appl. 24, 1863–1877 (2014).CAS 
Article 

Google Scholar 
19.Weber, B., Bowker, M., Zhang, Y. & Belnap, J. in Biological Soil Crusts: An Organizing Principle in Drylands (eds Weber, B., Büdel, B. & Belnap, J.) 479–498 (Springer, 2016).20.Reynolds, J. F. et al. Global desertification: building a science for dryland development. Science 316, 847–851 (2007).CAS 
Article 

Google Scholar 
21.Berdugo, M. et al. Global ecosystem thresholds driven by aridity. Science 367, 787–790 (2020).CAS 
Article 

Google Scholar 
22.Bestelmeyer, B. T. et al. Analysis of abrupt transitions in ecological systems. Ecosphere 2, e03360 (2011).Article 

Google Scholar 
23.Bestelmeyer, B. T. et al. Desertification, land use, and the transformation of global drylands. Front. Ecol. Environ. 13, 28–36 (2015).Article 

Google Scholar 
24.Beisner, B., Haydon, D. & Cuddington, K. Alternative stable states in ecology. Front. Ecol. Environ. 1, 376–382 (2003).25.Fukami, T. & Nakajima, M. Community assembly: alternative stable states or alternative transient states? Ecol. Lett. 14, 973–984 (2011).Article 

Google Scholar 
26.Belnap, J. & Büdel, B. in Biological Soil Crusts: An Organizing Principle in Drylands (eds Weber, B., Büdel, B. & Belnap, J.) 305–320 (Springer, 2016).27.Belnap, J. & Warren, S. D. Measuring restoration success: a lesson from Patton’s tank tracks. Ecol. Bull. 79, 33 (1998).28.Belnap, J. & Elderidge, D. in Biological Soil Crusts: Structure, Function and Management (eds Belnap, J. & Lange, O. L.) 363–383 (Springer, 2001).29.Turner, M. G. Disturbance and landscape dynamics in a changing world. Ecology 91, 2833–2849 (2010).Article 

Google Scholar 
30.Scheffer, M. & Carpenter, S. R. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol. 18, 648–656 (2003).Article 

Google Scholar 
31.Sala, O. E. & Lauenroth, W. K. Small rainfall events: an ecological role in semiarid regions. Oecologia 53, 301–304 (1982).32.Cayan, D. R. et al. Future dryness in the Southwest US and the hydrology of the early 21st century drought. Proc. Natl Acad. Sci. USA 107, 21271–21276 (2010).CAS 
Article 

Google Scholar 
33.Christensen, N. S., Wood, A. W., Nathalie, V., Lettenmaier, D. P. & Palmer, R. N. The effects of climate change on the hydrology and water resources of the Colorado river basin. Clim. Change 62, 337 (2004).Article 

Google Scholar 
34.Herrick, J. et al. Field soil aggregate stability kit for soil quality and rangeland health evaluations. Catena 44, 27–35 (2001).Article 

Google Scholar 
35.Escolar, C., Martínez, I., Bowker, M. A. & Maestre, F. T. Warming reduces the growth and diversity of biological soil crusts in a semi-arid environment: implications for ecosystem structure and functioning. Phil. Trans. R. Soc. B 367, 3087–3099 (2012).Article 

Google Scholar 
36.Scheffer, M. et al. Creating a safe operating space for iconic ecosystems: manage local stressors to promote resilience to global change. Science 347, 1317–1319 (2015).CAS 
Article 

Google Scholar 
37.Collins, S. L., Micheli, F. & Hartt, L. A method to determine rates and patterns of variability in ecological communities. Oikos 91, 285–293 (2000).Article 

Google Scholar 
38.Rillig, M. C. et al. The role of multiple global change factors in driving soil functions and microbial biodiversity. Science 366, 886–890 (2019).CAS 
Article 

Google Scholar 
39.IPCC. Climate Change 2014: Impacts, Adaptations, and Vulnerability (eds Field, C. B. et al.) (Cambridge Univ. Press, 2014).40.Mirzabaev, A. et al. in IPCC Special Report on Global Warming of 1.5 °C (eds Masson-Delmotte, V. et al.) Ch. 3 (WMO, 2018).41.Torres-Cruz, T. J. et al. Species-specific nitrogenase activity in lichen-dominated biological soil crusts from the Colorado Plateau, USA. Plant Soil 429, 113–125 (2018).CAS 
Article 

Google Scholar 
42.Omernik, J. M. & Griffith, G. E. Ecoregions of the conterminous United States: evolution of a hierarchical spatial framework. Environ. Manag. 54, 1249–1266 (2014).Article 

Google Scholar 
43.Kuske, C. R., Yeager, C. M., Johnson, S., Ticknor, L. O. & Belnap, J. Response and resilience of soil biocrust bacterial communities to chronic physical disturbance in arid shrublands. ISME J. 6, 886–897 (2011).Article 

Google Scholar 
44.Tucker, C. L., Ferrenberg, S. & Reed, S. C. Climatic sensitivity of dryland soil CO2 fluxes differs dramatically with biological soil crust successional state. Ecosystems 22, 15–32 (2018). https://doi.org/10.1007/s10021-018-0250-445.Cable, J. M. & Huxman, T. E. Precipitation pulse size effects on Sonoran Desert soil microbial crusts. Oecologia 141, 317–324 (2004).Article 

Google Scholar 
46.Karl, T. R., Knight, R. W. & Plummer, N. Trends in high-frequency climate variability in the twentieth century. Nature 377, 217–220 (1995).CAS 
Article 

Google Scholar 
47.Kunkel, K. E., Easterling, D. R., Redmond, K. & Hubbard, K. Temporal variations of extreme precipitation events in the United States: 1895–2000. Geophys. Res. Lett. 30, 1895–2000 (2003).Article 

Google Scholar 
48.Kim, J. A projection of the effects of the climate change induced by increased CO2 on extreme hydrologic events in the Western US. Clim. Change 68, 153–168 (2005).CAS 
Article 

Google Scholar 
49.Smith, S. J. et al. Climate change impacts for the conterminous USA: an integrated assessment part 1. Scenarios and context. Clim. Change 69, 7–25 (2005). https://doi.org/10.1007/1-4020-3876-3_250.Schwinning, S., Belnap, J., Bowling, D. R. & Ehleringer, J. R. Sensitivity of the Colorado Plateau to change: climate, ecosystems, and society. Ecol. Soc. 13, 28 (2008).51.Daly, C. et al. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol. 28, 2031–2064 (2008).Article 

Google Scholar 
52.Jonasson, S. The point intercept method for non-destructive estimation of biomass. Phytocoenologia 11, 385–388 (1983).Article 

Google Scholar 
53.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).54.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).55.Oksanen, A. J. et al. Vegan: Community Ecology Package. Rpackage version 2.5-7 https://CRAN.R-project.org/package=vegan (2020).56.Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46 (2008).Article 

Google Scholar 
57.Venables, W. & Ripley, B. Modern Applied Statistics with S. (Springer, 2002).58.Lenth, R., Singmann, H., Love, J., Buerkner, P. & Herve, M. Package ‘emmeans’ https://github.com/rvlenth/emmeans (2018).59.Signorell, A. DescTools: Tools for Descriptive Statistics (2021).60.Hallett, L. M. et al. codyn: an R package of community dynamics metrics. Methods Ecol. Evol. 7, 1146–1151 (2016).61.Wood, S. N. Generalized Additive Models: An Introduction with R 2nd edn 1–476 (CRC/Taylor & Francis, 2017).62.Bürkner, P. C. brms: an R package for Bayesian multilevel models using Stan. J. Stat. Softw. 80, 1–28 (2017).63.Gelman, A. & Rubin, D. B. Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 547–511 (1992).
Google Scholar 
64.Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat. Comput. 27, 1413–1432 (2017).Article 

Google Scholar 
65.Modrák, M., Barrett, M., Weber, F. & Coronado, E. bayesplot: Plotting for Bayesian Models. R package version 1.8.0 https://mc-stan.org/bayesplot/ (2021).66.Makowski, D., Ben-Shachar, M. & Lüdecke, D. bayestestR: describing effects and their uncertainty, existence and significance within the Bayesian framework. J. Open Source Softw. 4, 1541 (2019).Article 

Google Scholar 
67.Phillips, M. L., Howell, A., Lauria, C. M., Belnap, J. & Reed, S. C. Data and software code from two long-term experiments (1996–2011 and 2005–2018) at three sites on the Colorado Plateau of North America (US Geological Survey, 2021); https://doi.org/10.5066/P9RUN1TP More

Study areaThe North Sea is a continental sea connected to the Atlantic Ocean through the English Channel in the southwest and between northern Shetland along the 61° latitude parallel to Norway in the north (Fig. 1). It is bordered by Norway, Denmark, Germany, the Netherlands, Belgium, France and the UK, and has a surface of 570,000 km2. The North Sea has an average depth of 95 m, yet maximum depths of ca. 700 m are found in the Norwegian Trench. The maximum tidal amplitude of the North Sea can reach up to 8 m, average winter sea surface temperatures are ca. 6 °C and average summer temperatures reach ca. 17°C33. The English Channel encompasses the marine strait between the UK and France. It covers 75,000 km2, has an average depth of 63 m, a maximum depth of 174 m and can reach a maximum tidal amplitude up to 12 m. The average winter and summer sea surface temperatures in the English Channel are ca. 5 and 20 °C, respectively54.TaggingIn total, 320 silver eels were tagged with pop-off archival tags (Table 1; Supplementary Table S2). In Belgium, 238 eels were caught and tagged at a drainage system upstream of the Yser Estuary (hereafter referred to as the Belgian eels) in 2018–2020 via nets that were attached to gravitational discharge sluice gates (coordinates: 51.127 N, 2.761 E) in October, November and December (n2018 = 102, n2019 = 60 and n2020 = 76). In Germany, 82 eels were tagged in 2011 and 2012. In early December 2011, seven eels were caught at Lake Plön (coordinates: 54.137 N, 10.334 E) with fyke nets. During September, October and November 2012, eels were caught in the Rivers Eider (n = 30; coordinates: 54.190 N, 9.093 E) and Havel (n = 45; coordinates: 52.419 N, 12.571 E) with fyke and stow nets, respectively.Upon capture, the eels were anaesthetized with 0.3 ml/L clove oil (Belgium), 0.4 ml/L ethylene glycol monophenyl ether (Germany 2011) or 120 mg/L MS-222 (Germany 2012), and various morphometric characteristics were measured to identify the life stage55: total length (to the nearest mm), weight (to the nearest g), horizontal and vertical eye diameter (to the nearest 0.01 mm in Belgium and to the nearest 0.1 mm in Germany) and pectoral fin length (to the nearest 0.01 mm and 0.1 mm in Belgium and Germany, respectively). Given that their total body length was  > 450 mm, all eels were considered female55. According to the morphometrics, five Belgian eels could be considered in the premigratory stage (FIII); however, based on visual inspection, they were considered silver eels (i.e. silver-coloured abdomen, dark grey on the dorsal side, jaw hinge not proceeding beyond the eye, enlarged eyes and dark coloured pectoral fins). The other 315 eels identified as silver eels based on both morphometry and visual inspection (201 FIV stage and 114 FV stage).Eels weighing ≥ 550 g were externally fitted with a G5 PDST (CEFAS Technology Ltd, UK), which log temperature and pressure (providing information on depth). They were attached applying the three-point Westerberg attachment method56. Two tag types were used: one with a separate tag and pop-off mechanism (Germany) and one where both mechanisms were integrated (Belgium). The flotation collar of the PDSTs was painted bright red, contained contact information and a cash reward to stimulate retrieval by the general public (e.g. beach combers and fishermen). The seven eels caught in 2011 in Germany (minimum 1220 g) were fitted with PSATs (X-Tag, Microwave Telemetry Inc., USA), also using the Westerberg-method56. Like the PDSTs, the PSATs record temperature and pressure. After release, they drift to the surface and transmit the data to the user via the ARGOS satellite system (www.argos-system.org). For the specifications of the different tags, we refer to Supplementary Table S3.Upon recovery from the anaesthetic, eels tagged with PDSTs were released close to their capture locations in the rivers Eider (coordinates 2011: 54.381 N, 9.009 E; coordinates 2012: 54.379 N, 9.013 E), Elbe (coordinates 1: 53.793 N, 9.402 E; coordinates 2: 53.569 N, 9.700 E; coordinates 3: 53.396 N, 10.171 E) and Yser (coordinates: 51.135 N, 2.757 E) (Table 1). The seven eels captured for PSAT tagging in 2011 were held for several weeks in the Thünen Institute of Fisheries Ecology, then tagged and released the same day; others were tagged in the field.PreprocessingOnce downloaded, the temperature and pressure data obtained from the PDSTs was subsampled to 1-min (Belgian eels) or 2-min (German eels) intervals to reduce the datasets and improve geolocation calculation time; this discrepancy is due to the minimum logging rate of the tags (Supplementary Table S3). Linear regression was applied to correct for pressure sensor drift over time. Indeed, pressure values increased over time even if the tag was kept at atmospheric pressure level. The regression was applied between 15 min before release and the moment the tag popped off and reached the surface, since the tag was then considered at sea level and hence to be under zero pressure.The PSAT data were retrieved through the ARGOS satellite system as a subset with 15-min intervals and converted to values of pressure and temperature. Contemporaneous values of temperature and depth were not always transmitted due to the transmission method. As a consequence of the tag release programming, the transmission of the first position for one of the tags was only received five days after the tag reached the sea surface.GeolocationThe daily movements of each electronically tagged European eel were reconstructed using an adapted version of the tidal geolocation model of Pedersen et al.57. The geolocation model uses a novel Fokker–Planck based method to combine the tidal location method of Metcalfe and Arnold58 with a hidden Markov model (HMM), such that an individual’s daily location d is modelled conditionally on its previous location (d − 1), its inferred behavioural state ds, where behaviour is defined by a single diffusivity parameter (i.e. the maximum amount of movement permitted in a given day), and the observations made between d and d − 1. In this case, observations consisted of the recorded depth (m; D1, …, Dn) and temperature (°C; T1, …, Tn), where n is the number of measurements made per day (the HMM down-samples to 10-min intervals, hence 144 measurements per day), and any hydrostatic (tidal) data which are derived from the sinusoidal pressure cycle recorded in the depth data when a fish is at rest on the seafloor. In addition to bathymetry and tidal amplitude with phase, the model was developed to include sea surface temperature (SST), which can provide additional validation when fish are swimming at or near the surface (i.e. depth ≤ 20 m)59,60, and temperature at depth, which can provide additional validation when fish remain at depths well below the sea surface61,62.The model was run in three different configurations for each recovered dataset: (i) using the tidal location model only (as for Pedersen et al.57), hereafter termed TLM geolocation; (ii) using the TLM plus sea surface temperature (as for Wright et al.60), hereafter termed SST geolocation, and (iii) using temperature at the surface and sub-surface, hereafter termed 3D geolocation (Supplemental Fig. S3). The final trajectory output for the PDST Belgian eels and PSAT German eels was obtained via 3D geolocation, while SST geolocation was used for the PDST German eels. The reason for this discrepancy is that the German PDST eels stayed closer to the coast and in shallower water. Consequently, the 3D geolocation results were more prone to error due to coastal influences on water temperature. As a result, we used the SST geolocation method for these datasets to obtain more reliable results.Data for the model were derived from publicly available resources. Gridded global bathymetry data were obtained from the general bathymetric chart of the oceans (Gebco; British Oceanographic Data Centre, Liverpool, United Kingdom, 2009). Tidal constituents were obtained from the Oregon State University Tidal Prediction model, as described in Egbert and Erofeeva63. Sea surface temperature data were sourced from OSTIA64, while temperature at depth data were sourced from the operational Mercator global ocean analysis and forecast system65. These datasets were downloaded from the Copernicus Marine Environmental Model Service (CMEMS: documented here http://resources.marine.copernicus.eu/documents/PUM/CMEMS-GLO-PUM-001-024.pdf). Data were sourced so as to fit the spatial scale of the model (30°N to 80°N and from 110°W to 60°E) and coarsened to reduce model run-time by modifying the spatial grid to a 1/10th of a degree resolution. The output of the model is a nonparametric probability distribution of the geographical position from which a most probable location, for each day at liberty, and a most probable movement path can be estimated.Prior to running the model, a number of constraints and input parameters were defined to ensure that the model ran effectively. The recapture information was either set as (a) the latitude and longitude where the tag was recaptured, with a high confidence ( 200 m) and hence did not exhibit diel vertical migrations, the input estimates of longitude were based on a simple linear interpolation from release to estimated pop-up. However, for eels that did reach oceanic depths, the time of local noon was estimated (based on the timing of significant diel vertical migrations, as for Righton et al.16), and used to estimate longitude. Geolocation was conducted with MATLAB software66.Migration routesOnly datasets containing ≥ 100 km of net tracking distance were included for further analysis, leading to 54 datasets from the 96 retrieved tags and 320 tagged eels. The net tracking distance was identified as the distance along the reconstructed trajectory between the release of the tagged eel and the pop-off event. When an eel was ingested by a predator, leading to the tag tracking the predator rather than the eel, the data were excluded from the day the eel was predated. The 100 km cut-off point was arbitrarily chosen to select migration paths of sufficient length for further analysis (e.g. migration direction); tracks had a minimum deployment duration of 4 days.Migration speedTo exclude a size-effect, we first applied an independent two-sample t-test to confirm eel sizes (i.e. weight) did not differ between Belgian and German eels. The assumptions of normality (Shapiro–Wilk test), homogeneity of variances (F-test) and independence were met (weight measurements are individual-specific and therefore independent).Next, an independent two-sample t-test was conducted to test if the total migration speeds (i.e. the ground speed along the reconstructed trajectory between the release of the tagged eel and the pop-off or predation event) differed between Belgian and German eels. The assumptions were tested and met as described above.Finally, we tested if the daily migration speed (i.e. the ground speed along the reconstructed trajectory per day) differed according to the eel’s position (i.e. modelled latitude and longitude) via a linear mixed effects model. The tag IDs were implemented as a random effect to account for autocorrelation. Since the two-sample t-test showed a significant difference between Belgian and German total migration speeds, we performed a separate analysis on eels from both countries. Assumptions of normality, homogeneity of variances and independence were tested and met.The migration speed analyses were conducted in R (version 3.6.3)67. The packages ‘lme4’ and ‘nlme’ were used to conduct the linear mixed effects model.Ethical statementEels were tagged using approved protocols by trained and individually licensed scientists working under national project authority in accordance with institutional and national guides for the care and use of laboratory animals. These guidelines are consistent with Institutional Review Board/Institutional Animal Care and Use Committee guidelines. Tagging in Belgium was carried out in accordance with the Belgian national and regional regulations for animal welfare and treatment (Permit ID: EC INBO-011). Tagging in Germany followed German legislation concerning care and use of laboratory animals, and ethical permission for the experiments was given by the Ministry of Energy, Agriculture, the Environment, and Rural Areas of the federal state Schleswig–Holstein (reference numbers V312-72241.123-34 (90-8/11) and V311-7224.123.3 (93-6/12) for tagging in 2011 and 2012 respectively). More

Plant defence traitsWe compiled species level data for five plant traits: wood density (WD), leaf and stem spinescence, latex production, and leaf size, for tropical and extra-tropical South and Central American woody species (i.e., the Neotropical biogeographic realm). WD was obtained for 2577 species from ref. 44. We only used wood density data from Zanne et al.44, because this study used WD measured in stems, whereas most other studies with available data used WD measured in branches. Leaf size data were obtained for 2660 woody species from Wright et al.37. We did not include leaf size from herbaceous species because herbaceous and woody species are influence by different megafauna guilds, suggesting distinct mechanisms, and because this dataset37 only included data for 253 Neotropical herbaceous species. The presence or absence of stem (and/or branch) spines (mostly thorns, but also prickles) were obtained from Dantas and Pausas45 for Neotropical savanna and forest species (1004 species) and complemented with other literature sources for other ecoregions (listed in the supplementary materials) using the names of the species for which we had WD and Leaf Size data. Our final stem spines dataset included 2843 woody species. We also compiled data on the presence of latex in plant stems and leaves for all the species for which we had data on other traits (3160 species; references in the supplementary materials). Finally, we also compiled data on leaf spines. While we managed to find leaf spine data for a total of 2173 woody species, we found spinescence in leaves to be especially concentrated in the palm Family (Arecaceae; 198 out of 221 species with leaf spines). Moreover, out of the non-palm species, all but three species also presented stem spines, indicating that, for other taxa, leaf spines might be dependent on the presence of stem spines at the region (in palms, 51% have stem spines). Thus, we only used leaf spinescence data of palm species (694 species) from the global Palm Traits Database 1.046.For wood density and leaf size, we often had more than one trait value per species (1005 and 831 species with more than one trait value, respectively). Thus, we computed the species mean trait value. This rarely occurred for binary traits (spinescence and latex) and, when occurred, the maximum value was used (0 for absence and 1 for presence). This later decision was based in the assumption that omitting the presence of spines or latex is more likely than incorrectly reporting the presence when it is absent. Moreover, some of these traits can be plastic18.From species to ecoregionsWe searched for geographical distribution data (coordinates) from the Global Biodiversity Information Facility (GBIF) for all of the species in each species-trait dataset (Data available from GBIF using the following doi: WD: https://doi.org/10.15468/dl.3vua3x; Stem spines: https://doi.org/10.15468/dl.ar5ddj; Latex: https://doi.org/10.15468/dl.m8dzjd; Leaf spines: https://doi.org/10.15468/dl.vv8gw4; Leaf size: https://doi.org/10.15468/dl.k98nxc). For this search, we used tools provided by the “rgbif” package for R in which species names are updated to the most recent classification and the returned occurrences also include those associated with synonyms (i.e., the “backbone” method). We labelled the obtained geographical coordinates according to their ecoregion and biogeographical realm (following Dinerstein et al.47) and cropped out occurrences falling outside of the Neotropical realm. Since occurrence data was not available to all the species in our initial trait dataset, the number of species used to calculate ecoregion level means was reduced to 2110 species, for wood density, 2133, for leaf size, 2629, for stem spines, 2714, for latex, and 657, for leaf spines. A detailed evaluation of the representativity of this data in relation to ecoregion- and Neotropical- level patterns can be found in the Supplementary Methods. Based on the occurrence data and their ecoregion label, we built a species abundance (columns) by ecoregion (rows) matrix for each trait.We obtained ecoregion scale abundance-weighted means for continuous traits (WD and Leaf Size) by: (1) Multiplying species abundance in each grid cell of the ecoregion by the mean species value; (2) Summing up the row values; (3) dividing the resulting row sum by the total species abundance (row sum prior to trait multiplication), and (4) calculating the ecoregions’ means (across all of the grid cells). For Stem Spines and Latex (binary traits), we used a similar procedure, but the maximum (0 for absence and 1 for presence) value was used instead of the mean in step (1), and step (2) was directly used to calculate the number of presences (i.e., 1 s). Moreover, instead of the steps (3) and (4), we calculated the number of absences as the difference between the total abundance (row sums before trait multiplication) and the values obtained in step (2). This process resulted in weighted means for WD and stem spinescence for 173 ecoregions, and Leaf Size and Latex for 174, out of the 179 Neotropical ecoregions. For leaf spinescence, we used a similar approach, although, because of the fewer species, the abundance estimate from GBIF was less reliable. Thus, we transformed the ecoregion species abundance to presence/absence before multiplying the trait values (0/1 for absence/presence). We obtained leaf spinescence data for 159 out of the 179 Neotropical ecoregions. The species- and ecoregion- level data is provided in the Supplementary Data and in ref. 47.Historical megafauna distributionWe obtained data on historical distribution of megafauna species from the MegaPast2Future/PHYLACINE_1.2 dataset24, a dataset containing distribution maps (96.5 km of spatial resolution) and functional traits for mammal species of the last 130,000 years. From this dataset, we obtained the probable past distribution of extinct large mammal herbivore (hereafter, “megafauna”) species, if these species were still alive today (“Present Natural” scenario; see details below). The “Present Natural” distribution of extinct species in this database is based on the estimated historical distribution (i.e., preceding anthropogenic range modifications) of extant species that are known (from the fossil record) to have coexisted with the extinct species. In this approach, an extinct species is considered to have been present in a given grid cell if at least 50% of the extant species that were found coexisting with the extinct species in the fossil (and subfossil) record was predicted to have occurred in the same cell prior to anthropogenic range modifications24,48. This approach assumes that, since extant and extinct species coexisted in the same locations, they must have had similar ecological requirements. It also assumes that megafauna extinction had anthropogenic causes, instead of causes related to climate change49, which is largely accepted in the literature50.We extracted the “Present Natural” distribution of extinct mammal (coded “EP” for IUCN status; i.e., “extinct in prehistory”, meaning before 1500 CE) whose body mass was higher than 50 kg (megafauna), and for which at least 90% of their diet consisted of plants (i.e., strict herbivores). For each Ecoregion, we began by calculating two megafauna-related metrics: extinct megafauna species richness (Mrich) and their mean body mass (Mbm). For this, we cropped the distribution maps of the megafauna species (containing 1 for presence and 0 for absence of each species) to the Neotropical realm. To calculate Mrich, we (1) counted species presences within each of the grid cells in the global grid (i.e., calculated the cell’s megafauna richness); (2) assigned the corresponding ecoregion label to the resulting richness grid cells, subset the richness cell values corresponding to the Neotropical region; and (3) calculated the mean for each Neotropical ecoregion. For Mbm, we replaced the presences of the megafauna species in the initial raster object (grid cell map of each megafauna species) by their body masses and calculated the grid cell-level mean body mass, before calculating the ecoregion-level means. We also calculated megafauna density and secondary productivity based on allometric equations that relate these metrics to megafauna body mass. However, we did not used megafauna density and secondary productivity because they were strongly correlated to megafauna richness (Supplementary Fig. 3). More details on how these metrics were calculated can be found in the Supplementary Methods.We also obtained diet preference information from the literature for most megafauna species that occurred in the Neotropical region (details and references in the Supplementary Material). Based on these information, we calculated the richness of large browser (MBrich for megabrowser richness), grazer (MGrich for megagrazer richness), and mixed-feeder (MMfrich for mega mixed-feeder richness) species by sub setting the megafauna species by grid cell array before the richness calculation in order to select only species that were classified within the correspondent subgroup.Extant herbivore mammal distributionWe also compiled data on the distribution, body mass and diet of extant and recently extinct (i.e., extinct after 1500 CE) herbivore mammal species (for simplicity, called ‘extant’ species in this study). As with megafauna maps, the distributions used represented reconstructions for periods preceding anthropogenic reduction of extant herbivores ranges (“Present Natural” scenario), based on abiotic, biotic and geographic variables48, rather than the currently observed distribution. This scenario was used because modern anthropogenic range reductions are too recent to produce substantial geographic effects at this spatial scale. These data were obtained by sub setting the MegaPast2Future/PHYLACINE_1.2 dataset to exclude species that were coded “EP” for IUCN status and that were not strict herbivores (at least 90% of the diet constituting of plants). We subsequently associated diet information to these species using data from ref. 51 and excluded all species that did not feed mainly on aboveground vegetative plant tissues (i.e., species that fed mostly on fruits, seed, roots were excluded). This later filtering was because the number of herbivores that feed mostly on seed and fruit increase with decreasing size (and this dataset included small mammals). We subsequently calculated the same metrics as for the extinct megafauna species (except for the richness of mixed-feeders as our source for diets50 labelled species according to dominant feeding pattern). For this, we used the same approach described for extinct megafauna species. We did not use a size threshold for extant species because there were only 13 extant mammal herbivore species with over 50 kg in the Neotropical region, most of which were grazers (9 species; 4 species were mixed-feeders and none were browsers). Therefore, we relied on the mean body mass metric calculated for extant mammals to detect potential size-related effects.Climate, soil, fire, insularity, and hurricanesFor each Ecoregion, we obtained data on climate (mean annual precipitation and temperature, and rainfall seasonality) and soil (sand content, pH, and cation exchange capacity) variables. Climate data was obtained from WorldClim 2.1 (10 min spatial resolution) and was based on climate data from 1970 to 200052. Soil data were obtained from SoilGrids (5 km of spatial resolution)53, and consisted of mean values for two depths, 0.05 and 2 m. We calculated Ecoregion level means for all of the soil and climate variables after intersecting the climate and soil grid maps with the ecoregion map.We obtained the number (a proxy for frequency) and intensity of wildfires per ecoregion area using the MODIS active fire location product (MCD14ML)54. We only considered fires (i.e., hotspots) with detection confidence of 95% or higher occurring from November 2000 to December 2019 (both included). To ensure that only wildfires were considered, we associated each fire pixel with a land cover type (300 m of spatial resolution) from ref. 55 for a buffer area of 1000 m surrounding the fire pixel centroid. We excluded all of the fires occurring in areas in which more than 10% of the surrounding land cover pixels corresponded to agricultural, urban and water classes. We calculated the number of wildfires per ecoregion area by dividing the fire count of each Ecoregion by the ecoregion area, and multiplying the resulting value by the proportion of vegetated land cover pixels (same classes used to exclude fires in anthropogenic areas and water bodies above). Fire intensity was estimated as the average fire radiative power across all detected MODIS hotspots in the ecoregion. Ecoregions lacking large preserved vegetated areas (criteria above) were excluded from subsequent analyses.Using the ecoregion map, we also classified ecoregions into insular (1), when most of the ecoregion area was located in islands, vs. continental (0), otherwise. This was performed because island biogeography theory predicts that, in island, species richness should be low due to low colonization and high extinction rates. Insularity has also been shown to reduce megafauna body size (i.e., the island rule), even though the mechanisms are not fully understood56. We also compiled data on hurricane activity, as woody density was suggested to confer resistance against this disturbance57. We used data from 1990 to 2019 from the HURDAT2 dataset58, containing six-hourly information about the location of all of the known tropical and subtropical cyclones (0.1° latitude/longitude). We used the sum of hurricane occurrences per ecoregions divided by ecoregion area as an indicator of hurricane activity.Statistical analysesTo understand megafauna patterns, we began by fitting (multiple) regression models with habitat-related (fire, climate, soil) and geographical (insularity) variables as predictors. We expected that megafauna richness in general was higher under savanna conditions (arid nutrient-rich or mesic nutrient-poor environments with frequent fires)1,22. We also expected that megafauna richness and body mass were affected negatively by insularity (i.e., following the island biogeography theory and island rule). Before the analyses, we tested the correlations among all of the variables that would eventually be entered as predictors in the same model for both the megafauna and trait models (Supplementary Table 1), in order to avoid multicollinearity associated with highly correlated variables (here, r ≥ 0.60). Since mean annual precipitation and soil pH were strongly positively correlated (r = −0.78), for all of the analyses (including the analyses with functional traits, described below), model selection was performed separately for these two variables (i.e., two different model selection procedures, one containing each of the two variables among the initial set of predictors). We selected the best among the two resulting models as that with the lowest AIC (differences higher than two points in all of the cases). To make sure that no multicollinearity remained we also calculated the Variation Inflation Factor (VIF) for all of the predictor variables as 1/tolerance, where tolerance is calculated as 1 minus the R2 of all of the model regressing a predictive variable against all of the other predictors. In all of the models, VIF was 3.33 or smaller (i.e., a tolerance of 0.30 or higher), indicating absence of multicollinearity.Model simplification was carried interactively using stepwise (both forward and backward) searching for the model with the lowest AIC (using R’s “step” function) and subsequently retaining only the significant variables (p ≤ 0.05). We calculated the Pearson r statistics as a measure of effect size for the selected variables as well as the associated confidence intervals, using the packages “parameters” and “effectsize” for R. The average contribution of each predictor variable was also calculated, using the package “dominanceanalysis”, as the mean difference in R2 before and after removing the target variable from models containing all of the possible subset combinations of the selected predictor variables, including the full selected model.For testing whether the studied plant functional traits were related to our megafauna indicators, we fit linear models to WD and leaf size, and generalized linear models (GLM; binomial family) for spinescence and latescence, using ecoregion as the unit. For spinescence and latescence, we used the matrix containing the count of spiny/latex and non-spiny/non-latex plants (species abundance; for stem spines and latex) or number of species with or without spines (for leaf spines; see above) as response variables. The predictor variables included the animal indicators for extinct megafauna and extant herbivores, as well as climate, soil, and fire predictors (and, for WD, hurricane counts). Because total, as well as megagrazer, megabrowser, and mega mixed-feeder species richness were strongly positively correlated (Supplementary Table 1), we used the richness difference between grazers and browsers to evaluate the effect of diet (Supplementary Fig. 1). For consistency, we used the same diet variable for extant and extinct species. Since we did not identify strong correlations among extinct megafauna and extant herbivore indicators (Supplementary Table 1), these variables were all entered simultaneously in the same initial models. As with the analyses of the megafauna indices, we also used r as effect size and calculated the average predictor contribution in terms of R2 for these models. For the later, we used the MacFadden Pseudo-R2 in the GLM models as implemented in the “pscl” and “dominanceanalysis” packages for R, as this statistic is the most comparable with R2 from linear multiple regression (Maximum Likelihood and Cragg and Uhler’s Pseudo-R2 were also calculated for the logistic models), and adjusted R2 for continuous traits. Islands were not included in these models, as island plants were expected to respond differently due to the effects of insularity on animal species richness, precluding megafauna and extant mammal richness from being accurate proxies for consumer abundance. For stem spines, we always included a quadratic term to both megafauna and extant mammal herbivore body mass, as evidence suggest that medium-size herbivores (i.e., approximately 250 kg) are important selective drivers of this trait12. If a significant relationship with our herbivory indicators (both extant and extinct) were significant but not indicative of a selective effect by herbivores (for more defended plants), this relationship was discarded (along with related variables, such as diet); this happened only once, for leaf size, which increased with extant herbivore richness (Supplementary Table 8).For all of the general linear regression models, assumptions of normality, homoscesticity and lack of spatial autocorrelation in the residuals were checked using the Kolmogorov–Smirnov, Breusch–Pagan and Moran’s I tests, respectively. For the later, ecoregions were considered neighbours when they were adjacent and non-neighbour otherwise. In some cases, heteroscesticity was detected and, thus, the significance of the coefficients was tested using heteroskedasticity-consistent covariance matrix estimation. If one or more variable lost their significance they were stepwise removed from the final model, beginning by the least significant, until all remaining variables had a significant effect. Overdispersion in the generalized linear model was also detected and dealt with using overdispersed binomial logit models, as implemented in the “dispmod” package for R, in which weights are interactively calculated and used to maintain the residual deviance lower than the degrees of freedom. To confirm that the detected associations between megafauna indices and plant traits were robust, we also tested the coefficient significance using randomization of the plant species by ecoregion matrices (see Supplementary Methods for details).To test the prediction that Neotropical ecoregions could be broadly classified into the three hypothesised antiherbiomes, we used hierarchical clustering on principal component axes of the ecoregion by trait matrix (five plant traits, standardized to zero mean and unit variance). We selected the number of clusters associated with the highest loss of inertia (within group variability) when progressively increasing the number of clusters, using the R package “FactoMineR”. This procedure allowed the recognition of large regions characterised by specific patterns of defence strategies (‘antiherbiomes’). We subsequently tested for axes score, megafauna and environmental differences among the resulting antiherbiomes to verify whether and how trait, climate and soil patterns matched those described for African ecosystems, and to understand the megafaunal differences among the antiherbiomes. For these comparisons, we used Kruskal-Wallis and post-hoc pairwise Dunn tests, using the Benjamini & Hochberg59 (1995) correction of P-values for multiple comparisons in both cases, and exclusively included continental ecoregions. For spines, we used the proportion of spinescent plants/species (rather than the number of “yes” and “no” used on previous analyses) in the principal component analysis. Because palms were missing from 20 ecoregions, we completed the values for these ecoregions using predicted model probabilities. To better understand these associations between traits and the environmental and megafauna variables, we also regressed the PCA axes against the same predictors used for traits.We also developed a framework to identify forest ecoregions most likely to have experienced a biome shift after megafauna extinction using antiherbiome, biome and megafauna distribution data. Ecoregions likely to have experienced a savanna-to-forest shift since the Pleistocene are those that: (1) are currently forest-dominated; (2) are classified in antiherbiomes analogous to African arid nutrient-rich or mesic nutrient-poor savannas; and (3) were megafauna- and, especially, megagrazer- rich during the Pleistocene (richness equal or greater than the 0.75 quantile: 14 species for Mrich, and 3 for exclusively grazing species; MGrich). We validated the distribution of these areas with fossil evidence (22 sites) from the Last Glacial Maximum and mid-Holocene (see Supplementary Methods and Supplementary Table 9). For this, we also used information about the present dominant vegetation type in the fossil sites, extracted from the reference sources (see Supplementary Table 9), to segregate savanna-forest shifts from data coming from stable savanna patches within forest or long-term savanna regions. We also contrasted the predicted patterns with the present location of savanna patches within the Amazon Forest region from ref. 60.All statistical analyses and data handling were carried out in the R (v.4.0.2) environment, using the previously mentioned packages, in addition to FSA, gridExtra, grid, lattice, lmtest, latticeExtra, olsrr, raster, rgdal, rgeos, sandwich, spatialreg, spdep and vegan, using codes provided in ref. 47.Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Forest fires are disasters that cause extensive damage to the entire world in economic, ecological, and environmental ways. These fires can be caused by natural reasons, such as high temperatures that can create spontaneous combustion of dry fuel such as sawdust, leaves, lightning, etc., or by human activities, such as unextinguished campfires, arson, inappropriately burned debris, etc1. According to research, 90% of the world’s forest fire incidents have occurred as a result of the abovementioned human carelessness1. The increase in carbon dioxide levels in the atmosphere due to forest fires contributes to the greenhouse effect and climate change. Additionally, ash destroys much of the nutrients in the soil and can cause erosion, which may result in floods and landslides.
At earlier times, forest fires were detected using watchtowers, which were not efficient because they were based on human observations. In recent history and even the present day, several forest fire detection methods have been implemented, such as watchtowers, satellite image processing methods, optical sensors, and digital camera-based methods2, although there are many drawbacks, such as inefficiency, power consumption, latency, accuracy and implementation costs. To address these drawbacks, a forest fire detection system using wireless sensor networks is proposed in this paper.Wireless sensor networks (WSNs) are self-configured and infrastructure-free wireless networks that help monitor physical or environmental conditions and pass these data through the network to a designated location or sink where the data can be observed and analyzed3. Efficiency and low power consumption are the major advantages of a WSN. In the proposed detection system, wireless sensor nodes are deployed according to cellular architecture to cover the entire area with sensors to monitor temperature, relative humidity, light intensity level, and carbon monoxide (CO) level using a microcontroller, transceiver module, and power components. The power supply to the sensor node is provided using batteries as the primary power supply, and solar panels are used as the secondary power supply. These sensor nodes are specially designed with a spherical shape to withstand damage caused by environmental conditions as well as animals.The sensor readings for each parameter are checked with a preset threshold ratio and a ratio that is calculated continuously in the node in real time, and only the ratios that exceed the preset ratio are sent from the sensor node to the base station for further analytical processing. The network utilized for this transmission is in the architecture of tree topology considering facts such as low power consumption, reduced latency, less complexity, etc. Cluster heads are used in this network to gather data from several sensor nodes and pass them on to the base station or the gateway node. The gateway node is an interface that connects the network with the secondary analysis process.For the analysis process, a machine learning regression model was used along with threshold ratio analysis to enhance detection accuracy. For the training and testing process of the model, data were collected during the fire and no fire situations in different areas and under different climatic conditions. During the data collection process, 7000 data samples were collected, where a data sample included temperature, relative humidity, light intensity level, and CO level at a particular time. Eighty percent of the collected data were randomly used as training data for the model, and the remaining 20% were used as test data.If the outcome of the machine learning model indicates a fire in a specific area, a text message will be sent to the mobile phone numbers of the authorized officers in responsible units. As this process is designed with a minimum delay, the fire can be detected within the initial stage, and the responsible parties can take necessary actions in a shorter period, which will minimize the damage.Related workForest fire detection has been a focus of many researchers for the last decade because of increased forest fire case reports from all over the world due to severe damage to society and the environment. Many methods have been proposed to detect forest fires, such as camera-based systems, WSN-based systems, and machine learning application-based systems, with both positive and negative aspects and performance figures of detection. Due to the higher probability of accurate and early detection due to the use of multiple sensor sources and deployment of sensor nodes in areas not visible to satellites, wireless sensor networks have a more positive outlook, and they have become the more applicable technology in many fields4.Many researchers have focused on environmental parameters, such as air temperature, relative humidity, barometric pressure, sound, light intensity, soil moisture, and wind speed and direction, along with gases, such as CO, CO2, methane, H2, and hydrocarbons apart from smoke, to detect forest fire conditions by considering the variations in these parameters during a fire5,6, and sensors have been selected according to the range, sensitivity, power consumption, and cost7,8.As supplying power to a sensor node is a challenging task in forested areas, utilizing only battery options is difficult because they do not last long, and distributing power using a wire would require a higher cost to deploy in a large forest. Therefore, many researchers have proposed solar-powered systems as secondary power sources along with rechargeable batteries as the main power source4,6, while some researchers have proposed solar batteries because they last longer9. To reduce the power consumption of sensor nodes, techniques such as keeping selected components active while others are deactivated have been proposed10,11,12.Most WSN-based detection systems are centered around a base station due to the memory and processing limitations of the nodes. Important and partially processed data are transferred to the base station through wireless media for processing and enabling relevant actions, while the base station also acts as the gateway between the sensor nodes and the system user4,9.When constructing a WSN, communicating data among the relevant entities is the main objective, and star topology and mesh topology-based networks have been proposed in many papers because of the different attributes in their systems. A mesh topology was chosen over a star topology because of its ability to self-organize, self-configure, and automatically establish among nodes in a network13. As a smaller number of nodes involved for transmission results in minimum energy consumption, concepts based on cluster heads have been used14. To minimize the loss of energy and data packages during transfer, a cluster-tree network topology structure was proposed15. Considering the sensing range of a node, fault tolerance, and energy consumption, a paper has proposed applying the on-demand k-coverage technique that provides event detection using static nodes with variable sensing ranges. This technique utilizes the maximum detection performance with the minimum power consumption for an event16. A survey on rare event detection has mentioned many event detection strategies that deliver maximized detection capability, minimized detection delay and low energy consumption, such as duty cycle, component deactivation, overpopulation/node redundancy, collaboration, and energy harvesting17.To reduce deployment cost and power consumption, a paper proposes a novel localization scheme that divides the whole forest area into different grids and allocates them to respective zones with another 8 neighboring grids. One centroid node from those grids, which is called the initiator node, predicts whether the zone is highly active (HA), medium active (MA), and low active (LA). Here, HA zones send data continuously to the base node through the interior node, MA zones send data periodically, and LA zones do not transfer data in the status that manages power consumption effectively18. Another author proposed obtaining data from the sensors every 2 min if there is the potential of a forest fire or obtaining data every 15 min otherwise to reduce the energy wastage19.To place sensor nodes in the most effective configuration to detect fire conditions, a sensor node was proposed at three different heights to perform different parameter measurements, while some authors have suggested covering sensor nodes to avoid direct sunlight exposure and minimize the false alarm rate4,10. As the network connectivity of service providers in forest areas is not robust, communication techniques that use dedicated network paths such as LoRa, ZigBee, and XBee have been used as the communication infrastructure. When considering attributes such as transmission range, high security, low power consumption (LPWAN protocol), and other relevant configurations, most papers have suggested using the LoRa module for transmission13,19.Most papers have suggested having threshold value-based fire detection on a sensor node, and if the exceeded threshold has remained the same, then a sink determines the location and will send an alarm to the fire department11. Because of the environmental parameter variations according to the place and time, threshold values are configured by the user considering geographic situations, climatic changes, seasonal changes, etc. after sensors obtain the data from the surrounding10.A fusion information process was proposed, where information from multiple sources is considered in making the final decision, which is better than using those sources individually, and two algorithms based on the threshold ratio method and Dempster-Shafer theory were used4. To enhance the detection accuracy, machine learning applications have been proposed in many papers based on different machine learning approaches, such as support vector machine (SVM) classification14 and regression techniques, such as logistic regression. However, applying machine learning techniques to fire detection systems has many limitations, such as the limited amount of energy, the energy required for data processing, the short range of communication and limited computations, the complexity of ML algorithms when executing on sensor nodes, and the difficulty of being distributed on every sensor node20,21. More

1.Thompson, R. S., Parish, J. C., Maidment, S. C. R. & Barrett, P. M. Phylogeny of the ankylosaurian dinosaurs (Ornithischia: Thyreophora). J. Syst. Palaeontol. 10(2), 301–312 (2012).Article 

Google Scholar 
2.Arbour, V. M. & Currie, P. J. Ankylosaurid dinosaur tail clubs evolved through stepwise acquisition of key features. J. Anat. 227, 514–523. https://doi.org/10.1111/joa.12363 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
3.Brown, C. M. An exceptionally preserved armored dinosaur reveals the morphology and allometry of osteoderms and their horny epidermal coverings. PeerJ 5, e4066. https://doi.org/10.7717/peerj.4066 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
4.Butler, R. J. & Barrett, P. M. Palaeoenvironmental controls on the distribution of Cretaceous herbivorous dinosaurs. Sci. Nat. 95, 1027–1032 (2008).CAS 
Article 

Google Scholar 
5.Brown, C. M. et al. Dietary palaeoecology of an Early Cretaceous armoured dinosaur (Ornithischia; Nodosauridae) based on floral analysis of stomach contents. R. Soc. Open Sci. 7, 200305. https://doi.org/10.1098/rsos.200305 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
6.Bourke, J. M., Porter, W. R. & Witmer, L. M. Convoluted nasal passages function as efficient heat exchangers in ankylosaurs (Dinosauria: Ornithischia: Thyreophora). PLoS ONE 13, e0207381. https://doi.org/10.1371/journal.pone.0207381 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
7.Mallon, J. C. & Anderson, J. S. The functional and palaeoecological implications of tooth morphology and wear for the megaherbivorous dinosaurs from the dinosaur park formation (Upper Campanian) of Alberta, Canada. PLoS ONE 9(6), e98605. https://doi.org/10.1371/journal.pone.0098605 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
8.Ősi, A., Prondvai, E., Mallon, J. & Bodor, E. R. Diversity and convergences in the evolution of feeding adaptations in ankylosaurs (Dinosauria: Ornithischia). Hist. Biol. 29(4), 539–570. https://doi.org/10.1080/08912963.2016.1208194 (2017).Article 

Google Scholar 
9.Hayashi, S., Carpenter, K., Scheyer, T. M., Watanabe, M. & Suzuki, D. Function and evolution of ankylosaur dermal armor. Acta Palaeontol. Pol. 55(2), 213–228 (2010).Article 

Google Scholar 
10.Norman, D. B. Scelidosaurus harrisonii from the Early Jurassic of Dorset, England: Cranial anatomy. Zool. J. Linn. Soc. 188(1), 1–81. https://doi.org/10.1093/zoolinnean/zlz074 (2019).Article 

Google Scholar 
11.Galton, P. M. Endocranial casts of the plated dinosaur Stegosaurus (Upper Jurassic, Western USA): A complete undistorted cast and the original specimens of Othniel Charles Marsh. In The Armored Dinosaurs (ed. Carpenter, K.) 103–129 (Indiana University Press, 2001).
Google Scholar 
12.Galton, P. M. Skull bones and endocranial casts of stegosaurian dinosaur Kentrosaurus Hennig, 1915 from Upper Jurassic of Tanzania, East Africa. Geol. Palaeontol. 22, 123–143 (1988).
Google Scholar 
13.Kuzmin, I. et al. The braincase of Bissektipelta archibaldi—new insights into endocranial osteology, vasculature, and paleoneurobiology of ankylosaurian dinosaurs. Biol. Commun. 65(2), 85–156. https://doi.org/10.21638/spbu03.2020.201 (2020).Article 

Google Scholar 
14.Leahey, L. G., Molnar, R. E., Carpenter, K., Witmer, L. M. & Salisbury, S. W. Cranial osteology of the ankylosaurian dinosaur formerly known as Minmi sp. (Ornithischia: Thyreophora) from the Lower Cretaceous Allaru Mudstone of Richmond, Queensland, Australia. PeerJ 3, e1475. https://doi.org/10.7717/peerj.1475 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
15.Paulina-Carabajal, A., Lee, Y. N., Kobayashi, Y., Lee, H. J. & Currie, P. J. Neuroanatomy of the ankylosaurid dinosaurs Tarchia teresae and Talarurus plicatospineus from the Upper Cretaceous of Mongolia, with comments on endocranial variability among ankylosaurs. Palaeogeogr. Palaeoclimatol. Palaeoecol. 494, 135–146. https://doi.org/10.1016/j.palaeo.2017.11.030 (2018).Article 

Google Scholar 
16.Pereda-Suberbiola, X. & Galton, P. M. On the taxonomic status of the dinosaur Struthiosaurus austriacus Bunzel from the Late Cretaceous of Austria. C. R. Acad. Sci. Paris II 315, 1275–1280 (1992).
Google Scholar 
17.Pereda-Suberbiola, X. & Galton, P. M. Revision of the cranial features of the dinosaur Struthiosaurus austriacus Bunzel (Ornithischia: Ankylosauria) from the Late Cretaceous of Europe. Neues Jahrb. Geol. Palaontol. Abh. 191, 173–200 (1994).
Google Scholar 
18.Pereda-Suberbiola, X. & Galton, P. M. Reappraisal of the nodosaurid ankylosaur Struthiosaurus austriacus Bunzel from the Upper Cretaceous Gosau Beds of Austria. In The Armored Dinosaurs (ed. Carpenter, K.) 173–210 (Indiana University Press, 2001).
Google Scholar 
19.Pereda-Suberbiola, X. A revised census of European Late Cretaceous nodosaurids (Ornithischia: Ankylosauria): Last occurrence and possible extinction scenarios. Terra Nova 4, 641–648 (1992).ADS 
Article 

Google Scholar 
20.Pereda-Suberbiola, X. Ankylosaurian dinosaur remains from the Upper Cretaceous of Laño (Iberian Peninsula). Est. Mus. Cienc. Nat. de Álava 14(Número especial 1), 273–288 (1999).
Google Scholar 
21.Pereda-Suberbiola, X., Astibia, H. & Buffetaut, E. New remains of the armoured dinosaur Struthiosaurus from the Late Cretaceous of the Iberian peninsula (Lafio locality, Basque-Cantabric Basin). Bull. Soc. Géol. Fr. 166, 207–211 (1995).
Google Scholar 
22.Garcia, G. & Pereda-Suberbiola, X. A new species of Struthiosaurus (Dinosauria: Ankylosauria) from the upper cretaceous of Villeveyrac (Southern France). J. Vertbr. Paleontol 23(1), 156–165 (2003).Article 

Google Scholar 
23.Nopcsa, F. Dinosaurierreste aus Siebenbürgen V. Geologica Hungarica, Series Palaeontologica 4, 1–76 (1929).
Google Scholar 
24.Codrea, V. et al. More than just Nopcsa’s Transylvanian dinosaurs: A look outside the Haţeg Basin. Palaeogeogr. Palaeoclimatol. Palaeoecol. 293, 391–405. https://doi.org/10.1016/j.palaeo.2009.10.027 (2010).Article 

Google Scholar 
25.Ősi, A. & Prondvai, E. Sympatry of two ankylosaurs (Hungarosaurus and cf. Struthiosaurus) in the Santonian of Hungary. Cretac. Res. 44, 30–38. https://doi.org/10.1016/j.cretres.2013.03.006 (2013).Article 

Google Scholar 
26.Bunzel, E. Die Reptilfauna der Gosaformation in der Neuen Welt bei Wiener-Neustadt. Abhandlungen der Kaiserlich-Königlichen Geologischen Reichsanstalt. 5, 1–18 (1871).
Google Scholar 
27.Franzosa, J. & Rowe, T. Cranial endocast of the Cretaceous theropod dinosaur Acrocanthosaurus atokensis. J. Vertebr. Paleontol. 25, 859–864 (2005).Article 

Google Scholar 
28.Witmer, L. M. & Ridgely, R. C. New insights into the brain, braincase, and ear region of tyrannosaurs (Dinosauria, Theropoda), with implications for sensory organization and behavior. Anat. Rec. 292(9), 1266–1296. https://doi.org/10.1002/ar.20983 (2009).Article 

Google Scholar 
29.Schade, M., Rauhut, O. W. M. & Evers, S. W. Neuroanatomy of the spinosaurid Irritator challengeri (Dinosauria: Theropoda) indicates potential adaptations for piscivory. Sci. Rep. 10(9259), 1613–1616. https://doi.org/10.1038/s41598-020-66261 (2020).Article 

Google Scholar 
30.Witmer, L. M., Ridgely, R. C., Dufeau, D. L. & Semones, M. C. Using CT to peer into the past: 3D visualization of the brain and ear regions of birds, crocodiles, and nonavian dinosaurs. In Anatomical imaging: towards a new morphology (eds Endo, H. & Frey, R.) 67–87 (Springer, 2008).Chapter 

Google Scholar 
31.Evers, S. W. et al. Neurovascular anatomy of the protostegid turtle Rhinochelys pulchriceps and comparisons of membranous and endosseous labyrinth shape in an extant turtle. Zool. J. Linn. Soc. 187, 800–828. https://doi.org/10.1093/zoolinnean/zlz063 (2019).Article 

Google Scholar 
32.Ősi, A., Pereda Suberbiola, X. & Földes, T. Partial skull and endocranial cast of the ankylosaurian dinosaur Hungarosaurus from the Late Cretaceous of Hungary: Implications for locomotion. Palaeontol. Electro. 17(1), 18p (2014).
Google Scholar 
33.Sampson, S. D. & Witmer, L. M. Craniofacial anatomy of Majungasaurus crenatissimus (Theropoda: Abelisauridae) from the Late Cretaceous of Madagascar. J. Vertebr. Paleontol. 27, 32–104. https://doi.org/10.1671/0272-4634(2007)27[32:CAOMCT]2.0.CO;2 (2007).Article 

Google Scholar 
34.Walsh, S. A., Barrett, P. M., Milner, A. C., Manley, G. & Witmer, L. M. Inner ear anatomy is a proxy for deducing auditory capability and behaviour in reptiles and birds. Proc. R. Soc. B 276, 1355–1360 (2009).Article 

Google Scholar 
35.Porter, W. R. & Witmer, L. M. Vascular patterns in the heads of dinosaurs: Evidence for blood vessels, sites of thermal exchange, and their role in physiological thermoregulatory strategies. Anat. Rec. 303, 1075–1103. https://doi.org/10.1002/ar.24234 (2019).Article 

Google Scholar 
36.Benson, R. B. J., Starmer-Jones, E., Close, R. A. & Walsh, S. A. Comparative analysis of vestibular ecomorphology in birds. J. Anat. 231, 990–1018. https://doi.org/10.1111/joa.12726 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
37.Bronzati, M. et al. Deep evolutionary diversification of semicircular canals in archosaurs. Curr. Biol. 31(12), 2520-2529.e6. https://doi.org/10.1016/j.cub.2021.03.086 (2021).CAS 
Article 
PubMed 

Google Scholar 
38.Hanson, M., Hoffman, E. A., Norell, M. A. & Bhullar, B. S. The early origin of a birdlike inner ear and the evolution of dinosaurian movement and vocalization. Science 372(6542), 601–609. https://doi.org/10.1126/science.abb4305 (2021).ADS 
CAS 
Article 
PubMed 

Google Scholar 
39.Miyashita, T., Arbour, V. M., Witmer, L. M. & Currie, P. J. The internal cranial morphology of an armoured dinosaur Euoplocephalus corroborated by X-ray computed tomographic reconstruction. J. Anat. 219, 661–675. https://doi.org/10.1111/j.1469-7580.2011.01427.x (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
40.Paulina-Carabajal, A., Lee, Y. N. & Jacobs, L. L. Endocranial morphology of the primitive nodosaurid dinosaur Pawpawsaurus campbelli from the Early Cretaceous of North America. PLoS ONE 11, e0150845. https://doi.org/10.1371/journal.pone.0150845 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
41.Spoor, F. & Thewissen, J. G. M. Comparative and functional anatomy of balance in aquatic mammals. In Sensory Evolution on the Threshold, Adaptations in Secondarily Aquatic Vertebrates (eds Thewissen, J. G. M. & Nummela, S.) 65–81 (University of California Press, 2008).
Google Scholar 
42.Sobral, G. & Müller, J. Archosaurs and their kin: The ruling reptiles. In Evolution of the Vertebrate Ear (eds Clack, J. A. et al.) 285–326 (Springer, 2016).Chapter 

Google Scholar 
43.Arbour, V. M. Estimating impact forces of tail club strikes by ankylosaurid dinosaurs. PLoS ONE 4(8), e6738. https://doi.org/10.1371/journal.pone.0006738 (2009).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
44.Park, J. Y. et al. A new ankylosaurid skeleton from the Upper Cretaceous Baruungoyot Formation of Mongolia: Its implications for ankylosaurid postcranial evolution. Sci. Rep. 11, 4101. https://doi.org/10.1038/s41598-021-83568-4 (2021).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
45.Witmer, L. M., Chatterjee, S., Franzosa, J. & Rowe, T. Neuroanatomy of flying reptiles and implications for flight, posture and behaviour. Nature 425, 950–953. https://doi.org/10.1038/nature02048 (2003).ADS 
CAS 
Article 
PubMed 

Google Scholar 
46.Walsh, S. A. et al. Avian cerebellar floccular fossa size is not a proxy for flying ability in birds. PLoS ONE 8(6), e67176. https://doi.org/10.1371/journal.pone.0067176 (2013).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
47.Ferreira-Cardoso, S. et al. Floccular fossa size is not a reliable proxy of ecology and behaviour in vertebrates. Sci. Rep. 7(1), 2017. https://doi.org/10.1038/s41598-017-01981-0 (2005).CAS 
Article 

Google Scholar 
48.Ezcurra, M. D. et al. Enigmatic dinosaur precursors bridge the gap to the origin of Pterosauria. Nature 588(7838), 445–449. https://doi.org/10.1038/s41586-020-3011-4 (2020) (Epub 2020 Dec 9).ADS 
CAS 
Article 
PubMed 

Google Scholar 
49.Carpenter, K., Sanders, F., McWhinney, L. A. & Wood, L. Evidence for predator–prey relationships: examples for Allosaurus and Stegosaurus. In The Carnivorous Dinosaurs (ed. Carpenter, K.) 325–350 (Indiana University Press, 2005).
Google Scholar 
50.Mallison, H. CAD assessment of the posture and range of motion of Kentrosaurus aethiopicus HENNIG 1915. Swiss. J. Geosci. 103, 211–233. https://doi.org/10.1007/s00015-010-0024-2 (2010).Article 

Google Scholar 
51.Arbour, V. M., Lech-Hernes, N. L., Guldberg, T. E., Hurum, J. H. & Currie, P. J. An ankylosaurid dinosaur from Mongolia with in situ armour and keratinous scale impressions. Acta Palaeontol. Pol. 58(1), 55–64. https://doi.org/10.4202/app.2011.0081 (2013).Article 

Google Scholar 
52.Scheyer, T. M. & Sander, P. M. Histology of ankylosaur osteoderms: Implications for systematics and function. J. Vertebr. Paleontol. 24, 874–893 (2004).Article 

Google Scholar 
53.Carpenter, K. Ankylosauria. In The Complete Dinosaur 2nd edn (eds Brett-Surman, M. K. et al.) 505–526 (Indiana University Press, 2012).
Google Scholar 
54.Jerison, H. J. Evolution of the brain and intelligence 482 (Academic Press, 1973).
Google Scholar 
55.Marugan-Lobon, J., Chiappe, L. M. & Farke, A. A. The variability of inner ear orientation in saurischian dinosaurs: Testing the use of semicircular canals as a reference system for comparative anatomy. PeerJ 1, e124 (2013).Article 

Google Scholar 
56.Benoit, J. et al. A test of the lateral semicircular canal correlation to head posture, diet and other biological traits in “ungulate” mammals. Sci. Rep. 10, 19602. https://doi.org/10.1038/s41598-020-76757-0 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
57.Coombs, W. P. Jr. Forelimb muscles of the Ankylosauria (Reptilia, Ornithischia). J. Paleontol. 52, 642–657 (1978).
Google Scholar 
58.Maidment, S. C. R. & Barrett, P. M. Osteological correlates for quadrupedality in ornithischian dinosaurs. Acta Palaeontol. Pol. 59(1), 53–70 (2014).
Google Scholar  More

Symbiont diversity and distributionThe present study provides the first evidence of symbiosis in R. alisae, a species of nudibranchs. This is the most multiple symbiosis that have ever been recorded for marine invertebrates. While many organisms establish an exclusively one-on-one relationship with a single microbial species or microbes belonging to the same functional group5,12, there are also organisms that harbor multiple microbial species, in which symbiont–symbiont and host–symbiont interactions occur. Six phylotypes of chemoautotrophic bacteria were reported for mussel Idas sp. from a cold seep area11 and five extracellular symbionts for the gutless oligochaete worm Olavius algarvensis34. However, in these cases, symbioses involving bacteria and marine invertebrates are either endosymbiotic microbes co-occurring inside the host bacteriocytes5,11 or ectosymbiotic microbes associated with the external surfaces of the animals3,4,9,15,34, with the exception of scaly-foot snail from hydrothermal vents having partnerships simultaneously with epi- and endosymbiontic microbes35.Bacterial symbionts in R. alisae have appeared to be more diverse than was previously known for marine invertebrates. It is evident that the detected symbiont phylotypes differ greatly from all other known symbionts found in marine invertebrates. Labrenzia (Rodobacteriales) and Maritalea (Rhizobiales) have not been recorded as forming symbiotic associations with invertebrates or plants so far, although other members of the families Rodobacteriales and Rhizobiales are well known symbionts14. Strains of Bradyrhizobium, Burkholderia, Achromobacter, and Stenotrophomonas are reported as symbionts of plants, interacting with a vast majority of nodulating legume species and efficient in biological nitrogen fixation36. This may be important when considering the nature of these symbionts in the nudibranch. Symbioses between cyanobacteria and marine organisms are commonly found among marine plants, fungi, sponges, ascidians, corals, and protists37,38. Synechococcus, identified as dominant symbiont clones of R. alisae (Table S2), is a unicellular cyanobacterium common in the marine environment, providing a range of beneficial functions including photosynthesis, nitrogen fixation, UV protection, and production of defensive toxins8,9,37. Symbiotic interactions between actinobacteria and their host have been observed in insects, human, animals, and plants, where the bacteria provide the host with protection against pathogens and produce essential nutrients39. However, none of the members of the clade Actinobacteria recorded in R. alisae are known to live symbiotically.Arrangement of symbiotic associationDespite the high diversity of bacteria, they are well organized in the host. Dense clusters of rod-shaped bacteria, Labrenzia, Maritalea, Bradyrhizobium, Burcholderia, Achromobacter, and Stenotrophomonas, were found within host-derived vacuoles, referred to as bacteriocytes, inside epithelial cells of R. alisae (Fig. 3). Although such arrangement differs from that typical of bacteriocytes, which are usually considered as specialized cells of the hosts for harboring bacteria, it resembles that reported for scaly-food snail from hydrothermal vents, which harbor symbionts in the esophageal gland35. Bacteriocytes in the gastropod Lurifax vitreus found near hydrothermal vents also constitute a portion of the mantle epithelium; they have large vacuoles containing many live and dividing bacteria40. Each bacteriocyte was densely packed with certain symbionts, and the bacteriocytes were randomly distributed within the epithelium cells. A distinctly regular distribution pattern was observed in the gill epithelium of the mussel Bathymodiolus sp.: the thiotrophic symbionts occupy the apical region, and the methanotrophic symbionts are more abundant in the basal region of bacteriocytes4. In the mussel Idas sp., however, there is no spatial pattern of the six distinct bacterial phylotypes, and the symbionts are mixed within bacteriocytes11.Synechococcus dominated the cytoplasm of intestinal epithelium and, more rarely, epidermis cells, mainly as specialized cell type referred to as nitrogen-fixing heterocysts. They are visually similar to cyanobacteria from corals and sponges8,37.The phylogenetic diversity and the spatial organization of the symbiotic community in R. alisae were determined by the 16S rRNA analysis, which was consistent with the results of FISH and TEM. Unlike most symbioses of marine invertebrates when bacteria house specialized host cells5,11 or cover epidermis7,15, symbiotic association of R. alisae exhibited spatial partitioning between symbionts, which were unevenly distributed between the tissues (Table S2). It has been established that different members of the microbial community can complement each other in acquisition of various restrictive nutrients, confirming the importance of the functional diversity of symbionts41. Thus, Stenotrophomonas rhizophila and Bradyrhizobium build a beneficial association in the rhizosphere and can act synergistically on promoting growth and nutrient uptake of soybean36. Cyanobacteria can interact synergistically with beneficial members from the endophytic microbiome of rice seedlings42. The location of bacterium in the organism of R. alisae may, in fact, depend on the specific metabolic and ecological roles that the symbionts play, and also on the interaction with bacterium belonging to different physiological groups.Nature of symbiosisSymbiotic associations between microbes and invertebrates are acquired mainly in a nutrient-depleted environment where symbionts usually provide their hosts with essential nutrients and high-energy compounds1. In contrast to known symbioses between microbes and gutless invertebrates, which obtain nutrients exclusively from the bacteria, R. alisae, like most nudibranch species, is a sponge-eating predator. However, due to the lack of adipose tissue, sponges are distinguished by a low lipid content (0.4 to 1.5% of wet weight)43 and also by specific proteinaceous spongin fibers and chitin, a polysaccharide similar to cellulose that can be indigestible for some predators, which together indicate their low nutritional value. Furthermore, R. alisae feeds exclusively on the sponge O. pennata; therefore, in habitats with low prey availability, this nudibranch has to survive starvation while searching for sponge assemblages. We suppose that symbiotic bacteria of R. alisae contribute to the utilization of low-quality food, similarly to symbiotic bacteria from the genera Rhodobacter, Burkholderia, and Aeromonas associated with the detritivorous isopod Asellus aquaticus44.A fatty acid analysis, as a useful approach to clarifying the nature of symbiosis5,20,32, has confirmed the trophic interaction between symbionts and the nudibranch host (Table S2). Among the fatty acids of symbiotic bacteria in R. alisae, OBFA are a major acyl constituent of membranes in Stenotrophomonas45 and also in Actinobacteria, Arthrobacter, Iamia, Ilumatobacter, and Kocuria46. Cis-vaccenic acid is a major component of Maritalea30. Omega-cyclohexyl tridecanoic acid (cyclo19:0) is specific to Bradyrhizobium47, Burkholderia, and Achromobacter48. Linoleic acid is produced by cyanobacteria including marine species of Synecoccocus33; in nudibranch, it obviously serves as a precursor in the synthesis of arachidonic acid (20:4n-6), thus, providing additional evidence for the transfer of fatty acids from symbionts to the host. Mollusks are capable of converting linoleic acid to arachidonic acid, since they have enzymes required for its synthesis21. The presence of these bacteria-specific markers and the abundance of arachidonic acid confirm the metabolic role of symbionts in the nudibranch host.Among nutrients, biologically available nitrogen can be considered a restrictive nutrient for the sponge-eating R. alisae, which can be acquired with the help of nitrogen-fixing symbionts, also referred to as diazotrophs. R. alisae harbors Bradyrhizobium, Burkholderia, Achromobacter, and Stenotrophomonas that are efficient in biological nitrogen fixation previously found to be associated with nodulating legume species36. Symbiotic nitrogen fixers are known to be associated with a variety of marine invertebrates such as wood-boring bivalves, corals, sponges, sea urchins, tunicates, and polychaetes7,8,37. Moreover, the protection of the enzyme nitrogenase that catalyzes N2 fixation against oxygen is an important physiological requirement in bacteria such as symbiotic Bradyrhizobium, Burkholderia, Achromobacter, and Stenotrophomonas that are located in bacteriocytes and provide this protection. Synechococcus is known as a nitrogen-fixer37,49. It performs N2 fixation in heterocysts where nitrogenase is restricted under oxic conditions. Indeed, heterocysts of Synechococcus are abundant in the intestine cells of R. alisae (Fig. 5B–D).Nitrate assimilation is one of the major processes of nitrogen acquisition by many heterotrophic bacteria and cyanobacteria50,51. Symbionts of R. alisae can play an important role in the process of nitrate utilization through denitrification, dissimilatory nitrate reduction, and assimilatory nitrate reduction as a nitrogen source and synthesize it into organic nitrogen. The nitrate reducers, Labrenzia52, Stenotrophomonas53, Maritalea30, and Rhodobacteraceae29 are widely represented in R. alisae. Synechococcus also utilizes nitrate, nitrite, or ammonium for growth50. Thus, symbiotic bacteria may play a significant role in the N-budget of the nudibranch mollusk.The symbiotic bacteria of R. alisae, including Bradyrhizobium, Maritalea, Labrenzia, Burkholderia, Achromobacter, Stenotrophomonas, Arthrobacter, Iamia, Ilumatobacter, and Kocuria, are known as carboxydotrophic or carbon monoxide (CO) oxidizers54,55. Despite the toxicity of CO for multicellular organisms, numerous aerobic and anaerobic microorganisms can use CO as a source of energy and/or carbon for cell growth56. The marine worm Olavius algarvensis establishes symbiosis with chemosynthetic bacteria using CO, a substrate previously not known to play a role in symbiotic associations with marine invertebrates, as an energy source57. We do not rule out that the R. alisae symbionts also might exploit CO as carbon and energy source. Despite this, assumption may seem impossible taking in account the CO toxicity, but, since many invertebrates (mollusks, tube worm, etc.) use toxic sulfate, thiosulfate, and methane as an energy source1,15, this hypothesis is worth to be addressed.An important component of skeleton in marine sponges of the family Microcionidae, including O. pennata, is the structural polysaccharide chitin58. Some bacteria are capable of hydrolyzing chitin via the activity of chitinolytic enzymes and can utilize chitin as a source of carbon, nitrogen, and/or energy59. Chitinase activity was documented for strains of Labrenzia60, Burkholderia61, Arthrobacter62, Achromobacter63, Stenotrophomonas64, Alcaligenes65, and actinobacteria59 associated with R. alisae. Thus, these bacteria can work synergistically to digest chitin and spongin, contributing to feeding success of the host nudibranch which depends solely on low-quality, nitrogen- and carbon-deficient food available.Furthermore, direct evidence has confirmed that many bioactive compounds from invertebrates are produced by symbiotic microorganisms66,67. Many biologically active compounds including toxic and deterrent secretions have been identified in nudibranchs of the family Discodorididae68. Symbiotic bacteria may exhibit toxic activity to provide the host nudibranch with chemical defense against predators and environment. Bacteria, especially actinobacteria, living in a symbiotic relationship with R. alisae may help the host in defense, since nudibranch lack a shell, and secondary metabolites of bacteria can provide them with chemical defense against predators and environment, as has been reported for some marine invertebrates2,9,10.In complex associations, the integration and coexistence of symbionts depend on supplementary partnerships and mutual contribution to the host’s metabolism41. The most intensively studied cases are highly specialized associations, where both partners can only exist in close relationship with one another. The relatively high diversity of microbes in R. alisae complicates understanding the complex pattern of molecular and cellular interactions between the host and its symbionts. More