Ecology

Field observations
Since 1999, the FIA program has operated an extensive, nationally consistent forest inventory designed to monitor changes in forests across all lands in the US61. We used FIA data from 10 states in the continental western US (Washington, Oregon, California, Idaho, Montana, Utah, Nevada, Colorado, Arizona, and New Mexico) to quantify shifts in relative live tree density, excluding Wyoming due to a lack of repeated censuses (Fig. 1). This region spans a wide variety of climatic regimes and forest types, ranging from temperate rain forests of the coastal Pacific Northwest to pinyon-juniper woodlands of the interior southwest62. Although the spatial extent of the FIA plot network represents a large portion of the current range of all species examined in this study (Table 1), substantial portions of some species ranges (e.g., Douglas-fir) extend beyond the study region into Canada and/or Mexico and therefore were not fully addressed here.
The FIA program measures forest attributes on a network of permanent ground plots that are systematically distributed at a rate of ~1 plot per 2428 hectares across the US61. For trees, 12.7 cm DBH and larger, attributes (e.g., species, DBH, live/dead) are measured on a cluster of four 168 m2 subplots61. Trees 2.54–12.7 cm DBH are measured on a microplot (13.5 m2) contained within each subplot, and rare events such as very large trees are measured on an optional macroplot (1012 m2) surrounding each subplot61. In the event a major disturbance (i.e., >1 acre in size, resulting in mortality or damage to >25% of trees) has occurred between measurements on a plot, FIA field crews record the primary disturbance agent (e.g., fire) and estimated year of the event. In the western US, one-tenth of ground plots are measured each year, with remeasurements first occurring in 2011. Please see Data Availability for more information on forest inventory data accessibility.
Forest stability index
Allometric relationships between size and density of live trees make it difficult to interpret many indices of forest change19. Live tree density is expected to decline as trees grow in size, owing to increased individual demand for resources and growing space (i.e., competition)16,23. The expected magnitude of change in tree density, given some change in average tree size, varies considerably across forest communities63, site conditions64, and stand age classes23. Thus, we posit it is useful to contextualize observed changes in live tree density relative to those expected given shifts in average tree size within a stand. To this end, we developed the FSI, a measure of change in relative live tree density that can be applied in stands of any forest community and/or structural type.
To compute the FSI, we first develop a model of maximum size-density relationships for tree populations in our study system (Fig. 6). This model describes the theoretical maximum live tree density (({N}_{max }); in terms of tree number per unit area) attainable in stands as a function of their average tree size ((overline{S})) and will be used as a reference curve to determine the proportionate live tree density of observed stands (i.e., observed density with respect to theoretical maximum density). We use average tree basal area as an index of tree size (one, however, could also use biomass, volume, or other indices of tree size). For stand-type i, the general form of the maximum tree size-density relationship is given by

$${N}_{max }({bar{S}}_{i})={a}_{i}cdot {bar{S}}_{i}^{ {r}_{i}},$$
(1)

where a is a scaling factor that describes the maximum tree density at (bar{S}=1) and r is a negative exponent controlling the decay in maximum tree density with increasing average tree size. Such allometric size-density relationships (i.e., power functions) are widely accepted as quantitative law describing the behavior of even-aged plant populations under self-thinning conditions16,17, and have been used extensively to describe relative stand density in forests23,24. As detailed below, we allow both a and r to vary with stand-type i, as maximum size-density relationships have been shown to vary across forest communities and ecological settings63,65. Allowing a and r to vary by forest community type, for example, allows us to acknowledge that the maximum tree density attainable in a lodgepole pine stand is likely to differ from that of a pinyon-juniper stand with the same average tree size.
Fig. 6: Maximum size-density relationship for an example stand-type.

Individual points represent observed stand-level indices of tree density (N) and average tree size ((overline{S})). Maximum tree density (({N}_{max })) is modeled as a power function of average tree size within a stand. Here, we use quantile regression to estimate ({N}_{max }) as the 99th percentile of N conditional on (overline{S}). The resulting maximum size-density curve can then be used to compute the relative density of observed stands (RD), where relative density is defined as ratio of observed tree density (N) to maximum theoretical density (({N}_{max })), given (overline{S}). Source data are provided as a Source Data file.

Full size image

We next define an index of the relative density of a population of trees j (e.g., species, Pinus edulis) within a stand of type i (e.g., forest community type, pinyon/juniper woodland)

$${{rm{RD}}}_{ij}=sum _{h=1}frac{{N}_{hij}}{{N}_{max }({S}_{hi})},$$
(2)

where N is the density represented by tree h (in terms of tree number per unit area), and S is an index of individual-tree size (e.g., basal area, as used here). The denominator of Equation (2) represents the maximum tree density attainable in a stand of type i with average tree size equal to the size of tree h. We therefore express the relative density of a population j within stand-type i as a sum of the relative densities represented by individual trees within the stand. RD can be interpreted as the proportionate density, or stocking, of a population of trees within stand, where values range from 0 (population j is not present within a stand) to 1 (population j constitutes 100% of a stand and the stand is at maximum theoretical density given its size distribution). As we do in this study, one may apply any range of estimators to summarize the expected relative density of a population of trees j across a range of different stand-types (e.g., estimate the mean and variance of RDj across a broad region containing many different stand-types).
It is important to note that Equation (2) is approximately equal to a simpler method using aggregate indices (i.e., (frac{{sum }_{h = 1}{N}_{hij}}{{N}_{max }(overline{{S}_{i}})})) when tree size-distributions are normally distributed (even age-structures). However, the use of aggregate indices introduces class aggregation bias that results in overestimation of relative density in stands with non-normal size distributions (i.e., uneven age-structures), consistent with other indices of relative tree density66. In contrast, summing tree-level relative densities eliminates such bias and allows RD to accurately compare density conditions across stands in very different structural settings (e.g., even-aged plantation vs. irregularly structured old forest). Furthermore, the partitioning of relative density into tree-level densities allows RD to be accurately summarized within tree size-classes66. That is, it is possible to explicitly estimate the contribution of tree size-classes to overall stand density using RD.
For a given population j within stand-type i, we define the FSI as the average annual change in relative tree density observed between successive measurements of a stand

$${rm{FSI}}=frac{{{Delta }}{rm{RD}}}{{{Delta }}t},$$
(3)

where Δt is the number of years between successive measurement times t1 and t2 and ΔRD is the change in RD over Δt (i.e., ({{rm{RD}}}_{{t}_{2}}-{{rm{RD}}}_{{t}_{1}})). The FSI may also be expressed in units of percent change (%FSI), where average annual change in relative tree density is standardized by previous relative density

$$% {rm{FSI}}=frac{100cdot {rm{FSI}}}{{{rm{RD}}}_{{t}_{1}}}.$$
(4)

Here, stability is defined by zero net change in relative tree density over time (i.e., FSI equal to zero), but does not imply zero change in absolute tree density or tree size distributions. For example, a population exhibiting a decrease in absolute tree density (e.g., trees per unit area) may be considered stable if such decline is offset by a compensatory increase in average tree size. However, populations exhibiting expansion (i.e., ({{rm{RD}}}_{{t}_{1}} {,}{{rm{RD}}}_{{t}_{2}})) in relative tree density will be characterized by positive and negative FSI values, respectively.
Statistical analysis
We computed the FSI for all remeasured FIA plots in the western US (N = 24,229). We included plots on both public and private lands and considered all live stems (DBH ≥2.54 cm) in our analysis. As forest management can effect regional shifts in tree density, we excluded plots with evidence of recent (i.e., within 5 years of initial measurement) silvicultural treatment (e.g., harvesting, artificial regeneration, site preparation). All plot measurements occurred from 2001 to 2018, with an average remeasurement interval of 9.78 years (±0.005 years). For brevity, we restricted our analysis to consider the eight most abundant tree species in the western US. We identified the most abundant tree species using the rFIA R package60, defining abundance in terms of estimated total number of trees (DBH ≥ 2.54 cm) in the year 2018. We excluded species that exhibit non-tree growth habits (i.e., shrub, subshrub) across portions of the study region. All statistical analysis was conducted in Program R (4.0.0)67.
We developed a Bayesian quantile regression model to estimate maximum size-density relationships for stand-types observed within our study area. Here, we use TPH as an index of absolute tree density, average tree basal area ((overline{{rm{BA}}}); equivalent to tree basal area per hectare divided by TPH) as an index of average tree size, and forest community type to describe stand-types. We produced stand-level estimates of TPH and (overline{{rm{BA}}}) from the most recent measurements of FIA plots that (1) lack evidence of recent (within remeasurement period or preceding 5 years) disturbance and/or silvicultural treatment and (2) exhibit approximately normal tree diameter distributions (i.e., even-aged). Here we define an approximately normal tree diameter distribution as exhibiting Pearson’s moment coefficient of skewness between −1 and 1.
We transform the nonlinear size-density relationship to a linear function by taking the natural logarithm of TPH and (overline{{rm{BA}}}), and use a linear quantile mixed-effects model to estimate the 99th percentile of TPH conditional on (overline{{rm{BA}}}) (i.e., in log-log space) for all observed forest community types. We allowed both the model intercept and coefficient to vary across observed forest community types (i.e., random slope/intercept model), thereby acknowledging variation in the scaling factor (a) and exponent (r) of the maximum tree size-density relationship across stand-types. We place informative normal priors on the model intercept (μ = 7, σ = 1) and coefficient (μ = 0.8025, σ = 0.1) following the results of decades of previous research in maximum tree size-density relationships16,23,63,65.
The FIA program uses post-stratification to improve precision and reduce non-response bias in estimates of forest variables68, and we used these standard post-stratified estimators to estimate the mean and variance of the FSI for each species across their respective ranges within the study area (see Code Availability for all relevant code). Further, the FIA program uses an annual panel system to estimate current inventories and change, where inventory cycles consist of multiple panels, and individual panels are comprised of mutually exclusive subsets of ground plots measured in the same year within a region. Precision of point and change estimates can often be improved by combining annual panels within an inventory cycle (i.e., by augmenting current data with data collected previously). We used the simple moving average estimator implemented in the rFIA R package60 to compute estimates from a series of eight annual panels (i.e., sets of plots remeasured in the same year) ranging from 2011 to 2018. The simple moving average estimator combines information from annual panels with equal weight (i.e., irrespective of time since remeasurement), thereby allowing us to characterize long-term patterns in relative density shifts. We determine populations to be stable if the 95% confidence intervals for range-averaged FSI included zero. Alternatively, if confidence intervals of range-averaged FSI do not include zero, we determine the population to be expanding when the estimate is positive and declining when the estimate is negative.
To identify changes in species-size distributions, we used the simple moving average estimator to estimate the mean and variance of the FSI by species and size class across the range of each species within our study area. We assign individual trees to size-classes representing 10% quantiles of observed diameter distributions (i.e., diameter at 1.37 m above ground) of each species growing on one of seven site productivity classes (i.e., inherent capacity of a site to grow crops of industrial wood). That is, we allow size class definitions to vary among species and along a gradient of site productivity, thereby acknowledging intra-specific variation in diameter distributions arising from differences in growing conditions. The use of quantiles effectively standardizes absolute size distributions, simplifying both intra-specific and inter-specific comparison of trends in relative density shifts along species-size distributions.
We assessed geographic variation in species relative density shifts at two scales: ecoregion divisions and subsections69. Ecoregion divisions (shown for our study area in Fig. 1) are large geographic units that represent broad-scale patterns in precipitation and temperature across continents. Ecoregion subsections are subclasses of ecoregion divisions, differentiated by variation in climate, vegetation, terrain, and soils at much finer spatial scales than those represented by divisions. We again used the simple moving average estimator to estimate the mean and variance of the FSI by species within each areal unit (i.e., drawing from FIA plots within each areal unit to estimate mean and variance of the FSI). As a direct measure of changes in relative tree density, spatial variation in the FSI is indicative of spatial shifts in species distributions during the remeasurement interval (i.e., range expansion/contraction and/or within-range relative density shifts). That is, the distribution of populations shift toward regions increasing in relative density and away from regions decreasing in relative density during the temporal frame of sampling. We map estimates of the FSI for each areal unit to assess spatial patterns of changes in relative density and identify regions where widespread geographic shifts in species distributions may be underway.
We sought to quantify the average effect of forest disturbance processes on changes in the relative density of top tree species in the western US over the interval 2001–2018. To this end, we developed a hierarchical Bayesian model to determine the average severity and annual probability of disturbances (i.e., wildfire, insect outbreak, and disease) on sites where each species occurs. Average severity was modeled as

