Unbinned contigs expand known diversity in the global microbiome

Microbial life on Earth has deep evolutionary roots and is ubiquitous and abundant in extant ecosystems: prokaryotes (that is, Bacteria and Archaea) emerged more than 4 billion years ago¹, have been detected in every environment on the planet and are estimated to account for ~10³⁰ cells and 10¹⁶–10¹⁸ g of biomass^2,3. This long history of evolution across a comprehensive ecological range is reflected in a great accumulated phylogenetic diversity, vastly exceeding that of multicellular life⁴. Yet, it remains unclear just how diverse extant Bacteria and Archaea really are, with extrapolations from abundance distributions of individual samples predicting anywhere between 10⁶ and 10¹² prokaryotic species worldwide^5,6,7,8.

Only a fraction of this diversity has been accounted for in existing data. Using 16S ribosomal RNA (rRNA) amplicon sequences, previous studies estimated the number of species-level operational taxonomic units (OTUs) in publicly available datasets to be 35,500 in 2004⁹, 210,000 in 2014⁸, 109,000 in 2016¹⁰ and 740,000 in 2019⁷. Based on rarefaction curves, tracking the number of newly discovered types (species) as more data (samples) are added to the survey, these analyses concluded that, at least in some habitats, the rate of species discovery was slowing down.

Although 16S rRNA data continue to be more abundant in public repositories¹¹, integrated surveys of prokaryotic diversity increasingly rely on isolate genomes and metagenome-assembled genomes (MAGs) to overcome some of the limitations of amplicon-based datasets (Supplementary Discussion). In metagenomic data, clade-level diversity is often defined based on sequence similarity in a subset of near-universal taxonomic marker genes^12,13,14 or entire genomes¹⁵. For example, proGenomes3 (ref. ¹⁶) encompasses 41,171 species-level clusters of high-quality isolate genomes from the National Center for Biotechnology Information’s RefSeq and GenBank databases, whereas the Genome Taxonomy Database (GTDB) r226 delineates 143,614 species, combining isolate genomes and high-quality MAGs¹⁷. In addition, large-scale MAG catalogues continue to add to the census, reporting for example 4,644 species-level genomes in the human gut microbiome¹⁸, 8,300 species in the ocean¹⁹, 21,077 species in soils²⁰, between 100 and 13,000 species for the various habitats in the MGnify Genomes collection²¹ or 18,000 species across the habitats represented in the IMG/M database²². We recently developed SPIRE²³, which provides an integrated survey of microbial diversity across 99,146 manually annotated metagenomes from all sampled habitats on Earth, describing 107,068 species-level clusters, of which 92,134 are entirely based on MAGs (that is, lacking a cultivated representative).

Such genome-centric surveys have greatly expanded the scope of study for microbial diversity, ecology and evolution, in particular for uncultivated taxa: in GTDB r226, just 64 out of 189 phyla and 34,500 out of 143,614 species contain cultivated representatives. In other words, 66% of currently recognized phyla and 76% of species have been described solely based on MAGs. Nevertheless, MAG recovery remains limited in terms of both precision (due to residual binning artefacts²⁴) and recall (most bins remain incomplete with current automated binning tools²⁵). Indeed, with common workflows, most metagenomically assembled contigs remain unbinned: in SPIRE, just 2.36 terabase pairs (or 10%) of the total assembly are assigned to medium- or high-quality MAGs. These unbinned contigs circumscribe a space of discoverable microbial diversity (that is, the reservoir of recoverable genomes that remains untapped by current workflows). Recent studies estimated this genomically unrepresented diversity at around 135,000 species in 18,000 metagenomes²⁶ and 83,000 species in 249,000 metagenomes²⁷. Novel tools have emerged to make these genomes accessible: the recently developed Bin Chicken workflow of targeted co-assembly recovered genomes for 24,028 previously unrepresented species from just 800 sample groups²⁸.

In compiling this Resource, we conducted a survey of discoverable Bacterial and Archaeal diversity among a set of 92,187 metagenomes, sampled across Earth’s microbial habitats. We tracked the diversity represented in 120 Bacterial¹³ and 53 Archaeal¹⁴ marker genes originating from reference isolate genomes¹⁶, MAGs^17,23 and unbinned metagenomic contigs to quantify the number of lineages to be incrementally discovered in each set. Using habitat-stratified rarefactions, we tallied the amount of diversity that remains to be discovered in available sequence data and estimated the additional diversity expected to emerge as new metagenomes are added to the survey. We tracked lineages across marker gene trees to quantify how many deeper clades (that is, unrecognized phyla, classes, orders, families and genera) remain undiscovered across different habitats. We show that the size distributions of Bacterial and Archaeal clades (including those inferred from unbinned genes) follow power laws, in line with century-old hypotheses by Willis²⁹ and Yule³⁰, with ramifications for biodiversity theory and the study of microbial evolution. In other words, we ask how many and which microbial lineages are hiding in plain sight among unbinned contigs because they are missed by current toolkits, and how many more lineages we are poised to discover as we continue to sample Earth’s habitats.

Just 20–50% of discoverable species are captured by genomes

