Philippa Kaur

Castillo, C. P. et al. Agricultural Land Abandonment in the EU within 2015-2030. (Joint Research Centre (Seville site), 2018).van der Zanden, E. H., Verburg, P. H., Schulp, C. J. E. & Verkerk, P. J. Trade-offs of European agricultural abandonment. Land Use Policy 62, 290–301 (2017).Article 

Google Scholar 
Sun, Z. et al. Dietary change in high-income nations alone can lead to substantial double climate dividend. Nat. Food 3, 29–37 (2022).Article 
CAS 

Google Scholar 
Rasmussen, L. V. et al. Social-ecological outcomes of agricultural intensification. Nat. Sustainability 1, 275–282 (2018).Article 

Google Scholar 
Arneth, A. et al. Post-2020 biodiversity targets need to embrace climate change. Proc. Natl. Acad. Sci. 117, 30882–30891 (2020).Article 
CAS 

Google Scholar 
Seddon, N., Turner, B., Berry, P., Chausson, A. & Girardin, C. A. J. Grounding nature-based climate solutions in sound biodiversity science. Nat. Clim. Change 9, 84–87 (2019).Article 

Google Scholar 
Rey Benayas, J. M. & Bullock, J. M. Restoration of biodiversity and ecosystem services on agricultural land. Ecosystems 15, 883–899 (2012).Article 

Google Scholar 
Malhi, Y. et al. The role of large wild animals in climate change mitigation and adaptation. Current Biology 32, R181–R196 (2022).Article 
CAS 

Google Scholar 
Perino, A. et al. Rewilding complex ecosystems. Science 364, eaav5570 (2019).Article 
CAS 

Google Scholar 
Zhou, Y. et al. Limited increases in savanna carbon stocks over decades of fire suppression. Nature 603, 445–449 (2022).Article 
CAS 

Google Scholar 
Anderegg, W. R. L. et al. Climate-driven risks to the climate mitigation potential of forests. Science 368, eaaz7005 (2020).Article 
CAS 

Google Scholar 
Svenning, J.-C. Rewilding should be central to global restoration efforts. One Earth 3, 657–660 (2020).Article 

Google Scholar 
Dandy, N. & Wynne-Jones, S. Rewilding forestry. Forest Policy Econ. 109, 101996 (2019).Article 

Google Scholar 
Reino, L. et al. Does afforestation increase bird nest predation risk in surrounding farmland. Forest Ecol. Manag. 260, 1359–1366 (2010).Article 

Google Scholar 
Pausas, J. G. & Bond, W. J. Alternative biome states in terrestrial ecosystems. Trends Plant Sci. 25, 250–263 (2020).Article 
CAS 

Google Scholar 
Bakker, E. S. et al. Combining paleo-data and modern exclosure experiments to assess the impact of megafauna extinctions on woody vegetation. Proc. Natl. Acad. Sci. 113, 847–855 (2016).Article 
CAS 

Google Scholar 
Johnston, C. M. T. & Radeloff, V. C. Global mitigation potential of carbon stored in harvested wood products. Proceedings of the National Academy of Sciences 116, 14526–14531 (2019).Article 
CAS 

Google Scholar 
Shukla, P. R. et al. Climate Change and Land: an IPCC special report on climate change, desertification, landdegradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems. Intergovernmental Panel on Climate Change (2019).Hong, S. et al. Divergent responses of soil organic carbon to afforestation. Nat. Sustainability 3, 694–700 (2020).Article 

Google Scholar 
Bala, G. et al. Combined climate and carbon-cycle effects of large-scale deforestation. Proc. Natl. Acad. Sci. 104, 6550–6555 (2007).Article 
CAS 

Google Scholar 
Rohatyn, S., Yakir, D., Rotenberg, E. & Carmel, Y. Limited climate change mitigation potential through forestation of the vast dryland regions. Science 377, 1436–1439 (2022).Article 
CAS 

Google Scholar 
Beer, C., Zimov, N., Olofsson, J., Porada, P. & Zimov, S. Protection of permafrost soils from thawing by increasing herbivore density. Sci Rep-Uk 10, 4170 (2020).Article 
CAS 

Google Scholar 
Johnson, C. N. et al. Can trophic rewilding reduce the impact of fire in a more flammable world. Philos. Trans. R. Soc. B: Biological Sci. 373, 20170443 (2018).Article 

Google Scholar 
Kristensen, J. A., Svenning, J.-C., Georgiou, K. & Malhi, Y. Can large herbivores enhance ecosystem carbon persistence? Trends Ecol. Evol. https://doi.org/10.1016/j.tree.2021.09.006 (2021).Granado-Díaz, R., Villanueva, A. J. & Gómez-Limón, J. A. Willingness to accept for rewilding farmland in environmentally sensitive areas. Land Use Policy 116, 106052 (2022).Article 

Google Scholar 
Broughton, R. K. et al. Long-term woodland restoration on lowland farmland through passive rewilding. PloS one 16, e0252466 (2021).Article 
CAS 

Google Scholar 
Carver, S. et al. Guiding principles for rewilding. Conserv. Biology 35, 1882–1893 (2021).Article 

Google Scholar  More

The ecological baseline in TirezThe geology and the climate of the Tirez region favored the generation and maintenance of a type of hypersaline habitat characterized by extreme seasonality: the sulfate-chlorine waters, with sodium and magnesium cations, showed significant seasonal variations15. The alkaline pH, the low oxidant value for the redox potential of the water column and the highly reduced sediments imposed extreme conditions (see Table 1 and Supplementary Information for details). This extreme seasonality requires to define a valid representative ecological baseline to compare the ecology of the lagoon between 2002 and 2021 and, in this way, set the basis to proposing our model of ecological succession with increasing dryness as a “time-analog” for early Mars. Taxonomic data from 2002 is a snapshot of the community during one season, so we include in our discussion the results presented by Montoya et al. (2013) from a sample campaign carried out in 2005, because they16 analyzed both water and sediment and during both the wet and dry seasons.We consider here only the results obtained by Montoya et al.16 by gene cloning, since those obtained by isolation and sequencing are not comparable. At the level of large groups, no major seasonal differences were observed: Pseudomonadota, followed by Bacteroidetes, were the dominant phyla, in both water and sediments, and both in the dry and the rainy seasons; although Alphaproteobacteria was the dominant class in water, while Gammaproteobacteria was dominant in sediments (in both dry and rainy seasons). With respect to the archaeal domain, all the identified sequences were affiliated to Halobacteriales order, mainly Halorubrum (water) and Halobacterium (sediment), both within Halobacteriaceae family. We can consider these results presented in Montoya et al.16 as the “ecological baseline” for Tirez, however taken with a grain of salt, because only 43 bacterial and 35 archaeal sequences, including rainy and dry seasons and water and sediments, were considered for analysis.Prokaryotic diversity in 2002As can be expected for an extreme environment, the bacterial diversity detected in 2002 was low, although we cannot exclude the possibility that this may reflect the limitation of DNA sequencing techniques at the time. 59% of the obtained clones in the then-wet sediments corresponded to the Malaciobacter genus. Malaciobacter (previously named17 Arcobacter) is an aerotolerant Epsilonproteobacteria. Species within this genus are moderately halophilic, e.g., M. halophilus, capable to grow in up to 4% NaCl. Even though the role that Malaciobacter can play in the environment is not known, it seems to thrive in aquatic systems, like sewage, with a high organic matter content17: e.g., M. canalis, M. cloacae, or M. defluvii.After Malaciobacter-like clones, the next most numerous group belongs to the phylum Bacillota (27% of the sequenced clones; Table 2). Under stressful environmental conditions, members of the genus Virgibacillus produce endospores, a very useful property in an extreme and variable environment (ionic strength, temperature, light intensity), easy to compare with early Mars. Endospores facilitate species survival, allowing them to overcome drastic negative changes, like dry periods, and to germinate when the conditions are favorable again. The closest identified species was the halotolerant V. halodenitrificans, but with low homology, not far from other halotolerant (e.g., V. dokdonensis) or halophilic (e.g., V. marismortui) species within the same genus. The other Gram-positive clones belong to the order Clostridiales. These clones cluster in two taxonomic units related with the strictly anaerobic genus Tissierella.Despite the abundance of Pseudomonadota, their biodiversity was very low, reduced to only two genera within the Epsilon- and Delta-proteobacteria. Six sequences affiliated to Deltaproteobacteria, and clustered in one OTU (salB38, similarity 96.6% with Desulfotignum), were retrieved. Its presence in anaerobic media rich in sulfates (Table 1) seems reasonable. In fact, sulfate-reducing activity was detected using a specific enrichment assay.Finally, one taxon belonging to the phylum Spirochaetota (previously named Spirochaetes) was identified. The presence of Spirochaetota in this system is not strange because members of the genus Spirochaeta are very often found in mud and anaerobic marine environments rich in sulfates18. Moreover, the closest species to SalB63, although with a low similarity of 87%, was Spirochaeta bajacaliforniensis, a spirochete isolated19 from a microbial mat in Laguna Figueroa (Baja California), an extensive hypersaline lagoon with high gypsum content, very similar, although much bigger, than Tirez lagoon.The diversity within the domain Archaea was very low in 2002. The phylogenetic analysis of 96 clones indicate that they correspond to one specie belonging to the obligate halophile genus Methanohalophilus. Their high similarity (99.3%) with several species of Methanohalophilus, such as M. portucalensis (isolated from sediments of a solar saltern in Portugal), M. mahii (isolated from sediments of the Great Salt Lake), or M. halophilus (isolated20 from a cyanobacterial mat at Hamelin Pool, Australia), makes impossible its adscription to any particular species level. Methanohalophilus is strictly methylotrophic, which is consistent with this environment, given that the methylotrophic methanogenesis pathway, non-competitive at low-salt conditions, is predominant at high saline concentrations21. We further confirmed methanogenic activity in Tirez by the measurement of methane by gas chromatography in enrichment cultures.Prokaryotic diversity in context of other studies between 2002 and 2021It was challenging to establish a timeline for the succession of the populations involved, because the scarcity of data harvested and published so far from Tirez. However, combining our results with the few data available in Montoya et al.16 and Preston et al.22 on samplings carried out on 2005 and 2017, respectively, we can see a clear predominance of the phylum Pseudomonadota: Epsilonproteobacteria, i.e. Arcobacter-like, and Deltaproteobacteria, mainly sulfate-reducing bacteria (this work, sampling 2002), and Gammaproteobacteria16 when Tirez maintained a water film, to eventually a final predominance of Gammaproteobacteria, e.g. Chromatiales and Pseudomonadales, in the dry Tirez (this work, 2020 sampling). The Bacillales order has remained widely represented both in the wet and dry Tirez.Regarding the archaeal domain, the few references available (Refs.16,22; this work) confirm that the members of the Halobacteriaceae family are well adapted to both the humid and dry ecosystems of Tirez, being predominant in both conditions. Preston et al.22 found that the second most abundant group of archaea in the dry sediments of Tirez was the Methermicoccaceae family, within the Methanosarcinales order, Methanomicrobia class. Taking into account the results obtained in the dry Tirez (Preston et al.22; and this work, sample 2020), the methanogenic archaea have decreased drastically through time, probably due to salt stress and the competition with sulfate-reducing bacteria.Prokaryotic diversity in 2021From a metabolic point of view, most of the bacteria present today in the sediment are chemoorganotrophs, anaerobes, and halophilic or halotolerant. Scarce information is available about the predominant OTU, Candidate Division OP1. The OP1 division was one of the main bacterial phyla in a sulfur-rich sample in the deepest analyzed samples from the Red Sea sediments under brine pools23. In addition, the phylogenetically related Candidate division KB1 has been observed in deep-sea hypersaline anoxic basins at Orca Basin (Gulf of Mexico), and other hypersaline environments24. Eight of the nine genera identified show coverage greater than 1% of the sequences: i.e., Rubinisphaera, Halothiobacillus, Thiohalophilus, Anaerobacillus/Halolactibacillus, Halomonas, Halothermothrix, and Aliifodinibius are halophilic or halotolerant genera13,25.Regarding archaea, our analyses reveal archaeal groups that seem to thrive in sediments from extreme environments, e.g., marine brine pools/deep water anoxic basins or hypersaline lakes. The most abundant OTU, Thermoplasmata KTK4A, was found prominent and active in the sediment of Lake Strawbridge, a hypersaline lake in Western Australia26, and in soda-saline lakes in China27. The creation of a Candidatus Haloplasmatales, a novel order to include KTK4A-related Thermoplasmata, has been proposed27. On the other hand, both in the aforementioned soda-saline lakes in China27 and in a sulfur-rich section of the sediments from below the Red Sea brine pools23, retrieved sequences were assigned to Marine Benthic Groups B, D, and E. Finally, in the section of nitrogen-rich sediments from the aforementioned Red Sea brine pools, the unclassified lineage ST-12K10A represented the most abundant archaeal group. In the Tirez Lagoon sediment after desiccation, all Methanomicrobia readings belonged to this group.The significance and implications of an ecosystem characterized in 2021 by high diversity, high inequality, and lack of isolated representatives, resides in that Tirez is today an ecosystem in which many (most) of the species/OTUs present are dormant, and they do not play any metabolic role. Hence the high percentage of raretons, greater than 80% for both bacteria and archaea, which are actually present in the lagoon but with only one or two copies each. Only those species adapted to the conditions imposed by the extreme environment are able to actually thrive, and consequently only a few species carry out all the metabolic activity. We conclude that the microbiota in Tirez today represents an ecosystem with a high resilience capacity in the face of environmental changes that may occur.We want to clearly highlight that the technique available in 2002 to study the microbiota of the Tirez lagoon only allowed to obtain a low-resolution image, but that was the state-of-the-art procedure at the time, and the Tirez lagoon cannot be sampled again with the conditions back in 2002, which no longer exist and are not expected to return. Although we have kept in storage several samples of water and sediment from the 2002 Tirez lagoon, it is reasonable to assume that those laboratory microcosms would have chemically and microbiologically changed during the last 20 years, and as such no longer represent reliable replicas of the original lagoon, so we cannot use them for the purposes of this work. Therefore, we are aware that any comparisons of the 2002 laboratory results with the much more robust results obtained by Illumina in 2021 need to be taken with a grain of salt. With all the precautions required, in a high-level, first-order comparison, the most noticeable difference between 2002 and 2021 is a drastic change in the microbial Tirez population. Only some OTUs within Bacillales (Virgibacillus/Anaerobacillus), sulfate-reducing Deltaproteobacteria (Desulfotignum/Desulfobacteraceae-Desulfovibrio), and Spirochaetes are shared among the 2002 and 2021 samples. This comparison is enough for the purposes of this work, as we are interested in the evolution of the lagoon system as a whole to establish a “time-analog” with the wet-to-dry transition on early Mars, and not in the particular outcome of each and every OTU in Tirez. With the results at hand, we conclude that, since 2002, the lacustrine microbiota has shifted to one more adapted to the extreme conditions in the dry sediments, derived from the gradual and persistent desiccation concluding ca. 7 years ago (i.e., completely desiccated in 2015), such as lack of light, absence of oxygen, and lack of water availability. This shift has likely been triggered because organisms that were originally in the lagoon but at low abundance in 2002 became dominant as they were better adapted to desiccation, and because the incoming of new microorganisms transported by birds or wind28.Lipid biomarkers analysis of the desiccated lake sedimentsThe analysis of cell membrane-derived lipid compounds on the dry lake sediments at present allow to provide another perspective of the microbial communities inhabiting the Tirez lagoon, by contributing additional information about the ecosystem and depositional environment. It is important to note that, analyzing only the 2021 lake sediments, we cannot differentiate between lipidic biomarkers of the microorganisms inhabiting Tirez in 2002 and before from those left behind by the microorganisms living in the dried sediments today. Instead, the analyses of lipid biomarkers provide clues about the different microorganisms that have populated Tirez through time, including both older communities inhabiting the former aqueous system and also younger communities better adapted to the present dry conditions. Thus, the lipid biomarkers analysis can be considered as a time-integrative record of the microbial community inhabiting Tirez during the last decades.Based on the molecular distribution of lipid biomarkers, the presence of gram-positive bacteria was inferred from the relative abundance of the monounsaturated alkanoic acid C18:1[ω9], or iso/anteiso pairs of alkanoic acids from 12 to 17 carbons29 with dominance of i/a-C15:0 and i/a-C17:0 (Fig. 3B). In contrast, generally ubiquitous alkanoic acids such as C16:1[ω7], C18:1[ω7], or C18:2[ω6] suggested a provenance rather related to gram-negative bacteria30. The combined detection of the i/a-C15:0 and i/a-C17:0 acids, with dominance of the iso over the anteiso congeners, together with other biomarkers such as the mid-chain branched 10Me16:0, the monounsaturated C17:1, or the cyclopropyl Cy17:0 and Cy19:0 acids, may be associated with a community of SRB31 in today´s dry sediments of Tirez. Specifically, most of those alkanoic acids have been found in a variety of Deltaproteobacteria and/or Bacteroidota (previously named Bacteroidetes). The presence of archaea was deduced from the detection of prominent peaks of archaeol in the polar fraction32 (Fig. 3C), as well as squalene and relatives (dihydrosqualene and tetrahydrosqualene) in the apolar fraction33 (Fig. 3A). Squalene and a variety of unsaturated derivatives are present in the neutral lipid fractions of many archaea with high abundances in saline lakes34. The relative abundance of autotrophs over heterotrophs35 can be estimated by the ratio of the autotrophically-related pristane and phytane over the both autotrophically- and heterotrophically-produced n-C17 and n-C18 alkanes ([Pr + Ph]/[n-C17 + n-C18]). A ratio of 0.56 in the Tirez sediments suggest the presence of a relevant proportion of heterotrophs in the ancient lacustrine system.Furthermore, the lipid biomarkers analysis was able to detect compounds specific of additional microbial sources, such as cyanobacteria36 (n-C17, C17:1, or 7Me-C15 and 7Me-C17), microalgae and/or diatoms (phytosterols37; or C20:5, and C22:6 alkanoic acids30), and other photoautotrophs (phytol and potentially degradative compounds such as pristane and phytane31). A relatively higher preservation of the cell-membrane remnants (i.e., lipids) compared to the DNA-composing nucleic acids may contribute to explain the lack of detection of cyanobacteria, diatoms and microalgae, and other phototrophs by DNA analysis (a deficit in our results shared with Montoya et al.16, and Preston et al.22). Although abundant in higher plants38, sterols such as those detected here (i.e., the sterols campesterol, stigmasterol, and β-sitosterol, as well as ergosterol) are also major sterols in some microalgal classes37 (such as Bacillariophyceae, Chrysophyceae, Euglenophyceae, Eustigmatophyceae, Raphidophyceae, Xanthophyceae, and Chlorophyceae), cyanobacteria (β-sitosterol), and fungi (ergosterol39).The carbon isotopic composition of lipid biomarkers provides a rapid screening of the carbon metabolism in a system, by recognizing the principal carbon fixation pathways used by autotrophs. The range of δ13C values measured in the Tirez sediments (from − 33.9 to − 16.1‰) denotes a mixed use of different carbon assimilation pathways, involving mostly the reductive pentose phosphate (a.k.a. Calvin–Benson–Bassham or just Calvin) cycle (from − 19 to − 30‰), and in lesser extent the reductive acetyl-CoA (a.k.a. Wood–Ljungdahl) pathway (from − 28 to − 44‰), and/or the reverse tricarboxylic acid (rTCA) cycle (from − 12 to − 21‰).The lipids synthesized by microorganisms using the Calvin or reductive acetyl-CoA pathway are typically depleted relative to the bulk biomass, particularly those produced via de latter pathway. In the dry Tirez sediments, the majority of the lipid compounds are more depleted in 13C than the bulk biomass (Fig. 4). In particular, the branched alkane DiMeC18 (Fig. 4A) and the SRB-indicative 10Me16:0 acid (Fig. 4B) showed the most depleted δ13C values and suggested the use of the reductive acetyl-CoA pathway. The rest of lipid compounds showed isotopic signatures (from − 16.1 to − 31.4‰) compatible with the prevalence of the Calvin pathway. These values may directly reflect the autotrophic activity of microorganisms fixing carbon via the Calvin cycle or heterotrophic activity of microorganisms growing on their remnants. Thus, the saturated and linear alkyl chains of lipids (i.e., n-alkanes, n-alkanoic acids, and n-alkanols) showing the most negative δ13C values (e.g., alkanes n-C17 and C17:1; or acid C18:1[ω7]) reflect prokaryotic sources of Calvin-users autotrophs (e.g., cyanobacteria or purple sulfur bacteria), while the rest of compounds with slightly less negative δ13C values instead stem from the autotrophic activity of eukaryotes also users of the Calvin cycle (unsaturated fatty acids and sterols) or from the metabolism of heterotrophs such as SRB (iso/anteiso-, other branched, and cyclopropyl fatty acids) and haloarchaea (isoprenoids, phytanol, and archaeol). All in all, the compound-specific isotope composition of the dry sediments in the today´s Tirez lagoon may indirectly reflect the prevailing autotrophic mechanisms in the present lacustrine system of Tirez, by showing isotopic signatures of secondary lipids similar to their carbon source40.In addition, the use of a number of lipid molecular ratios or proxies allow further characterization of the lacustrine ecosystem and depositional environment. For example, the average chain length of the n-alkanes (24.1) suggests a relevant presence of eukaryotic biomass in the lacustrine sediments, since long-chained alkanes ( > C20) are known to originate from epicuticular leaf waxes in higher plants41. Highlighting the relevance of eukaryotes and their ecological roles is one of the major contributions of this work, because previous studies on the microbial ecology of hypersaline environments have been focused primarily on prokaryotes42.The proportion of odd n-alkanes of high molecular chain (i.e., n-C27, n-C29, and n-C31) over even n-alkanes of low molecular chain (i.e., n-C15, n-C17, and n-C19) provides an estimate of the relative abundance of terrigenous over aqueous biomass43, which in Tirez is TAR = 1.8. The Paq index may also be used to differentiate the proportion of terrigenous versus aquatic (emergent and submerged) plant biomass44. A Paq of 0.3 in the Tirez sediments from 2021 supported the relative abundance of land plants. Finally, the depositional environment in the lacustrine system of Tirez may be also characterized analyzing the ratio of pristane over phytane (Pr/Ph), which is higher than 1 when phytol degrades to pristane under oxic conditions45. Assuming that both isoprenoids in the Tirez sediments derived from phytol31, according to their similarly depleted δ13C (Fig. 4A), we can conclude that the sediments in the Tirez lagoon were deposited under predominantly oxic conditions (i.e., Pr/Ph ratio of 1.1).In summary, the lipid biomarkers study revealed useful information about the depositional environment and lacustrine ecosystem, including the presence of active or past autotrophic metabolisms involving prokaryotes (e.g., cyanobacteria and purple sulfur bacteria) and eukaryotes (plants, diatoms and other microalgae), as well as heterotrophic metabolisms of likely SRB and haloarchaea growing on Calvin-users exudates. These results are quite in agreement with the microbial community previously reported16 in sediments from the wet and dry seasons: abundant Gammaproteobacteria and Alphaproteobacteria, together with Algae and Cyanobacteria, dinoflagellates and filamentous fungi, Bacillota, Actinomycetota (previously named Actinomycetes), and a halophilic sulfate-reducing Deltaproteobacteria.Tirez as the first astrobiological “time-analog” for early MarsEarly Mars most likely had a diversity of environments in terms of pH, redox conditions, geochemistry, temperature, and so on. Field research in terrestrial analog environments contribute to understand the habitability of this diversity of environments on Mars in the past, because terrestrial analogues are places on Earth characterized by environmental, mineralogical, geomorphological, or geochemical conditions similar to those observed on present or past Mars9. Therefore, so far analogs have been referred to terrestrial locations closely similar to any of the geochemical environments that have been inferred on Mars, i.e., they are “site-analogs” that represent snapshots in time: one specific environmental condition at a very specific place and a very specific time. Because of this, each individual field analog site cannot be considered an adequate representation of the changing martian environmental conditions through time. Here we introduce the concept of astrobiological “time-analog”, referred to terrestrial analogs that may help understand environmental transitions and the related possible ecological successions on early Mars. In this sense, they should be “time-resolved analogs”: dynamic analog environments where we can analyze changes over time. To the best of our knowledge, this is the first study that looks at the environmental microbiology of a Mars astrobiological analog site over a significant and long period of change, and try to understand the ecological successions to put them in the context of martian environmental evolution.As Mars lost most of its surface water at the end of the Hesperian5,9,12, this wet-to-dry global transition can be considered the major environmental perturbation in the geological history of Mars, and therefore merits to be the first one to be assigned a “time-analog” for its better understanding and characterization. The drying of Mars was probably a stepwise process, characterized by multiple transitions between drier and wetter environments12,47, and therefore the seasonal fluctuations and eventual full desiccation of Tirez represent a suitable analog to better understand possible ecological transitions during the global desiccation of most of the Mars’s surface before the Amazonian (beginning 3.2 Ga).To introduce Tirez as the first Mars astrobiological “time-analog” of the wet-to-dry transition on early Mars, the objective of this study was threefold: first, we wanted to identify the dominant prokaryotic microorganisms in the active Tirez lagoon 20 years ago, a unique hypersaline ecosystem with an ionic composition different from that of marine environments, and therefore potentially analogous to ancient saline lacustrine environments on Mars during the Noachian and into the Hesperian46,47. Our results provide a preliminary basis to hypothesize how the microbial communities on the Noachian Mars could have developed in salty environments with dramatically fluctuating water availability. The requirement to deal with important variations in ionic strength and water availability, involving at times the complete evaporation of the water, could have represented additional constraints48 for microorganisms on early Mars.The second objective of this investigation was the identification of the microbial community inhabiting the desiccated Tirez sediments today, after all the water was lost, as a potential analog to desiccated basins on Mars at the end of the Hesperian1,3,4,47. Our results suggest that hypothetical early microbial communities on early Mars, living with relative abundance of liquid water during the Noachian, would have been forced to adapt to increasingly desiccating surface environments, characterized by extreme conditions derived from the persistent dryness and lack of water availability. Our investigation in Tirez suggest that hypothetical microorganisms at the end of the Hesperian would have needed to evolve strategies similar to those of microorganisms on Earth adapted to living at very low water activity49, to thrive in the progressively desiccating sediments.And the third objective of this investigation was the identification of the lipidic biomarkers left behind by the microbial communities in Tirez, as a guide to searching and identifying the potential leftovers of a hypothetical ancient biosphere on Mars. Lipids (i.e., fatty acids and other biosynthesized hydrocarbons) are structural components of cell membranes bearing recognized higher resistance to degradation relative to other biomolecules, thus with potential to reconstruct paleobiology in a broader temporal scale than more labile molecules50. Our results reinforce the notion that lipidic biomarkers should be preferred targets in the search for extinct and/or extant life on Mars precisely because they are so recalcitrant. More