$${y}_{jk} sim {rm{normal}}({alpha }_{j}+sum _{l}{beta }_{jl}cdot {x}_{lk}, {varsigma }_{j}^{2}),$$
(5)

where yjk is the FSI of species j on plot k, αj is a species-specific intercept, βjl is a species-specific coefficient corresponding to the binary variable xlk that takes the value of 1 if disturbance l occurred within plot k measurement interval and 0 otherwise. The intercept and regression coefficients each received an uninformative normal prior distribution. The species-specific residual standard deviation ςj received a uninformative uniform prior distribution70.
On average, disturbance will occur at the midpoint of plot remeasurement periods, assuming temporal stationarity in disturbance probability over the study period. As plots in this study are remeasured on 10-year intervals, we assume that tree populations have, on average, 5 years to respond to any disturbance event. Hence, our definition of disturbance severity, βjl’s, cannot be interpreted as the immediate change in relative tree density resulting from disturbance. Rather, disturbance severity (as defined here) includes the immediate effects of disturbance, as well as 5 years of change in relative tree density prior to and following disturbance (where disturbance is assumed to be functionally instantaneous).
Annual probability of disturbance l on plot k was modeled as

$${x}_{lk} sim {rm{binomial}}({{Delta }}{t}_{k},{psi }_{jl}),$$
(6)

where Δtk is the number of years between successive measurements of plot k, viewed here as the number of binomial “trials,” and ψjl is the species-specific probability for disturbance which was assigned a beta(1,1) prior distribution. Hence, annual probability of disturbance is assumed to vary by species j and by disturbance type l.
We estimate the mean effect of forest disturbance processes on changes in species-specific relative tree density by multiplying the posterior distributions of βjl and ψjl. That is, we multiply species-specific disturbance severity by disturbance probability to yield an estimate of the mean change in relative density caused by disturbance over the study period. We then standardize these values across species by dividing by the average relative density of each species at the beginning of the study period. Thus, standardized values can be interpreted as the annual proportionate change in the relative tree density of each species resulting from disturbance over the period 2001–2018.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

Six sponge species, R. odorabile, C. matthewsi, C. foliascens, S. flabelliformis, I. ramosa and C. orientalis (a bioeroding sponge), were selected for metagenomic sequencing (7 ± 0.5 Gbp) as these species represent dominant habitat forming taxa on tropical and temperate Australian reefs and exhibit high intraspecies similarity in their microbiomes. In addition, previously published microbial MAGs from I. ramosa and Aplysina aerophoba were analysed [8, 12], including 62 additional unpublished MAGs from A. aerophoba. The recovered MAGs, averaging 86 ± 12% completeness and 2 ± 2% contamination, made up 72 ± 21% relative abundance of their respective communities (by read mapping) on average and spanned the vast majority of microbial lineages typically seen in marine sponges [45] (Fig. S1 and Table S1), including the bacterial phyla Proteobacteria (331 MAGs), Chloroflexota (242), Actinobacteriota (155), Acidobacteriota (97), Gemmatimonadota (60), Latescibacterota (44; including lineages Anck6, PAUC34 and SAUL), Cyanobacteria (43), Bacteroidota (38), Poribacteria (35), Dadabacteria (22; including SBR1093), Nitrospirota (22), Planctomycetota (15), UBP10 (14), Bdellovibrionota (13), Patescibacteria (9; includes Candidate Phylum Radiation), Spirochaetota (8), Nitrospinota (7), Myxococcota (4), Entotheonella (2) and the archaeal class Nitrososphaeria (21; phylum Crenarchaeota), hereafter referred to by their historical name “Thaumarchaeota” for name recognition. Mapping of the metagenomic reads to the recovered MAGs showed that the communities had high intraspecies similarity across replicates, consistent with previous 16S rRNA gene-based analyses (Fig. S1). In general, taxa present in A. aerophoba, C. foliascens, C. orientalis and S. flabelliformis appeared unique to those sponge species, with only one dominant lineage present in C. orientalis (order Parvibaculales). In contrast, several Actinobacteriota, Acidobacteriota and Cyanobacteria populations were shared across C. matthewsi, R. odorabile and I. ramosa. Further, members of the Thaumarchaeota were detected in all sponge species and were particularly abundant in S. flabelliformis at 12 ± 4% relative abundance (Fig. S1). Addition of these sponge MAGs to genome trees comprising all publicly available sponge symbionts (N = 1188 MAGs) resulted in a phylogenetic gain of 44 and 75% for Bacteria and Archaea, respectively, reflecting substantial novel genomic diversity (Fig. 1).
Comparative genomic analysis of the sponge-derived MAGs provided unique insights into the distribution of metabolic pathways across sponge symbiont taxa. For example, microbial oxidation of ammonia benefits the sponge host by preventing ammonia from accumulating to toxic levels [46], a process thought to be mediated by both symbiotic Bacteria and Archaea (i.e. Thaumarchaeota) [33]. Prior identification of ammonia oxidisers has been based on functional inference from phylogeny (16S rRNA gene amplicon surveys) [47] or homology to specific Pfams (metagenomes) [33]. However, the CuMMO gene family is diverse, encompassing functionally distinct relatives that include amoA, particulate methane monooxygenases and hydrocarbon monooxygenases that cannot be distinguished by homology alone [35]. We used GraftM [32] to recover CuMMO genes from the sponge MAGs and their metagenomic assemblies, as well as previously sequenced metagenomic assemblies from six additional sponge microbiomes where bacterial amoA gene sequences had been identified [33]. Phylogenetic analysis of the recovered CuMMO genes showed that all archaeal homologues came from Thaumarchaeota and fell within the archaeal amoA clade. In contrast, bacterial CuMMO sequences were identified exclusively in MAGs from the phylum UBP10 (formerly unclassified Deltaproteobacteria) and from an unknown taxonomic group in the previous metagenomic assemblies [33]. All recovered bacterial and taxonomically unidentified CuMMO placed within the Deltaproteobacteria/Actinobacteria hmo clade, indicating these genes are specific for hydrocarbons rather than ammonia (Fig. S2). The finding that Thaumarchaeota are the only microbes within any of the surveyed sponge species capable of oxidising ammonia, and their ubiquity across sponges, suggests they are a keystone species for this process.
To further investigate the distribution of functions within the sponge microbiome, a set of highly complete ( >85%) sponge symbiont MAGs were grouped by principal components analysis based on their KEGG and Pfam annotations, as well as orthologous clusters that reflected all gene content. Similar analysis conducted on 37 MAGs from the sponge Aplysina aerophoba suggested the presence of functional guilds, with MAGs from disparate microbial phyla carrying out similar metabolic processes [12] (e.g. carnitine catabolism). Here, we find that MAGs clustered predominately by microbial taxonomy (phylum) rather than function in all three analyses (Fig. S3). While functional guilds could not be identified based on analysis of total genome content, this does not preclude the existence of such guilds based on more specific metabolic pathways.
To identify pathways enriched within the sponge microbiome, sponge-associated MAGs with >85% completeness (N = 798) were compared with a set of coral reef and coastal seawater MAGs (N = 86), 31 derived from published datasets [31] and 55 from this study (Table S1). Seawater MAGs with >85% genome completeness (93 ± 4% completeness and 2 ± 2% contamination; Table S1) spanned the bacterial phyla Proteobacteria (48 MAGs), Bacteroidota (13), Planctomycetota (5), Myxococcota (5), Gemmatimonadota (3), Marinisomatota (3), Actinobacteriota (3), Verrucomicrobiota (2), Cyanobacteriota (2), Bdellovibrionota (1) and the archaeal phylum Nanoarchaeota (1). Comparative analysis revealed that sponge symbionts were enriched in metabolic pathways for carbohydrate metabolism, defence against infection by MGE, amino acid synthesis, eukaryote-like gene repeat proteins (ELRs) and cell–cell attachment (Tables S2–S4).
Genes belonging to GH and carbohydrate esterase (CE) families (Table S2) acting on starch (GH77), arabinose (CAZY families GH127 and GH51), fucose (GH95 and GH29) and xylan polymers (CE7 and CE15), were enriched in sponge-associated lineages, likely reflecting the hosts critical role in catabolising dissolved organic matter (DOM) present in reef seawater (Fig. 2). Microbial GHs from the GH77 family target starch, the main sugar storage compound in marine algae [48], whereas GHs from families 51 and 127 are known to act on plant arabinosaccharides, such as the hydroxyproline-linked arabinosaccharides found in algal extensin glycoproteins [49, 50]. GH127 enzymes are also required for microbial degradation of carrageenan, a complex heteropolysaccharide produced by red algae [51]. Members of the fucosidase GH95 and GH29 enzyme families are known to degrade fucoidan, a complex fucosaccharide prominent in brown algae [50, 52]. Notably, arabino- and fucopolysaccharides also make up a significant proportion of coral mucus, a major component of DOM in coral reefs that sponges have been shown to utilise [53, 54]. Supporting this observation, isotopic investigation of the fate of coral mucus and algal polysaccharides in sponges showed that the microbiome participates in metabolism of these compounds, particularly in sponges with high microbial abundance and diversity [4, 5]. Enzymes from the CE families 15 and 7 have been primarily characterised in terrestrial plants where they act as glucuronyl esterases and acetyl-xylan esterases, degrading lignocellulose and removing acetyl groups from hemicellulose [55] (e.g. xylans). Characterisation of CE15 and CE7 from marine microbes is rare, though activity on xylans, which are a structural component of marine algae, has previously been demonstrated [55,56,57].
Fig. 2: Phylogenetic tree showing the distribution of glycosyl hydrolases and esterases across MAGs with >85% completeness (N = 884).

Values represent the copy number of each gene per MAG. Internal branches of the tree are coloured by phylum, while the outer strip is coloured by class. Both are listed clockwise in the order in which they appear. Seawater MAGs are denoted by grey labels with red text.

Full size image