We tracked the incremental discovery of Bacterial and Archaeal diversity among 92,187 metagenomic samples in SPIRE²³ along rarefaction curves using species-level clusters of assembled taxonomic marker genes (Fig. 1a). We estimated the total discoverable diversity in the studied dataset to be ~705,000 Bacterial and ~27,000 Archaeal species (Supplementary Table 1 and Fig. 1b); rarefaction trajectories (Extended Data Fig. 1) and overall species count estimates (Extended Data Fig. 2) were remarkably consistent across the different independently queried marker genes. Just ~20,000 Bacterial and ~800 Archaeal species comprised cultivated representatives in proGenomes3; an additional ~25,000 and ~1,700 species (for a total representation of 6.4 and 8.8%, respectively) mapped to GTDB r220; and a further ~80,500 and ~4,200 species were contributed by SPIRE MAGs, for a total species-level diversity representation within genome-based datasets of 17.8% among Bacteria and 24.6% among Archaea. In other words, up to ~75–80% of species-level clusters were not captured by genomes of cultivated organisms or MAGs.

These unbinned species were over-represented among uncultivated phyla (recognized solely based on MAGs in the GTDB, with no cultivated isolates; Cohen’s d = 0.33; Wilcoxon P = 0.01; Extended Data Fig. 3a) and, compared with species with a cultivated representative, were associated with higher genomic GC content (R² = 0.15; P = 1.6 × 10⁻¹⁶), but not average genome size (P = 0.338), coding density (P = 0.418) or estimated clade size (P = 0.889) in a multiple linear regression (Extended Data Fig. 3b). Species-level cluster sizes followed a long-tailed distribution, driven by a large fraction of singleton clusters (containing just one sequence; Extended Data Fig. 4). Although this is consistent with expectations of the so-called rare biosphere of lineages with low prevalence and abundance³¹, a subset of singleton clusters may also arise from assembly artefacts or spurious sequences. When conservatively considering only non-singleton marker gene clusters (containing ≥2 sequences from independent sources), diversity estimates were accordingly reduced to a total of ~249,000 Bacterial and ~12,800 Archaeal species (Fig. 1a,b and Supplementary Table 1); this can be considered a lower bound (Supplementary Discussion), yet still corresponded to an increase of 98% relative to the discoverable species set represented by genomes. Indeed, although around three in four unbinned marker gene clusters were singletons, these accounted for <20% of unbinned genes, with the remainder mapping to reference genomes, MAGs or larger unbinned gene clusters (Extended Data Fig. 4). Taken together, this suggests that the vast majority of unbinned genes represent a real biological signal and probably belong to low-quality (incomplete or contaminated) MAGs, enriched in uncultivated clades.

At the same time, marker genes for around half of the species defined by (isolate) genomes in proGenomes3 and GTDB r220 were not represented in metagenomic assemblies (total reference dataset sizes indicated by dotted lines in Fig. 1a), and around 40% were singletons (Extended Data Fig. 4). This suggests that species with cultivated representatives and those pinpointed for isolation are themselves often rare, occurring at low prevalence and/or low abundance in environmental data. Thus, the discovery gap goes both ways: isolate genomes and MAGs captured only 20–50% of discoverable species in metagenomic contigs, but metagenomic assemblies likewise failed to capture around half of the diversity represented by isolates.

Species discovery remains in full swing across Earth’s habitats

The number of discoverable species varied greatly between individual habitats (Fig. 1b), with moderate correlation to sampling effort (ρ_Spearman = 0.51 for Bacteria and 0.09 for Archaea). Estimated Archaeal species counts were generally two to three orders of magnitude below those for Bacteria in host-associated and anthropogenic habitats (including some environments with almost no detectable Archaea), roughly within 1 log in aquatic habitats and soils and almost on par in extreme environments. Genomic representation of discoverable species likewise varied across habitats, but was generally within one order of magnitude of the total.

We next tested whether species discovery showed any signs of slowing down by calculating species discovery coefficients (termed α) based on individual, habitat-stratified rarefaction curves (Fig. 1c and Methods). We observed that genome-based discovery of Archaeal species was fully saturated (α ≤ 0) or approaching saturation (low positive α) in most host-associated and anthropogenic habitats. Specifically, all Archaeal species represented in proGenomes3 and GTDB r220 that can be metagenomically discovered using current approaches have already been accounted for in habitats such as human infant gut, oral or skin samples; adding more samples to the survey is not expected to increase metagenomic representation. Similarly, anthropogenic environments (built environments and wastewater), wetlands and extreme environments (hot springs and hydrothermal vents) appear saturated in terms of the metagenomic discovery of Archaeal reference species, whereas even MAG-based discovery showed clear signs of slowing down. In contrast, Archaeal species discovery coefficients among unbinned contigs were substantially higher, indicating that the existing large discovery gap between binned and unbinned contigs will continue to widen.