Sampling sites and sample preparationCulms of three- to four-year-old of Bambusa vulgaris Schrad. ex J.C., B. maculata Widjaja, and Gigantochloa atter (Hassk) Kurz ex Munro were collected from naturally bamboo forests at four sites in Lombok Island, Indonesia23. The culm age was estimated based on some morphological features (the presence of culm sheath, color, and sound created by tapping with fingers) checked by an experienced bamboo farmer. Figure 1 shows the map of sampling sites and climatic conditions of the sites. Ten individual culms in each species at each site were collected from different clumps and cut 20 cm above the ground (Fig. 2). A total of 120 culms (three species × four sites × 10 individual culms from 10 individual clumps) were collected in the present study (Fig. 2). To determine the longitudinal variations of the anatomical characteristics and mechanical properties, the internode section was collected at 2-m intervals from 2 to 8 m above the ground; a total of 480 internode sections. (120 culms × four heights) were obtained from three species (Fig. 2). The collection of bamboo culms was permitted by Indonesian Institute of Science (Reference no. B-206/SKIKH/KS.02.04/X/2020) and complied with relevant guidelines and regulations of Indonesian CITES Management Authority, Ministry of Environment and Forestry, Indonesia. In addition, the voucher specimen was deposited at the Herbarium Lesser Sunda, University of Mataram, Indonesia under the voucher number of DSR01, 02, and 03 (specimens were identified by Mr. Niechi Valentino). Table 1 shows the culm diameter at 1.3 m above the ground, total culm height, and mean value of culm thickness at four positions23.Figure 1Locations and climate conditions of sampling sites in the present study23. Note: Site I, Tempos (8°41′59″ S, 116°8′40″ E); Site II, Kabul (8°47′21″ S, 116°10′21″ E); Site III, Keruak (8°45′45″ S, 116°28′54″ E); Site IV, Genggelang (8°23′16″ S, 116°15′35″ E). *, mean annual precipitation. The value in the bracket is the mean annual temperature. Climate data were provided from Nusa Tenggara River Basin Management I, Indonesia. Mean monthly temperature and precipitation were calculated by averaging monthly values from 2016 to 2018. Bars indicate the mean values of precipitation. Circles indicate the mean values of temperature. The graph was originally created by R27 (version 4.0.3, https://www.R-project.org/).Full size imageFigure 2Photographs of the clumps in three bamboo species (a–c) and schematic diagrams of experimental procedures (d). Note: a, B. vulgaris; b, B. maculata; c, G. atter. The specimens of fiber area measurement and mechanical properties have the whole culm thickness (including the cortex and inner part of the culm) in the radial direction.Full size imageTable 1 Mean values and standard deviations of growth characteristics in three bamboo species at each site23.Full size tableAnatomical characteristicsThe internode sections were split into two parts: the strips (10 mm in the longitudinal direction) and the small blocks (10 [T] mm by 10 [L] mm by culm thickness in the radial direction) (Fig. 2). The strips and small blocks were the samples for measuring fiber length and fiber area, respectively. In the present study, the fiber area was defined as the sheaths area around the vascular bundles24.To determine the fiber length, small sticks (not including the cortex and the most inner part of the culm) were obtained from the strips with a razor blade (Fig. 2). Randomly selected sticks from each height position (without separation of collected positions of the samples within the radial direction of the culm in a height) were macerated with Schultze’s solution (100 mL of 35% nitric acid containing 6 g potassium chloride) at 70 °C for two hours. The length of 50 fibers was measured in each sample with a digital caliper (CD-15CX, Mitutoyo, Kawasaki, Japan) on a microprojector (V-12B, Nikon, Tokyo, Japan).To measure the fiber area, one block was taken at each height position on each individual culm (Fig. 2). The transverse sections of the blocks were polished with sandpaper sheet (#180, 3 M Japan, Tokyo, Japan), and then their images were captured using a microscope digital camera (DS-2210, Sato Shouji Inc., Kawasaki, Japan) attached to a stereo microscope (SZX12, Olympus, Tokyo, Japan). The fiber area was determined by ImageJ25 (version 1.53e). Binarized images were prepared by ImageJ to distinguish as clearly as possible between the vascular bundle and the background (Fig. 3). The darker area of binarized images in Fig. 3 was identified as fiber sheaths. The fiber area was calculated as follows:$$FAleft( % right) , = A_{fs} /A_{c} times {1}00$$
(1)
where FA = fiber area (%), Afs = the transverse-sectional area of fiber sheath in bamboo culm (mm2), and Ac = the transverse-sectional area of bamboo culm (mm2).Figure 3The photomicrographs of transverse section in B. vulgaris (a and d), B. maculata (b and e), and G. atter (c and f). Note: a, b and c, original image; d, e and f, binarized image processed by ImageJ25 (version 1.53e, https://imagej.nih.gov/ij/). The darker area in photomicrographs (d, e and f) is fiber sheath area.Full size imageMechanical propertiesThe following mechanical properties of culm were measured: bending properties (MOE and MOR), CS, and tensile properties (TM and TS). A total of 480 specimens (one specimen × four heights in an individual × ten individuals × three species × four sites) without node were obtained in each property (Fig. 2).The strips (10 [T] mm × 200 [L] mm × varied culm thickness in the radial direction) were prepared as the specimens for the static bending test (Fig. 2). The static bending test was conducted using a universal testing machine (MSC 5/500–2, Tokyo Testing Machine, Tokyo, Japan). A load was applied to the center of the specimen on the outer cortex surface with 180 mm span and 3 mm min−1 load speed. Due to larger thickness (exceeded 12.9 mm = 180 mm of span / 14) in the radial direction, the span / depth ratio in some specimens was less than 14, indicating that MOR in some specimens might be underestimated due to the occurrence of the shearing strength26. Of 480 specimens, the large culm thickness exceeded 12.9 mm was total 19 specimens from B. vulgaris species collected at 2 m height position from different sites (Site I = four specimens, Site II = six specimens, Site III = four specimens, and Site IV = five specimens). However, all these 19 specimens were broken at the tension side of the specimens during static bending test, which was the normal breaking forms of bending specimens with span / depth ratio less than 14.The load and deflection were recorded with a personal computer, and then MOE and MOR were calculated by the following formulae:$$MOE , left( {GPa} right) , = Delta Pl^{3} / , 4Delta Ybh^{3} , times 10^{ – 3}$$
(2)
$$MOR , left( {MPa} right) , = , 3Pl/ , 2bh^{2}$$
(3)
where ΔP = difference between upper and lower proportional limit within the range of elasticity (N), l = length of the span (mm), ∆Y = deflection due to ∆P (mm), b = width of the specimen (mm), h = height of the specimen (mm), and P = maximum load (N).The compressive test specimen (10 [T] mm × 20 [L] mm × culm thickness in the radial direction) was also prepared (Fig. 2). The test was conducted using a universal testing machine (RTF-2350, A&D, Tokyo, Japan) with a load speed of 0.3 mm min−1. The compressive strength parallel to grain (CS) was calculated by the following formula:$${text{CS }}left( {{text{MPa}}} right) , = P/A_{0}$$
(4)
where P = maximum load (N), and A0 = the cross-sectional area of the specimen (mm2).The tensile tests were conducted using bone-shaped specimens (Fig. 2). The specimen length was 230 (L) mm with a 20 (T) mm width of the specimen grip. The cross-sectional area of the specimen was 2 mm in the tangential direction by culm thickness in the radial direction. A strain gage type extensometer (SG25-10A, A&D, Tokyo, Japan) was used to detect the elongation in the test specimen. The specimen grip sections were attached to small boards (75 mm in length × 40 mm in width × 5 mm in thickness) and then were clamped between the metal grip of a universal testing machine (RTC-2410, A&D, Tokyo, Japan). The tensile load was applied at 1 mm min−1. The tensile strength (TS) and Young’s modulus (TM) were calculated by the following formulae:$${text{TS }}left( {{text{MPa}}} right) , = P/A_{0}$$
(5)
$${text{TM }}left( {{text{GPa}}} right) , = Delta Pl/A_{0} Delta l times {1}0^{{ – 3}}$$
(6)
where P = maximum load (N), A0 = the cross-sectional area of the specimen (mm2), ∆P = difference between upper and lower proportional limit within the range of elasticity (N), l = gauge length (mm), and ∆l = elongation of the original gauge length (mm).The moisture content and air-dry density of each specimen were measured after each mechanical testing by the oven-dry method. The moisture content and air-dry density of the specimen at testing were listed in Table S1.Statistical analysisThe statistical analyses were conducted using R software (version 4.0.3)27. To evaluate the longitudinal variations of the measured properties in each species, the y-intercept, linear, and nonlinear mixed-effects models with each measured property value as a responsible variable, the height position as a fixed effect, and site and individual culm as random effects were developed by the “lmer” function in “lme4” packages28 and the “nlme” function in the “nlme” package29. The following four full models were developed and compared:Model I (y-intercept model):$$Y_{ijk} = alpha_{{1}} + Site_{{{1}k}} + Culm_{{{1}jk}} + e_{ijk}$$
(7)
Model II (linear model):$$Y_{ijk} = , (beta_{0} + Site_{0k} + Culm_{0jk} )X_{ijk} + beta_{{1}} + Site_{{{1}k}} + Culm_{{{1}jk}} + e_{ijk}$$
(8)
Model III (logarithmic model):$$Y_{ijk} = , (gamma_{0} + Site_{0k} + Culm_{0jk} ){text{ ln }}left( {X_{ijk} } right) + gamma_{{1}} + Site_{{{1}k}} + Culm_{{{1}jk}} + e_{ijk}$$
(9)
Model IV (quadratic model):$$begin{gathered} Y_{ijk} = , (zeta_{0} + Site_{0k} + Culm_{0jk} )X_{ijk}^{{2}} + , (zeta_{{1}} + Site_{{{1}k}} + Culm_{{{1}jk}} )X_{ijk} hfill \ + zeta_{{2}} + Site_{{{2}k}} + Culm_{{{2}jk}} + e_{ijk} hfill \ end{gathered}$$
(10)

where Yijk is measured property at the ith height position from the jth individual culm within the kth site, Xijk is the ith height position from the jth individual culm within the kth site, α1, β0, β1, γ0, γ1, ζ0, ζ1, and ζ2 are the fixed effects, Site0k, Site1k, and Site2k are the random effect at the site level, Culm0jk, Culm1jk, and Culm2jk are the random effects at the individual culm level, and eijk is residual. Total 36 derived models (three y-intercept models, 15 linear models, nine logarithmic models, and nine quadratic models) were developed. The model selection was conducted using the Akaike information criterion30. The model with the minimum AIC value was regarded as the most parsimonious model among developed models. In addition, the differences in AIC (ΔAIC) ≤ 2 indicate no significant differences between models, and a simpler model with fewer parameters is preferred31. To evaluate the longitudinal variation, estimated values of each property was calculated at 0.1 m interval from 2.0 to 8.0 m above the ground using fixed-effect parameters of the selected models. Mean value and standard deviation were obtained from the estimated values from 2.0 to 8.0 m in each property. In addition, the coefficient of variation was also calculated from the mean value and standard deviation. The longitudinal variation patterns were classified into four types (Types A to D) based on the model selection (Fig. 4). Although model II to IV was selected, longitudinal variation with the coefficient of variation less than 3.0% was regarded as stable (Type A in Fig. 4).Figure 4Classification of longitudinal variation of bamboo culm property. Note: Lines or curves indicate formulae with fixed-effect parameters in the selected mixed-effect model for explaining longitudinal variation (Tables 3, 4, 5). Coefficient of variation calculated from mean values and standard deviation from 2 to 8 m above the ground estimated by fixed-effect parameters values less than 3.0% is regard as stable variation (Type A), even in selected model is Model II to IV.Full size imageGeographic variations in each bamboo property were estimated by evaluating the variance component of sites and culms as random effects by using the intercept-only linear mixed-effects model. The full model is described as follows:$$Y_{ijk} = mu + Site_{k} + Culm_{jk} + e_{ijk}$$
(11)
where Yijk is the bamboo property at the ith height position of the jth individual culm within kth site, μ is the model intercept or grand mean, Sitek is the random effect of the kth site, Culmjk is random effect of jth individual culm within kth site, and eijk is the residual. The contribution of each level of variation was calculated as a percentage of the total random variation in the best model32,33. More

Global Forest Resources Assessment 2020 (FAO, 2020).Taubert, F. et al. Global patterns of tropical forest fragmentation. Nature 554, 519–522 (2018).ADS 
CAS 

Google Scholar 
Vancutsem, C. et al. Long-term (1990–2019) monitoring of forest cover changes in the humid tropics. Sci. Adv. 7, eabe1603 (2021).ADS 

Google Scholar 
Foley, J. A. et al. Amazonia revealed: Forest degradation and loss of ecosystem goods and services in the Amazon Basin. Front. Ecol. Environ. 5, 25–32 (2007).
Google Scholar 
Barlow, J. et al. The future of hyperdiverse tropical ecosystems. Nature 559, 517–526 (2018).ADS 
CAS 

Google Scholar 
Brandon, K. Ecosystem services from tropical forests: Review of current science. SSRN J. https://doi.org/10.2139/ssrn.2622749 (2014). Article 

Google Scholar 
Indarto, J. & Mutaqin, D. J. An overview of theoretical and empirical studies on deforestation. MPRA. Paper No. 70178 (2016).Geist, H. J. & Lambin, E. F. Proximate causes and underlying driving forces of tropical deforestation: Tropical forests are disappearing as the result of many pressures, both local and regional, acting in various combinations in different geographical locations. Bioscience 52, 143–150 (2002).
Google Scholar 
Angelsen, A. & Kaimowitz, D. Rethinking the causes of deforestation: Lessons from economic models. World Bank Res. Obs. 14, 73–98 (1999).CAS 

Google Scholar 
Contreras-Hermosilla, A. The Underlying Causes of Forest Decline (Center for International Forestry Research, 2000).
Google Scholar 
Turner, B. L. et al. Two types of global environmental change: Definitional and spatial-scale issues in their human dimensions. Glob. Environ. Change 1, 14–22 (1990).
Google Scholar 
Meyer, W. B. & Turner, B. L. Human population growth and global land-use/cover change. Ann. Rev. Ecol. Syst. 2, 39–61 (1992).
Google Scholar 
Miyamoto, M., Mohd Parid, M., Noor Aini, Z. & Michinaka, T. Proximate and underlying causes of forest cover change in Peninsular Malaysia. For. Policy Econ. 44, 18–25 (2014).
Google Scholar 
Lim, C. L., Prescott, G. W., De Alban, J. D. T., Ziegler, A. D. & Webb, E. L. Untangling the proximate causes and underlying drivers of deforestation and forest degradation in Myanmar. Conserv. Biol. 31, 1362–1372 (2017).
Google Scholar 
Carodenuto, S. et al. A methodological framework for assessing agents, proximate drivers and underlying causes of deforestation: Field test results from southern cameroon. Forests 6, 203–224 (2015).
Google Scholar 
Curtis, P. G., Slay, C. M., Harris, N. L., Tyukavina, A. & Hansen, M. C. Classifying drivers of global forest loss. Science 361, 1108–1111 (2018).ADS 
CAS 

Google Scholar 
Hosonuma, N. et al. An assessment of deforestation and forest degradation drivers in developing countries. Environ. Res. Lett. 7, 044009 (2012).ADS 

Google Scholar 
Köthke, M., Leischner, B. & Elsasser, P. Uniform global deforestation patterns—An empirical analysis. For. Policy Econ. 28, 23–37 (2013).
Google Scholar 
Busch, J. & Ferretti-Gallon, K. What drives deforestation and what stops it? A meta-analysis. Rev. Environ. Econ. Policy 11, 3–23 (2017).
Google Scholar 
Ferrer Velasco, R. F., Köthke, M., Lippe, M. & Günter, S. Scale and context dependency of deforestation drivers: Insights from spatial econometrics in the tropics. PLoS One 15, e0226830 (2020).CAS 

Google Scholar 
Lambin, E. F. et al. Effectiveness and synergies of policy instruments for land use governance in tropical regions. Glob. Environ. Change 28, 129–140 (2014).
Google Scholar 
Börner, J., Schulz, D., Wunder, S. & Pfaff, A. The effectiveness of forest conservation policies and programs. Ann. Rev. Resour. Econ. 12, 45–64 (2020).
Google Scholar 
Bemelmans-Videc, M.-L., Rist, R. C. & Vedung, E. Carrots, Sticks & Sermons: Policy Instruments and their Evaluation (Transaction Publishers, 1998).
Google Scholar 
Seymour, F. & Harris, N. L. Reducing tropical deforestation. Science 365, 756–757 (2019).ADS 
CAS 

Google Scholar 
Lambin, E. F. et al. The role of supply-chain initiatives in reducing deforestation. Nat. Clim. Change 8, 109–116 (2018).ADS 

Google Scholar 
Wolff, S. & Schweinle, J. Effectiveness and economic viability of forest certification: A systematic review. Forests 13, 798 (2022).
Google Scholar 
Müller, R., Pistorius, T., Rohde, S., Gerold, G. & Pacheco, P. Policy options to reduce deforestation based on a systematic analysis of drivers and agents in lowland Bolivia. Land Use Policy 30, 895–907 (2013).
Google Scholar 
Tegegne, Y. T., Lindner, M., Fobissie, K. & Kanninen, M. Evolution of drivers of deforestation and forest degradation in the Congo Basin forests: Exploring possible policy options to address forest loss. Land Use Policy 51, 312–324 (2016).
Google Scholar 
Hoffmann, C., García Márquez, J. R. & Krueger, T. A local perspective on drivers and measures to slow deforestation in the Andean-Amazonian foothills of Colombia. Land Use Policy 77, 379–391 (2018).
Google Scholar 
Henders, S., Ostwald, M., Verendel, V. & Ibisch, P. Do national strategies under the UN biodiversity and climate conventions address agricultural commodity consumption as deforestation driver?. Land Use Policy 70, 580–590 (2018).
Google Scholar 
Salvini, G. et al. How countries link REDD+ interventions to drivers in their readiness plans: implications for monitoring systems. Environ. Res. Lett. 9, 074004 (2014).ADS 

Google Scholar 
Bos, A. B. et al. Integrated assessment of deforestation drivers and their alignment with subnational climate change mitigation efforts. Environ. Sci. Policy 114, 352–365 (2020).CAS 

Google Scholar 
Fritz, S. et al. A continental assessment of the drivers of tropical deforestation with a focus on protected areas. Front. Conserv. Sci. https://doi.org/10.3389/fcosc.2022.830248 (2022).Article 

Google Scholar 
Lawrence, D. & Vandecar, K. Effects of tropical deforestation on climate and agriculture. Nat. Clim. Change 5, 27–36 (2015).ADS 

Google Scholar 
Fedele, G., Locatelli, B., Djoudi, H. & Colloff, M. J. Reducing risks by transforming landscapes: Cross-scale effects of land-use changes on ecosystem services. PLoS One 13, e0195895 (2018).
Google Scholar 
Yackulic, C. B. et al. Biophysical and socioeconomic factors associated with forest transitions at multiple spatial and temporal scales. Ecol. Soc. https://doi.org/10.5751/ES-04275-160315 (2011).Article 

Google Scholar 
Loran, C., Ginzler, C. & Bürgi, M. Evaluating forest transition based on a multi-scale approach: Forest area dynamics in Switzerland 1850–2000. Reg. Environ. Change 16, 1807–1818 (2016).
Google Scholar 
Moonen, P. C. et al. Actor-based identification of deforestation drivers paves the road to effective REDD+in DR Congo. Land Use Policy 58, 123–132 (2016).
Google Scholar 
Strassburg, B. The tragedy of the tropics: A dynamic, cross-scale analysis of deforestation incentives. Working Paper—Centre for Social and Economic Research on the Global Environment No. 07-02 (2007).López-Carr, D. et al. Space versus place in complex human–natural systems: Spatial and multi-level models of tropical land use and cover change (LUCC) in Guatemala. Ecol. Model. 229, 64–75 (2012).
Google Scholar 
Hoang, N. T. & Kanemoto, K. Mapping the deforestation footprint of nations reveals growing threat to tropical forests. Nat. Ecol. Evol. 5, 845–853 (2021).
Google Scholar 
Pendrill, F. et al. Agricultural and forestry trade drives large share of tropical deforestation emissions. Glob. Environ. Change 56, 1–10 (2019).
Google Scholar 
Ferrer Velasco, R. et al. Towards accurate mapping of forest in tropical landscapes: A comparison of datasets on how forest transition matters. Remote Sens. Environ. 274, 112997 (2022).ADS 

Google Scholar 
Jayathilake, H. M., Prescott, G. W., Carrasco, L. R., Rao, M. & Symes, W. S. Drivers of deforestation and degradation for 28 tropical conservation landscapes. Ambio 50, 215–228 (2021).
Google Scholar 
Minang, P. A. et al. REDD+Readiness progress across countries: Time for reconsideration. Clim. Policy 14, 685–708 (2014).
Google Scholar 
Current pledges | Bonn challenge. https://www.bonnchallenge.org/pledges. Accessed: 15th August 2022.Nansikombi, H. et al. Can de facto governance influence deforestation drivers in the Zambian Miombo?. For. Policy Econ. 120, 102309 (2020).
Google Scholar 
Sullivan, A., York, A., White, D., Hall, S. & Yabiku, S. D. Jure versus de facto institutions: Trust, information, and collective efforts to manage the invasive mile-a-minute weed (Mikania micrantha). Int. J. Commons 11, 171–199 (2017).
Google Scholar 
Busch, J. & Amarjargal, O. Authority of second-tier governments to reduce deforestation in 30 tropical countries. Front. For. Glob. Change https://doi.org/10.3389/ffgc.2020.00001 (2020).Article 

Google Scholar 
Sandström, C., Eckerberg, K. & Raitio, K. Studying conflicts, proposing solutions—Towards multi-level approaches to the analyses of forest conflicts. For. Policy Econ. 33, 123–127 (2013).
Google Scholar 
Hoogstra-Klein, M. A., Permadi, D. B. & Yasmi, Y. The value of cultural theory for participatory processes in natural resource management. For. Policy Econ. 20, 99–106 (2012).
Google Scholar 
de Jong, W., Ruiz, S. & Becker, M. Conflicts and communal forest management in northern Bolivia. For. Policy Econ. 8, 447–457 (2006).
Google Scholar 
Eckerberg, K. & Sandström, C. Forest conflicts: A growing research field. For. Policy Econ. 33, 3–7 (2013).
Google Scholar 
Sierra, R., Calva, O. & Guevara, A. La Deforestación en el Ecuador, 1990–2018. Factores promotores y tendencias recientes, 216 (2021).Wasserstrom, R. & Southgate, D. Deforestation, agrarian reform and oil development in Ecuador, 1964–1994. Nat. Resour. 04, 31 (2013).
Google Scholar 
Wiebe, P. C., Zhunusova, E., Lippe, M., Ferrer Velasco, R. & Günter, S. What is the contribution of forest-related income to rural livelihood strategies in the Philippines’ remaining forested landscapes?. For. Policy Econ. 135, 102658 (2022).
Google Scholar 
Le, H. D., Smith, C. & Herbohn, J. What drives the success of reforestation projects in tropical developing countries? The case of the Philippines. Glob. Environ. Change 24, 334–348 (2014).
Google Scholar 
Carandang, A. P. et al. Analysis of key drivers of deforestation and forest degradation in the Philippines. Deutsche Gesellschaft für Internationale Zusammenarbeit (GIZ) (2013).Phiri, D., Morgenroth, J. & Xu, C. Four decades of land cover and forest connectivity study in Zambia—An object-based image analysis approach. Int. J. Appl. Earth Obs. Geoinf. 79, 97–109 (2019).ADS 

Google Scholar 
Nansikombi, H., Fischer, R., Kabwe, G. & Günter, S. Exploring patterns of forest governance quality: Insights from forest frontier communities in Zambia’s Miombo ecoregion. Land Use Policy 99, 104866 (2020).
Google Scholar 
Zhang, H., Wang, P. & Wood, J. Does institutional quality matter for the nexus between environmental quality and economic growth?: A tropics perspective. In Business, Industry, and Trade in the Tropics (eds Wood, J. et al.) (Routledge, 2022).
Google Scholar 
Reed, J., Van Vianen, J., Deakin, E. L., Barlow, J. & Sunderland, T. Integrated landscape approaches to managing social and environmental issues in the tropics: Learning from the past to guide the future. Glob. Change Biol. 22, 2540–2554 (2016).ADS 

Google Scholar 
Fischer, R. et al. Interplay of governance elements and their effects on deforestation in tropical landscapes: Quantitative insights from Ecuador. World Dev. 148, 105665 (2021).
Google Scholar 
Torres, B., Vasco, C., Günter, S. & Knoke, T. Determinants of agricultural diversification in a hotspot area: Evidence from colonist and indigenous communities in the Sumaco biosphere reserve Ecuadorian Amazon. Sustainability 10, 1432 (2018).
Google Scholar 
Ojeda Luna, T., Zhunusova, E., Günter, S. & Dieter, M. Measuring forest and agricultural income in the Ecuadorian lowland rainforest frontiers: Do deforestation and conservation strategies matter?. For. Policy Econ. 111, 102034 (2020).
Google Scholar 
Kazungu, M. et al. Effects of household-level attributes and agricultural land-use on deforestation patterns along a forest transition gradient in the Miombo landscapes Zambia. Ecol. Econ. 186, 107070 (2021).
Google Scholar 
Kleemann, J. et al. Deforestation in continental ecuador with a focus on protected areas. Land 11, 268 (2022).
Google Scholar 
Mulenga, M. M. & Roos, A. Assessing the awareness and adoptability of pellet cookstoves for low-income households in Lusaka, Zambia. J. Energy South. Afr. 32, 52–61 (2021).
Google Scholar 
Eguiguren, P., Ojeda Luna, T., Torres, B., Lippe, M. & Günter, S. Ecosystem service multifunctionality: Decline and recovery pathways in the amazon and chocó lowland rainforests. Sustainability 12, 7786 (2020).CAS 

Google Scholar 
Vasco, C., Torres, B., Pacheco, P. & Griess, V. The socioeconomic determinants of legal and illegal smallholder logging: Evidence from the Ecuadorian Amazon. For. Policy Econ. 78, 133–140 (2017).
Google Scholar 
van der Ploeg, J., van Weerd, M., Masipiqueña, A. B. & Persoon, G. A. Illegal logging in the Northern Sierra Madre Natural Park, the Philippines. Conserv. Soc. 9, 202–215 (2011).
Google Scholar 
Liu, D. S., Iverson, L. R. & Brown, S. Rates and patterns of deforestation in the Philippines: Application of geographic information system analysis. For. Ecol. Manag. 57, 1–16 (1993).
Google Scholar 
Boquet, Y. Environmental challenges in the Philippines. In The Philippine Archipelago (ed. Boquet, Y.) 779–829 (Springer International Publishing, 2017). https://doi.org/10.1007/978-3-319-51926-5_22.Chapter 

Google Scholar 
MAGAP. ATPA: Reconversión Agro productiva Sostenible en la Amazonía Ecuatoriana (2014).Jones, K. W. et al. Forest conservation incentives and deforestation in the Ecuadorian Amazon. Environ. Conserv. 44, 56–65 (2017).
Google Scholar 
Lindsey, P. A. et al. Underperformance of African protected area networks and the case for new conservation models: Insights from Zambia. PLoS One 9, e94109 (2014).ADS 

Google Scholar 
Fischer, R. et al. Effectiveness of policy instrument mixes for forest conservation in the tropics – a stakeholder perspective from Ecuador, the Philippines and Zambia. Land Use Policy https://doi.org/10.1016/j.landusepol.2023.106546 (2022).Article 

Google Scholar 
Gurney, G. G. et al. Biodiversity needs every tool in the box: Use OECMs. Nature 595, 646–649 (2021).ADS 
CAS 

Google Scholar 
Maxwell, S. L. et al. Area-based conservation in the twenty-first century. Nature 586, 217–227 (2020).ADS 
CAS 

Google Scholar 
Priebe, J. et al. Transformative change in context—Stakeholders’ understandings of leverage at the forest–climate nexus. Sustain. Sci. 17, 1921–1938 (2022).
Google Scholar 
Höhl, M. et al. Forest landscape restoration—What generates failure and success?. Forests 11, 938 (2020).
Google Scholar 
Köthke, M., Ahimbisibwe, V. & Lippe, M. The evidence base on the environmental, economic and social outcomes of agroforestry is patchy—An evidence review map. Front. Environ. Sci. https://doi.org/10.3389/fenvs.2022.925477 (2022).Article 

Google Scholar 
Fischer, R., Giessen, L. & Günter, S. Governance effects on deforestation in the tropics: A review of the evidence. Environ. Sci. Policy 105, 84–101 (2020).
Google Scholar 
Bare, M., Kauffman, C. & Miller, D. C. Assessing the impact of international conservation aid on deforestation in sub-Saharan Africa. Environ. Res. Lett. 10, 125010 (2015).ADS 

Google Scholar 
Vuohelainen, A. J., Coad, L., Marthews, T. R., Malhi, Y. & Killeen, T. J. The effectiveness of contrasting protected areas in preventing deforestation in Madre de Dios. Peru. Environ. Manag. 50, 645–663 (2012).ADS 

Google Scholar 
Hull, V. & Liu, J. Telecoupling: A new frontier for global sustainability. Ecol. Soc. 23, 41 (2018).
Google Scholar 
Aitchison, J. The statistical analysis of compositional data. J. Roy. Stat. Soc. 44, 139–160 (1982).MathSciNet 
MATH 

Google Scholar 
Norman, G. Likert scales, levels of measurement and the “laws” of statistics. Adv. Health Sci. Educ. 15, 625–632 (2010).
Google Scholar 
Day, M., Gumbo, D., Moombe, K. B., Wijaya, A. & Sunderland, T. Zambia Country Profile: Monitoring, Reporting and Verification for REDD+ Vol. 113 (CIFOR, 2014).
Google Scholar 
Piotrowski, M. Nearing the tipping point. Drivers of Deforestation in the Amazon Region (2019).Sarker, P. K., Fischer, R., Tamayo, F., Navarrete, B. T. & Günter, S. Analyzing forest policy mixes based on the coherence of policies and the consistency of legislative policy instruments: A case study from Ecuador. For. Policy Econ. 144, 102838 (2022).
Google Scholar 
Likert, R. A technique for the measurement of attitudes. Arch. Psychol. 22(140), 55–55 (1932).
Google Scholar 
Altinsoy, M. et al. Ambulatory ECG monitoring for syncope and collapse in United States, Europe, and Japan: The patients’ viewpoint. J. Arrhythm. 37, 1023–1030 (2021).
Google Scholar 
R Core Team. R: A language and environment for statistical computing. https://www.R-project.org/. (R Foundation for Statistical Computing, Vienna, Austria, 2022).Kassambara, A. rstatix: Pipe-friendly framework for basic statistical tests. R package version 0.7.0 (2021).Kassambara, A. & Mundt, F. factoextra: Extract and visualize the results of multivariate data analyses. R package version 1.0.7 (2020).Komsta, L. & Novometsky, F. moments: Moments, cumulants, skewness, kurtosis and related tests. R package version 0.14.1 (2022).Zhang, Y., Zhou, M. & Shao, Y. mvnormalTest: Powerful tests for multivariate normality. R package version 1.0.0 (2020).Kassambara, A. ggpubr: ‘ggplot2’ Based Publication Ready Plots. R package version 0.4.0 (2020).Wickham, H. et al. Welcome to the Tidyverse. JOSS 4, 1686 (2019).ADS 

Google Scholar 
Bache, S. M. & Wickham, H. magrittr: A Forward-Pipe Operator for R. R package version 2.0.3 (2022).Ushey, K., Allaire, J., Wickham, H. & Ritchie, G. rstudioapi: Safely Access the RStudio API. R package version 0.13 (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer International Publishing, 2016). https://doi.org/10.1007/978-3-319-24277-4.Book 
MATH 

Google Scholar 
Wilkins, D. treemapify: Draw Treemaps in ‘ggplot2’. R package version 2.5.5 (2021).Shapiro, S. S. & Wilk, M. B. An analysis of variance test for normality (complete samples). Biometrika 52, 591–611 (1965).MathSciNet 
MATH 

Google Scholar 
Mardia, K. V. Measures of multivariate skewness and kurtosis with applications. Biometrika 57, 519–530 (1970).MathSciNet 
MATH 

Google Scholar 
Kruskal, W. H. & Wallis, W. A. Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621 (1952).MATH 

Google Scholar 
Dunn, O. J. Multiple comparisons using rank sums. Technometrics 6, 241–252 (1964).
Google Scholar 
Conover, W. J. & Iman, R. L. Rank transformations as a bridge between parametric and nonparametric statistics. Am. Stat. 35, 124–129 (1981).MATH 

Google Scholar 
Student,. The probable error of a mean. Biometrika 6, 1–25 (1908).MATH 

Google Scholar 
Tukey, J. W. Comparing individual means in the analysis of variance. Biometrics 5, 99–114 (1949).MathSciNet 
CAS 

Google Scholar 
Jolliffe, I. T. Principal Component Analysis (Springer, 2002).MATH 

Google Scholar  More

Guidetti, R. & Bertolani, R. B. Tardigrade taxonomy: An updated check list of the taxa and a list of characters for their identification. Zootaxa 845, 1–46. https://doi.org/10.11646/zootaxa.845.1.1 (2005).Article 

Google Scholar 
Degma, P. & Guidetti, R. Notes to the current checklist of Tardigrada. Zootaxa 1579, 41–53. https://doi.org/10.11646/zootaxa.1579.1.2 (2007)Article 

Google Scholar 
Vicente, F. & Bertolani, R. Considerations on the taxonomy of the phylum Tardigrada. Zootaxa 3626, 245–248. https://doi.org/10.11646/zootaxa.3626.2.2 (2013).Article 

Google Scholar 
Degma, P. & Guidetti, R. Actual checklist of Tardigrada species. (Version 41: Edition: 16-05-2022). (2009–2022).Ramazzotti, G. & Maucci, W. Il phylum Tardigrada. III edizione riveduta e aggiornata. Mem. Ist. Ital. Idrobiol. 41, 1–1012 (1983).
Google Scholar 
Beasley, C. W. The phylum Tardigrada. in English Translation P. 3rd edn (eds Ramazzotti, G. & Maucci, W.) 1–1014 (Abilene, USA, 1995).Nelson, D. R., Guidetti, R., Rebecchi, L., Kaczmarek, Ł. & McInnes, S. Phylum Tardigrada. in Thorp and Covich’s Freshwater Invertebrates 505–522 (Elsevier, 2020). https://doi.org/10.1016/B978-0-12-804225-0.00015-0.Da Cunha, A. X. & do Nascimento-Ribeiro, F. A fauna de Tardígrados da Ilha da Madeira. Mem. Estud. Mus. Zool. Univ. Coimbra 1–24 (1962).Fontoura, P., Pilato, G. & Lisi, O. Tardigrada from Santo Antão Island (Archipelago of Cape Verde, West Africa) with the description of a new species. Zootaxa 2838, 30–40. https://doi.org/10.11646/zootaxa.2838.1.2 (2011).Article 

Google Scholar 
Gąsiorek, P., Vončina, K. & Michalczyk, Ł. Echiniscus testudo (Doyère, 1840) in New Zealand: Anthropogenic dispersal or evidence for the ‘Everything is Everywhere’ hypothesis?. N. Z. J. Zool. 46, 174–181. https://doi.org/10.1080/03014223.2018.1503607 (2019).Article 

Google Scholar 
Guidetti, R., Schill, R. O., Bertolani, R., Dandekar, T. & Wolf, M. New molecular data for tardigrade phylogeny, with the erection of Paramacrobiotus gen. nov. J. Zool. Syst. Evol. 47, 315–321. https://doi.org/10.1111/j.1439-0469.2009.00526.x (2009).Article 

Google Scholar 
Kaczmarek, Ł, Gawlak, M., Bartels, P. J., Nelson, D. R. & Roszkowska, M. Revision of the genus Paramacrobiotus Guidetti et al., 2009 with the description of a new species, re-descriptions and a key. Ann. Zool. 67, 627–656. https://doi.org/10.3161/00034541ANZ2017.67.4.001 (2017).Article 

Google Scholar 
Marley, N. J. et al. A clarification for the subgenera of Paramacrobiotus Guidetti, Schill, Bertolani, Dandekar and Wolf, 2009, with respect to the International Code of Zoological Nomenclature. Zootaxa 4407, 130–134. https://doi.org/10.11646/zootaxa.4407.1.9 (2018).Article 
CAS 

Google Scholar 
Guidetti, R., Cesari, M., Bertolani, R., Altiero, T. & Rebecchi, L. High diversity in species, reproductive modes and distribution within the Paramacrobiotus richtersi complex (Eutardigrada, Macrobiotidae). Zool. Lett. 5, 1–28. https://doi.org/10.1186/s40851-018-0113-z (2019).Article 

Google Scholar 
Stec, D., Krzywański, Ł, Zawierucha, K. & Michalczyk, Ł. Untangling systematics of the Paramacrobiotus areolatus species complex by an integrative redescription of the nominal species for the group, with multilocus phylogeny and species delineation in the genus Paramacrobiotus. Zool. J. Linn. Soc. 188, 694–716. https://doi.org/10.1093/zoolinnean/zlz163 (2020).Article 

Google Scholar 
Murray, J. Scottish Tardigrada, a review of our present knowledge. Ann. Scot. Nat. Hist. 78, 88–95 (1911).
Google Scholar 
Murray, J. XXV.—Arctic Tardigrada, collected by Wm. S. Bruce. Trans. R. Soc. Edinb. 45, 669–681 (1907).Article 

Google Scholar 
Ramazzotti, G. Tre nouve specie di Tardigradi ed altre specie poco comuni. Atti Soc. Nat. Milano 10, 284–291 (1956).
Google Scholar 
Schill, R. O., Förster, F., Dandekar, T. & Wolf, M. Using compensatory base change analysis of internal transcribed spacer 2 secondary structures to identify three new species in Paramacrobiotus (Tardigrada). Org. Divers. Evol. 10, 287–296. https://doi.org/10.1007/s13127-010-0025-z (2010).Article 

Google Scholar 
Kaczmarek, Ł et al. Integrative description of bisexual Paramacrobiotus experimentalis sp. Nov. (Macrobiotidae) from republic of Madagascar (Africa) with microbiome analysis. Mol. Phylogenet. Evol. 145, 106730. https://doi.org/10.1016/j.ympev.2019.106730 (2020).Article 

Google Scholar 
Bertolani, R. Partenogenesi geografica triploide in un Tardigrado (Macrobiotus richtersi). Rend. Acc. Naz. Lincei. Ser. 8, 487–489 (1971).
Google Scholar 
Bertolani, R. Sex ratio and geographic parthenogenesis in Macrobioutus (Tardigrada). Experientia 28, 94–95. https://doi.org/10.1007/BF01928285 (1972).Article 

Google Scholar 
Bertolani, R. L. partenogenesi nei Tardigradi. Boll. Zool. 39, 577–581. https://doi.org/10.1080/11250007209431414 (1972).Article 

Google Scholar 
Bertolani, R. Cytology and Reproductive Mechanisms in Tardigrades. I. 93–114 (East Tennesse State University Press, Johnson City, 1982).
Google Scholar 
Lemloh, M., Brümmer, F. & Schill, R. O. Life-history traits of the bisexual tardigrades Paramacrobiotus tonollii and Macrobiotus sapiens. J. Zool. Syst. Evol. Res. 49, 58–61. https://doi.org/10.1111/j.1439-0469.2010.00599.x (2011).Article 

Google Scholar 
Guil, N. & Giribet, G. A comprehensive molecular phylogeny of tardigrades-adding genes and taxa to a poorly resolved phylum-level phylogeny. Cladistics 28, 21–49. https://doi.org/10.1111/j.1096-0031.2011.00364.x (2012).Article 

Google Scholar 
Kosztyła, P. et al. Experimental taxonomy confirms the environmental stability of morphometric traits in a taxonomically challenging group of microinvertebrates. Zool. J. Linn. Soc. 178, 765–775. https://doi.org/10.1111/zoj.12409 (2016).Article 

Google Scholar 
Kaczmarek, Ł et al. New records of Antarctic Tardigrada with comments on iterpopulation variability of the Paramacrobiotus fairbanksi Schill, Förster, Dandekar and Wolf, 2010. Diversity 12, 108. https://doi.org/10.3390/d12030108 (2020).Article 

Google Scholar 
Stec, D., Vecchi, M., Calhim, S. & Michalczyk, Ł. New multilocus phylogeny reorganises the family Macrobiotidae (Eutardigrada) and unveils complex morphological evolution of the Macrobiotus hufelandi group. Mol. Phylogenet. Evol. 160, 106987. https://doi.org/10.1016/j.ympev.2020.106987 (2021).Article 

Google Scholar 
Stec, D., Smolak, R., Kaczmarek, Ł & Michalczyk, Ł. An integrative description of Macrobiotus paulinae sp. Nov. (Tardigrada: Eutardigrada: Macrobiotidae: hufelandi group) from Kenya. Zootaxa 4052, 501–526. https://doi.org/10.11646/zootaxa.4052.5.1 (2015).Article 

Google Scholar 
Bryce, D. On some moss-dwelling Cathypnadae; with descriptions of five new species. Sci. Gossip Lond. 28, 271–275 (1892).
Google Scholar 
Casquet, J., Thebaud, C. & Gillespie, R. G. Chelex without boiling, a rapid and easy technique to obtain stable amplifiable DNA from small amounts of ethanol-stored spiders. Mol. Ecol. Resour. 12(1), 136–141. https://doi.org/10.1111/j.1755-0998.2011.03073.x (2012).Article 
CAS 

Google Scholar 
Stec, D., Kristensen, R. M. & Michalczyk, Ł. An integrative description of Minibiotus ioculator sp. nov. from the Republic of South Africa with notes on Minibiotus pentannulatus Londoño et al., 2017 (Tardigrada: Macrobiotidae). Zool. Anz. 286, 117–134. https://doi.org/10.1016/j.jcz.2020.03.007 (2020).Article 

Google Scholar 
Stec, D., Zawierucha, K. & Michalczyk, Ł. An integrative description of Ramazzottius subanomalus (Biserov, 1985 (Tardigrada) from Poland. Zootaxa 4300, 403–420. https://doi.org/10.11646/zootaxa.4300.3.4 (2017).Article 

Google Scholar 
Mironov, S. V., Dabert, J. & Dabert, M. A new feather mite species of the genus Proctophyllodes Robin, 1877 (Astigmata: Proctophyllodidae) from the Long-tailed Tit Aegithalos caudatus (Passeriformes: Aegithalidae)—Morphological description with DNA barcode data. Zootaxa 3253, 54–61. https://doi.org/10.11646/zootaxa.3253.1.2 (2012).Article 

Google Scholar 
White, T. J., Bruns, T., Lee, S. & Taylor, J. PCR Protocols: A Guide to Methods and Application 315–322. https://doi.org/10.1016/B978-0-12-372180-8.50042-1 (Academic Press, 1990).Book 

Google Scholar 
Folmer, O., Black, M., Hoeh, W., Lutz, R. & Vrijenhoek, R. Phylogenetic uncertainty. Mol. Mar. Biol. Biotechnol. 3, 294–299 (1994).CAS 

Google Scholar 
Vecchi, M. & Stec, D. Integrative descriptions of two new Macrobiotus species (Tardigrada, Eutardigrada, Macrobiotidae) from Mississippi (USA) and Crete (Greece). ZSE 97, 281–306. https://doi.org/10.3897/zse.97.65280 (2021).Article 

Google Scholar 
Thulin, G. Über die phylogenie und das system der. Hereditas 11, 207–266. https://doi.org/10.1111/j.1601-5223.1928.tb02488.x (1928).Article 

Google Scholar 
Stec, D. Mesobiotus datanlanicus sp. nov., a new tardigrade species (Macrobiotidae: Mesobiotus harmsworthi group) from Lâm Đồng Province in Vietnam. Zootaxa 4679, 164–180. https://doi.org/10.11646/zootaxa.4679.1.10 (2019).Article 

Google Scholar 
Katoh, K. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. NAR 30, 3059–3066. https://doi.org/10.1093/nar/gkf436 (2002).Article 
CAS 

Google Scholar 
Katoh, K. & Toh, H. Recent developments in the MAFFT multiple sequence alignment program. Brief. Bioinform. 9, 286–298. https://doi.org/10.1093/bib/bbn013 (2008).Article 
CAS 

Google Scholar 
Vaidya, G., Lohman, D. J. & Meier, R. SequenceMatrix: Concatenation software for the fast assembly of multi-gene datasets with character set and codon information. Cladistics 27, 171–180. https://doi.org/10.1111/j.1096-0031.2010.00329.x (2011).Article 

Google Scholar 
Lanfear, R., Calcott, B., Ho, S. Y. & Guindon, S. PartitionFinder: Combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29(6), 1695–1701. https://doi.org/10.1093/molbev/mss020 (2012).Article 
CAS 

Google Scholar 
Xia, X., Xie, Z., Salemi, M., Chen, L. & Wang, Y. An index of substitution saturation and its application. Mol. Phylogenet. Evol. 26, 1–7. https://doi.org/10.1016/S1055-7903(02)00326-3 (2003).Article 
CAS 

Google Scholar 
Xia, X. & Lemey, P. Assessing substitution saturation with DAMBE. In The Phylogenetic Handbook (eds Lemey, P. et al.) 615–630. https://doi.org/10.1017/CBO9780511819049.022 (Cambridge University Press, 2009).Chapter 

Google Scholar 
Ronquist, F. & Huelsenbeck, J. P. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. https://doi.org/10.1093/bioinformatics/btg180 (2003).Article 
CAS 

Google Scholar 
Rambaut, A., Suchard, M. A., Xie, D. & Drummond, A. J. Tracer v1. 6. 2014. (2015).Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. https://doi.org/10.1093/bioinformatics/btu033 (2014).Article 
CAS 

Google Scholar 
Bandelt, H., Forster, P. & Röhl, A. Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. https://doi.org/10.1093/oxfordjournals.molbev.a026036 (1999).Article 
CAS 

Google Scholar 
Ehrenberg, C. G. Beitrag zur Bestimmung des Stationären Mikroskopischen Lebens in bis 20,000 Fuss Alpenhöhe. (1859).Guil, N. & Guidetti, R. A new species of Tardigrada (Eutardigrada: Macrobiotidae) from Iberian Peninsula and Canary Islands (Spain). Zootaxa 889, 1–11. https://doi.org/10.11646/zootaxa.889.1.1 (2005).Article 

Google Scholar 
Plate, L. H. Beiträge zur Naturgeschichte der Tardigraden. Zool. Jahrb. Abteilung Anat. Ontog. Tiere 3, 487–550. https://doi.org/10.5962/bhl.part.1265 (1889).Article 

Google Scholar 
Kaczmarek, Ł, Kayastha, P., Roszkowska, M., Gawlak, M. & Mioduchowska, M. Integrative redescription of the Minibiotus intermedius (Plate, 1888)—The type species of the genus Minibiotus R.O. Schuster, 1980. Diversity 14, 356. https://doi.org/10.3390/d14050356 (2022).Article 
CAS 

Google Scholar 
Londoño, R., Daza, A., Lisi, O. & Quiroga, S. New species of waterbear Minibiotus pentannulatus (Tardigrada: Macrobiotidae) from Colombia. Rev. Mex. Biodivers. 88, 807–814. https://doi.org/10.1016/j.rmb.2017.10.040 (2017).Article 

Google Scholar 
Vecchi, M. et al. Macrobiotus naginae sp. nov., a new Xerophilous Tardigrade species from Rokua Sand Dunes (Finland). Zool. Stud. 61, e22 (2022).
Google Scholar 
Stec, D., Dudziak, M. & Michalczyk, Ł. Integrative descriptions of two new Macrobiotidae species (Tardigrada: Eutardigrada: Macrobiotoidea) from French Guiana and Malaysian Borneo. Zool. Stud. 59, e23 (2020).
Google Scholar 
Stec, D., Roszkowska, M., Kaczmarek, Ł & Michalczyk, Ł. Paramacrobiotus lachowskae, a new species of Tardigrada from Colombia (Eutardigrada: Parachela: Macrobiotidae). N. Z. J. Zool. 45, 43–60. https://doi.org/10.1080/03014223.2017.1354896 (2018).Article 

Google Scholar 
Sugiura, K., Matsumoto, M. & Kunieda, T. Description of a model tardigrade Paramacrobiotus metropolitanus sp. nov. (Eutardigrada) from Japan with a summary of its life history, reproduction and genomics. Zootaxa 5134, 92–112. https://doi.org/10.11646/zootaxa.5134.1.4 (2022).Article 

Google Scholar 
Tumanov, D. V. Three new species of Macrobiotus (Eutardigrada, Macrobiotidae, tenuis-group) from Tien Shan (Kirghizia) and Spitsbergen. J. Limnol. 66, 40. https://doi.org/10.4081/jlimnol.2007.s1.40 (2007).Article 

Google Scholar 
Zawierucha, K., Kolicka, M. & Kaczmarek, Ł. Re-description of the Arctic tardigrade Tenuibiotus voronkovi (Tumanov, 2007 (Eutardigrada; Macrobiotidea), with the first molecular data for the genus. Zootaxa 4196, 498. https://doi.org/10.11646/zootaxa.4196.4.2 (2016).Article 

Google Scholar 
Stec, D., Tumanov, D. T. & Kristensen, R. M. Integrative taxonomy identifies two new tardigrade species (Eutardigrada: Macrobiotidae) from Greenland. EJT 614, 1–40. https://doi.org/10.5852/ejt.2020.614 (2020).Article 

Google Scholar 
Fontaneto, D., Flot, J.-F. & Tang, C. Q. Guidelines for DNA taxonomy, with a focus on the meiofauna. Mar. Biodiv. 45, 433–451. https://doi.org/10.1007/s12526-015-0319-7 (2015).Article 

Google Scholar 
Zhang, J., Kapli, P., Pavlidis, P. & Stamatakis, A. A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29, 2869–2876. https://doi.org/10.1093/bioinformatics/btt499 (2013).Article 
CAS 

Google Scholar 
Puillandre, N., Brouillet, S. & Achaz, G. ASAP: Assemble species by automatic partitioning. Mol. Ecol. Resour. 21, 609–620. https://doi.org/10.1111/1755-0998.13281 (2021).Article 

Google Scholar 
Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T. & Calcott, B. PartitionFinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol. Biol. Evol. 34(3), 772–773. https://doi.org/10.1093/molbev/msw260 (2017).Article 
CAS 

Google Scholar 
Roszkowska, M., Stec, D., Gawlak, M. & Kaczmarek, Ł. An integrative description of a new tardigrade species Mesobiotus romani sp. nov. (Macrobiotidae: harmsworthi group) from the Ecuadorian Pacific coast. Zootaxa 4450, 550–564. https://doi.org/10.11646/zootaxa.4450.5.2 (2018).Article 

Google Scholar 
Pilato, G. & Binda, M. G. Definition of families, subfamilies, genera and subgenera of the Eutardigrada, and keys to their identification. Zootaxa 2404, 1–54. https://doi.org/10.11646/zootaxa.2404.1.1 (2010).Article 

Google Scholar 
Kaczmarek, Ł & Michalczyk, Ł. The Macrobiotus hufelandi group (Tardigrada) revisited. Zootaxa 4363, 101–123. https://doi.org/10.11646/zootaxa.4363.1.4 (2017).Article 

Google Scholar 
Michalczyk, Ł & Kaczmarek, Ł. A description of the new tardigrade Macrobiotus reinhardti (Eutardigrada: Macrobiotidae, harmsworthi group) with some remarks on the oral cavity armature within the genus Macrobiotus Schultze. Zootaxa 331, 1–24. https://doi.org/10.11646/zootaxa.331.1.1 (2003).Article 

Google Scholar 
Kaczmarek, Ł, Cytan, J., Zawierucha, K., Diduszko, D. & Michalczyk, Ł. Tardigrades from Peru (South America), with descriptions of three new species of Parachela. Zootaxa 3790, 357–379. https://doi.org/10.11646/zootaxa.3790.2.5 (2014).Article 

Google Scholar 
Kiosya, Y., Pogwizd, J., Matsko, Y., Vecchi, M. & Stec, D. Phylogenetic position of two Macrobiotus species with a revisional note on Macrobiotus sottilei Pilato, Kiosya, Lisi & Sabella, 2012 (Tardigrada: Eutardigrada: Macrobiotidae). Zootaxa 4933, 113–135. https://doi.org/10.11646/zootaxa.4933.1.5 (2021).Article 

Google Scholar 
Pilato, G. Analisi di nuovi caratteri nello studio degli Eutardigradi. Animalia 8, 51–57 (1981).
Google Scholar 
Michalczyk, Ł & Kaczmarek, Ł. The Tardigrada Register: a comprehensive online data repository for tardigrade taxonomy. J. Limnol. 72, e22. https://doi.org/10.4081/jlimnol.2013.s1.e22 (2013).Article 

Google Scholar 
Bertolani, R. et al. Phylogeny of Eutardigrada: New molecular data and their morphological support lead to the identification of new evolutionary lineages. Mol. Phylogenet. Evol. 76, 110–126. https://doi.org/10.1016/j.ympev.2014.03.006 (2014).Article 

Google Scholar 
Perry, E., Miller, W. R. & Kaczmarek, Ł. Recommended abbreviations for the names of genera of the phylum Tardigrada. Zootaxa 4608, 145. https://doi.org/10.11646/zootaxa.4608.1.8 (2019).Article 

Google Scholar 
Degma, P., Michalczyk, Ł & Kaczmarek, Ł. Macrobiotus derkai, a new species of Tardigrada (Eutardigrada, Macrobiotidae, huziori group) from the Colombian Andes (South America). Zootaxa 1731, 1–23. https://doi.org/10.11646/zootaxa.1731.1.1 (2008).Article 

Google Scholar 
Kaczmarek, Ł, Michalczyk, Ł & Diduszko, D. Some tardigrades from Siberia (Russia, Baikal region) with a description of Macrobiotus garynahi sp. nov. (Eutardigrada: Macrobiotidae: richtersi group). Zootaxa 1053, 35–45. https://doi.org/10.11646/zootaxa.1053.1.3 (2005).Article 

Google Scholar 
Michalczyk, Ł & Kaczmarek, Ł. Macrobiotus huziori, a new species of Tardigrada (Eutardigrada: Macrobiotidae) from Costa Rica (Central America). Zootaxa 1169, 47–59. https://doi.org/10.11646/zootaxa.1169.1.3 (2006).Article 

Google Scholar 
Michalczyk, L. & Kaczmarek, L. A new species Macrobiotus magdalenae (Tardigrada: Eutardigrada: Macrobiotidae, richtersi group) from Costa Rican rain forest (Central America). N. Z. J. Zool. 33, 189–196. https://doi.org/10.1080/03014223.2006.9518444 (2006).Article 

Google Scholar 
Michalczyk, Ł, Kaczmarek, Ł & Węglarska, B. Macrobiotus sklodowskae sp. nov. (Tardigrada: Eutardigrada: Macrobiotidae, richtersi group) from Cyprus. Zootaxa 1371, 45–56. https://doi.org/10.11646/zootaxa.1371.1.4 (2006).Article 

Google Scholar 
Tumanov, D. V. Notes on the Tardigrada of Thailand, with a description of Macrobiotus alekseevi sp. nov. (Eutardigrada, Macrobiotidae). Zootaxa 999, 1–6. https://doi.org/10.11646/zootaxa.999.1.1 (2005).Article 

Google Scholar 
Doyère, M. Memoire sur les tardigrades. Ann. Sci. Nat Zool. Ser. 2, 269–362 (1840).
Google Scholar 
Richters, F. Tardigrada. In Handbuch der Zoologie Vol. 3 (eds Kükenthal, W. & Krumbach, T.) 58–61 (Walter de Gruyter & Co., Berlin and Leipzig, 1926).
Google Scholar 
Stec, D., Cancellario, T. & Fontaneto, D. Diversification rates in Tardigrada indicate a decreasing tempo of lineage splitting regardless of reproductive mode. Org. Divers. Evol. 22(4), 965–974. https://doi.org/10.1007/s13127-022-00578-4 (2022).Article 

Google Scholar 
Dellicour, S. & Flot, J.-F. The hitchhiker’s guide to single-locus species delimitation. Mol. Ecol. Resour. 18, 1234–1246. https://doi.org/10.1111/1755-0998.12908 (2018).Article 

Google Scholar 
Magoga, G., Fontaneto, D. & Montagna, M. Factors affecting the efficiency of molecular species delimitation in a species-rich insect family. Mol. Ecol. Resour. 21, 1475–1489. https://doi.org/10.1111/1755-0998.13352 (2021).Article 
CAS 

Google Scholar 
Kaczmarek, Ł, Michalczyk, Ł & McInnes, S. J. Annotated zoogeography of non-marine Tardigrada. Part I: Central America. Zootaxa 3763, 1–62. https://doi.org/10.11646/zootaxa.3763.1.1 (2014).Article 

Google Scholar 
Kaczmarek, Ł, Michalczyk, Ł & McInnes, S. J. Annotated zoogeography of non-marine Tardigrada. Part II: South America. Zootaxa 3923, 1–107. https://doi.org/10.11646/zootaxa.3923.1.1 (2015).Article 

Google Scholar 
Kaczmarek, Ł, Michalczyk, Ł & Mcinnes, S. J. Annotated zoogeography of non-marine Tardigrada. Part III: North America and Greenland. Zootaxa 4203, 1–249. https://doi.org/10.11646/zootaxa.4203.1.1 (2016).Article 

Google Scholar 
Mcinnes, S. J., Michalczyk, Ł & Kaczmarek, Ł. Annotated zoogeography of non-marine Tardigrada. Part IV: Africa. Zootaxa 4284, 1. https://doi.org/10.11646/zootaxa.4284.1.1 (2017).Article 

Google Scholar 
Michalczyk, Ł, Kaczmarek, Ł & Mcinnes, S. J. Annotated zoogeography of non-marine Tardigrada. Part V: Australasia. Zootaxa 5107, 1–119. https://doi.org/10.11646/zootaxa.5107.1.1 (2022).Article 

Google Scholar 
Pilato, G., Claxton, S. & Binda, M. G. Tardigrades from Australia. III. Echiniscus marcusi and Macrobiotus peteri, new species of tardigrades from New South Wales. Animalia 16, 43–48 (1989).
Google Scholar 
Pilato, G., Binda, M. G. & Lisi, O. Eutardigrada from New Zealand, with descriptions of two new species. N. Z. J. Zool. 33, 49–63. https://doi.org/10.1080/03014223.2006.9518430 (2006).Article 

Google Scholar 
Bartels, P. J., Pilato, G., Lisi, O. & Nelson, D. R. Macrobiotus (Eutardigrada, Macrobiotidae) from the Great Smoky Mountains National Park, Tennessee/North Carolina, USA (North America): Two new species and six new records. Zootaxa 2022, 45–57. https://doi.org/10.11646/zootaxa.2022.1.4 (2009).Article 

Google Scholar 
Binda, M. G., Pilato, G., Moncada, E. & Napolitano, A. Some tardigrades from Central Africa with the description of two new species: Macrobiotus ragonesei and M. priviterae (Eutardigrada Macrobiotidae). Trop. Zool. 14, 233–242. https://doi.org/10.1080/03946975.2001.10531155 (2001).Article 

Google Scholar 
Pilato, G., Binda, M. G. & Lissi, O. Notes on tardigrades of the Seychelles with the description of two new species. Ital. J. Zool. 71, 171–178 (2004).Article 

Google Scholar 
Pilato, G., Binda, M. G. & Lisi, O. Three new species of eutardigrades from the Seychelles. N. Z. J. Zool. 33, 39–48. https://doi.org/10.1080/03014223.2006.9518429 (2006).Article 

Google Scholar 
Pilato, G., Binda, M. G. & Lisi, O. Notes on tardigrades of the Seychelles with the description of three new species. Ital. J. Zool. 71, 171–178. https://doi.org/10.1080/11250000409356569 (2004).Article 

Google Scholar 
Pilato, G., Binda, M. G. & Catanzaro, R. Remarks on some tardigrades of the African fauna with the description of three new species of Macrobiotus Schultze 1834. Trop. Zool. 4, 167–178. https://doi.org/10.1080/03946975.1991.10539487 (1991).Article 

Google Scholar 
Maucci, W. & Durante Pasa, M. V. Tardigradi muscicoli delle Isole Andamane. Boll. Mus. Civ. St. Nat. Verona 7, 281–291 (1980).
Google Scholar 
Iharos, G. Neuere Daten zur Kenntnis der Tardigraden-Fauna von Neuguinea. Opusc. Zool. Bp. 11, 65–73 (1973).
Google Scholar 
Binda, M. G. & Pilato, G. Macrobiotus savai and Macrobiotus humilis, two new species of tardigrades from Sri Lanka. Boll. Accad. Gioenia Sci. Nat. Catania 34, 101–111 (2001).
Google Scholar 
Pilato, G. Macrobiotus centesimus, new species of eutardigrade from the South America. Boll. Accad. Gioenia Sci. Nat. Catania 33, 97–101 (2000).
Google Scholar 
Daza, A., Caicedo, M., Lisi, O. & Quiroga, S. New records of tardigrades from Colombia with the description of Paramacrobiotus sagani sp. nov. and Doryphoribius rosanae sp. nov. Zootaxa 4362, 29–50. https://doi.org/10.11646/zootaxa.4362.1.2 (2017).Article 

Google Scholar 
Claps, M. C. & Rossi, G. C. Tardígrados de Uruguay, con descripción de dos nuevas especies (Echiniscidae, Macrobiotidae). Iheringia Sér. Zool. 83, 17–22 (1997).
Google Scholar 
Iharos, G. Neue tardigraden-arten aus ungarn (neuere beitrage zur kenntnis der tardigraden-fauna ungarns. 6.). Acta Zool. Acad. Sci. Hung. 12(1–2), 111 (1966).
Google Scholar 
Pilato, G., Kiosya, Y., Lisi, O. & Sabella, G. New records of Eutardigrada from Belarus with the description of three new species. Zootaxa 3179, 39–60. https://doi.org/10.11646/zootaxa.3179.1.2 (2012).Article 

Google Scholar 
Pasa, D. & Maucci, W. Moss Tardigrada from the Scandinavian Peninsula. in Second International Symposium on Tardigrada, Vol. 79(25). 47–85 (1979).Lisi, O., Binda, M. G. & Pilato, G. Eremobiotus ginevrae sp. nov. and Paramacrobiotus pius sp. nov., two new species of Eutardigrada. Zootaxa 4103, 344–360. https://doi.org/10.11646/zootaxa.4103.4.3 (2016).Article 

Google Scholar 
Biserov, V. I. Macrobiotus lorenae sp. n., a new species of Tardigrada (Eutardigrada Macrobiotidae) from the Russian Far East. Arthr Sel. 5, 145–149 (1996).
Google Scholar 
Biserov, V. I. Tardigrades of the Caucasus with a taxonomic analysis of genus Ramazzottius. Zool. Anz. 236, 139–159 (1997).
Google Scholar 
Morek, W. et al. Redescription of Milnesium alpigenum Ehrenberg, 1853 (Tardigrada: Apochela) and a description of Milnesium inceptum sp. nov., a tardigrade laboratory model. Zootaxa 4586(1), 35. https://doi.org/10.11646/zootaxa.4586.1.2 (2019).Article 

Google Scholar 
Morek, W., Surmacz, B., López-López, A. & Michalczyk, Ł. “Everything is not everywhere”: Time-calibrated phylogeography of the genus Milnesium (Tardigrada). Mol. Ecol. 30, 3590–3609. https://doi.org/10.1111/mec.15951 (2021).Article 

Google Scholar 
Mogle, M. J., Kimball, S. A., Miller, W. R. & McKown, R. D. Evidence of avian-mediated long-distance dispersal in American tardigrades. PeerJ 6, e5035. https://doi.org/10.7717/peerj.5035 (2018).Article 

Google Scholar 
Vuori, T., Calhim, S. & Vecchi, M. A lift in snail’s gut provides an efficient colonization route for tardigrades. Ecology 103, e3702. https://doi.org/10.1002/ecy.3702 (2022).Article 

Google Scholar 
Książkiewicz, Z. & Roszkowska, M. Experimental evidence for snails dispersing tardigrades based on Milnesium inceptum and Cepaea nemoralis species. Sci. Rep. 12(4421), 1–10. https://doi.org/10.1038/s41598-022-08265-2 (2022).Article 
ADS 
CAS 

Google Scholar  More

Study areaIn Colombia, the Amazon region represents 42.3% of the territory with an estimated area of 483,164 km2. In this area, 14% is dominated by agricultural lands, secondary vegetation and fragmented forests. Currently, 86% of the area corresponds to natural areas in a good state of conservation, where forests are the dominant coverage6. In the northwest area, the region borders the Andean Cordillera and Orinoquía to the north. The political-administrative division includes the departments Amazonas, Caquetá, Guainía, Guaviare, Putumayo and Vaupés, and part of the departments Cauca, Meta, Nariño and Vichada. The human population is estimated at ~ 1.4 million, with a density of 2.5 inhab/km2. Internal conflict and poverty make this region one of the most important population dynamics in the country in terms of displacement36. The geographical location of the study area and the spatial pattern of the loss of forests that occurred between 2002 and 2016 are shown in Fig. 1.Figure 1Study area. Colombian Amazon and location of Amazonian tropical forests that were lost between 2002 and 2016. (Maps were generated using software ArcGis 10.7.1 https://www.esri.com).Full size imageLand cover maps and variables for change analysisThematic land cover maps used in this research were produced by the Colombian Amazon Land Cover Monitoring System (SIMCOBA) of the Amazon Institute for Scientific Research SINCHI (https://siatac.co/simcoba/). SIMCOBA has prepared land cover maps for the periods 2002, 2007, 2012, 2014, 2016 and 2018. Three of the land cover maps prepared were used in this study: 2002, 2016 and 2018 a scale of 1:100,00033. The maps were generated from the visual interpretation of a mosaic of Landsat 5 Thematic Mapper (TM) and Landsat 8 Operational Land Imager (OLI) images, using the PIAO technique (Photo Interprétation Assistée par Ordinateur). The classification categories of the land cover maps were based on the Corine land cover methodology adapted for Colombia37.The SIMCOBA system calculates the annual rates of Amazon forest loss (forest loss/ha/annual) by comparing the cover maps of the last two periods and subtracting from the previous map those forests that are no longer present in the most current map (Fig. 3). This process only considers the forests loss and the permanent forests. New forests due to natural regeneration or restoration are omitted in the calculations6.To facilitate the interpretation of changes and cover transitions, the classification categories of the maps were re-categorized into 7 types: “Amazon forests”, “floodplain forests”, “fragmented forests and secondary vegetation”, “grasslands and shrublands”, “water bodies and wetlands”, “pastures and crops” and “urban and artificialized cover”. The land cover maps were resampled at a resolution of 60 m × 60 m to facilitate the computational analysis of the explanatory model, the simulations of the scenarios, and to keep the detailed spatial resolution of the coverage and explanatory variables16.A geospatial database was created with a set of variables for the cover changes to create an explanatory model for each transition. Driving factors of change are grouped into the following variables: (1) accessibility, (2) climate, (3) landscape features, (4) production practices and environmental degradation, (5) landscape management, (6) socioeconomy, and (7) soil characteristics. We considered 41 explanatory variables (see supplementary information Table S1).Accessibility variables such as roads and navigable rivers were obtained from the geodatabase at a scale of 1:100,000 of the Agustín Codazzi Geographical Institute of Colombia (IGAC). Bioclimatic temperature data were obtained from Worldclim v1.438. Cover variables (e.g., patch sizes Amazon forests and distance to pastures and crops) were created using the software ArcGis (v.10.7.1)39 from the 2002 land cover map to understand which drivers were more influential in the dynamics of land-use changes since 2002 that resulted in the distribution of land cover in 2016.Degradation variables, such as advances of the agricultural frontier, were obtained from the Territorial Environmental Information System of the Colombian Amazon (SIAT-AC)40; livestock density data came from the Colombian Agricultural Institute (ICA); the fire density were processed from MODIS and VIIRS images (https://siatac.co/puntos-de-calor/); and the location of mining titles was obtained from the National Mining Agency.The information on the landscape features and socioeconomic variables was obtained from different sources: (1) the limit of the protected natural areas was provided by the National System of Protected Areas (SINAP)41, (2) the Amazon Forest Reserve areas (Second Law of 1959) were obtained from the Ministry of Environment and Sustainable Development (MADS), (3) the location of the indigenous reservations was provided by the Ministry of the Interior, and (4) the limits of the areas of Indigenous Reservations and the legal status of the Amazonian region were obtained from the SINCHI cartographic database40.Socioeconomic information was spatialized from data from the National Administrative Department of Statistics (DANE). Soil-type data were obtained from IGAG, and topographic and altitudinal variables were derived from a DEM at 100 m resolution from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER V003) sensor42. All explanatory variables were resampled at a resolution of 60 m.Patterns of land cover changes and transitionsThe transformation patterns of territory are mainly defined by human intentions and the activities that these groups plan to develop after making the land cover changes, as well as the dynamics of vegetation regeneration43. In this study, these changes in the study area were obtained and analyzed employing the Land Change Modeller (LCM) module of TerrSet34 and using the land cover maps for 2002 and 2016 as input information (Fig. 2).Figure 2(Source: Open Data—SINCHI Institute https://datos.siatac.co/pages/coberturas) (Maps were generated using software ArcGis 10.7.1 ).Land cover maps 2002, 2016 and 2018, produced by the Colombian Amazon Land Cover Monitoring System (SIMCOBA) of the Amazonian Research Institute SINCHIFull size imageTo represent dynamics and changes in the vegetation during the study period, a total of 14 transitions of greater importance in terms of area were considered (transitions with an area  More

Loreau, M. et al. Biodiversity and ecosystem functioning: current knowledge and future challenges. Science 294, 804–808 (2001).Article 
ADS 
CAS 

Google Scholar 
Pimm, S. L., Russell, G. J., Gittleman, J. L. & Brooks, T. M. The future of biodiversity. Science 269, 347–350 (1995).Article 
ADS 
CAS 

Google Scholar 
Newbold, T. et al. Global effects of land use on local terrestrial biodiversity. Nature 520, 45–50. https://doi.org/10.1038/nature14324 (2015).Article 
ADS 
CAS 

Google Scholar 
Cardoso, P. et al. Scientists’ warning to humanity on insect extinctions. Biol. Conserv. 242, 108426. https://doi.org/10.1016/j.biocon.2020.108426 (2020).Article 

Google Scholar 
Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE 12, e0185809. https://doi.org/10.1371/journal.pone.0185809 (2017).Article 
CAS 

Google Scholar 
Seibold, S. et al. Arthropod decline in grasslands and forests is associated with landscape-level drivers. Nature 574, 671–674. https://doi.org/10.1038/s41586-019-1684-3 (2019).Article 
ADS 
CAS 

Google Scholar 
Sánchez-Bayo, F. & Wyckhuys, K. A. G. Worldwide decline of the entomofauna: A review of its drivers. Biol. Cons. 232, 8–27. https://doi.org/10.1016/j.biocon.2019.01.020 (2019).Article 

Google Scholar 
Yang, L. H. & Gratton, C. Insects as drivers of ecosystem processes. Curr. Opin. Insect Sci. 2, 26–32. https://doi.org/10.1016/j.cois.2014.06.004 (2014).Article 

Google Scholar 
Bowler, D. E., Heldbjerg, H., Fox, A. D., de Jong, M. & Böhning-Gaese, K. Long-term declines of European insectivorous bird populations and potential causes. Conserv. Biol. 33, 1120–1130. https://doi.org/10.1111/cobi.13307 (2019).Article 

Google Scholar 
Biesmeijer, J. C. et al. Parallel declines in pollinators and insect-pollinated plants in Britain and the Netherlands. Science 313, 351–354. https://doi.org/10.1126/science.1127863 (2006).Article 
ADS 
CAS 

Google Scholar 
Tscharntke, T., Klein, A. M., Kruess, A., Steffan-Dewenter, I. & Thies, C. Landscape perspectives on agricultural intensification and biodiversity – ecosystem service management. Ecol. Lett. 8, 857–874. https://doi.org/10.1111/j.1461-0248.2005.00782.x (2005).Article 

Google Scholar 
Scherber, C. et al. Bottom-up effects of plant diversity on multitrophic interactions in a biodiversity experiment. Nature 468, 553–556. https://doi.org/10.1038/nature09492 (2010).Article 
ADS 
CAS 

Google Scholar 
Siemann, E., Tilman, D. & Haarstad, J. Insect species diversity, abundance and body size relationships. Nature 380, 704–706. https://doi.org/10.1038/380704a0 (1996).Article 
ADS 
CAS 

Google Scholar 
Borer, E. T., Seabloom, E. W. & Tilman, D. Plant diversity controls arthropod biomass and temporal stability. Ecol. Lett. 15, 1457–1464. https://doi.org/10.1111/ele.12006 (2012).Article 

Google Scholar 
Ebeling, A. et al. Plant diversity effects on arthropods and arthropod-dependent ecosystem functions in a biodiversity experiment. Basic Appl. Ecol. 26, 50–63. https://doi.org/10.1016/j.baae.2017.09.014 (2018).Article 

Google Scholar 
Ebeling, A. et al. Plant diversity induces shifts in the functional structure and diversity across trophic levels. Oikos 127, 208–219. https://doi.org/10.1111/oik.04210 (2018).Article 

Google Scholar 
Ebeling, A. et al. Plant diversity impacts decomposition and herbivory via changes in aboveground arthropods. PLoS ONE 9, e106529. https://doi.org/10.1371/journal.pone.0106529 (2014).Article 
ADS 
CAS 

Google Scholar 
Marquard, E. et al. Plant species richness and functional composition drive overyielding in a six-year grassland experiment. Ecology 90, 3290–3302 (2009).Article 

Google Scholar 
Tilman, D. et al. Diversity and productivity in a long-term grassland experiment. Science 294, 843–845. https://doi.org/10.1126/science.1060391 (2001).Article 
ADS 
CAS 

Google Scholar 
Simons, N. K. et al. Resource-mediated indirect effects of grassland management on arthropod diversity. PLoS ONE 9, e107033. https://doi.org/10.1371/journal.pone.0107033 (2014).Article 
ADS 
CAS 

Google Scholar 
Wardle, D. A., Nicholson, K. S., Bonner, K. I. & Yeates, G. W. Effects of agricultural intensification on soil-associated arthropod population dynamics, community structure, diversity and temporal variability over a seven-year period. Soil Biol. Biochem. 31, 1691–1706 (1999).Article 
CAS 

Google Scholar 
Luff, M. L. & Rushton, S. P. The ground beetle and spider fauna of managed and unimproved upland pasture. Agr. Ecosyst. Environ. 25, 195–206 (1989).Article 

Google Scholar 
Dennis, P., Young, M. R., Howard, C. L. & Gordon, I. J. The response of epigeal beetles (Col, Carabidae, Staphylinidae) to varied grazing regimes on upland Nardus stricta grasslands. J. Appl. Ecol. 34, 433–443 (1997).Article 

Google Scholar 
Murdoch, W. W., Evans, F. C. & Peterson, C. H. Diversity and pattern in plants and insects. Ecology 53, 819–829 (1972).Article 

Google Scholar 
Siemann, E., Tilman, D., Haarstad, J. & Ritchie, M. Experimental tests of the dependence of arthropod diversity on plant diversity. Am. Nat. 152, 738–750 (1998).Article 
CAS 

Google Scholar 
Joern, A. & Laws, A. N. Ecological mechanisms underlying arthropod species diversity in grasslands. Annu. Rev. Entomol. 58, 19–36. https://doi.org/10.1146/annurev-ento-120811-153540 (2013).Article 
CAS 

Google Scholar 
Hunter, M. D. & Price, P. W. Playing chutes and ladders: Heterogeneity and relative roles of bottom-up and top-down forces in natural communities. Ecology 73, 724–732 (1992).Article 

Google Scholar 
Knops, J. M. H. et al. Effects of plant species richness on invasion dynamics, disease outbreaks, insect abundances and diversity. Ecol. Lett. 2, 286–293 (1999).Article 
CAS 

Google Scholar 
Thiele, H. U. Carabid beetles in their environment. A study on habitat selection by adaptions in physiology and behaviour. (Springer- Verlag, 1977).Harvey, J. A., van der Putten, W. H., Turin, H., Wagenaar, R. & Bezemer, T. M. Effects of changes in plant species richness and community traits on carabid assemblages and feeding guilds. Agr. Ecosyst. Environ. 127, 100–106 (2008).Article 

Google Scholar 
Luff, M. L. Use of Carabids as environmental indicators in grasslands and cereals. Ann. Zool. Fenn. 33, 185–195 (1996).
Google Scholar 
Kotze, D. J. et al. Forty years of carabid beetle research in Europe—from taxonomy, biology, ecology and population studies to bioindication, habitat assessment and conservation. ZooKeys https://doi.org/10.3897/zookeys.100.1523 (2011).Article 

Google Scholar 
Barnes, A. D. et al. Biodiversity enhances the multitrophic control of arthropod herbivory. Sci. Adv. 6, eabb6603. https://doi.org/10.1126/sciadv.abb6603 (2020).Article 
ADS 

Google Scholar 
Bianchi, F. J. J. A., Booij, C. J. H. & Tscharntke, T. Sustainable pest regulation in agricultural landscapes: A review on landscape composition, biodiversity and natural pest control. Proc. R. Soc. B: Biol. Sci. 273, 1715–1727. https://doi.org/10.1098/rspb.2006.3530 (2006).Article 
CAS 

Google Scholar 
Lövei, G. L. & Magura, T. Ground beetle (Coleoptera: Carabidae) diversity is higher in narrow hedges composed of a native compared to non-native trees in a Danish agricultural landscape. Insect Conserv. Divers. 10, 141–150. https://doi.org/10.1111/icad.12210 (2017).Article 

Google Scholar 
Loreau, M. Consumers as maximizers of matter and energy flow in ecosystems. Am. Nat. 145, 22–42. https://doi.org/10.1086/285726 (1995).Article 

Google Scholar 
Mielke, L. et al. Nematode grazing increases the allocation of plant-derived carbon to soil bacteria and saprophytic fungi, and activates bacterial species of the rhizosphere. Pedobiologia 90, 150787. https://doi.org/10.1016/j.pedobi.2021.150787 (2022).Article 

Google Scholar 
Holland, J. M. & Luff, M. L. The effects of agricultural practices on Carabidae in temperate agroecosystems. Integr. Pest Manag. Rev. 5, 109–129. https://doi.org/10.1023/A:1009619309424 (2000).Article 

Google Scholar 
Roscher, C. et al. The role of biodiversity for element cycling and trophic interactions: An experimental approach in a grassland community. Basic Appl. Ecol. 5, 107–121 (2004).Article 

Google Scholar 
Weisser, W. W. et al. Biodiversity effects on ecosystem functioning in a 15-year grassland experiment: Patterns, mechanisms, and open questions. Basic Appl. Ecol. https://doi.org/10.1016/j.baae.2017.06.002 (2017).Article 

Google Scholar 
Freude, H., Harde, K. W. & Lohse, G. A. Die Käfer Mitteleuropas Bd.1–11. (Goecke & Evers, 1965–83).Koch, K. Die Käfer Mitteleuropas. Ökologie Bd.1–6. (Goecke & Evers, 1989–95).R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2021).Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
Schmid, B., Baruffol, M., Wang, Z. & Niklaus, P. A. A guide to analyzing biodiversity experiments. J. Plant Ecol. 10, 91–110. https://doi.org/10.1093/jpe/rtw107 (2017).Article 

Google Scholar 
Zuur, A., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. Mixed effects models and extensions in ecology with R. (Springer, 2009).Oksanen, J. et al. vegan: Community Ecology Package v. 2.6–2 (2022).Lenth, R. et al., emmeans: Estimated Marginal Means, aka Least-Squares Means v. 1.8.1-1 (2022).Lovei, G. L. & Sunderland, K. D. Ecology and behavior of ground beetles (Coleoptera: Carabidae). Annu. Rev. Entomol. 41, 231–256 (1996).Article 
CAS 

Google Scholar 
Ravenek, J. M. et al. Long-term study of root biomass in a biodiversity experiment reveals shifts in diversity effects over time. Oikos 123, 1528–1536. https://doi.org/10.1111/oik.01502 (2014).Article 

Google Scholar 
Root, R. Organization of a plant -arthropod association in simple and diverse habitats: the fauna of collards (Brassica oleracea). Ecol. Monogr. 43, 95–124 (1973).Article 

Google Scholar 
Duelli, P. & Obrist, M. K. Regional biodiversity in an agricultural landscape: The contribution of seminatural habitat islands. Basic Appl. Ecol. 4, 129–138 (2003).Article 

Google Scholar 
Perner, J. & Malt, S. Assessment of changing agricultural land use: Response of vegetation, ground-dwelling spiders and beetles to the conversion of arable land into grassland. Agr. Ecosyst. Environ. 98, 169–181 (2003).Article 

Google Scholar 
Purtauf, T., Dauber, J. & Wolters, V. Carabid communities in the spatio-temporal mosaic of a rural landscape. Landsc. Urban Plan. 67, 185–193 (2004).Article 

Google Scholar 
Eisenhauer, N. et al. Biotic interactions, community assembly, and eco-evolutionary dynamics as drivers of long-term biodiversity–ecosystem functioning relationships. Res. Ideas Outcomes https://doi.org/10.3897/rio.5.e47042 (2019).Article 

Google Scholar 
Guerrero-Ramirez, N. R. et al. Diversity-dependent temporal divergence of ecosystem functioning in experimental ecosystems. Nat. Ecol. Evol. 1, 1639–1642. https://doi.org/10.1038/s41559-017-0325-1 (2017).Article 

Google Scholar 
Reich, P. B. et al. Impacts of biodiversity loss escalate through time as redundancy fades. Science 336, 589–592. https://doi.org/10.1126/science.1217909 (2012).Article 
ADS 
CAS 

Google Scholar 
Isbell, F. I., Polley, H. W. & Wilsey, B. J. Biodiversity, productivity and the temporal stability of productivity: Patterns and processes. Ecol. Lett. 12, 443–451. https://doi.org/10.1111/j.1461-0248.2009.01299.x (2009).Article 

Google Scholar 
Blake, S., Foster, G. N., Fisher, G. E. J. & Ligertwood, G. L. Effects of management practices on the carabid faunas of newly established wildflower meadows in southern Scotland. Ann. Zool. Fenn. 33, 139–147 (1996).
Google Scholar 
Boetzl, F. A., Krimmer, E., Krauss, J. & Steffan-Dewenter, I. Agri-environmental schemes promote ground-dwelling predators in adjacent oilseed rape fields: Diversity, species traits and distance-decay functions. J. Appl. Ecol. 56, 10–20. https://doi.org/10.1111/1365-2664.13162 (2019).Article 

Google Scholar 
Knapp, M., Seidl, M., Knappová, J., Macek, M. & Saska, P. Temporal changes in the spatial distribution of carabid beetles around arable field-woodlot boundaries. Sci. Rep. 9, 8967. https://doi.org/10.1038/s41598-019-45378-7 (2019).Article 
ADS 
CAS 

Google Scholar  More

Our objective in developing the Eco-ISEA3H database37 was to compile a coordinated, global set of tabular data, characterizing environmental conditions and the geographic distributions of large mammalian species. The database was built on the ISEA3H DGGS, a multi-resolution system of global grids, each grid dividing the Earth’s surface into discrete, equal-area hexagonal cells. These cells constitute areal units of observation, uniformly resampling data provided in different coordinate reference systems, spatial resolutions, geographic data models, and file formats. We included data at six consecutive ISEA3H resolutions, in which cell centroid spacing ranges from 29 kilometers to approximately 450 kilometers.Eco-ISEA3H themes and variables were derived from 17 geospatial data sources, and represent 3,033 features to be used for ML-based predictive modeling. Source datasets were published in raster or vector format, data models built on fundamentally different representations of spatial phenomena. Raster datasets comprise regular arrays of pixels, each pixel holding a value, while vector datasets comprise point, line, and polygon features, each feature defined by one or more (x, y) coordinate pairs and attributed with one or more values. Our task was to integrate these disparate source datasets, resampling and summarizing the values of raster pixels and vector features via the discrete, equal-area cells of the ISEA3H global grid system. The hexagonal cells on which the Eco-ISEA3H database37 is built thus serve as unifying observational units for SDM and ecometric analysis and modeling.From the statistical and ML perspective, each areal observational unit is characterized by (1) a set of environmental variables, representing climatic conditions, soil and near-surface lithology, land cover, and physical geography; and (2) a set of occurrence variables, representing the present and estimated natural distributions of large mammalian species. Predictive modeling tasks for statistical and ML modeling can be formulated in two directions: predicting species’ occurrences as a function of climatic and other environmental conditions (as in SDM studies), or predicting climatic and other environmental conditions as a function of species’ occurrences and functional traits (as in ecometric studies).Spatial units of observationTo study continuous spatial phenomena over a region of interest, it is often necessary to divide the region into a number of discrete, areal observational units, which may be used in statistical summaries and/or modeling. Machine learning methods for ecometric and species distribution modeling require discrete observational units, each characterized by two sets of variables, one describing environmental conditions, the other species’ geographic distributions. A major question in data representation concerns the form of these units; defining discrete spatial units of observation constitutes a well-known problem in geography, termed the modifiable areal unit problem (MAUP)38. As we change the size of proposed observational units, or change the boundaries between units while holding unit areas constant, measures of interest within these units – and derived summary statistics and model parameters – may differ; these are termed the “scale” and “zone” effects, respectively38.Our objective in utilizing the ISEA3H DGGS34 was to implement a robust spatial division of the Earth’s surface. The grid cells of the DGGS discretize the Earth’s sphere, forming, at each DGGS resolution, a global set of areal observational units with which to sample and summarize source datasets. To be optimally effective in the observation, simulation, and visualization of spatial phenomena, such a grid must meet certain structural criteria. We propose, modifying the Goodchild Criteria39, the DGGS grid must contain (1) contiguous, (2) equivalent observational units, (3) minimizing intra-unit variability, (4) having uniform topology with neighboring units, and (5) being visually effective, facilitating interpretation and communication. Each criterion will be discussed in detail; further, we will argue the ISEA3H DGGS selected for this study satisfies these criteria.Contiguity & congruencyWe suggest that a regular tiling maximally satisfies the criteria of (1) contiguity and (2) equivalence. A tiling is simply a set of shapes which cover a plane without gaps or overlaps40. A regular tiling is one of a class of tilings in which the tiles – our observational units – are highly equal; such tilings are monohedral, and composed of congruent, regular (equiangular and equilateral) polygons. Thus, regular tilings are also highly symmetrical, being vertex-, edge-, tile-, and flag-transitive. Three regular polygons may be used to create a regular tiling: the equilateral triangle, the square, and the regular hexagon40.With this suggestion, we follow common convention; in ecology, grids of square (or rectangular) cells are most often utilized, motivated in part by the use of raster datasets41, made of rectilinear rows and columns of pixels. However, it should be noted that while the square cells of these grids are equal in the coordinate reference system in which they are defined, such cells are rarely congruent, or indeed even square, on the Earth’s surface. The properties of the ISEA projection selected for this DGGS – area preservation, and relatively low angular distortion – serve to retain considerable congruency when inversely projecting grid cells to the spherical surface of the Earth.CompactnessTo accurately represent the spatially continuous phenomena of the Earth system, the grid cells of a DGGS – the areal observational units used in summarizing, modeling, and visualizing – must effectively discretize these phenomena. Thus, the DGGS must be structured such that (3) intra-unit variability is minimized, and inter-unit variability is maximized. In this way, patterns of variation among units more accurately represent patterns of variation inherent in the phenomena.Intra-unit variability may be minimized, in expectation, by compact observational units. Tobler’s oft-cited first law of geography serves as explanation: “everything is related to everything else, but near things are more related than distant things”42. Thus, compact units, in which all portions of the interior are nearer each other, are expected to contain less interior variability than elongated units, in which portions of the interior may be more distant. Given these properties, compact units are optimal in the context of DGGS development, elongated units in the context of efficient ecological sampling.Regular hexagons are the most compact of the three polygons – the equilateral triangle, square, and regular hexagon – admitting regular tilings. This compactness may be expressed in several related and complementary ways. First, of any equal-area tiling, regular hexagons have the minimum possible ratio of perimeter to area43. In minimizing perimeter length per unit area, regular hexagons are thus the most circle-like of the polygons admitting equal-area tilings. Relatedly, regular hexagonal packing is the highest-density arrangement of equal-area circles on a plane44.Finally, a regular hexagonal lattice optimally quantizes a plane; of the polygons admitting regular tilings, regular hexagons minimize the mean squared distance of any point to the nearest polygon centroid45. This distance, or “dimensionless second moment,” quantifies the more qualitative notion of interior nearness discussed in relation to Tobler’s Law.TopologyIn addition to maximally satisfying the (3) compactness criterion, regular hexagons have a topological advantage over equilateral triangles and squares. Of these three regular polygons, hexagons have the simplest relationship with neighbors in a tiling or grid, each (4) uniformly sharing an edge with the six adjacent hexagons forming its first-order neighborhood. Triangles and squares, in contrast, share only a single vertex with three or four neighbors, respectively, and an edge with three or four neighbors, complicating the definition of neighborhood in these grids.It follows that hexagonal topology has greater angular resolution than edge-based triangular or square topologies; movement may be simulated between cells in six directions, rather than in three or four, respectively. These properties – neighborhood simplicity and angular resolution – were confirmed by Golay46, in the context of pattern transformation operations on two-dimensional arrays. Further, these properties likely account for the widespread use of hexagonal grids in strategy board games, since these grids were introduced in the early 1960s47.Differing grid topologies affect the results of ecological models simulating dispersal. White and Kiester48, for example, found the topology of the network of communities in a neutral community ecology model – in which simulated communities had hexagonal neighborhoods, or von Neumann, Moore, or Margolus neighborhoods – affected modeled species abundances and diversities, but in complex ways, which differed given different model parameter values. (Note that the four neighbors with which a square cell shares an edge are termed its rook, or von Neumann neighborhood, and these plus the four neighbors with which it shares a single vertex its queen, or Moore neighborhood.)VisualizationFinally, in addition to these gains in representational accuracy, (5) hexagonal tilings are more visually effective than square tilings. Whether used in cartography or other two-dimensional data visualization, tilings inevitably create visual lines, artifacts of the lattice of shared edges between tiles49. Given our “sense of gravitational balance,” Carr et al.49 argue the horizontal and vertical lines of square tilings strongly distract the human eye, obscuring data-driven patterns in a dataset so visualized. The non-orthogonal lines of hexagonal tilings, however, feature less prominently, and thus distract less from patterns of interest49.Note that this is not an issue of aesthetics only: maps are often essential tools in scientific reasoning and communication, and effective visualization is important. Indeed, Carr et al.49 suggest this visual advantage makes a stronger case for hexagonal tilings than the representational advantages discussed previously.DGGS sampling workflowsThe set of scripted workflows developed to incorporate spatial datasets into the Eco-ISEA3H database37 utilize published spatial libraries and packages for Python and R, and include several validation steps, intended to verify the integrity of source datasets and the fidelity of the transfer to the DGGS. Workflows developed for raster datasets are presented in Fig. 1, and workflows for vector datasets in Fig. 2.Fig. 1Workflow developed to incorporate raster datasets into the ISEA3H DGGS.Full size imageFig. 2Workflow developed to incorporate vector datasets into the ISEA3H DGGS.Full size imageTo begin, one general principle guides each workflow: each source dataset is processed in its native coordinate reference system. In all cases, a representation of the DGGS is developed in the coordinate reference system of the source dataset, and used in summarizing that dataset. The guiding premise here is that the spatial dataset is as the authors intended it in the coordinate reference system in which it is published and distributed.This is especially relevant for vector polygon datasets. Consider, for example, certain species’ range polygons published by the IUCN Red List50; these polygons are defined only roughly, having relatively few, widely spaced vertices, connected by arcs many hundreds of kilometers in length. These arcs are “straight” in the plate carrée projection with which the dataset’s WGS84 latitude/longitude coordinates are visualized by default. If vertex coordinates were projected into another coordinate reference system, the arcs would be similarly “straight” in this new system, and thus potentially trace very different paths across the Earth’s surface. Absent information to the contrary, we assume the arcs are as intended in the reference system in which the data are distributed.The spatial structure of raster datasets depends similarly on each dataset’s coordinate reference system; rasters are made of rows and columns of pixels, rectilinear and orthogonal only in the raster’s native coordinate reference system. We assume raster pixels are “atomic” units, each indivisible and representative of the area it natively covers. Thus, we query the DGGS at each pixel’s centroid, and assign the pixel wholly to the coincident DGGS cell.Raster dataset processingIf necessary, source raster datasets were first converted to the GeoTIFF file format, so that the files were readable in the open-source GIS software used later in the processing workflow. GeoTIFF files are simply Tag Image File Format (TIFF) image files with embedded georeferencing information, describing the dataset’s spatial extent and coordinate reference system. Hierarchical Data Format Release 4 (HDF4) files were converted to GeoTIFF format using the Geospatial Data Abstraction Library (GDAL) translate utility51.Next, raster tiles containing ISEA3H hexagon identification (HID) indexing numbers were generated; these integer HIDs uniquely identify each cell at a given ISEA3H resolution. A set of HID raster tiles was required for each source raster dataset, for each ISEA3H resolution, because (1) GeoTIFF rasters are able to hold only a single value at each pixel; and (2) HIDs sequentially number cells at a given ISEA3H resolution, from 1 to the number of cells present at that resolution. Thus, HIDs are not unique between resolutions; HID 84, for example, identifies a cell at each ISEA3H resolution 2 and higher.The HID raster tiles generated for a source raster dataset matched that dataset’s grid resolution, extent, and coordinate reference system precisely; thus, there was a one-to-one correlation between the pixels of the HID raster tiles and the source raster dataset tiles. For each tile, pixel centroid coordinates were passed to the dggridR package52 for R, which returned the ISEA3H cell identification number for that location. In this way, the pixels of the source raster were treated as indivisible units, assigned wholly to a particular HID on the basis of each pixel’s centroid. HID rasters were written in GeoTIFF format using the raster package53 for R.In equal-area projected coordinate reference systems, simple counts of the number of raster pixels assigned to each HID were sufficient to determine each ISEA3H cell’s total area. In all other cases – for example, for raster datasets using the World Geodetic System 1984 (WGS84) coordinate reference system – raster tiles containing pixel areas were generated. These areas were calculated by passing each pixel’s corner coordinates to the GeographicLib library54 for Python.Finally, source raster dataset tiles, HID raster tiles, and area raster tiles (for source rasters using non-authalic coordinate reference systems) were superimposed to generate summary tabular files, describing the features of the source raster dataset by ISEA3H cell. The specifics of this process, which utilized functions of the raster package53 for R, depended on whether the source raster contained discrete, categorical values, or continuous, real-numbered values.Discrete themesFor each source raster dataset containing discrete pixel values, one or more of the following summary statistics were calculated. While the centroid attribute requires a simple point sample, the fraction and mode attributes are area-integrated, and involve a multiple-step sampling process. For rasters using an authalic coordinate reference system, the raster package’s crosstab function53 was used to generate a contingency table for each tile; applied to source raster and HID raster tiles, the function tallied the number of pixels of each class coincident with each HID, for each tile. These tile-specific tables were then summed, to obtain total counts of pixels of each class within each HID.For rasters using a non-authalic coordinate reference system, area raster tiles were required as well. For each tile, a vector of classes present in the source raster was assembled. For each of these classes in turn, a mask raster tile was generated, retaining pixels belonging to the class, and screening pixels belonging to all other classes. This mask was applied to the area raster tile, and retained pixels were summed within each HID using the raster package’s zonal function53. Thus, a contingency table was compiled for each raster tile, containing the area of each class within each HID. Finally, these tile-specific tables were summed, to obtain the total area of each class within each HID.

Centroid. The centroid attribute records the categorical value occurring at each ISEA3H cell’s centroid. Where the source raster dataset contains a null value at a centroid, the cell is assigned a flag signifying no value is available.

Fraction. The fraction attributes record the proportion of each ISEA3H cell’s area covered by each categorical value. For example, the Köppen-Geiger climate classification system, as implemented by Beck et al.55, includes 30 classes, listed in Table 4. Thus, each ISEA3H cell has an associated set of 30 fraction attributes for this dataset, recording the proportions of the cell’s area covered by the 30 categorical values, from tropical rainforest (Af) to polar tundra (ET).

Mode. The mode attribute records the categorical value covering the greatest proportion of each ISEA3H cell’s area. For example, if an ISEA3H cell had a fraction value of 0.4 for some hypothetical categorical value A, 0.3 for B, and 0.3 for C, it would be assigned a mode value of A. A mode attribute is specified for cells in which the sum of the fraction attributes is greater than or equal to 0.2; where fraction attributes total less than 0.2, a flag signifying no value is assigned.

Continuous variablesFor each source raster dataset containing continuous pixel values, one or more of the following summary statistics were calculated. Again, the centroid attribute requires only a simple point sample, while the mean attribute is area-integrated, requiring area raster tiles for source rasters using a non-authalic coordinate reference system.

Centroid. The centroid attribute records the continuous value occurring at each ISEA3H cell’s centroid. Where the source raster dataset contains a null value at a centroid, the cell is assigned a flag signifying no value is available.

Mean. The mean attribute records the area-weighted arithmetic mean of the continuous values of raster pixels within each ISEA3H cell. For raster datasets in authalic coordinate reference systems, the area-weighted mean is equivalent to the simple mean of the values of raster pixels within each cell; however, in all other cases, pixel values are weighted by pixel areas per the equation below, in which wi and xi indicate the area and value, respectively, of each pixel i within an ISEA3H cell containing n pixels.

$$overline{x}=frac{{sum }_{i=1}^{n}{w}_{i}{x}_{i}}{{sum }_{i=1}^{n}{w}_{i}}$$For each tile, source raster values and area values were multiplied, pixel by pixel, using the raster package’s * arithmetic operator53. The resulting product raster tile, as well as the area raster tile, were then summed within each HID using the raster package’s zonal function53. Finally, these tile-specific tables were summed, to obtain both the numerator (summed product values) and denominator (summed area values) for the above equation, for each HID.Vector dataset processingSource vector datasets incorporated into the Eco-ISEA3H database37 contain polygon features, discrete areas assigned a categorical value. A dataset may (1) contain polygons of several different classes; for example, the vector shapefile published by Olson et al.56 contains ecoregion polygons, each assigned to one of several biogeographic realms. Alternatively, a dataset may (2) represent a single class, with polygons indicating class presence; for example, the shapefiles published by the IUCN Red List50 each represent a species’ geographic range, with polygons indicating regions the species is present. In both cases, the summary statistics discussed in reference to raster datasets containing discrete values may be calculated.Prior to inclusion in the Eco-ISEA3H database37, source vector datasets were preprocessed. To simplify the geographic representation of the class(es) of interest – that is, to remove unnecessary polygon boundaries – dataset polygons were dissolved, either on the class attribute in case (1), or globally in case (2), using the QGIS open-source desktop GIS application. The geodesic areas of dissolved polygons were then calculated using the GeographicLib library54. Finally, the geometries of dissolved polygons were checked for conformance with the OGC Simple Feature Access standard57 using the Shapely library58 for Python, ensuring these features served as valid input in the processing workflow to follow.The intersection of source dataset polygons and ISEA3H cell polygons is central to the vector processing workflow. Source polygons result from the preliminary simplification and verification steps just discussed; cell polygons result from polygonizing a set of HID raster tiles for the ISEA3H resolution of interest. The polygonizing procedure utilized the open-source GDAL command-line tools polygonize and ogrmerge51, as well as the GeographicLib54 and Shapely58 libraries. Polygonizing HID raster tiles of the appropriate coordinate reference system (specifically, the system matching that of the source polygon dataset) ensured HID polygon boundaries displayed both proper geodesic curvature and the shape distortion induced by the ISEA map projection.Intersection is a set-theoretic operation, returning polygons representing each coincident class/HID combination. The operation was implemented via the Shapely library58, and the geodesic areas of intersected polygons were calculated via the GeographicLib library54. Note that the scripted intersection tools developed for the Eco-ISEA3H database37 allow limiting the ISEA3H cells included in a single tool run, to break the processing of large datasets into manageable pieces. Runs may be limited to a user-specified range of HIDs. Additionally, if cells at the next coarser or finer ISEA3H resolution have been intersected with the source dataset, cells retained by the operation may be used as a spatial index; a list of coincident HIDs at the ISEA3H resolution of interest may be generated, and used to limit tool runs.An output shapefile is written, containing intersected polygons attributed with the geodesic area, the HID, and in case (1), the source class. Next, an additional verification of the geometries of these intersected polygons is performed. Each intersected polygon is superimposed over the original ISEA3H cell polygon having the same HID. If intersected polygons have too few vertices to be valid, or are not contained by the original cell polygon from which each was derived, these polygons are flagged for review and revision. This step was implemented to catch geometry errors observed early in the development of the Eco-ISEA3H intersection tools.Finally, the geodesic areas of intersected polygons are totaled, and the total area of each class within each HID is calculated. Dividing by the geodesic areas of the original ISEA3H cell polygons, these class totals are expressed as fractions of each cell’s total area. In two final verification steps, (1) the total intersected area of each class, across all HIDs, is compared to the area of the same class in the source dataset; and (2) class fraction values are confirmed to be less than or equal to unity within each HID. Deviations are flagged for review and revision.Data sources & themesThe Eco-ISEA3H database37 incorporates 17 source datasets, characterizing the Earth’s climate, geology, land cover, and physical geography, as well as human population density and the geographic ranges of nearly 900 large mammalian species. Data sources are listed in Table 1. We first present a brief overview of these sources, and describe sources and themes in greater detail in the following sections.Table 1 Source datasets and themes included in the Eco-ISEA3H database37. Each dataset is described by full and abbreviated name, source, spatial resolution (for datasets published/distributed at more than one resolution), version, and scenario. Each theme is described by full and abbreviated name and type (whether it contains discrete, categorical values or continuous, real-valued variables).Full size tableClimate is characterized primarily by temperature- and precipitation-based averages and extremes, summarized over the past 50 to 70 years, and forecasted for 40 to 60 years in the future under the RCP 8.5 climate change scenario59; data sources include WorldClim30,31, ENVIREM60, and the ETCCDI extremes indices derived by Sillmann et al.61,62 from ERA-4063 and CCSM464. Additionally, present climate is classified via the Köppen-Geiger climate classification system, from GLOH2O55. Geological data include soil types, from the Digital Soil Map of the World (DSMW)65; near-surface rock types, from the Global Lithological Map (GLiM)66; and sedimentary basin types67. Human geography is quantified by human population density, from the Gridded Population of the World (GPW)68. Land cover is described by the International Geosphere-Biosphere Programme (IGBP) cover classification scheme, from MCD12Q169; and by percent tree, non-tree, and non-vegetated cover, from MOD44B70. The Earth’s physical geography is characterized by continental and island landmasses, from Natural Earth; lakes and wetlands, from the Global Lakes and Wetlands Database (GLWD)71; biogeographic realms56; and terrestrial topography and ocean bathymetry, from ENVIREM60 and SRTM30_PLUS72. Finally, distributional data include the present and estimated natural ranges of large mammalian species, from the IUCN Red List50 and the Phylogenetic Atlas of Mammal Macroecology (PHYLACINE)73,74.Climate

ENVIREM. The ENVIREM (ENVIronmental Rasters for Ecological Modeling) dataset60 contains 16 climatic variables derived from WorldClim v1.4 monthly temperature and precipitation30, and extraterrestrial radiation. These are intended to compliment the WorldClim v1.4 bioclimatic variables30, capturing additional environmental features directly relevant to floral and faunal physiology and ecology60. Source rasters at 30 arc-second resolution were summarized by area-weighted mean at ISEA3H resolutions 8 and 9. Variable codes, descriptions, and units are listed in Table 2. Title and Bemmels60, and references therein, provide full definitions and calculation methods for these variables.Table 2 Codes, descriptions, and units for the 16 ENVIREM climatic variables, from Title and Bemmels60.Full size table

ETCCDI. A comprehensive set of 27 climate extremes indices was defined by the Expert Team on Climate Change Detection and Indices (ETCCDI); these generally capture “moderate” extremes, having recurrence intervals of a year or shorter, and are based on observed/simulated daily temperature and precipitation61,62. Sillmann et al.61,62 derive these indices from results of a number of global climate models and atmospheric reanalyses, several of which were incorporated in the Eco-ISEA3H database37. Given the relatively low-resolution grids used in modeling and reanalysis, these source rasters were interpolated to ISEA3H cell centroids by inverse (geodesic) distance weighting (IDW). Variable codes, descriptions, and units are listed in Table 3. Sillmann et al.61 provide full definitions and calculation methods for these indices.Table 3 Codes, descriptions, and units for the 27 ETCCDI climate extremes indices, from Sillmann et al.61,62.Full size table

The Eco-ISEA3H database37 includes ETCCDI variables based on results of the ERA-40 reanalysis63, produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). The reanalysis combines past meteorological observations with a weather forecasting model, producing a global representation of the state of the atmosphere for each reanalysis time step, usually a six-hour interval63. These were averaged for the period 1958 to 2001, the 44 full years for which the ERA-40 reanalysis was conducted, and were interpolated to ISEA3H resolutions 5 to 9.Additionally, the database includes ETCCDI variables based on results of the Community Climate System Model v4 (CCSM4), a global climate model developed for CMIP564. These were averaged for the period 1950 to 2000, to match the approximate period covered by WorldClim v1.4, and for the period 2061 to 2080, to match the final interval for which CCSM4 model results were downscaled/debiased using WorldClim v1.430. Variables were interpolated to ISEA3H resolution 9.ETCCDI variables for this latter period represent conditions under Representative Concentration Pathway (RCP) 8.5, the RCP resulting in the highest radiative forcing (8.5 W/m2) by 210059. This scenario was selected such that future conditions maximally different from the present might be considered; in RCP 8.5, rapid population growth, and relatively slow growth in per capita income and technological development, lead to high energy demand without associated climate mitigation policies, resulting in greenhouse gas emissions and atmospheric concentrations increasing significantly in the coming decades59.

Köppen-Geiger Climate Classification. As implemented by Beck et al.55, the Köppen-Geiger system classifies the Earth’s terrestrial climates into five primary classes, and further into 30 subclasses, based on a set of threshold criteria referencing monthly mean temperature and precipitation. These climate classes are ecologically significant, as regions within each class support floral communities sharing common characteristics. Beck et al.55 utilize four climatic datasets, including WorldClim v1.x and v2.x, adjusted to the period 1980 to 2016, to define the present-day classes incorporated in the Eco-ISEA3H database37. The source raster at 30 arc-second resolution was summarized by fraction and mode at ISEA3H resolution 9. Variable codes and descriptions are listed in Table 4.Table 4 Codes and descriptions for the 30 Köppen-Geiger climate classes, from Beck et al.55.Full size table

WorldClim v1.4. The first-generation WorldClim dataset30 contains four monthly themes, each with 12 variables, characterizing monthly temperature and precipitation; additionally, it contains 19 bioclimatic variables, derived from the monthly variables, capturing biologically relevant seasonal and annual features of the climate system. These bioclimatic variables, first developed for the BIOCLIM species distribution modeling (SDM) package75, are used extensively in SDM studies; a recent synthesis found most were included in more than 1,000 published MaxEnt SDMs (of 2,040 reviewed)76.

WorldClim monthly temperature and precipitation rasters are interpolated from weather station observations averaged for the approximate period 1950 to 2000. The interpolation was done using thin plate smoothing splines, with latitude, longitude, and elevation as predictor variables30. These rasters characterize present-day climate, and further served as an observational baseline with which the predictions of CMIP5 global climate models were downscaled and bias-corrected.The 19 bioclimatic variables, for both present-day and future conditions (the latter averaged for the period 2061 to 2080, from the CCSM4 RCP 8.5 simulation), were incorporated into the Eco-ISEA3H database37; source rasters at 30 arc-second resolution were summarized by area-weighted mean at ISEA3H resolution 9. Variable codes, descriptions, and units are listed in Table 5. O’Donnell and Ignizio77 provide full definitions and calculation methods for these variables.Table 5 Codes, descriptions, and units for the 19 WorldClim bioclimatic variables, from v1.430 and v2.031.Full size table

WorldClim v2.0. The second-generation WorldClim dataset31 contains seven monthly themes, each with 12 variables, characterizing monthly temperature, precipitation, solar radiation, wind speed, and vapor pressure; additionally, it contains the standard set of 19 bioclimatic variables, derived from monthly temperature and precipitation.

As in the first-generation dataset, monthly rasters were interpolated from weather station observations, averaged here for the approximate period 1970 to 200031. Again, thin plate smoothing splines were used in the interpolation, but with additional covariates included for one or more interpolated features: distance to coast, computed extraterrestrial radiation, and three satellite-derived observations – cloud cover, and maximum and minimum land surface temperature, from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument.The 12 source rasters for each of the seven monthly themes, at 30 arc-second resolution, were summarized by centroid at ISEA3H resolutions 5 to 10. Additionally, the 19 source bioclimatic rasters, at 30 arc-second resolution, were summarized by centroid at ISEA3H resolutions 5 to 10, and by area-weighted mean at ISEA3H resolutions 6 to 9. Codes, descriptions, and units for the bioclimatic variables are listed in Table 5.Geol10ogy

DSMW. The Digital Soil Map of the World (DSMW)65 describes the geographic distribution and physical and chemical properties of the world’s soils. The DSMW was digitized from the FAO-UNESCO Soil Map of the World, printed at 1:5,000,000 scale. Each digitized mapping unit is assigned a number of soil attributes; here we classify units via the DOMSOI attribute, the dominant soil or land unit code. The DSMW includes 117 soils in 26 major soil groupings, as well as six other land units, for a total of 123 DOMSOI classes. The source vector dataset was summarized by fraction and mode at ISEA3H resolutions 5, 6, and 9. Variable codes and descriptions are listed in Table 6.Table 6 Codes and descriptions for the 123 DSMW soil and land units, from the FAO65.Full size table

GLiM. The Global Lithological Map (GLiM)66 represents the rock and unconsolidated sediments at or near the Earth’s terrestrial surface; this geological material is a source of geochemical flux to the Earth’s soils, biosphere, and hydrosphere. Hartmann and Moosdorf66 compiled the map and accompanying database from 92 regional geological maps and 318 literature sources. Rock was classified into 16 first-level lithological classes; 12 second-level and 14 third-level subclasses further describe specific mineralogical and physical properties.

The source vector dataset was summarized by centroid at ISEA3H resolution 9. Variable codes and descriptions are listed in Table 7. The attribute assigned each ISEA3H cell takes the form xxyyzz; underscore characters (_) in the yy and/or zz slots indicate the second- and/or third-level subclasses were undefined.Table 7 Codes and descriptions for the 16 first-level, 12 second-level, and 14 third-level GLiM lithological classes, from Hartmann and Moosdorf66.Full size table

Sedimentary Basins. Sedimentary basins are areas of subsidence in the Earth’s crust, in which sediments eroded from uplands are deposited and potentially preserved for a million or more years67, thus entering the planet’s long-term geological record. Nyberg and Howell67 delineate active sedimentary basins, covering both the Earth’s terrestrial surface and marine areas over continental crust. The authors operationally defined basins as low-relief areas containing Quaternary Period sediments, and further classified the basins by tectonic setting, identifying backarc, forearc, foreland, extensional, intracratonic, passive margin, and strike-slip basins on the basis of published literature and geological maps67.

Terrestrial basins were incorporated in the Eco-ISEA3H database37. Note that no terrestrial backarc basins were delineated. The source vector dataset was summarized by fraction and mode at ISEA3H resolution 9.Human geography

GPW. Human population density is one of several measures of human presence and activity which together define the human “footprint,” associated with profound, adverse effects on natural systems78. Given this pervasive impact, data characterizing degree of human influence are used as predictors in some ecological models, including SDMs28. The Gridded Population of the World (GPW)68 density dataset represents the global distribution of human population density, developed using census records, population registers, and the administrative boundaries of approximately 13.5 million national and subnational units. Density, measured by population count per square kilometer, was estimated every five years, from 2000 to 2020, inclusive. The source raster dataset for each year, at 30 arc-second resolution, was summarized by area-weighted mean at ISEA3H resolutions 6 to 9.

Land cover

MCD12Q1. The Moderate Resolution Imaging Spectroradiometer (MODIS) land cover type (MCD12Q1) dataset69 describes land cover globally, via six different classification schemes. The Eco-ISEA3H database37 includes land cover classified via the International Geosphere-Biosphere Programme (IGBP) scheme, initially developed for the DISCover dataset79; the IGBP scheme includes 16 land cover classes, 13 natural and three anthropogenically modified. The MCD12Q1 dataset is derived from reflectance data collected by the MODIS instruments aboard the Terra and Aqua satellites; the two instruments observe the entirety of the Earth’s surface every one to two days, recording reflectance in 36 spectral bands.

MCD12Q1 land cover is estimated annually. For each year, reflectance time-series data are smoothed and gap-filled via smoothing splines; derived spectro-temporal features are used as input to a random forest classifier; and output land cover classifications are post-processed, to incorporate prior knowledge and reduce inter-annual variability69. The source raster dataset for 2001 and 2014 to 2018, inclusive, at approximately 500 meter resolution, was summarized by centroid, fraction, and mode at ISEA3H resolutions 5 to 10. Variable codes and descriptions are listed in Table 8.Table 8 Codes and descriptions for the 16 IGBP land cover classes, from Friedl and Sulla-Menashe69.Full size tableMOD44B. The MODIS vegetation continuous fields (VCF) dataset (MOD44B)70 describes global land cover quantitatively, as fractions of three cover components: tree canopy, non-tree canopy, and non-vegetated, barren cover. Note that canopy cover, as defined here, indicates the area over which light is intercepted; this differs from crown cover, which indicates the area covered by a plant’s crown regardless of light interception/penetration. The MOD44B dataset is derived from reflectance data collected by the MODIS instrument aboard the Terra satellite; for each annual VCF estimate, reflectance time-series data are used as input to a bagged ensemble of linear regression trees70. The source raster dataset for 2018, at approximately 250 meter resolution, was summarized by area-weighted mean at ISEA3H resolution 9.Physical geography

Biogeographic Realms. As defined by Olson et al.56, the eight terrestrial biogeographic realms are the broadest divisions of the Earth’s terrestrial flora and fauna; these may be further subdivided into biomes and ecoregions, the latter containing distinct natural communities. Olson et al.56 developed this hierarchical system primarily for global and regional conservation planning. Realm, biome, and ecoregion delineations are based on expert knowledge, contributed by more than 1,000 scientists working in relevant fields; these divisions thus incorporate knowledge of endemic taxa, unique species assemblages, and local geological and biogeographical history56. Realms were included in the Eco-ISEA3H database37 to provide a high-level classification of the Earth’s biogeography, from a source frequently cited in the scientific literature. The source vector dataset was summarized by fraction and mode at ISEA3H resolutions 5 to 9. Variable codes and descriptions are listed in Table 9.Table 9 Codes and descriptions for the eight biogeographic realms, from Olson et al.56.Full size table

ENVIREM. In addition to the climatic variables discussed previously, the ENVIREM dataset60 contains two topographic variables, derived from SRTM30_PLUS. These two indices characterize terrain roughness, a measure of variability in local elevation; and topographic wetness, a function of slope and upgradient contributing area. Source rasters at 30 arc-second resolution were summarized by area-weighted mean at ISEA3H resolutions 8 and 9. Variable codes, descriptions, and units are listed in Table 10.Table 10 Codes, descriptions, and units for the two ENVIREM topographic variables, from Title and Bemmels60.Full size table

GLWD. The Global Lakes and Wetlands Database (GLWD)71, Level 3, represents the maximum extent of lakes, reservoirs, rivers, and a number of wetland types, comprising 12 waterbody classes in total. Lehner and Döll71 compiled the three levels of the GLWD by combining seven source map and attribute datasets, and suggest Level 3 may be useful as input in global hydrologic and climatic modeling. The source raster dataset at 30 arc-second resolution was summarized by fraction and mode at ISEA3H resolution 9. Variable codes and descriptions are listed in Table 11.Table 11 Codes and descriptions for the 12 GLWD waterbody classes, from Lehner and Döll71.Full size table

Natural Earth. Natural Earth is a public-domain collection of raster and vector datasets developed for production cartography. Three vector themes describing physical geography were incorporated: Land, which includes continents and major islands; Islands, which includes additional minor islands; and Lakes, which includes lakes and reservoirs. Source vector datasets at 1:10,000,000 scale were summarized by fraction at ISEA3H resolutions 5 to 9. Further, fractions for a Terra theme were calculated, by adding per-cell Land and Islands, and subtracting Lakes. The Terra theme may be thresholded (for example, at a fraction value ≥0.5) to identify terrestrial ISEA3H cells, excluding cells covered primarily by ocean or freshwater habitat.

SRTM30_PLUS. The SRTM30_PLUS dataset72 is a global digital elevation model (DEM), representing the Earth’s terrestrial topography and ocean bathymetry. A number of elevation sources were incorporated in developing the DEM; terrestrial topography was derived from the Shuttle Radar Topography Mission (SRTM) at latitudes between ±60°, from GTOPO30 in the Arctic, and from GLAS/ICESat in the Antarctic. Ocean bathymetry was derived from satellite radar altimetry, calibrated on 298 million corrected ship-based depth soundings, gathered from several sounding sources72. The source raster dataset at 30 arc-second resolution was summarized by area-weighted mean at ISEA3H resolutions 6 to 10.

Species rangesFrom the Red List and the Phylogenetic Atlas, the geographic ranges of species belonging to four mammalian orders were sampled: Artiodactyla (even-toed ungulates), Perissodactyla (odd-toed ungulates), Primates, and Proboscidea (elephants). These species are primarily large-bodied herbivores, and as such are frequently the subject of dental ecometrics research; for example, averaged dental traits of communities of these mammals have been used to predict measures of local precipitation, at both global3 and regional11 scales.

IUCN Red List. The International Union for Conservation of Nature’s (IUCN) Red List of Threatened Species50 comprises global assessments of the conservation status of nearly 150,000 floral, faunal, and fungal species. The Red List includes expert-delineated geographic ranges for most of these species, including most extant mammalian species. For each species, portions of the range for which the species’ presence was coded extant, and for which its origin was coded native or reintroduced, were sampled. Source vector datasets were summarized by fraction at ISEA3H resolutions 8 to 9 (Artiodactyla and Perissodactyla), 9 (Primates), and 7 to 9 (Proboscidea).

PHYLACINE. The Phylogenetic Atlas of Mammal Macroecology (PHYLACINE)73,74 includes trait, phylogeny, and geographic range data for all mammalian species known from the last interglacial period (approximately 130,000 years ago) to the present, both extant and recently extinct. PHYLACINE includes species’ ranges under two scenarios, both of which were incorporated: present-day ranges, from the IUCN v2016.3; and “present-natural” ranges, for which each species’ present-day range was modified to estimate its distribution under current climatic conditions, but absent anthropogenic pressure. This included, among eight modification categories, reconnecting fragmented ranges, by filling suitable intervening habitat; and expanding ranges reduced by human activity, by filling suitable adjacent habitat. Present-natural range modifications are documented for each species in PHYLACINE’s metadata, and intended to mitigate human impact on the results of macroecological analysis and modeling. Source rasters at approximately 100 kilometer resolution were summarized by centroid at ISEA3H resolution 9. More