GHs acting on sialic acids (GH33) and glycosaminoglycans (GH88) were also enriched in the sponge-associated MAGs and may act on compounds found within sponge tissue [13] (Fig. 2). In contrast, no genes for the degradation of collagen (collagenases), one of the main structural components of the sponge skeleton were identified. Sialic acid-linked residues are found in the sponge mesohyl [58], and although the impact of cleavage on the host is unknown, analogy can be made to other symbioses. For example, sialidases are common in the commensal bacteria present in the human gut where they are used to cleave and metabolise the sialic acid-containing mucins lining the gut wall [59]. Increased sialidase activity is associated with gut dysbiosis and inflammation [60] and careful control of sialidase-containing commensals is therefore necessary to maintain gut homoeostasis [59]. As glycosaminoglycans are also part of sponge tissue [13, 61], the same may apply to microorganisms encoding GH88 family enzymes. However, these genes are also implicated in the degradation of external sugar compounds, such as ulvans, a major sugar storage compound found in green algae that can make up to 30% of their dry weight [62]. Thus, the ecological role of GH88 family enzymes within the sponge microbiome requires further investigation.
Enrichment of GHs and CEs was largely restricted to the Poribacteria, Latescibacteria (class UBA2968), Spirochaetota, Chloroflexota (classes UBA2235 and Anaerolineae, but not Dehalococcoides) and Acidobacteriota (class Acidobacteriae). These findings corroborate previous targeted genomic characterisations of the Chloroflexota and Poribacteria [13, 14] but show that they are part of a larger set of polysaccharide-degrading lineages. Identification of disparately related microbial taxa across several sponge lineages (Figs. 1 and 2) that encode similar pathways for polysaccharide degradation, and therefore occupy a similar ecological niche, supports the existence of functional guilds within the sponge microbiome when viewed at the level of individual pathways. Given the fundamental role of marine sponges in recycling coral reef DOM, studies targeting these specific guilds are needed to quantify their contribution to reef DOM transformation.
Because sponges filter and retain biomass from an extensive range of reef taxa (eukaryotic algae, bacteria, archaea, etc), they are exposed to a greatly expanded variety of MGEs from these organisms, including viruses, transposable elements and plasmids [33, 63]. For this reason, sponge-associated microorganisms likely require a diverse toolbox of molecular mechanisms for resisting infection. Both RM and CRISPR systems are capable of recognising and cleaving MGEs as part of the bacterial immune repertoire. RM systems are part of the innate immune system of bacteria and archaea and are encoded by a single (Type II) or multiple proteins (Type I, III and IV) that recognise and cleave foreign DNA based on a defined target sequence. In contrast, CRISPR systems are part of the adaptive immune system of some bacteria and archaea and encode a target sequence derived from the genome of a previous infective agent that is used by a CRISPR-associated protein (CAS) to identify and cleave foreign DNA. RM (Fig. S4) and CAS (Fig. S5) genes were both enriched (Table S3) in the sponge-associated MAGs and relatively evenly distributed across taxa, with the exception of the Planctomycetota and Verrucomicrobiota, where they were largely absent. As these MAGS average 93 +/− 5% completeness, this result is not likely due to genome incompleteness. This finding contrasts with comparative investigations of Planctomycetota genomes from other environments [64] and additional research is required to ascertain the mechanisms used by sponge-associated Planctomycetota and Verrucomicrobiota to avoid infection. Although Type III RM genes were enriched in sponge MAGs, they were also present in all seawater MAGs. In contrast, Types I and II RM genes were present almost exclusively in the sponge-associated MAGs. In conjunction with an enrichment in CRISPR systems, this expanded repertoire of defence systems likely reflects the increased burden from MGEs associated with the hosts role in filtering and concentrating diverse sources of reef biomass. Supporting this hypothesis, metagenomic surveys of sponge-associated viruses revealed a more diverse viral population than what could be recovered from the surrounding seawater [63]. Further, we found that genes encoding toxin-antitoxin systems, which are present on MGEs, such as plasmids, were also enriched in sponge-associated MAGs. These observations suggest that RM and CRISPR systems are important features of microbe-sponge symbiosis, allowing the symbionts to colonise and persist within their host by avoiding viral infection or being overtaken by MGEs.
Pathways for the synthesis of amino acids were also enriched in the sponge microbiome. The inability of animals to produce several essential amino acids has been proposed as a primary reason that they harbor microbial symbionts [65,66,67,68] and it has long been thought that sponges acquire at least some of their essential amino acids from their microbiome [69, 70]. Further, gene-centric characterisation of the Xestospongia muta and R. odorabile microbiomes revealed pathways to synthesise and transport essential amino acids [33, 70]. However, these same amino acid pathways are also used catabolically by the microorganisms, and transporters could simply be importing amino acids into the microbial cell. Further, as sponges are almost constantly filter feeding, essential amino acids could be acquired through consumption of microorganisms present in seawater. Comparison of sponge MAGs with those from seawater revealed enrichment of specific pathways for the synthesis of lysine, arginine, histidine, threonine, valine and isoleucine (Table S4). However, visualisation of the distribution of these genes revealed that almost all MAGs in both sponges and seawater produce all amino acids, though specific lineages may use different pathways to achieve this (Fig. S6). The enrichment observed in the sponge MAGs was therefore ascribed to differences in pathway completeness between sponge-associated and seawater microbes, rather than an enhanced ability of sponge symbionts to produce any specific amino acid. In contrast, compounds, such as taurine, carnitine and creatine have also been proposed as important host-derived carbon sources for symbionts [69], but pathways for their catabolism were enriched in seawater rather than sponge-associated MAGs. While these findings do not invalidate the possibility that microbial communities play a role in amino acid provisioning to the host or that they utilise host-derived taurine, carnitine, or creatine, they suggest that these are not key processes mediating microbe-sponge symbiosis.
To form stable symbioses, bacteria must persist within the sponge tissue and avoid phagocytosis by host cells. Microbial proteins containing ELR motifs have been identified in a range of animal and plant-associated microbes and are thought to modulate the host’s intracellular processes to facilitate stable symbiotic associations [71, 72]. For example, ELR-containing proteins from sponge-associated microbes have been shown to confer the ability to evade host phagocytosis when experimentally expressed in E. coli [10, 73]. ELR-containing proteins from the ankyrin (ARP), leucine-rich, tetratricopeptide and HEAT repeat families were enriched in the sponge-associated MAGs. In contrast, WD40 repeats were not found to be enriched but are included here as they have previously been reported as abundant in Poribacteria and symbionts of other marine animals [13, 31]. Most ELRs were present across all taxa but were much more prevalent in specific lineages (Fig. 3). For example, sponge-associated Poribacteria, Latescibacterota and Acidobacteriota encoded a high proportion of all ELR types, while other lineages, such as the Gemmatimonadota (average 0.25% coding genes per sponge-associated MAG versus 0.09% in seawater MAGs), Verrucomicrobiota (2%), Deinococcota (0.85%), Acidobacteriota (0.20%; specifically class Luteitaleia at 0.55%) and Dadabacteria from C. orientalis (0.62%) encoded a comparatively high percentage of ARPs and Nitrospirota encoded a high percentage of HEAT_2 family proteins (0.55% versus 0.05% in seawater MAGs) relative to other taxa. In contrast, ELR abundances were substantially lower, or absent, in the Actinobacteriota, the class Bacteroidia within the phylum Bacteroidota, and the Thaumarchaeota, suggesting these microorganisms utilise alternative mechanisms to maintain their stable associations with the host.
Fig. 3: Phylogenetic tree showing the distribution of eukaryote-like repeat proteins—ankyrin (ARP), leucin-rich (LRR), tetratricopeptide (TPR), HEAT and WD40—across MAGS with >85% completeness (N = 884).

Values represent the percentage of coding genes per MAG devoted to each gene class. Internal branches of the tree are coloured by phylum, while the outer strip is coloured by class, and both are listed clockwise in the order in which they appear. MAGs from seawater are denoted by grey labels with red text.

Full size image

The mechanisms by which ELRs interact with sponge cells remains largely unknown, although microbes in other host systems are known to deliver ELR-containing effector proteins into host cells via needle-like secretion systems (types III, IV and V) or extracellular contractile injection systems [74, 75], where they interact with the cellular machinery of the host to modify its behaviour. In sponges, it is also possible that ELRs could be secreted into the extracellular space by type I or II secretion systems. Interestingly, although most sponge MAGs encoded eukaryote-like proteins (Fig. 3), few lineages encoded the necessary genes to form secretion systems (Fig. S7). It is therefore unlikely that ELRs are introduced to the sponge host via traditional secretion pathways used in other animal-symbiont systems.
Maintaining stable association with the sponge may also require mechanisms for attachment to the host tissue. For example, cadherin domains are Ca2+-dependent cell–cell adhesion proteins that are abundant in eukaryotes and have been found to serve the same function in bacteria [76]. Similarly, fibronectin III domains mediate cell adhesion in eukaryotes, but also occur in bacteria where they play various roles in carbohydrate binding and biofilm formation [77, 78]. In addition, some bacterial pathogens utilise fibronectin-binding proteins to gain entry into host tissue by binding to host fibronectin [77, 78]. Genes containing cadherin domains were enriched in the sponge-associated MAGs and were identified in most bacterial lineages, but were notably absent in the Cyanobacteriota and Verrucomicrobiota (Fig. 4). Genes containing fibronectin III domains and those for fibronectin-binding proteins were also enriched in sponge-associated MAGs and were distributed across most lineages, though were particularly abundant in the Actinobacteriota and Chloroflexota. However, although fibronectin III-containing genes were taxonomically widespread, those encoding fibronectin-binding proteins were restricted to the phyla Poribacteria, Gemmatimonadota, Latescibacterota, Cyanobacteriota, class Anaerolineae within the Chloroflexota (but not Dehalococcoidia), class Rhodothermia within the Bacteroidota, Spirochaetota, Nitrospirota and the archaeal phylum Thaumarchaeota. Interestingly, the taxonomic distribution of these genes shares significant overlap with lineages encoding the genes for sponge sialic acid and glyosaminoglycans degradation, suggesting that attachment to the host may be necessary for utilisation of these carbohydrates (Fig. 2). However, as the host, bacterial, and archaeal components of the sponge holobiont have fibronectin III domains, symbionts encoding fibronectin-binding proteins may use these to adhere to the host tissue or potentially to form biofilms (bacteria–bacteria attachment). In either case, the enrichment and wide distribution of cadherins, fibronectins and fibronectin-binding proteins in the sponge MAGs suggests that cell–cell adhesion is critical for successful establishment in the sponge niche.
Fig. 4: Phylogenetic tree showing the distribution of cadherins, fibronectins and fibronectin-binding proteins across MAGS with >85% completeness (N = 884).

Values represent the copy number of each gene per MAG. Internal branches of the tree are coloured by phylum while the outer strip is coloured by class. Both are listed clockwise in the order in which they appear. Seawater MAGs are denoted by grey labels with red text.

Full size image

Distribution of genes encoding ELRs, polysaccharide-degrading enzymes (GHs and CEs), cadherins, fibronectins, RMs and CRISPRs across distantly related taxa suggests that they were either acquired from a common ancestor or that they represent more recent LGT events, potentially mediated by MGEs, which are enriched in sponge-associated microbial communities [69]. Here, we identify 4963 LGTs from five sponges for which sufficient sequence data were available ( >100 Mbp total sequence length across all MAGs), as well as 136 LGTs from seawater MAGs, averaging 1.64 and 0.52 LGTs per Mbp sequences, respectively (Fig. 5 and Table S5). Sequence similarity of LGTs from MAGs within a sponge species was higher than between sponge species, indicating relatively recent gene transfers (Fig. S8). A higher frequency (Fig. S9) and lower genetic divergence of LGTs among MAGs derived from the same sponge species likely results from the close physical distance between members of each microbiome, as has been observed in other host-symbiont systems [79]. The identification of lateral transfers between microbes from different sponge species may highlight the horizontal acquisition of these microbes or that a recent ancestor inhabited the same host. Notably, LGTs included a subset of genes that were enriched within the sponge-associated MAGs, such as GH33 (sialidases) and CE7 (acetyl-xylan esterases), attachment proteins (cadherins and fibronectin III), RM and CAS proteins, and members of all ELR families other than WD40 (Figs. 6 and S10). The observation that a significant number of sponge-enriched genes were laterally transferred between disparate microbial lineages suggests that the processes they mediate provide a strong selective advantage within the sponge niche, though further research is required to validate these findings.
Fig. 5: Visualisation of LGTs detected within the MAGs for the five sponges passing the cumulative MAG length criteria ( >100 Mbp).

The inner strip is coloured by phylum while the outer strip is coloured by host sponges. Bands connect donors and recipients, with their colour corresponding to the donors and the width correlating to the number of LGTs.

Full size image