Bacterial species discovery coefficients were generally higher, indicating that, with current approaches, Bacterial diversity will continue to be discovered at higher rates than for Archaea. The human mouth, urogenital tract and adult gut, as well as porcine and bovine gut environments, are the only tested habitats in which the metagenomic discovery of isolate-represented diversity is approaching saturation. Yet, even in these habitats, as across all others, MAG-based species discovery is expected to continue at high rates, outpaced only by species discovery among unbinned contigs (α ≥ 0.5 across most habitats). Soils, wetlands, the rhizosphere, freshwater habitats and the gastrointestinal tracts of non-mammalian, non-avian animals stood out as particular hotspots of untapped diversity, with virtually no slowdown in species discovery (α ≥ 0.8). Overall, the marked difference in species discovery coefficients between reference genome sets, MAGs and unbinned contigs suggests a widening gap between discovered (genomically represented) and discoverable (unbinned) diversity, as the latter will continue to outpace the former.

Beyond species discovery within individual habitats, we explored how much each habitat contributed to the total survey. We tracked habitat-stratified species accumulation curves, quantifying how much diversity each sequentially added habitat contributed incrementally. This revealed striking discrepancies between sampling effort and observed unique diversity (Extended Data Fig. 5). Although human-associated, gastrointestinal and built-environment metagenomes represented almost three-quarters of sampling effort, they only accounted for 38 and 8% of discoverable Bacterial and Archaeal species, respectively. In contrast, soils and aquatic habitats—in particular wetlands—remain major reservoirs of unexplored microbial diversity (Supplementary Discussion).

Uncaptured deeply branching lineages may enrich the Tree of Life

We next explored how these trends translate to deeper taxonomic levels. We inferred gene phylogenies for each of the considered 53 Archaeal and 122 Bacterial markers, estimated the number of represented genus-, family-, order-, class- and phylum-level clades by cutting each tree at level-specific relative evolutionary divergence (RED)³² cutoffs, and summarized the resulting clade counts (Methods). As shown in Fig. 2 for Archaeal marker RNA polymerase subunit E (rpoE), unbinned sequences added considerable tip-level (that is, approximately species-level) diversity for most recognized clades and often also provided relevant ecological context. For example, for the family Nitrosopumilaceae (within the Thermoproteota phylum³³) that encompasses several known ammonia-oxidizing lineages, we observed a substantial expansion of covered diversity among marine genera (for example, Nitrosopelagicus or Nitrosopumilus) while linking other subclades to wetland-, groundwater-, soil- and plant-associated habitats, supported across several phylogenies.

Marker gene phylogenies recovered previously recognized deep-branching clades as mostly homogenous groups and recapitulated common hypotheses regarding their branching order³⁴. Our estimates of reference phylum-level clades (containing sequences from proGenomes3 or GTDB) exceeded the number of phyla recognized in GTDB r220 based on full species trees for both Archaea (estimate: 36 ± 12; GTDB r220: 19) and Bacteria (238 ± 82; GTDB: 175), in part because our RED partitioning algorithm tended to split large clades (such as the TACK superphylum, comprising Thaumarchaeota (Nitrososphaerota), Aigarchaeota, Crenarchaeota and Korarchaeota; represented as the single phylum Thermoproteota in the GTDB) or groups predicted to have skewed evolutionary rates (such as DPANN Archaea; named for the phyla Diapherotrites, Parvarchaeota, Aenigmarchaeota, Nanoarchaeota and Nanohaloarchaeota) into multiple phyla. All explored phylogenies contained varying numbers of deep (phylum-level) branches encompassing only sequences in unbinned contigs or SPIRE MAGs. For example, the Archaeal rpoE tree in Fig. 2 contains an unclassified phylum-level branch comprising 758 sequences from different habitats, organized into 109 tips (gene clusters) and five predicted family-level clades. Overall, we estimated a median of ten putative phylum-level clades in Archaea (4 ± 3 from SPIRE MAGs and 6 ± 8 among unbinned contigs), corresponding to an increase of ~28% over estimated reference clades, and 145 clades in Bacteria (28 ± 29 among MAGs and 117 ± 137 in unbinned contigs), corresponding to an increase over reference estimates of ~61%. Towards shallower taxonomic levels, our estimates of reference clade counts approximated those in GTDB r220 (with the exception of the species level; see above), whereas estimates on discoverable clades among SPIRE MAGs and unbinned contigs strongly increased with increasing taxonomic resolution (Fig. 3a and Table 1).

The recovery of predicted clades was heterogeneous across habitats (Fig. 3b). Gut and non-intestinal human samples covered only a small fraction of predicted clades, at both deep (phylum-level) and shallower (genus-level) resolution, despite comprising the majority of the sequencing effort. The majority of predicted higher-level clades were observable in soil, rhizosphere, wastewater and aquatic environments. Hydrothermal vents stood out as particular hotspots of Archaeal diversity: accounting for just 296 metagenomes in our dataset (0.3% of the total), they contained representatives of four-fifths of reference and two-thirds of MAG-based and unbinned predicted Archaeal phylum-level clades, representing the widest coverage of deep lineages observed for any tested habitat.

Prokaryotic diversity follows Willis’ law and Yule curves

In 1922, Willis observed that taxonomic clade size distributions of several plant and animal groups follow power laws: the frequency of genera decreases with the exponent of genus size²⁹, resulting in hollow curves dominated by singleton clades (containing just one species), with a heavy tail of very large clades (disproportionally species-rich genera). Beginning with Yule’s seminal work explaining this observation by postulating preferential attachment of novel species to existing genera, proportionally to genus size³⁰, researchers have sought to identify the evolutionary mechanisms and processes that give rise to extant biodiversity^{35,36,37,38,39}.

We hypothesized that prokaryotic diversity follows the same empirical laws and reasoned that our dataset was uniquely suited to test this, as our clade-level groups were data driven (inferred based on RED values), taxonomy agnostic (defined from marker phylogenies) and comprehensive (sampled across Earth’s habitats). We found that Archaeal and Bacterial diversity indeed follow Willis’ law: clade frequencies decrease with increasing clade size in a power law relationship (Fig. 4), along the full tested range of more than five orders of magnitude of clade counts. Willis’ law held across taxonomic scales, from species within genera to species within phyla (Fig. 4a), and also for subclades at higher taxonomic levels, from genera within families to classes within phyla (Fig. 4b). Fitting naive power law equations for each individual marker gene, we computed a marker-gene-specific Willis coefficient ω from the slope in a log–log plot (Methods). We observed a clear trend of decreasing ω (implying a heavier tail; that is, stronger bias towards large clades) across taxonomic levels for species counts (from a ω value of ~1.4 for species within genera to a ω value of ~0.1 for species within phyla), and more similar ω values ranging from 1.4–1.7 for subclade-within-clade counts (Extended Data Fig. 6a).

Although naive power laws provided good overall fits, we observed that clade size distributions were disproportionately heavy tailed on the right (large clades were over-represented) while also deviating towards the left (lower frequencies for small clades), thereby showing characteristics originally used by Yule to deduce underlying processes³⁰. We therefore fit the data to Yule–Simon distributions (Fig. 4a,b). Estimates of the Yule–Simon parameter ρ were consistent across marker genes and between Bacteria and Archaea (Extended Data Fig. 6b): ρ decreased with increasing taxonomic depth for species counts (ρ = ~1.0 for genera to ρ = ~0.4 for phyla), with an inverse trend for subclade counts (ρ = ~0.85 for genera within families to ρ = ~1.50 for classes within phyla). ρ can be thought of as a rich-get-richer coefficient that describes preferential attachment in a Yule process: low ρ values indicate a strong preferential attachment of new subclades to large existing clades. Therefore, our results indicate that new species tend to arise in larger existing clades, with a skew that is strongest in phyla and weakest for genera (that is, the dominance of species-rich phyla, such as Bacillota or Pseudomonadota, is more pronounced than that of species-rich genera, such as Prochlorococcus). The inverse was observed for subclades within clades: new genera arise with a stronger bias towards existing large families than new classes towards existing large phyla. For orders within classes and classes within phyla, we observed a disparity between Archaea and Bacteria that was only in part attributable to sampling noise (fewer Archaeal clades leading to noisier fits). This may indicate that an Archaeal order or class is not fully equivalent to Bacterial orders or classes in terms of phylogenetic depth, or that Archaeal and Bacterial diversity at these levels is organized differently.

Diversity patterns vary across habitats and phylogeny

We explored how clade size distributions are shaped by environment and phylogeny by refitting curves for individual habitats (Fig. 4c) and phyla (Fig. 4d). Yule–Simon coefficients were generally lower in well-sampled habitats (ρ_Spearman = −0.50 for Bacteria and −0.25 for Archaea), with clear differences between broader environments (Fig. 4c). Intestinal and human non-intestinal habitats showed the lowest ρ values, arguably because they are narrowly defined and correspond to just one or few host species. In contrast, higher coefficients were observed for aquatic environments (in particular hot springs and groundwater) and soils, but also non-mammalian guts, indicating that diversity is less dominated by few large clades and that newly discovered species in these environments are more likely to belong to undescribed deeper clades, whereas they are expected to belong to (known) large clades in the intestine.

Yule–Simon coefficients were strongly associated with total sampled species richness within recognized phyla for both Bacteria (ρ_Spearman = −0.71) and Archaea (ρ_Spearman = −0.79; Fig. 4d): well-sampled phyla such as Pseudomonadota or Bacillota_A had low ρ coefficients, whereas smaller phyla were less dominated by large clades. Two notable exceptions were Methanobacteriota, which deviated towards lower ρ (fewer species sampled, but highly concentrated), and Patescibacteria, with an estimated ~19,500 species distributed into more evenly sized clades (ρ = 1.33). Yet, nearly all tested prokaryotic phyla had higher Yule–Simon coefficients than eukaryotic kingdoms and even viruses, inferred based on reference taxonomies of the Global Biodiversity Information Facility (GBIF; Extended Data Fig. 7). Moreover, domain-wide ρ estimates based on our dataset were substantially lower than those based on GBIF Bacterial and Archaeal taxonomies, GTDB r226 reference taxonomies⁴⁰ or 16S rRNA OTUs from the Microbe Atlas Project^11,41. This suggests that although eukaryotic diversity (in particular among plants and fungi) is more concentrated into few large clades than for Bacteria and Archaea, the organization of biodiversity follows similar and consistent patterns across the entire Tree of Life.