Fig. 6: Visualisation of gene flow among microbial phyla for gene families enriched in sponge-associated MAGs.

The inner ring and band connecting donor and recipient is coloured by protein family of the gene being transferred, with the width of the band correlating to the number of LGTs. Recipient MAGs are shown in grey. The outer ring is coloured by microbial phylum. Representation of RM and CAS gene LGTs can be found in Fig. S10.

Full size image

Sponges are important constituents of coral reef ecosystems because of their critical role in DOM cycling and retention via the sponge-loop. Despite their importance, functional characterisation of sponge symbiont communities has been restricted to just a few lineages of interest, potentially biasing our view of sponge symbiosis. Here we present a comprehensive characterisation of sponge symbiont MAGs spanning the complete range of taxa found in marine sponges (Fig. 7), most of which were previously uncharacterised. We revealed enrichment in glycolytic enzymes (GHs and CEs) reflecting specific functional guilds capable of aiding the sponge in the degradation of reef DOM. Further, we identified several ELRs, CRISPRs and RMs that likely facilitate stable association with the sponge host, showing the specificity of ELR types with individual microbial lineages. We also clarified the role of Thaumarchaeota as a keystone taxon for ammonia oxidation across sponge species and showed that processes previously thought to be important, such as amino acid provisioning and taurine, creatine and carnitine metabolism are unlikely to be central mechanisms mediating sponge-microbe symbiosis. Many of the enriched genes are laterally transferred between microbial lineages, suggesting that LGT plays an important role in conferring a selective advantage to specific sponge-associated microorganisms. Taken together, these data illustrate how evolutionary processes have distributed and partitioned ecological functions across specific sponge symbiont lineages, allowing them to occupy or share specific niches and live symbiotically with their sponge hosts.
Fig. 7: Schematic overview of microbial interactions with the host as inferred from the functional potential encoded by the sponge-associated microbial MAGs.

Fbn fibronectin, cdh cadherins, RM restriction-modification systems, CAS CRISPR-associated proteins, ELP eukaryotic-like repeat proteins, CE7 carbohydrate esterase family 7, GH33 glycosyl hydrolase family 33.

Full size image More

This announcement sent shockwaves around the globe, largely lauding it as an important step in the right direction4. China’s decision is also unprecedented at several levels, which could result in profound and far-reaching impacts for both humans and wildlife.
First, this decision was initiated and adopted by China’s highest legislature — the Standing Committee of the National People’s Congress — with the explicit endorsement of President Xi1,2. In contrast, China’s response to the SARS outbreak in 2003 — a short-lived ban on the trade and consumption of palm civets5 — was initiated by various government agencies at lower legislative levels. In responding to the COVID-19 pandemic, the political will in China has never been stronger or more overt across multiple levels of government6. Within a few months, all 31 provinces in China have published provincial legislation on wildlife farming and consumption. Perhaps more importantly, from May to July, the People’s Congress standing committees at the national and provincial levels conducted nation-wide evaluations on the effectiveness of these policies and their enforcement. The committees concluded that policies were generally well implemented, but there was room for improvement on certain aspects, including finding alternative livelihoods for affected wildlife traders, continuing to revise the protected species lists, and addressing loopholes in wildlife trade monitoring and habitat conservation. In terms of positive outcomes, joint actions and special operations from the government have closed 12,000 wildlife-related businesses, intensified monitoring efforts to include over four million e-commerce platforms, and removed 990,000 online sources of information associated with wildlife trade6.
Second, China’s current decision includes a series of new legislations to build on the achievements of current actions by enhancing the regulation of wildlife farms and markets (Box 1). The revision of the country’s Wildlife Protection Law is expected to bring about long-term and systematic changes to wildlife conservation. Additionally, China is also revising its List of Protected Animals. Species threatened by consumption, such as the pangolin and yellow-breasted bunting, are being promoted to the highest protection level (Class I Wildlife species)7. Furthermore, China’s Ministry of Agriculture and Rural Affairs published an updated Catalogue of Animal Genetic Resource in May 2020. Among wild animals, only species in this catalogue can be farmed or consumed8. There are 64 species of wild animals that are being farmed for consumption, but are not yet included in this catalogue for various reasons (for example, to reduce the risk of sourcing animals from the wild). Nevertheless, the Ministry of Forestry and Grassland has categorized them into two groups: the farming of 45 species (for example, bamboo rat and civet cat) is due to be banned by the end of 2020, and the remaining 19 species (for example, several species of snakes) are allowed to be farmed for non-consumption uses9. Furthermore, the disbursement of government financial compensation to the farmers affected by these new legislations, amounting to over a billion US dollars, is expected to be completed by the end of 2020.
Third, the current ban is likely to galvanize rapid and widespread knock-on actions and impacts. For example, Guangdong has already banned wild vertebrate animals as pets10. The consumption of dogs and cats is banned in the city of Shenzhen11. Pangolin scales are removed as a key ingredient in traditional Chinese medicine, although it is still included as an ingredient in patent medicines in the 2020 Chinese Pharmacopoeia12,13. The pre-COVID wildlife management system in China, especially in its wildlife farming industry, has long been criticized by conservationists to be disordered and outdated. The lack of incentives and capacities had been a major barrier to change for government agencies in the country. The COVID-19 pandemic, at great human and economic costs, has mainstreamed the discourse of wildlife conservation for human well-being, clarified legislations on what species can be farmed, and provided a policy framework for systematic and enforceable wildlife management and conservation. These actions are exactly what scientists have long called for to minimize the risk of zoonotic disease transmission and outbreaks in the future14. More

1.
Ferguson BA, Dreisbach TA, Parks CG, Filip GM, Schmitt CL. Coarse-scale population structure of pathogenic Armillaria species in a mixed-conifer forest in the Blue Mountains of northeast Oregon. Can J Res. 2003;33:612–23.
Article  Google Scholar 
2.
Fricker MD, Heaton LLM, Jones NS, Boddy L. The mycelium as a network. Microbiol Spectr. 2017;5:1–32.
Google Scholar 

3.
Treseder KK, Lennon JT. Fungal traits that drive ecosystem dynamics on land. Microbiol Mol Biol Rev. 2015;79:243–62.
CAS  PubMed  PubMed Central  Article  Google Scholar 

4.
Smith SE, Read DJ. Mycorrhizal symbiosis. 3rd ed. London, UK: Academic Press; 2008.
Google Scholar 

5.
Maron JL, Marler M, Klironomos JN, Cleveland CC. Soil fungal pathogens and the relationship between plant diversity and productivity. Ecol Lett. 2011;14:36–41.
PubMed  Article  Google Scholar 

6.
García-Guzmán G, Heil M. Life histories of hosts and pathogens predict patterns in tropical fungal plant diseases. New Phytol. 2013.

7.
Lin X, Alspaugh JA, Liu H, Harris S. Fungal morphogenesis. Cold Spring Harb Perspect Med. 2014;5:a019679.
PubMed  Article  CAS  PubMed Central  Google Scholar 

8.
Griffin DH. Fungal physiology. 2nd ed. New York, NY: John Wiley and Sons; 1994.
Google Scholar 

9.
Trewavas A. What is plant behaviour? Plant, Cell Environ. 2009;32:606–16.
Article  Google Scholar 

10.
Karban R. Plant behaviour and communication. Ecol Lett. 2008;11:727–39.
PubMed  Article  Google Scholar 

11.
de Kroon H, Mommer L. Root foraging theory put to the test. Trends Ecol Evol. 2006;21:113–6.
PubMed  Article  Google Scholar 

12.
Novoplansky A. Developmental plasticity in plants: implications of non-cognitive behavior. Evol Ecol. 2002;16:177–88.
Article  Google Scholar 

13.
Lovett Doust L. Population dynamics and local specialization in a clonal perennial (Ranunculus repens): II. The dynamics of leaves, and a reciprocal transplant-replant experiment. J Ecol. 1981;69:757–68.
Article  Google Scholar 

14.
Saiz H, Bittebiere A-K, Benot M-L, Jung V, Mony C. Understanding clonal plant competition for space over time: a fine-scale spatial approach based on experimental communities. J Veg Sci. 2016;27:759–70.
Article  Google Scholar 

15.
Andrews JH. Comparative ecology of microorganisms and macroorganisms. New York, NY: Springer; 1991.
Google Scholar 

16.
Carlile MJ. The success of the hypha and mycelium. In: Gow NAR, Gadd GM, editors. The growing fungus. London: Chapman & Hall; 1995. pp. 3–20.

17.
Boddy L. Saprotrophic cord-forming fungi: meeting the challenge of heterogeneous environments. Mycologia. 1999;91:13.
Article  Google Scholar 

18.
Bielčik M, Aguilar-Trigueros CA, Lakovic M, Jeltsch F, Rillig MC. The role of active movement in fungal ecology and community assembly. Mov Ecol. 2019;7:36.
PubMed  PubMed Central  Article  Google Scholar 

19.
Ritz K, Young IM. Interactions between soil structure and fungi. Mycologist. 2004;18:52–9.
Article  Google Scholar 

20.
Harris K, Young IM, Gilligan CA, Otten W, Ritz K. Effect of bulk density on the spatial organisation of the fungus Rhizoctonia solani in soil. FEMS Microbiol Ecol. 2003;44:45–56.
CAS  PubMed  Article  PubMed Central  Google Scholar 

21.
Otten W, Harris K, Young IM, Ritz K, Gilligan CA. Preferential spread of the pathogenic fungus Rhizoctonia solani through structured soil. Soil Biol Biochem. 2004;36:203–10.
CAS  Article  Google Scholar 

22.
Burges A, Nicholas DP. Use of soil sections in studying the amount of fungal hyphae in soil. Soil Sci. 1961;92:25–9.
Article  Google Scholar 

23.
Dechesne A, Wang G, Gülez G, Or D, Smets BF. Hydration-controlled bacterial motility and dispersal on surfaces. Proc Natl Acad Sci USA. 2010;107:14369–72.
CAS  PubMed  Article  PubMed Central  Google Scholar 

24.
Wolfaardt GM, Hendry MJ, Birkham T, Bressel A, Gardner MN, Sousa AJ, et al. Microbial response to environmental gradients in a ceramic-based diffusion system. Biotechnol Bioeng. 2008;100:141–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

25.
Otten W, Pajor R, Schmidt S, Baveye PC, Hague R, Falconer RE. Combining X-ray CT and 3D printing technology to produce microcosms with replicable, complex pore geometries. Soil Biol Biochem. 2012;51:53–5.
CAS  Article  Google Scholar 

26.
Deng J, Orner EP, Chau JF, Anderson EM, Kadilak AL, Rubinstein RL, et al. Synergistic effects of soil microstructure and bacterial EPS on drying rate in emulated soil micromodels. Soil Biol Biochem. 2015;83:116–24.
CAS  Article  Google Scholar 

27.
Whitesides GM. The origins and the future of microfluidics. Nature. 2006;442:368–73.
CAS  PubMed  PubMed Central  Article  Google Scholar 

28.
Stanley CE, Grossmann G, Casadevall i Solvas X, DeMello AJ. Soil-on-a-chip: microfluidic platforms for environmental organismal studies. Lab Chip. 2016;16:228–41.
CAS  PubMed  Article  PubMed Central  Google Scholar 

29.
Held M, Kaspar O, Edwards C, Nicolau DV. Intracellular mechanisms of fungal space searching in microenvironments. Proc Natl Acad Sci USA. 2019;116:13543–52.
CAS  PubMed  Article  PubMed Central  Google Scholar 

30.
Soufan R, Delaunay Y, Gonod LV, Shor LM, Garnier P, Otten W, et al. Pore-scale monitoring of the effect of microarchitecture on fungal growth in a two-dimensional soil-like micromodel. Front Environ Sci. 2018;6:68.
Article  Google Scholar 