There are multiple layers of discovery in the description of prokaryotic lineages: the isolation of organisms for experimental access and formal taxonomic approval remains standard^42,43,44,45, yet the delineation of genomes directly from metagenomes has enabled the characterization of a greatly extended range of clades without cultivated representatives, and even variants in (amplified) marker gene sequences can be informative of further diversity. In this context, our study addresses four basic questions: (1) how many lineages can be discovered in currently available data; (2) how much of this diversity is missed by current genome-centric approaches; (3) should we expect discovery to slow down as more data come in; and (4) how is prokaryotic diversity organized across habitats and phylogeny?

Our results extend previous observations^18,20,23,46 that isolates greatly underestimate tractable genomic diversity in metagenomes. Yet, strikingly, this discovery gap goes both ways: only half of the species encompassed by genomes of cultivated representatives (and, by extension, reference MAGs) were detected in our metagenomic dataset. These levels of detection are higher than previous reports²⁶, probably because the present dataset is more comprehensive, although saturating discovery curves for several environments indicate that further detections are not expected when adding more samples (at least when assuming consistent sampling bias). This tallies with the idea that cultivation efforts and metagenomic sampling often target different (micro-)environments and that cultivability does not necessarily reflect abundance in natural communities⁴⁷.

The vast majority of discoverable diversity in our dataset extended beyond all genome-centric approaches. As a lower bound, defined as independent detection in at least two distinct samples, one species remains unbinned for every species discovered among cultivated organisms (proGenomes3), reference MAGs (GTDB r220) or dataset-specific MAGs (SPIRE). When considering all detected sequence variants, estimates go up to four discoverable species per genomically defined species, or 580,000 uncaptured Bacterial species and 21,000 uncaptured Archaeal species, exceeding recent reports of 83,000 species based on raw read mapping to conserved marker gene windows using SingleM and sandpiper²⁷ or 135,000 species based on metagenomic assembly²⁶.

At broader taxonomic levels, our data likewise suggest a wealth of discoverable clades: we estimate that fewer than one in four Bacterial genera and just three in five Bacterial phyla that can be inferred from assemblies are currently recognized in reference databases. These are estimates across many phylogenies so that issues affecting individual gene trees (such as paralogues or uneven evolutionary rates^48,49,50) are expected to mostly even out. Nevertheless, the proper characterization of lineages requires genomic context, or at least phylogenies of multiple concatenated informative markers combined with complex models of evolution^51,52, beyond the scope of the present study. Although these numbers carry increasing uncertainties at broader taxonomic levels, this puts the total (currently discoverable) Bacterial phyla in the range of 300–400, compared with 30–40 among Archaea, although we note that RED-based clade definitions at phylum and class levels were broader in Archaea than in Bacteria. Indeed, this ratio of roughly one order of magnitude difference between Bacterial and Archaeal clade counts held across taxonomic levels, supporting previous observations that Archaea are much less diverse than Bacteria^7,53,54.

We found that Bacterial and Archaeal diversity are organized along common fundamental patterns. Empirical scaling laws are prevalent in microbiology, from biogeographic taxa–area relationships^55,56 to macroecological patterns⁶. Our work provides comprehensive evidence that prokaryotic clade size distributions likewise follow power laws, supporting century-old hypotheses by Willis²⁹ and Yule³⁰. Power laws are a hallmark of fractal geometries, and the remarkable consistency of inferred coefficients—not only between Bacteria and Archaea, but also across eukaryotic kingdoms and viruses, and even for subclades within clades at different depths—implies that consistent (and possibly predictable) evolutionary forces mould the Tree of Life into recurring patterns from root(s) to tips (that is, biodiversity is indeed fundamentally fractal)^35,57,58,59. Our data strongly support Yule’s hypothesis^30,39 that one central pattern is preferential attachment: the rich get richer as new clades arise preferentially within already large clades. The effect varies across habitats and phylogeny and can also be interpreted in light of clade discovery: assuming random sampling, far more species need to be discovered in intestinal environments or species-rich phyla such as Pseudomonadota (with low Yule–Simon coefficients) than in hot springs or among Patescibacteria (higher coefficients) to identify a previously undescribed deeper clade.

We caution against extrapolating total richness estimates from our rarefaction curves, as future data may break out from previous sampling bias (Supplementary Discussion). Indeed, our inferred rarefaction curves and species discovery coefficients are conservative lower-bound estimates, as they assume continued sampling from the same underlying distribution, whereas in fact the vast majority of the planet remains a metagenomic terra incognita, and overlooked endemic clades in unsampled sites may disproportionately contribute to global biodiversity^60,61. Despite these caveats, soils, aquatic habitats and non-mammalian guts stood out as key reservoirs of discoverable diversity and hotspots primed for future discovery, although our findings suggest that the discovery of undescribed diversity will continue unmitigated for almost all considered habitats as data continue to accrue near-exponentially in public repositories. In the meantime, even existing data may still hold countless surprises, as vast prokaryotic diversity remains hidden in plain sight.

Data sources and pre-processing