31.
Schmieder SS, Stanley CE, Rzepiela A, van Swaay D, Sabotič J, Nørrelykke SF, et al. Bidirectional propagation of signals and nutrients in fungal networks via specialized hyphae. Curr Biol. 2019;29:217–28.e4.
CAS  PubMed  Article  Google Scholar 

32.
Aleklett K, Kiers ET, Ohlsson P, Shimizu TS, Caldas VE, Hammer EC. Build your own soil: exploring microfluidics to create microbial habitat structures. ISME J. 2018;12:312–9.
PubMed  Article  Google Scholar 

33.
Veresoglou SD, Wang D, Andrade-Linares DR, Hempel S, Rillig MC. Fungal decision to exploit or explore depends on growth rate. Micro Ecol. 2018;75:289–92.
Article  Google Scholar 

34.
Lange M, Smith, Alexander H. The coprinus ephemerus group. Mycologia. 1953;45:747–80.
Article  Google Scholar 

35.
Hernández-Rodríguez M, Oria-de-Rueda JA, Martín-Pinto P. Post-fire fungal succession in a Mediterranean ecosystem dominated by Cistus ladanifer L. Ecol Manag. 2013;289:48–57.
Article  Google Scholar 

36.
Hughes KW, Petersen RH. Transatlantic disjunction in fleshy fungi III: Gymnopus confluens. MycoKeys. 2015;9:37–63.
Article  Google Scholar 

37.
Garnier-Delcourt M, Reckinger C, Tholl M-T, Turk J. Notes mycologiques luxembourgeoises. IV. Bull Soc Nat Luxemb. 2011;112:39–50.
Google Scholar 

38.
Maynard DS, Bradford MA, Covey KR, Lindner D, Glaeser J, Talbert DA, et al. Consistent trade-offs in fungal trait expression across broad spatial scales. Nat Microbiol. 2019;4:846–53.
CAS  PubMed  Article  Google Scholar 

39.
Held M, Edwards C, Nicolau DV. Probing the growth dynamics of Neurospora crassa with microfluidic structures. Fungal Biol. 2011;115:493–505.
PubMed  Article  Google Scholar 

40.
Dowson CG, Rayner ADM, Boddy L. Spatial dynamics and interactions of the woodland fairy ring fungus, Clitocybe nebularis. N Phytol. 1989;111:699–705.
Article  Google Scholar 

41.
Fukasawa Y, Savoury M, Boddy L. Ecological memory and relocation decisions in fungal mycelial networks: responses to quantity and location of new resources. ISME J. 2020;14:380–8.
PubMed  Article  PubMed Central  Google Scholar 

42.
Held M, Binz M, Edwards C, Nicolau DV. Dynamic behaviour of fungi in microfluidics: a comparative study. In: Proc. SPIE 7182, Imaging, Manipulation, and Analysis of Biomolecules, Cells, and Tissues VII. 2009. pp. 718213. https://doi.org/10.1117/12.822464.

43.
Hanson KL, Nicolau DV, Filipponi L, Wang L, Lee AP, Nicolau DV. Fungi use efficient algorithms for the exploration of microfluidic networks. Small. 2006;2:1212–20.
CAS  PubMed  Article  PubMed Central  Google Scholar 

44.
Falconer RE, Houston AN, Otten W, Baveye PC. Emergent behavior of soil fungal dynamics: influence of soil architecture and water distribution. Soil Sci. 2012;177:111–9.
CAS  Article  Google Scholar 

45.
Schmidt MWI, Torn MS, Abiven S, Dittmar T, Guggenberger G, Janssens IA. et al. Persistence of soil organic matter as an ecosystem property. Nature. 2011;478:49–56.
CAS  PubMed  Article  PubMed Central  Google Scholar 

46.
Nezhad AS, Packirisamy M, Bhat R, Geitmann A. In vitro study of oscillatory growth dynamics of camellia pollen tubes in microfluidic environment. IEEE Trans Biomed Eng. 2013;60:3185–93.
PubMed  Article  PubMed Central  Google Scholar 

47.
Tayagui A, Sun Y, Collings D, Garrill A, Nock V. An elastomeric micropillar platform for the study of protrusive forces in hyphal invasion. Lab Chip. 2017;17:3643–53.
CAS  PubMed  Article  PubMed Central  Google Scholar 

48.
Brand A, Gow NA. Mechanisms of hypha orientation of fungi. Curr Opin Microbiol. 2009;12:350–7.
CAS  PubMed  PubMed Central  Article  Google Scholar 

49.
Pantidou M, Watling R, Gonou Z. Mycelial characters, anamorphs, and teleomorphs in genera and species of various families of Agaricales in culture. Mycotaxon. 1983;17:409–32.
Google Scholar 

50.
Boekhout T, Stalpers J, Verduin SJW, Rademaker J, Noordeloos ME. Experimental taxonomic studies in Psilocybe sect. Psilocybe. Mycol Res. 2002;106:1251–61.
Article  Google Scholar 

51.
Valenzuela E, Garnica S. Pseudohelicomyces, a new anamorph of Psilocybe. Mycol Res. 2000;104:738–41.
Article  Google Scholar 

52.
Gadd GM, Bahri-Esfahani J, Li Q, Rhee YJ, Wei Z, Fomina M, et al. Oxalate production by fungi: Significance in geomycology, biodeterioration and bioremediation. Fungal Biol Rev. 2014;28:36–55.
Article  Google Scholar  More

SST data
I used the National Oceanic and Atmospheric Administration (NOAA) daily optimum interpolation SST gridded dataset V2.0 to identify MHWs in the period 1 January 1982 to 31 December 20184,28. The dataset is a blend of observations from satellites, ships and buoys and includes bias adjustment of satellite and ship observations to compensate for platform differences and sensor biases. Remotely sensed SSTs were obtained through the Advanced Very High Resolution Radiometer and interpolated daily onto a 0.25° × 0.25° spatial grid globally. Data were provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at https://www.esrl.noaa.gov/psd/ accessed in January 2019.
I also obtained data from simulation models implemented in the fifth phase of the Coupled Model Intercomparison Project (CMIP5). I used the first ensemble member r1i1p1 of 12 coupled Earth System Models that allowed the analysis of variation in daily SSTs in all scenarios: CNMR-CM5, GFDL-CM3, GFDL-ESM2G, GFDLESM2M, IPSL-CM5A-LR, IPSL-CM5A-MR, MPI-ESM-LR, MPI-ESM-MR, MIROC5, MIROC-ESM, MRI-CGCM3, MIROC-ESM-CHEM. A climatology (the statistical properties of the timeseries, including the mean, variance, seasonal cycle and quantiles; see section “Identifying marine heatwaves” below for derivation) was obtained from historical simulations over the period 1861–2005 and used to identify MHWs for the historical scenario and for simulated SSTs over the period 2006–2100 following high and low emission scenarios (RCP 8.5 and RCP 2.6, respectively; RCP, representative concentration pathway). All the analyses on simulated data were implemented on a multi-model ensemble obtained by averaging the twelve models, unless otherwise indicated. For comparisons with the simulated SSTs, the satellite 0.25° × 0.25° data were regridded by averaging daily onto a regular 1° × 1° grid. The climatology for the satellite MHWs was derived from the whole observational period (1982–2018).
Whether a fixed climatology is appropriate instead of using shifting baselines to define MHWs is a matter of debate13,40. Here, the historical scenario provides a common reference to gauge shifts in the spatiotemporal dynamics of projected MHWs under high (RCP 8.5) and low (RCP 2.6) emission scenarios.
Identifying marine heatwaves
I identified MHWs from daily observed and simulated SST timeseries within each 1° × 1° cell following Hobday et al.27, who define a MHW as an anomalously warm water event with daily SSTs exceeding the seasonally varying 90th percentile (climatological threshold) for at least 5 consecutive days. The climatological mean and threshold were computed for each calendar day within a 11-day window centered on the focal day across all years within the climatological period. The mean and threshold were further smoothed by applying a 31-day moving average. Two events with a break of less than 3 days were considered the same MHW. I then derived characteristic metrics of MHWs, including duration, intensity and frequency and linked them to network properties (see below, Network analysis). Only SST timeseries with less than 10% of missing data were used in the analysis. I used the R package heatwaveR to identify marine heatwaves from SSTs41.
Topological data analysis and the Mapper workflow
Topological Data Analysis (TDA) is a collection of statistical methods based on topology, the field of mathematics that deals with the study of shapes, to find structure in complex datasets24. The Mapper algorithm is one tool of TDA that allows reducing high-dimensional data into a combinatorial object that encapsulates the original topological and geometric information of the data, such that points close to each other are more similar than distant points. The combinatorial object, also called a shape graph, is indeed a network with nodes and edges. The statistical properties of the TDA-based Mapper algorithm and how it relates to other non-linear dimensionality reduction techniques have been discussed in Ref.22. Here, I briefly summarize the five key steps involved in a Mapper analysis (Fig. 1). The first step of MAPPER consists of collapsing a raster stack of spatiotemporal data of MHWs into a binary 2D matrix where rows are timeframes (days) and columns are 1° × 1° cells arranged sequentially to represent the occurrence of MHWs across the global ocean. The first column of the matrix corresponds to the upper-left pixel of the raster centered at 89.5°N and − 180°W and the subsequent 364 columns represent adjacent pixels within the same latitude. Column 366 is centered at 88.5°N and − 180°W and so on, with the last column of the matrix corresponding to the lower-right pixel at − 89.5°S and 180°E. Although this scheme would result in matrices with 64,800 columns (360 × 180), I used reduced matrices in computations by excluding pixels on land or where missing SST values prevented the identification of MHWs. The final size of the matrices used in the analysis is (rows × columns) 13,514 × 42,365 for observed SSTs and 52,960 × 41,968, 34,675 × 41,074 and 34,675 × 41,482 for simulated SSTs under the historical, RCP 2.6 and RCP 8.5 scenarios, respectively.
The second step of Mapper involves dimensionality reduction or filtering. I used the Uniform Manifold Approximation and Projection dimensionality reduction (UMAP) algorithm to perform nonlinear dimensionality reduction42. This algorithm is similar to t-distributed Stochastic Neighbor Embedding (tSNE), which is widely used in machine learning43. The advantage of UMAP is that it has superior run time performance compared to tSNE, while retaining the ability to preserve the local structure of the original high-dimensional space after projection into the low-dimensional space.
The third step of Mapper consists of dividing the output range generated by the filtering process into overlapping bins. The number of bins and the amount of overlap are determined by the resolution (R) and gain (G) parameters, respectively. I used an optimization procedure to objectivity identify the combination of parameters R and G that best localized timeframes with similar cumulative intensity of MHWs nearby in the network (see “Parameter search and sensitivity analysis”). This procedure selected R = 24 and G = 45 for observed SSTs and for the RCP 8.5 scenario, R = 22 and G = 45 for the historical scenario and R = 12 and G = 25 for the RCP 2.6 scenario.
The fourth step of Mapper consists of partial clustering of timeframes within bins. Although Mapper is flexible and can accommodate different clustering methods and distance functions, I employed single-linkage clustering with Euclidean distance25. It is worth noting that this approach does not involve averaging of timeframes within clusters, so the original information is preserved in a compressed representation of the data.
The fifth and final step involves the generation of the network graph from the low dimensional compressed representation of the data. Clusters become nodes in the network and nodes become connected if they share one or more timeframes. I implemented the TDA-based Mapper algorithm using a parallelized version of function mapper2D in the R package TDAmapper44.
Network analysis
I employed two widely used measures of network topology, modularity and node degree, to compare the structure of the four MHW networks. Modularity describes the strength of division of a network into communities—i.e. cohesive groups of nodes that have dense connections among them and that are only sparsely connected with nodes in other groups. High modularity indicates the presence of distinct regimes of spatiotemporal dynamics of MHWs. As a second measure of network structure, I used mean node degree—where the degree of a node is the number of edges that are incident to that node. High mean node degree indicates that many nodes share one or more timeframes and depicts similar spatiotemporal patterns of MHWs within those nodes. In contrast, low node degree indicates the occurrence of many isolated nodes with few timeframes in common and more isolated MHWs. I computed modularity and node degree with functions ‘modularity’ and ‘degree’ in the R package igraph45.
To provide significance tests for the observed measures of network topology and to evaluate if they originated simply from non-stationarity properties in the original data, I run two null models based on surrogate data for each of the four networks. Surrogate data can be obtained through the Fourier transform of the original timeseries, shuffling the phases and applying the inverse transform to generate the surrogate series46. Phase randomization preserves the power spectrum, autocorrelation function and other linear properties of the data, but not the amplitude distribution. To address this potential drawback, I generated surrogate timeseries via the Theiler’s Amplitude Adjusted Fourier Transform (AAFT) using function ‘surrogate’ in the R package fractal, which also preserves the amplitude distribution of the original timeseries47. Using this approach, I applied two schemes of randomization—one employing a random sequence for each timeseries and one employing the same sequence for all timeseries. Randomizing using a fixed sequence for all timeseries (constant phase) randomizes the nonlinear properties of the data while preserving linear properties, such as the linear cross-correlation function. The randomization scheme based on random sequences (random phase) also disrupts linear relationships in the data. A significant departure of the observed statistic from the null model under constant phase randomization allows rejecting the null hypothesis that the observed time series is a monotonic nonlinear transformation of a Gaussian process. A significant departure from the null model under random phase allows rejecting also the null hypothesis that the original data come from a linear Gaussian process. To assess significance, 1000 randomizations were performed for each network under each scheme of random and constant phase and a two-tailed test was performed at α = 0.025 to account for multiple testing (Bonferroni correction).
To quantify temporal transitions of MHWs, I estimated node degree of the temporal connectivity matrix (TCM) obtained from each network. The degree for each node in the TCM was estimated by counting the number of non-zero edges connected to that node22. Temporal fluctuations in node degree were benchmarked against the confidence intervals of the random phase null model, estimated as twice the standard deviation of the null distribution. A Generalized Additive Model (GAM) smooth function was fitted to node degree data to visualize temporal trends. The timing of collapse of node degree for the historical scenario was estimated as the year when the smooth curve intersected the upper confidence limit of the null model. To provide a measure of uncertainty, I obtained analogous estimates of the year of collapse by repeating the whole analysis for each of the twelve ESMs separately and computing the median and the bootstrap standard error (n = 1000) of these estimates. A similar analysis was done on the RCP 8.5 data to estimate the year when node degree increased again and diverged significantly from the null distribution. To determine the duration of the different period of connectivity identified in the TCMs, I used the change point algorithm implemented in function cpt.mean of package changepoint48.
Parameter search and sensitivity analysis
To objectively identify parameters R and G (resolution and gain) as part of the binning process in the Mapper algorithm, I used an optimization procedure that best localized timeframes with similar patterns of MHWs nearby in the network. Localization can be done for any of the properties of MHWs. I used cumulative intensity as the localization criterion since it was a good proxy for other properties of MHWs, such as duration (r = 0.88, p  More

1.
Garrison, W. M., Morrison, D. C., Hamilton, J. G., Benson, A. A. & Calvin, M. Reduction of carbon dioxide in aqueous solutions by ionizing radiation. Science 114, 416–418 (1951).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 
2.
Draganić, Z. D., Draganić, I. G. & Borovičanin, M. The radiation chemistry of aqueous solutions of hydrogen cyanide in the megarad dose range. Radiat. Res. 66, 42–53 (1976).
PubMed  Article  ADS  PubMed Central  Google Scholar 

3.
Bar-Nun, A. & Hartman, H. Synthesis of organic compounds from carbon monoxide and water by UV photolysis. Origins Life 9, 93–101 (1978).
CAS  Article  ADS  Google Scholar 

4.
Miller, S. L. & Urey, H. C. Organic compound synthesis on the primitive earth. Science 130, 245–251 (1959).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

5.
Pasek, M. A., Dworkin, J. P. & Lauretta, D. S. A radical pathway for organic phosphorylation during schreibersite corrosion with implications for the origin of life. Geochim. Cosmochim. Acta 71, 1721–1736 (2007).
CAS  Article  ADS  Google Scholar 

6.
Lim, R. W. J. & Fahrenbach, A. C. Radicals in prebiotic chemistry. Pure Appl. Chem. 92, 1971–1986 (2020).
CAS  Article  Google Scholar 

7.
Studer, A. & Curran, D. P. Catalysis of radical reactions: A radical chemistry perspective. Angew. Chem. Int. Ed. 55, 58–102 (2016).
CAS  Article  Google Scholar 

8.
Shock, E. L. et al. Quantifying inorganic sources of geochemical energy in hydrothermal ecosystems, Yellowstone National Park, USA. Geochim. Cosmochim. Acta 74, 4005–4043 (2010).
CAS  Article  ADS  Google Scholar 

9.
Bím, D., Maldonado-Domínguez, M., Rulíšek, L. & Srnec, M. Beyond the classical thermodynamic contributions to hydrogen atom abstraction reactivity. Proc. Natl. Acad. Sci. USA 115, E10287–E10294 (2018).
PubMed  Article  CAS  PubMed Central  Google Scholar 

10.
Mayer, J. M. Hydrogen atom abstraction by metal–oxo complexes: Understanding the analogy with organic radical reactions. Acc. Chem. Res. 31, 441–450 (1998).
CAS  Article  Google Scholar 

11.
Gutowski, M. & Kowalczyk, S. A study of free radical chemistry: Their role and pathophysiological significance. Acta Biochim. Pol. 60, 1–16 (2013).
CAS  PubMed  Article  PubMed Central  Google Scholar 

12.
Moran, J. & Rauscher, S. Energy and self-organization at the origin of metabolism. Commun. Chem. (in rev.).

13.
Nghe, P. et al. Prebiotic network evolution: Six key parameters. Mol. BioSyst. 11, 3206–3217 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

14.
Jolley, C. & Douglas, T. Topological biosignatures: Large-scale structure of chemical networks from biology and astrochemistry. Astrobiology 12, 29–39 (2012).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

15.
Solé, R. V. & Munteanu, A. The large-scale organization of chemical reaction networks in astrophysics. Europhys. Lett. 68, 170–176 (2004).
Article  ADS  CAS  Google Scholar 

16.
Shenhav, B., Solomon, A., Lancet, D. & Kafri, R. in Transactions on Computational Systems Biology I (ed. Priami, C.) 14–27 (Springer, Berlin, 2005).
Google Scholar 

17.
Brown, J. H. et al. The fractal nature of nature: Power laws, ecological complexity and biodiversity. Philos. Trans. R. Soc. Lond. B 357, 619–626 (2002).
Article  Google Scholar 

18.
Walker, S. I. & Mathis, C. in Prebiotic Chemistry and Chemical Evolution of Nucleic Acids (ed. Menor-Salvár, C.) 263–291 (Springer, Berlin, 2018).
Google Scholar 

19.
Hordijk, W., Hein, J. & Steel, M. Autocatalytic sets and the origin of life. Entropy 12, 1733–1742 (2010).
CAS  Article  ADS  Google Scholar 

20.
Albert, R. Scale-free networks in cell biology. J. Cell. Sci. 118, 4947–4957 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

21.
Liu, R., Mao, G. & Zhang, N. Research of chemical elements and chemical bonds from the view of complex network. Found. Chem. 21, 193–206 (2019).
CAS  Article  Google Scholar 

22.
Estrada, E. The complex networks of earth minerals and chemical elements. MATCH Commun. Math. Comput. Chem. 59, 605–624 (2008).
MathSciNet  CAS  MATH  Google Scholar 

23.
Fricker, M. D., Boddy, L., Nakagaki, T. & Bebber, D. P. In Adaptive Biological Networks (eds. Gross, T. & Sayama, H.) 51–70 (Springer, Berlin, 2009).
Google Scholar 

24.
Nicolis, G. Chemical chaos and self-organization. J. Phys. Condens. Matter 2, SA47–SA62 (1990).
CAS  Article  ADS  Google Scholar 

25.
Pérez-Mercader, J. In Astrobiology (eds. Horneck, G. & Baumstark-Khan, C.) 337–360 (Springer, Berlin, 2002).
Google Scholar 

26.
Li, W. Expansion-modification systems: A model for spatial 1/f spectra. Phys. Rev. A 43, 5240–5260 (1991).
MathSciNet  CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

27.
Albert, R. & Barabási, A.-L. Topology of evolving networks: Local events and universality. Phys. Rev. Lett. 85, 5234–5237 (2000).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

28.
Barabási, A.-L. & Albert, R. Emergence of scaling in random networks. Science 286, 509–512 (1999).
MathSciNet  PubMed  MATH  Article  ADS  PubMed Central  Google Scholar 

29.
Marković, D. & Gros, C. Power laws and self-organized criticality in theory and nature. Phys. Rep. 536, 41–74 (2014).
MathSciNet  Article  ADS  Google Scholar 

30.
Adler, R., Feldman, R. & Taqqu, M. (eds.) A Practical Guide to Heavy Tails: Statistical Techniques and Applications (Springer, Berlin, 1998).
Google Scholar 

31.
Patten, B. C. & Higashi, M. Modified cycling index for ecological applications. Ecol. Modell. 25, 69–83 (1984).
Article  Google Scholar 

32.
Essington, T. E. & Carpenter, S. R. Nutrient cycling in lakes and streams: Insights from a comparative analysis. Ecosystems 3, 131–143 (2000).
CAS  Article  Google Scholar 

33.
Christian, R. R. & Thomas, C. R. Network analysis of nitrogen inputs and cycling in the Neuse River estuary, North Carolina, USA. Estuaries 26, 815–828 (2003).
CAS  Article  Google Scholar 

34.
Allesina, S. & Ulanowicz, R. E. Cycling in ecological networks: Finn’s index revisited. Comput. Biol. Chem. 28, 227–233 (2004).
CAS  PubMed  MATH  Article  PubMed Central  Google Scholar 

35.
Loreau, M. Material cycling and the stability of ecosystems. Am. Nat. 143, 508–513 (1994).
Article  Google Scholar 

36.
DeAngelis, D. L. et al. Nutrient dynamics and food-web stability. Annu. Rev. Ecol. Syst. 20, 71–95 (1989).
Article  Google Scholar 

37.
Artzy-Randrup, Y. & Stone, L. Connectivity, cycles, and persistence thresholds in metapopulation networks. PLoS Comput. Biol. 6, e1000876 (2010).
MathSciNet  PubMed  PubMed Central  Article  ADS  CAS  Google Scholar 

38.
Newsholme, E. A. & Crabtree, B. Substrate cycles in metabolic regulation and in heat generation. Biochem. Soc. Symp. 41, 61–109 (1976).
CAS  Google Scholar 

39.
Kritz, M. V., dos Santos, M. T., Urrutia, S. & Schwartz, J.-M. Organising metabolic networks: Cycles in flux distributions. J. Theor. Biol. 265, 250–260 (2010).
MathSciNet  PubMed  MATH  Article  PubMed Central  Google Scholar 

40.
Valentine, J. W. & May, C. L. Hierarchies in biology and paleontology. Paleobiology 22, 23–33 (1996).
Article  Google Scholar 

41.
McShea, D. W. The hierarchical structure of organisms: A scale and documentation of a trend in the maximum. Paleobiology 27, 405–423 (2001).
Article  Google Scholar 

42.
Trebilco, R., Baum, J. K., Salomon, A. K. & Dulvy, N. K. Ecosystem ecology: Size-based constraints on the pyramids of life. Trends Ecol. Evol. 28, 423–431 (2013).
PubMed  Article  PubMed Central  Google Scholar 

43.
Lindeman, R. L. The trophic-dynamic aspect of ecology. Bull. Math. Biol. 53, 167–191 (1991).
Article  Google Scholar 

44.
Kleidon, A. & Lorenz, R. D. (eds.) Non-equilibrium Thermodynamics and the Production of Entropy: Life, Earth, and Beyond (Springer, Berlin, 2005).
Google Scholar 

45.
Goldenfeld, N. & Woese, C. Life is physics: Evolution as a collective phenomenon far from equilibrium. Annu. Rev. Condens. Matter Phys. 2, 375–399 (2011).
CAS  Article  ADS  Google Scholar 

46.
Braakman, R. & Smith, E. The compositional and evolutionary logic of metabolism. Phys. Biol. 10, 011001 (2013).
PubMed  Article  ADS  CAS  PubMed Central  Google Scholar 

47.
Ji, S. Molecular Theory of the Living Cell: Concepts, Molecular Mechanisms, and Biomedical Applications (Springer, Berlin, 2012).
Google Scholar 

48.
Yi, R. et al. A continuous reaction network that produces RNA precursors. Proc. Natl. Acad. Sci. USA 117, 13267–13274 (2020).
CAS  PubMed  Article  PubMed Central  Google Scholar 

49.
Yi, R., Hongo, Y., Yoda, I., Adam, Z. R. & Fahrenbach, A. C. Radiolytic synthesis of cyanogen chloride, cyanamide and simple sugar precursors. ChemistrySelect 3, 10169–10174 (2018).
CAS  Article  Google Scholar 

50.
Ritson, D. & Sutherland, J. D. Prebiotic synthesis of simple sugars by photoredox systems chemistry. Nat. Chem. 4, 895–899 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

51.
Ferus, M. et al. High energy radical chemistry formation of HCN-rich atmospheres on early Earth. Sci. Rep. 7, 6275 (2017).
PubMed  PubMed Central  Article  ADS  CAS  Google Scholar 

52.
Getoff, N. Significance of solvated electrons (eaq−) as promoters of life on Earth. In Vivo 28, 61–66 (2014).
CAS  PubMed  PubMed Central  Google Scholar 

53.
Negrón-Mendoza, A., Draganić, Z. D., Navarro-González, R. & Draganić, I. G. Aldehydes, ketones, and carboxylic acids formed radiolytically in aqueous solutions of cyanides and simple nitriles. Radiat. Res. 95, 248–261 (1983).
Article  ADS  Google Scholar 

54.
Adam, Z. R. et al. Estimating the capacity for production of formamide by radioactive minerals on the prebiotic Earth. Sci. Rep. 8, 265 (2018).
PubMed  PubMed Central  Article  ADS  CAS  Google Scholar 

55.
Bedau, M. A. et al. Open problems in artificial life. Artif. Life 6, 363–376 (2000).
CAS  PubMed  Article  PubMed Central  Google Scholar 

56.
Grassberger, P. in Information Dynamics NATO ASI Series (Series B: Physics) (eds. Atmanspacher, H. & Scheingraber, H.) 15–33 (Springer, Berlin, 1991).
Google Scholar 

57.
Kaneko, K. Chaos as a source of complexity and diversity in evolution. Artif. Life 1, 163–177 (1993).
Article  Google Scholar 

58.
Buhl, D. & Ponnamperuma, C. Interstellar molecules and the origin of life. Sp. Life Sci. 3, 157–164 (1971).
CAS  ADS  Google Scholar 

59.
Airapetian, V. S., Glocer, A., Gronoff, G., Hébrard, E. & Danchi, W. Prebiotic chemistry and atmospheric warming of early Earth by an active young Sun. Nat. Geosci. 9, 452–455 (2016).
CAS  Article  ADS  Google Scholar 

60.
Paranicas, C., Cooper, J. F., Garrett, H. B., Johnson, R. E. & Sturner, S. J. in Europa (eds. Pappalardo, R. T. et al.) 529–544 (University of Arizona Press, Tucson, 2009).
Google Scholar 

61.
Takano, Y., Masuda, H., Kaneko, T. & Kobayashi, K. Formation of amino acids from possible interstellar media by γ-rays and UV irradiation. Chem. Lett. 31, 986–987 (2002).
Article  Google Scholar 

62.
Powner, M. W., Gerland, B. & Sutherland, J. D. Synthesis of activated pyrimidine ribonucleotides in prebiotically plausible conditions. Nature 459, 239–242 (2009).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

63.
Clauset, A., Shalizi, C. R. & Newman, M. E. J. Power-law distributions in empirical data. SIAM Rev. 51, 661–703 (2009).
MathSciNet  MATH  Article  ADS  Google Scholar 

64.
Grohe, M. in Proceedings of the 39th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems 1–16 (Portland, OR, USA, 2020).

65.
Grover, A. & Leskovec, J. in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining 855–864 (San Francisco, CA, USA, 2016).

66.
Palumbo, E. et al. in The Semantic Web: European Semantic Web Conference Vol. 11155, 117–120 (Springer, Crete, Greece, 2018).

67.
Kim, M., Baek, S. H. & Song, M. Relation extraction for biological pathway construction using node2vec. BMC Bioinform. 19, 206 (2018).
Article  CAS  Google Scholar 

68.
Shen, Z., Chen, F., Yang, L. & Wu, J. Node2vec representation for clustering journals and as a possible measure of diversity. J. Data Inf. Sci. 4, 79–92 (2019).
Google Scholar 

69.
Barabási, A.-L. & Oltvai, Z. N. Network biology: Understanding the cell’s functional organization. Nat. Rev. Genet. 5, 101–113 (2004).
PubMed  Article  CAS  PubMed Central  Google Scholar 

70.
Jeong, H., Tombor, B., Albert, R., Oltvai, Z. N. & Barabási, A.-L. The large-scale organization of metabolic networks. Nature 407, 651–654 (2000).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

71.
Ritson, D. J. & Sutherland, J. D. Synthesis of aldehydic ribonucleotide and amino acid precursors by photoredox chemistry. Angew. Chem. Int. Ed. 52, 5845–5847 (2013).
CAS  Article  Google Scholar 

72.
Fahrenbach, A. C. et al. Common and potentially prebiotic origin for precursors of nucleotide synthesis and activation. J. Am. Chem. Soc. 139, 8780–8783 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

73.
Muñoz, M. A. Colloquium: Criticality and dynamical scaling in living systems. Rev. Mod. Phys. 90, 031001 (2018).
MathSciNet  Article  ADS  Google Scholar 

74.
Langton, C. G. Computation at the edge of chaos: Phase transitions and emergent computation. Phys. D Nonlinear Phenom. 42, 12–37 (1990).
MathSciNet  Article  ADS  Google Scholar 

75.
Bak, P., Tang, C. & Wiesenfeld, K. Self-organized criticality: An explanation of the 1/f noise. Phys. Rev. Lett. 59, 381–384 (1987).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

76.
Gaveau, B., Moreau, M. & Toth, J. Scenarios for self-organized criticality in dynamical systems. Open Syst. Inf. Dyn. 7, 297–308 (2000).
MathSciNet  MATH  Article  Google Scholar 

77.
Bak, P. & Paczuski, M. Complexity, contingency, and criticality. Proc. Natl. Acad. Sci. USA 92, 6689–6696 (1995).
CAS  PubMed  Article  ADS  PubMed Central  Google Scholar 

78.
Hoffmann, H. & Payton, D. W. Optimization by self-organized criticality. Sci. Rep. 8, 2358 (2018).
PubMed  PubMed Central  Article  ADS  CAS  Google Scholar 

79.
Lovecchio, E., Allegrini, P., Geneston, E., West, B. J. & Grigolini, P. From self-organized to extended criticality. Front. Physiol. 3, 98 (2012).
PubMed  PubMed Central  Article  Google Scholar 

80.
Lima-Mendez, G. & van Helden, J. The powerful law of the power law and other myths in network biology. Mol. BioSyst. 5, 1482–1493 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

81.
Broido, A. D. & Clauset, A. Scale-free networks are rare. Nat. Commun. 10, 1017 (2019).
Article  ADS  CAS  Google Scholar 

82.
Stumpf, M. P. H. & Porter, M. A. Critical truths about power laws. Science 335, 665–666 (2012).
MathSciNet  CAS  PubMed  MATH  Article  ADS  PubMed Central  Google Scholar 

83.
Mitzenmacher, M. A brief history of generative models for power law and lognormal distributions. Internet Math. 1, 226–251 (2003).
MathSciNet  MATH  Article  Google Scholar 

84.
Glassman, I., Yetter, R. A. & Glumac, N. G. Combustion 41–69 (Elsevier, New York, 2015).
Google Scholar 

85.
Gleiss, P. M., Stadler, P. F., Wagner, A. & Fell, D. A. Relevant cycles in chemical reaction networks. Adv. Complex Syst. 4, 207–226 (2001).
MathSciNet  MATH  Article  Google Scholar 

86.
Dančík, V., Basu, A. & Clemons, P. in Systems Biology (eds. Prokop, A. & Csukas, B.) 129–178 (Springer, Berlin, 2013).
Google Scholar 

87.
Patten, B. C., Higashi, M. & Burns, T. P. Trophic dynamics in ecosystem networks: Significance of cycles and storage. Ecol. Modell. 51, 1–28 (1990).
Article  Google Scholar 

88.
Orgel, L. E. The implausibility of metabolic cycles on the prebiotic Earth. PLoS Biol. 6, e18 (2008).
PubMed  PubMed Central  Article  CAS  Google Scholar 

89.
Monks, P. S. Gas-phase radical chemistry in the troposphere. Chem. Soc. Rev. 34, 376–395 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

90.
Platt, U. et al. in Tropospheric Chemistry: Results of the German Tropospheric Chemistry Programme (eds. Seiler, W. et al.) 359–394 (Springer, Berlin, 2002).
Google Scholar 

91.
Vasas, V., Fernando, C., Santos, M., Kauffman, S. & Szathmáry, E. Evolution before genes. Biol. Direct 7, 1 (2012).
PubMed  PubMed Central  Article  Google Scholar 

92.
Robertson, M. P. & Joyce, G. F. The origins of the RNA world. Cold Spring Harb. Perspect. Biol. 4, a003608 (2012).
PubMed  PubMed Central  Article  CAS  Google Scholar 

93.
Damer, B. & Deamer, D. The hot spring hypothesis for an origin of life. Astrobiology 20, 429–452 (2020).
PubMed  PubMed Central  Article  ADS  Google Scholar 

94.
Martin, W., Baross, J., Kelley, D. & Russell, M. J. Hydrothermal vents and the origin of life. Nat. Rev. Microbiol. 6, 805–814 (2008).
CAS  PubMed  Article  PubMed Central  Google Scholar 

95.
Soloveichik, D., Cook, M., Winfree, E. & Bruck, J. Computation with finite stochastic chemical reaction networks. Nat. Comput. 7, 615–633 (2008).
MathSciNet  MATH  Article  Google Scholar 

96.
Bastian, M., Heymann, S. & Jacomy, M. in Proceedings of the Third International AAAI Conference on Weblogs and Social Media 361–362 (2009).

97.
Alstott, J., Bullmore, E. & Plenz, D. powerlaw: A Python package for analysis of heavy-tailed distributions. PLoS One 9, e85777 (2014).
PubMed  PubMed Central  Article  ADS  CAS  Google Scholar  More

Affiliations

Department of Biology, McGill University, Montreal, Quebec, Canada
Brian Leung & Anna L. Hargreaves

Bieler School of Environment, McGill University, Montreal, Quebec, Canada
Brian Leung

Department of Biological Sciences, Simon Fraser University, Burnaby, British Columbia, Canada
Dan A. Greenberg

School of Biology and Ecology, University of Maine, Orono, ME, USA
Brian McGill

Mitchell Center for Sustainability Solutions, University of Maine, Orono, ME, USA
Brian McGill

Centre for Biological Diversity, University of St Andrews, St Andrews, UK
Maria Dornelas

Indicators and Assessments Unit, Institute of Zoology, Zoological Society of London, London, UK
Robin Freeman

Authors
Brian Leung

Anna L. Hargreaves

Dan A. Greenberg

Brian McGill

Maria Dornelas

Robin Freeman

Corresponding author
Correspondence to Brian Leung. More

1.
Huwyler, F., Kappeli, J., Serafimova, K., Swanson, E. & Tobin, J. Conservation Finance: Moving Beyond Donor Funding Toward an Investor-driven Approach (WWF, Credit Suisse and McKinsey & Company, 2014); http://go.nature.com/2Ka5Y2u
2.
Deutz, A. et al. Financing Nature: Closing the Global Biodiversity Financing Gap: Full Report (Paulson Institute, Nature Conservancy and Cornell Atkinson Center for Sustainability, 2020).

3.
Halpern, B. et al. Gaps and mismatches between global conservation priorities and spending. Conserv. Biol. 20, 56–64 (2006).
Article  Google Scholar 

4.
James, A., Gaston, K. J. & BalmfordA. Can we afford to conserve biodiversity? BioScience 51, 43–52 (2001).
Article  Google Scholar 

5.
McCarthy, D. et al. Financial costs of meeting global biodiversity conservation targets: current spending and unmet needs. Science 338, 946–949 (2012).
CAS  Article  Google Scholar 

6.
Nature’s Dangerous Decline ‘Unprecedented’; Species Extinction Rates ‘Accelerating’ (IPBES, 2019); http://go.nature.com/2V4ZBN9

7.
The Global Risks Report 2020 (WEF, 2020); https://go.nature.com/3ahNfg8

8.
IUCN Views on the Preparation, Scope and Content of the Post-2020 Global Biodiversity Framework (IUCN, 2018); https://go.nature.com/2WlW3ti

9.
Biodiversity: Finance and the Economic and Business Case for Action (OECD, 2019); https://go.nature.com/3h0F9Kc

10.
Parker, C. & Cranford, M. The Little Biodiversity Finance Book. A Guide to Proactive Investment in Natural Capital (Global Canopy Program, 2010); https://go.nature.com/3mwyxUJ

11.
Coad, L. et al. Widespread shortfalls in protected area resourcing undermine efforts to conserve biodiversity. Front. Ecol. Environ. 17, 259–264 (2019).
Article  Google Scholar 

12.
Kearney, S. G. et al. Estimating the benefit of well-managed protected areas for threatened species conservation. ORYX 54, 276–284 (2020).
Article  Google Scholar 

13.
Waldron, A. et al. Protecting 30% of the Planet for Nature: Costs, Benefits and Economic Implications (IIASA, 2020); https://go.nature.com/387GkDq

14.
Stepping, K. M. K. & Meijer, K. S. The challenges of assessing the effectiveness of biodiversity-related development aid. Trop. Conserv. Sci. https://doi.org/10.1177/1940082918770995 (2018).

15.
Waldron, A. et al. Targeting global conservation funding to limit immediate biodiversity declines. Proc. Natl Acad. Sci. USA 110, 12144–12148 (2018).
Article  Google Scholar 

16.
Gallo-Cajiao, E. et al. Crowdfunding biodiversity conservation. Conserv. Biol. 32, 1426–1435 (2018).
Article  Google Scholar 

17.
Parker, C., Cranford, M., Oakes, N. & Leggett, M. The Little Biodiversity Finance Book 3rd edn (Global Canopy Programme, 2012).

18.
Arlaud, M. et al. in Towards a Sustainable Bioeconomy: Principles, Challenges and Perspectives (eds Filho, W. L. et al.) Ch. 5 (Springer, 2018); https://doi.org/10.1007/978-3-319-73028-8_5

19.
Rawat, U. S. & Agarwal, N. K. Biodiversity: concept, threats and conservation. Environ. Conserv. J. 16, 19–28 (2015).
Article  Google Scholar 

20.
Gorobets, A. Wild fauna conservation: IUCN-CITES match is required. Ecol. Indic. 112, 106091 (2020).
Article  Google Scholar 

21.
Rodrigues, A. S. L. et al. The value of the IUCN Red List for conservation. Trends Ecol. Evol. 21, 71–76 (2006).
Article  Google Scholar 

22.
Rao, M., Naro-Maciel, E. & Sterling, E. Protected Areas and Biodiversity Conservation II: Management and Effectiveness (Network of Conservation Educators and Practitioners, 2009).

23.
Adams, V. M., Iacona, G. D. & Possingham, H. P. Weighing the benefits of expanding protected areas versus managing existing ones. Nat. Sustain. 2, 404–411 (2019).
Article  Google Scholar 

24.
BIOFIN The Biodiversity Finance Initiative Workbook 2018 (United Nations Development Programme, 2018).

25.
Costanza, R. et al. The value of the world’s ecosystem services and natural capital. Nature 387, 253–260 (1997).
CAS  Article  Google Scholar 

26.
Costanza, R. et al. Changes in the global value of ecosystem services. Glob. Environ. Change 26, 152–158 (2014).
Article  Google Scholar 

27.
Naidoo, R. et al. Global mapping of ecosystem services and conservation priorities. Proc. Natl Acad. Sci. USA 105, 9495–9500 (2008).
CAS  Article  Google Scholar 

28.
Turner, W. et al. Global conservation of biodiversity and ecosystem services. BioScience 57, 868–873 (2007).
Article  Google Scholar 

29.
Balmford, A. et al. Economic reasons for conserving wild nature. Science 297, 950–953 (2002).
CAS  Article  Google Scholar 

30.
Hily, E. et al. Assessing the cost-effectiveness of a biodiversity conservation policy: a bio-econometric analysis of Natura 2000 contracts in forests. Ecol. Econ. 119, 197-208 (2015).

31.
Ferraro, P. J., McIntosh, C. & Ospina, M. The effectiveness of the US endangered special act: an econometric analysis using matching methods. J. Environ. Econ. Manag. 54, 245–261 (2007).
Article  Google Scholar 

32.
Waldron, A. et al. Targeting global conservation funding to limit immediate biodiversity declines. Proc. Natl Acad. Sci. USA 110, 12144–12148 (2013).
CAS  Article  Google Scholar 

33.
Waldron, A. et al. Reductions in global biodiversity loss predicted from conservation spending. Nature 551, 364–367 (2017).
CAS  Article  Google Scholar 

34.
Richerzhagen, C. et al. Why We Need More and Better Biodiversity Aid Briefing Paper 13 (German Development Institute, 2016); https://go.nature.com/2K0S9Dz

35.
Myers, N., Mittermeier, R. A., Mittermeier, C. G., Da Fonseca, G. A. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).