We used three principal data sources in our analysis: (1) proGenomes3 (ref. ¹⁶), a quality-filtered subset of 907,388 isolate genomes sourced from the National Center for Biotechnology Information’s RefSeq and GenBank databases⁶², organized into 41,171 species-level clusters based on a set of 40 SpecI marker genes¹²; (2) GTDB r220 (April 2024¹⁷), a curated set of 596,859 isolate genomes and MAGs, clustered into 107,235 Bacterial and 5,869 Archaeal species-level clusters sharing ≥95% whole-genome average nucleotide identity; and (3) SPIRE version 1.1 (ref. ²³), encompassing a curated set of 92,187 shotgun metagenomic samples assembled into ~23.3 terabase pairs of contigs, a subset of which were binned into 1,158,468 medium- or high-quality MAGs, organized into 107,078 species-level clusters that partly overlap with those in proGenomes3.

Genomes from all three datasets were downloaded and taxonomically (re-)classified against GTDB r220 using GTDB-Tk version 2.4.0 (ref. ⁶³), and consensus taxonomies for species-level clusters were inferred based on adjusted majority votes, as described previously²³. Compared with SPIRE version 1.0, 6,959 metagenomic samples were excluded for the present analysis based on data type and provenance (for example, excluding samples explicitly enriched for viruses) or insufficient assembly size. For the remaining 92,187 samples, habitat annotations and contextual data were updated by manually curating against an extended microntology version 0.3.0 encompassing 103 terms and categories (ref. ⁶⁴; reference set and updated annotations available via https://spire.embl.de/downloads and Supplementary Table 2) and further organized into 32 higher-level categories (Supplementary Table 1).

Estimates of Bacterial and Archaeal diversity based on taxonomic marker genes

The majority of analyses presented in the main text are based on 168 near-universal taxonomic marker genes, as established by the GTDB: 120 Bacterial (bac120; ref. ¹³) and 53 Archaeal (arc53; ref. ¹⁴) markers, with an overlap of five genes. We downloaded profile hidden Markov models for these marker sets as part of the GTDB-Tk version 2.4.0 database and used the HMMER version 3.4 (ref. ⁶⁵) hmmsearch routine to identify and extract marker gene sequences among predicted open reading frames. For the final filtered set of marker genes (see below), this yielded on average 904,654 ± 98,879 sequences per Bacterial marker in proGenomes3, 93,198 ± 4,733 among GTDB r220 species representatives, 840,023 ± 133,904 in medium- and high-quality SPIRE MAGs and 3,051,465 ± 1,290,024 among unbinned SPIRE contigs (that is, those not assigned to a MAG passing genome quality filters in SPIRE version 1). For Archaea, we obtained 2,451 ± 150 sequences from proGenomes3, 4,898 ± 410 from GTDB r220 species representatives, 15,504 ± 2,593 from SPIRE MAGs and 107,419 ± 75,823 from unbinned SPIRE contigs (Supplementary Table 3 provides a per-gene overview).

For each marker gene, we clustered extracted sequences using the MMseqs2 cascading clustering routine⁶⁶ with a cutoff of 96.5% sequence similarity (see below) and otherwise default parameters. To convert the number of marker gene sequence clusters into a corresponding number of species-level clusters, we estimated conversion factors as follows. First, we generated marker gene cluster discovery curves via iterative logarithmic rarefaction (that is, we downsampled the number of considered gene sequences along a logarithmic scale (10, 20, …, 100 sequences; 200, 300, …, 1,000 sequences and so on) with ten iterations at each step). At each rarefaction point, we recorded the number of discovered marker gene sequence clusters and the number of represented species or species-level genome clusters in proGenomes3, GTDB r220 and among SPIRE MAGs, considering each data source individually. We then used linear regression models of the type number_of_species ~ number_of_gene_clusters with a forced intercept at 0 along these rarefactions to estimate gene-cluster-to-species conversion factors (that is, the number of newly discovered species per newly discovered marker gene cluster). Extended Data Fig. 8 shows the resulting fits for ten randomly selected marker genes, each for Archaea and Bacteria. Based on benchmarks of marker gene sequence similarity cutoffs ranging from 95.0 to 99.5%, we based further analyses on 96.5% clusters as these showed very robust linear fits (with standard errors in the range of 10⁻³) with conversion factors closest to identity (that is, roughly one species discovered per marker gene cluster discovered), and high consistency across the different species-level reference clusterings in the underlying datasets (based on 40 SpecI marker genes in proGenomes3, 95% whole-genome average nucleotide identity in GTDB and a combination of both approaches among SPIRE MAGs). Finally, we filtered the set of considered marker genes by excluding redundant markers (part of both the arc53 and bac120 sets and therefore prone to cross-mappings) and removing markers for which predicted species numbers deviated by more than 20% from the domain-level median (in particular for some Archaeal markers, detected sequence counts were inflated, possibly due to the detection of Bacterial paralogues). Sequence counts per data source for the remaining set of 28 Archaeal and 102 Bacterial marker genes are available in Supplementary Table 3.

Inference of global and habitat-resolved species discovery curves

For each of 32 broadly defined habitat categories (ranging from 267–19,659 samples per group; Supplementary Table 1), as well as for the combined set of all 92,187 metagenomes in SPIRE version 1.1, we generated Bacterial and Archaeal species discovery (or rarefaction) curves. We iteratively subsampled the number of considered metagenomes along a logarithmic rarefaction scale, with five random permutations per step. At each rarefaction point and for each considered marker gene, we inferred the number of discovered species based on the number of discovered marker gene clusters, using gene-specific conversion factors, as described above. Rarefaction permutations were averaged per marker gene and then summarized per marker gene set (bac120 for Bacteria and arc53 for Archaea) as median predicted species counts at each rarefaction step across marker genes (Fig. 1). Extended Data Fig. 1 shows relative errors (coefficient of variation; that is, the standard deviation divided by the mean) at each rarefaction step for each marker gene; consistently, the coefficient of variation decreased to <5% when including 1,000–10,000 samples for Bacteria (that is, when using ~1–10% of the total data), with drops at slightly higher rarefaction steps for Archaea.

In an independent approach, we performed an incremental rarefaction, stratified by habitats: we sequentially added samples from each of the 32 broadly defined habitat categories, starting with adult, healthy human gut samples, then randomized the discovery order within each habitat block and tracked the globally discovered species (Extended Data Fig. 5 and Supplementary Figs. 1 and 2). In other words, we tracked the contribution of each habitat to overall discovered species diversity, beyond the diversity discovered in previously considered habitats.

Counts of discovered marker gene clusters (and, by proxy, inferred discovered species) were hierarchically stratified by data source to account for overlap between datasets as follows: all clusters containing at least one sequence originating from a genome in proGenomes3 were labelled as proGenomes3, irrespective of the origin of other sequences within the same cluster (data series marked in dark blue in main Fig. 1 and Extended Data Fig. 5); clusters containing sequences from GTDB r220 (and any other data source except proGenomes3) were labelled as GTDB (light blue in Fig. 1 and Extended Data Fig. 5); clusters containing sequences from SPIRE MAGs (but neither from proGenomes3 nor GTDB) were labelled as SPIRE MAG (orange data series); clusters containing only sequences from unbinned SPIRE contigs were labelled as unbinned, non-singleton if they contained at least two sequences from different contigs (dark yellow data series) and unbinned, all otherwise (light yellow data series). Thus, clusters labelled as GTDB represent sequence diversity contained in GTDB r220 beyond that already represented in proGenomes3, and a substantial subset of proGenomes3-labelled clusters also contained sequences originating from GTDB, SPIRE MAGs or unbinned contigs. Similarly, SPIRE MAG clusters corresponded to diversity not already covered by proGenomes3 and GTDB r220, and unbinned clusters encompassed diversity not represented in any genome or MAG in the dataset. Moreover, by design we only tracked marker gene clusters in metagenomic samples, meaning that a non-negligible subset of genome-based clusters from proGenomes3 and GTDB r220 were not represented in rarefaction runs and calculations, as corresponding sequences were not detected in any of the 92,187 considered metagenomic assemblies (main text).

Calculation of species discovery coefficients

We quantified the rate at which species-level clusters are discovered in each habitat (and globally, across all habitats) using the equation

where S is the number of newly discovered species per N samples added to the survey; k is a proportionality constant; and γ is a saturation coefficient. An analogous formula is commonly used to describe a (microbial) pangenome’s openness^67,68 (that is, the degree to which genes are discovered as more genomes from the same species are considered). Moreover, the approach is conceptually and mathematically related to the Heridan–Heaps law in linguistics, which describes the number of distinct words in a document as a function of that document’s length⁶⁹.

We fit equation (1) to each habitat-specific and global rarefaction curve described above, stratified by data source, resulting in estimates for γ for each considered marker gene; Extended Data Fig. 9 shows the example fits for a randomly selected subset of marker genes. For more intuitive interpretability, we calculated a species discovery coefficient α as

and summarized values across marker genes within each domain-specific set as median values. Thus defined, α scales on [-∞, 1], although only mildly negative values are expected to be observed in practice. If α ≤ 0, species discovery in a habitat is saturated, meaning that adding more samples of the same type is not expected to add species to the survey, analogous to a closed microbial pangenome where additional strains do not add genes. Values for α in [0, 1] correspond to unsaturated species discovery curves where additional samples continue to add species to the survey (analogous to open pangenomes). Lower α values indicate a more pronounced flattening off in the rarefaction curve, indicating a more pronounced slowdown in species discovery, and higher α values indicate a less pronounced decrease in the rate of species discovery. For α → 1, species discovery is fully unsaturated, meaning that each newly added sample adds species to the survey, with no discernible flattening off in the species discovery curve.

Phylogenetic inference and analyses

Phylogenetic trees for each marker gene were constructed as follows. We first realigned cluster representative amino acid sequences against the respective marker gene hidden Markov models using the hmmalign routine from HMMER version 3.4 (ref. ⁶⁵), then we trimmed the resulting pseudo-multiple sequence alignments to informative columns using clipkit version 1.4.1 (ref. ⁷⁰) in kpi-gappy mode and removed sequences with >70% gaps among the remaining columns. Based on the resulting alignments, we inferred phylogenetic trees using FastTree 2 version 2.1.11 (ref. ⁷¹) under the Whelan–Goldman model⁷².

We calculated RED values³² for each node in the resulting trees using castor package version 1.8.3 in R⁷³. To infer marker-gene-specific RED cutoffs at the phylum, class, order, family and genus levels, we subset phylogenies to sequences originating from fully taxonomically classified genomes in proGenomes3, GTDB r220 or among SPIRE MAGs, cut the trees at incremental RED values (with a RED tolerance of ±0.1) and compared the resulting partitions (that is, tip sets) versus assigned taxonomic labels, quantifying partition consistency as adjusted mutual information⁷⁴. These steps were performed iteratively on trees re-rooted using each recognized phylum in turn as an outgroup and then summarized across alternative root positions (ignoring the current outgroup) to obtain average RED values per internal node, in a workflow adapted from ref. ³².

After confirming that the resulting RED cutoffs at peak adjusted mutual information (that is, best capturing the distribution of established taxonomic labels) well approximated those reported in ref. ³² across individual Archaeal marker genes, we used GTDB’s reference cutoffs (with a ±0.1 tolerance interval) for Bacterial marker gene trees and cut the full trees at the respective levels. The number of resulting clusters (that is, nodes where the tree was cut) and their composition (that is, each node’s descendant tips) were then used as estimates of clade-level groups at each taxonomic level, for each marker gene. We manually selected genes that showed the most consistent and robust profiles across iterations (and rootings) and whose phylogenies were not ostensibly impacted by paralogues, bringing the final sets of considered genes to 29 Archaeal and 53 Bacterial markers. Marker gene phylogenies were visualized using the ggtree and ggtreeExtra packages in R⁷⁵.

Marker gene trees were also used for the analyses underlying cumulative marker gene distributions, with the results presented in Extended Data Fig. 4. Internal nodes corresponding to phyla recognized in the GTDB were identified as the most recent common ancestors of tips representing marker gene clusters containing genes from reference genomes or MAGs classified to the focal phylum. Phylum-level subtrees cut at the respective most recent common ancestor nodes were then profiled for the ratio of tips representing only unbinned genes to tips representing genomes (that is, binned genes), as well as for average genomic features (GC content, genome size and coding density) and the estimated total number of represented species. To alleviate the effects of genome misclassification and errors in phylogenetic inference, analyses were repeated five times for each phylum, each time using a random subset of 80% of genomes to identify the most recent common ancestor node.

Fitting taxonomic clade size distributions

Four datasets were used to fit clade sizes to previously suggested distributions: (1) the above-defined clade-level groups (based on RED cutoffs), for each Archaeal and Bacterial marker gene; (2) GTDB r226 Archaeal and Bacterial reference taxonomies (accessed via https://gtdb.ecogenomic.org/); (3) the GBIF backbone taxonomy for eukaryotes, Bacteria, Archaea and viruses⁷⁶, accessed using the taxizedb R package⁷⁷; and (4) the Microbe Atlas Project¹¹ 16S rRNA reference database version 3.0 (ref. ⁴¹), using 97% OTUs to represent species and 90% OTUs to represent genus. For each of these taxonomies, we first fit Willis curves²⁹ relating the frequency of clades to their size (that is, the number of subclades they contained) as a regular power law:

where C is the number of higher clades (for example, genus) containing S subclades (for example, species); a is a proportionality constant; and ω is the scaling coefficient, referred to as the Willis coefficient here for simplicity. In other words, equation (3) relates the frequency of a clade to its size, with positive ω values indicating that large clades (for example, genera containing many species) are exponentially less common than small clades (for example, genera containing few or only one species). Otherwise, given the similarity of equation (3) to equation (1) above, the interpretation of ω is closely related to that of the above-defined species discovery coefficient α.

We moreover fit clade sizes to the distribution first suggested by Yule³⁰ specifically for this purpose and later formalized by Simon⁷⁸:

where f(C_S) is the non-negative frequency of clades of size S (that is, containing S subclades); Β is the beta function; and ρ > 0 is the Yule–Simon shape parameter. The distribution results from a Yule process where newly arising subclades preferentially attach to larger existing clades, proportionally to the existing clades’ sizes, sometimes colloquially referred to as a rich-get-richer process. For sufficiently large counts C_S, the frequency f(C_S) approximates a power law as in equation (3), with coefficient ω at ∼ρ + 1. Parameter ρ indicates the shape of the distribution and the strength of preferential attachment, with a large ρ value indicating a sharper drop in the distribution (less pronounced preferential attachment, a less tail-heavy distribution and higher dominance by small clades), whereas a smaller ρ value indicates more tail-heavy distributions and stronger preferential attachment effects (large clades are larger; small clades are fewer). In other words, ρ can be thought of as a rich-get-richer coefficient, with low ρ values indicating a stronger dominance of fewer large clades in the clade size distribution.

Habitat-stratified analyses (Fig. 4c) were conducted by recomputing clade size distributions using only genes assembled from samples of the focal habitat. For phylum-resolved analyses (Fig. 4d), subtrees corresponding to recognized phyla were iteratively extracted from full marker gene phylogenies (see above) and clade size distributions were recomputed. In both cases, estimated ρ values were summarized as the median across marker genes.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.