36.
Karousakis, K. Evaluating the Effectiveness of Policy Instruments for Biodiversity: Impact Evaluation, Cost-effectiveness Analysis and Other Approaches Environment Working Paper No.141 (OECD, 2018).

37.
Isaza, C., Bofill, W. & Cabrera, H. Cost-effective species conservation: an application to Huemul (Hippocamelus bisulcus) in Chile. Environ. Dev. Econ. 12, 535–551 (2007).
Article  Google Scholar 

38.
Alix-Garcia, J. M., Shapiro, E. N. & Sims, K. R. Forest conservation and slippage: evidence from Mexico’s national payments for ecosystem services program. Land Econ. 88, 613–638 (2012).
Article  Google Scholar 

39.
Bare, M. Assessing the impact of international conservation aid on deforestation in sub-Saharan Africa. Environ. Res. Lett. 10, 125010 (2015).
Article  Google Scholar 

40.
Ferraro, P. J. et al. More strictly protected areas are not necessarily more protective: evidence from Bolivia, Costa Rica, Indonesia, and Thailand. Environ. Res. Lett. 8, 025011 (2013).
Article  Google Scholar 

41.
Lindsey, P. A. et al. More than $1 billion needed annually to secure Africa’s protected areas with lions. Proc. Natl Acad. Sci. USA 115, E10788–E10796 (2018).
CAS  Article  Google Scholar 

42.
Bonham, C. et al. Conservation trust funds, protected area management effectiveness and conservation outcomes: lessons from the global conservation fund. Parks 20, 89–100 (2014).
Article  Google Scholar 

43.
Hein, Lars et al. Progress in natural capital accounting for ecosystems. Science 367, 514–515 (2020).
CAS  Article  Google Scholar 

44.
Natural Capital Accounting and Valuing Ecosystem Services Project (UN, 2019); http://go.nature.com/2K2jsxn

45.
Ecosystem Valuation and Natural Capital Accounting (Gaborone Declaration for Sustainability in Africa, 2012); http://www.gaboronedeclaration.com/nca

46.
Climate Public Expenditure and Institutional Review (CPEIR) (UNDP, 2015); https://go.nature.com/2K0C7tp

47.
BIOFIN Workbook: Mobilising Resources for Biodiversity and Sustainable Development (UND, 2016); https://go.nature.com/3p1PDMb

48.
Shieh, G. Effect size, statistical power, and sample size for assessing interactions between categorical and continuous variables. Br. J. Math. Stat. Psychol. 72, 136–154 (2019).
Article  Google Scholar 

49.
Leon, A. C. & Heo, M. Sample sizes required to detect interactions between two binary fixed-effects in a mixed-effects linear regression model. Comput. Stat. Data Anal. 53, 603–608 (2009).
Article  Google Scholar 

50.
Marques, A. et al. Increasing impacts of land use on biodiversity and carbon sequestration driven by population and economic growth. Nat. Ecol. Evol. 3, 628–637 (2019).
Article  Google Scholar 

51.
Tilman, D. et al. Future threats to biodiversity and pathways to their prevention. Nature 546, 73–81 (2017).
CAS  Article  Google Scholar 

52.
Luther, D. A. et al. Determinants of bird conservation—action implementation and associated population trends of threatened species. Conserv. Biol. 30, 1338–1346 (2016).
Article  Google Scholar 

53.
Hoffmann, M. et al. The impact of conservation on the status of the world’s vertebrates. Science 330, 1503–1509 (2010).
CAS  Article  Google Scholar 

54.
Brooks, T. M. et al. Analysing biodiversity and conservation knowledge products to support regional environmental assessments. Sci. Data 3, I60007 (2016).
Article  Google Scholar 

55.
Keith, D. A. et al. Scientific foundations for an IUCN Red List of ecosystems. PLoS ONE 8, e62111 (2013).
CAS  Article  Google Scholar 

56.
Kaufmann, D., Kraay, A. & Mastruzzi, M. The worldwide governance indicators: methodology and analytical issues. Hague J. Rule Law 3, 220–246 (2011).
Article  Google Scholar 

57.
Akaike, H. Information Theory and an Extension of the Maximum Likelihood Principle (Academiai Kiado, 1973).

58.
Bozdogan, H. Model selection and Akaike’s Information Criterion (AIC): the general theory and its analytical extensions. Psychometrika 52, 345–370 (1987).
Article  Google Scholar 

59.
Angrist, J. D. & Pischke, J.-S. Mostly Harmless Econometrics: An Empiricist’s Companion (Princeton Univ. Press, 2009); http://go.nature.com/3r5t6zA

60.
Wooldridge, J. M. Econometric Analysis of Cross Section and Panel Data 2nd edn (MIT Press, 2010). More