Ecology

Filtering threshold for handling spurious sequencesWe first used bacterial communities of known composition (simplified communities) to assess the occurrence of spurious taxa and to determine at which relative abundances they begin to appear. To propose a cutoff that is potentially applicable to different 16S rRNA gene amplicon studies, we included reference data obtained with different variable regions and sequencing pipelines and originating from both in vitro an in vivo communities varying in number and type of species (max. 58) (Tables 1 and 2). To determine a filtering threshold that allowed exclusion of most spurious taxa, we recorded the relative abundance of the first spurious OTU occurring in each of the reference community datasets (Fig. 2a). Median values of approx. 0.12% relative abundance were observed (Fig. 2b). Besides one outlier in the mock communities (0.44% relative abundance), all values were below 0.25% relative abundance.Fig. 2: Determination of filtering thresholds using artificial communities of known composition in vitro (mock; n = 9 different types; 21 replicates in total) and in mice (gnotobiotes; n = 4 different communities; 28 mice in total).a Example of the occurrence of all molecular species detected without filtering in the gut of a gnotobiotic mouse [49]. The arrow indicates the position of the first spurious molecular species, all following taxa being considered as having a high risk of being spurious (light gray bars in the enlarged inset). b Distribution of the relative abundances of first occurring spurious molecular species (as shown in panel a) across all mock communities and samples from gnotobiotes. The orange dashes on the y-axis indicate the consensus threshold of 0.25% relative abundance, above which no spurious taxa occurred with the exception of one outlier in a mock community at a relative abundance of 0.44%. c Comparison of various standard filtering cutoffs (see explanations in the text) in terms of spurious taxa (i.e., those molecular species not matching sequences of the known species contained in the artificial communities). d Corresponding percentages of positive hits retained by the different filtering strategies, with positive hits being defined as the reference sequences found in the respective amplicon datasets. e Percentage of spurious taxa and positive hits in the same reference communities using the DADA2 pipeline for analysis based on amplicon sequence variants (ASVs) [6]. f Effect of filtering thresholds at increments of 0.05% relative abundance on the detection of spurious taxa and positive hits in all mock and gnotobiotic datasets for OTUs (upper panel) and ASVs (lower panel). Lines correspond to mean values; ribbons represent standard deviations.Full size imageWithout any filtering, sequence clustering generated an average of 508 ± 355 OTUs (min. 52; max. 1081) per mock community (10–58 target species in theory) and 105 ± 50 OTUs (min, 55; max. 215) per gnotobiotic community (4–12 target species in theory). Up to 87% of these OTUs were spurious (i.e., they did not match the expected classification of species contained in the corresponding artificial community) (Fig. 2c). On average, the proportion of spurious OTUs in both the mock communities and samples from gnotobiotic mice was slightly lower after removing singletons, although this did not reach statistical significance (50.8 vs. 64.3%, p = 0.227; 57.5% vs. 65.7%; p = 0.70, pairwise comparison by t-test, including Benjamini–Hochberg correction following ANOVA). Interestingly, the proportion of spurious molecular species was higher in gnotobiotic mice independent of filtering (p  0.50) (Fig. 2d). Note that the diversity of reference communities in the gnotobiotic mice was relatively low (4–12 members; Table 2), resulting in a marked drop in the percentage of positive hit (8–25%) when even just one true member is excluded after filtering because of its low relative abundance (which is an expectable event considering a classical, exponentially decreasing distribution of species occurrence in gut environments).We next employed the widely used ASV analysis approach to confirm the aforementioned results. Processing of the same simplified communities generated a total number of 42 ± 25 ASVs (min. 16; max. 98) for mock communities (10–58 target species) and 14 ± 8 ASVs (min. 4; max. 25) for gnotobiotes (4–12 target species). Altogether, a marked decrease in spurious taxa was observed compared with OTU clustering, with an average of 8.6 ± 11.8 and 4.4 ± 6.4% spurious sequences for mock and gnotobiotic communities, respectively (comparison of purple box plots in Fig. 2e, top panels, and Fig. 2c). Of note, the DADA2 pipeline used for the ASV approach does not infer sequence variants that are only supported by a single read (singletons) due to a lack of confidence in their existence relative to sequencing errors. Consequently, data corresponding to “no filtering” with the OTU-based approach were not generated. On average, the first spurious ASV occurred at a relative abundance of 0.10 ± 0.32%. By applying the cutoff of 0.25% relative abundance, spurious sequences were completely removed (except for three outlying samples), albeit with a slight drop in positive hits for both mock and gnotobiotic communities (Fig. 2e).To obtain a more comprehensive view on how filtering thresholds affect the detection of spurious taxa, all datasets (mock and gnotobiotic mice) were processed using a range of relative abundance filtering thresholds (from 0 to 0.5% at increments of 0.05%) after either OTU- or ASV-based processing of raw sequence reads (Fig. 2f). These data indicate that filtering thresholds between 0.1 and 0.3% are appropriate to reduce the occurrence of spurious taxa to 600 of the 678 spurious OTUs occurred in fewer than five of the ten sequencing runs tested, with approximately 450 of them occurring in only one run (Fig. 3c). This observation indicates that the majority of spurious taxa are sporadic cross-contaminations rather than generalist artifacts across sequencing runs, suggesting that fully independent technical replicates would improve data quality. Although most of the spurious taxa were characterized by relative abundances between 0.25 and 2% in the IMNGS-amplicon datasets tested, they represented very dominant populations in a few samples (Fig. 3d).Fig. 3: Origin and occurrence of spurious taxa.a Taxonomic profile and ecological distribution. Inner ring: SILVA-based classification of all non-redundant spurious molecular species at the phylum and family level. Outer colored ring: sample type characterized by the highest prevalence for the given taxon. Outer bars: corresponding highest prevalence values. Only samples with relative abundances >0.25% for any given OTU were counted as positive for prevalence calculation. The total numbers of samples considered were: human, 46,153; soil, 29,864; freshwater, 13,977; mouse, 10,409; marine, 8478. b Distribution of the spurious taxa across sample types. The exclusivity of each OTU for any given sample type was assessed using a Z-test: those assumed to be non-specific for any given sample type appear in red (p 0.25% in at least one replicate were kept). Richness was calculated using ampvis2 [29]. Applying the 0.25% cutoff decreased the number of observed ASVs from 408 ± 71 to 139 ± 5 and, more importantly, the IQR from 101 to 7 (Fig. 6b). Unweighted UniFrac distances within and between runs as calculated using ampvis2 were also compared before and after filtering. Sequences were aligned using MAFFT [30] and phylogeny was inferred using FastTree. Whilst the community makeup in the soil sample varied substantially between sequencing runs without additional filtering, the 0.25% cutoff reduced this variation to the level observed within runs without filtering (Fig. 6c). Replicates within a run were very similar after applying the 0.25% cutoff. Altogether, these data serve as an independent confirmation that stringent filtering delivers more stable values obtained for the exact same sample sequenced in replicates across several sequencing runs. More

2,809 draft MAGs from the rumen ecosystemWe amassed 3.2 terabase pairs (Tbp) of data from 346 publicly available and 66 new rumen metagenome datasets (Supplementary Table 1). The metagenomes were from cattle (312 samples, 2.1 Tbp), sheep (75 samples, 888.4 gigabase pairs (Gbp)), moose (9 samples, 108.8 Gbp), deer (8 samples, 62.9 Gbp), and bison (8 samples, 52.3 Gbp). Metagenomes were assembled independently to reduce the influence of strain variation and improve the recovery of closely related genomes18,19. Following refinement, dereplication, and filtering of resulting population genomes, we identified 2,809 nonredundant MAGs satisfying the following criteria: dRep20 genome quality score ≥60, ≥75% complete, ≤10% contamination, N50 ≥5 kbp, and ≤500 contigs.The median estimated completeness and contamination of the MAGs were 89.7% and 0.9%, respectively (Fig. 1a and Supplementary Data 1). Further, recovered MAGs had a median genome size of 2.2 Mbp, a median of 131 contigs, and a median N50 of 28.3 kbp (Fig. 1b). The proposed minimum information about a MAG (MIMAG) specifies high-quality draft genomes to have an estimated ≥90% completeness, ≤5% contamination, at least 18 tRNAs, and contain 23S, 16S, and 5S rRNA genes21. It remains challenging to reconstruct rRNA genes from short metagenomic reads due to the high sequence similarity of rRNA genes in closely related species. As a result, despite high estimated completeness and low contamination rates, only 20 MAGs meet the MIMAG standards for a high-quality draft genome. We identified a 16S rRNA gene in 197 of the MAGs. The remaining MAGs are characterized as medium-quality MAGs under the MIMAG standards.Fig. 1: Genomic properties of 2,809 rumen MAGs.a CheckM completeness and contamination estimates for the 2,809 population genomes recovered from rumen metagenomes. The size of the point on the scatter plot corresponds to the dRep genome quality score, where Quality = Completeness − (5 ⋅ Contamination) + (Contamination ⋅ (Strain Heterogeneity/100)) + 0.5 ⋅ (({mathrm{log}},)(N50). The reported MAGs meet the following minimum criteria: genome quality score ≥60, ≥75% complete, ≤10% contamination, N50 ≥5 kbp, and ≥500 contigs. b The frequency distribution of the number of contigs and genome sizes of reconstructed MAGs.Full size imageThe majority of bacterial MAGs belonged to phyla Firmicutes or Bacteroidota (2,326; Fig. 2a and Supplementary Data 1). Additionally, we assembled 12 bacterial genomes from the superphylum Patescibacteria. At lower taxonomic ranks, Lachnospiraceae (415) and Prevotella (398) were the dominant family and genus identified among the assembled bacterial genomes. The most prevalent archaeal family and genus were Methanobacteriaceae (45) and Methanobrevibacter (35), respectively (Fig. 2b). The recovered MAGs represent several new taxonomic lineages, as four genomes could not be classified at the rank of order, 16 at the rank of family, and 243 at the genus rank.Fig. 2: Phylogenetic relationships and coverage patterns of near-complete MAGs.a Phylogenomic analysis of 1,163 near-complete (≥90% complete, ≤5% contamination, and N50 ≥15 kbp) bacterial MAGs and (b) 20 near-complete archaeal MAGs inferred from the concatenation of phylogenetically informative proteins. Layers below the genomic trees designate bacterial phylum or archaeal genus based on GTDB taxonomic assignments, genomic size (0–5 Mbp), and the mean number of bases with ≥1× coverage in a rumen metagenomic dataset (layer color indicates the ruminant the data was collected from). The mean number of bases with ≥1× coverage was used as input for hierarchical clustering of rumen metagenomic datasets based on Euclidean distance and Ward linkage. The bacterial and archaeal phylogenetic trees are provided as Supplementary Data 6 and Supplementary Data 7, respectively.Full size imageSpecies-level overlap between reference genomes, the Hungate1000 Collection, and rumen MAGsTo further characterize the assembled genomes, we compared the MAGs to other rumen-specific genome collections, specifically genomes generated from the Hungate1000 project3 and MAGs identified from the Stewart et al. studies4,5. We clustered genomes based on approximate species-level thresholds (≥95% ANI) and calculated the intersection between MAGs in the current study and the Hungate1000 Collection (410 genomes)3, MAGs from Stewart et al. (4,941 genomes)4,5, and a dereplicated genome collection from the GTDB (22,441 genomes, see Methods)22, which includes reference isolate genomes and some environmental MAGs23. It should be noted that we used the raw data from the first of the Stewart et al. studies4 (Supplementary Table 1), but with different assembly and binning approaches. Approximately one-third of the MAGs (1,007) did not exhibit ≥95% ANI with a genome in the GTDB, Stewart et al. MAGs, or the Hungate1000 isolates (Fig. 3a). When considering the pairwise intersections between the datasets, 98 (3.5%), 933 (33.2%), and 1,438 (51.2%) of the MAGs in the current study had ≥95% ANI with a genome in the Hungate1000 Collection3, GTDB22, and Stewart et al.4,5, respectively. One hundred twenty-one (29.5%), 552 (2.5%), and 3,125 (63.2%) of the genomes from the Hungate1000 Collection3, GTDB22, and Stewart et al.4,5 displayed ≥95% ANI with a MAG from the current study. Together, these results indicate that we recovered a majority of previous rumen genomic diversity with additional lineages not previously identified in other major rumen genomic collections.Fig. 3: Genomes sharing ≥95% ANI between databases and the characterization of rumen-specific 95% ANI clusters.a The approximate number of species overlapping amongst rumen-specific and reference genomic datasets. Genomes demonstrating ≥95% ANI were considered to be shared between two datasets. Presented are a subset of intersections in which a MAG from the current study was the query genome. b The number of genomes comprising each of the 3,541 95% ANI clusters generated from 8,160 rumen microbial genomes in the current study, the Hungate1000 Collection3, and Stewart et al. studies4, 5. c Rarefaction analysis based on subsampling 95% ANI clusters at steps of 500 genomes indicates the 8,160 genomes from recently published rumen genomic collections still only represent a fraction of expected microbial species diversity in the rumen ecosystem. Phylogenomic relationships of the 1,781 near-complete bacterial (d) and 35 near-complete archaeal (e) representative genomes with the highest dRep genome quality score from the 3,541 95% ANI clusters generated from 8,160 rumen-specific genomes. Near-complete genomes were defined as being ≥90% complete, having ≤5% contamination, and contig N50 ≥15 kbp. Layers surrounding the genomic trees indicate the bacterial phyla or archaeal genera and the log normalized number of genomes from each rumen genomic collection belonging to the same 95% ANI cluster. The bacterial and archaeal phylogenetic trees are provided as Supplementary Data 8 and Supplementary Data 9, respectively.Full size imageWe applied an additional clustering approach to identify the approximate number of species represented by the rumen-specific genomes assembled in this study, in the Hungate1000 Collection3, and Stewart et al.4,5. A 95% ANI threshold yielded 3,541 clusters from the combination of the datasets (Supplementary Data 2). Of the 3,541 clusters, 2,024 contained a MAG from the current study, and 1,135 were composed exclusively of MAGs from the current study. In comparison, 2,175 and 286 clusters were comprised of genomes from Stewart et al.4,5 and the Hungate1000 Collection3, respectively. The majority of 95% ANI clusters (2,166) are only comprised of a single genome (Fig. 3b). Furthermore, a rarefaction curve suggests the 8,160 genomes from the genomic collections analyzed here only represent a fraction of the estimated microbial species diversity in the rumen (Fig. 3c). The genome with the best dRep score from each cluster was used to generate a phylogenetic tree highlighting the species diversity within each rumen genomic collection and represents the vast diversity of rumen bacterial (Fig. 3d) and archaeal (Fig. 3e) genomes published to date.As stated previously, the median genome size of reconstructed MAGs was 2.2 Mbp, smaller than the median size of genomes from the Hungate1000 project (3.1 Mbp)3. To provide an assessment at a finer resolution, genome sizes of MAGs and Hungate1000 genomes3 belonging to the same 95% ANI cluster were compared (Supplementary Fig. 1). Adjusted sizes of MAGs and Hungate1000 genomes that are ≥95% complete displayed a regression coefficient of 0.96 with a slope of 0.86, indicating the binning process likely did not lead to extensive losses and systematic biases in the reconstructed genomes. Instead, it further highlights that current culturing approaches have not brought large portions of rumen microbial diversity into culture and putatively supports previous findings from the human gut that revealed genome-reduction in uncultured bacteria24.Rumen metagenome classification rates using reference and rumen-specific genomesUtilizing an approach similar to Stewart et al.4,5, we investigated the influence of MAGs on rates of metagenomic read classification. The baseline for read classification was the standard Kraken database containing bacterial, archaeal, fungal, and protozoal RefSeq genomes25. Each rumen-specific dataset was incrementally added to the Kraken RefSeq genomic database in the following order to build new databases: the Hungate1000 Collection3, MAGs from Stewart et al.4,5, and MAGs from the current study. Each individual and collective database was used for classification of sample reads that underpinned metagenomic binning and from a rumen metagenomic dataset not used in the reconstruction of MAGs26. MAGs from the current work classified more reads from deer, moose, and sheep metagenomes, while the more numerous MAGs from Stewart et al.4,5 classified more reads from bison and cattle metagenomes (Supplementary Fig. 2a). The addition of MAGs improves classification relative to databases primarily based on cultured isolates, like the Hungate1000 Collection3 (Supplementary Fig. 2b). Using the combination of all reference and rumen-specific genomes, the median classification rate on an independent set of cattle metagenomes was 62.6%.Phylogenetic characterization of biosynthetic gene clustersMicrobial genome mining is a powerful tool for natural product discovery. We sought to explore the extent of secondary metabolite diversity coded by the MAGs in the current study, the Hungate1000 Collection3, and Stewart et al. MAGs4,5. We identified 14,814 BGCs encoded by the 8,160 rumen-specific genomes using antiSMASH27 (Fig. 4a and Supplementary Data 3). The majority of BGCs were NRPS (5,346), followed by aryl polyenes (2,800), sactipeptides (2,126), and bacteriocins (1,943). Only a few PKS were identified (75). Firmicutes harbored the vast majority of clusters for NRPS, sactipeptide, lantipeptide, lassopeptide, and bacteriocin synthesis (Fig. 4b). At lower taxonomic ranks, DTU089 (979), Bacteroidaceae (934), and Lachnospiraceae (923) coded for the bulk of NRPS gene clusters. Moreover, Acidaminococcaceae genomes contained 21.2% of identified bacteriocins and Ruminococcus spp. possessed the bulk of sactipeptides and lantipeptides. Archaea were predicted to code 737 BGCs, including an average of 3.8 NRPS gene clusters per genome (Fig. 4a).Fig. 4: Characterization of BGCs from 8,160 rumen genomes and MAGs.a Number and types of BGCs identified from select phyla in genomes from the Hungate1000 Collection3, Stewart et al. studies4, 5, and the current study. b Phylogenomic analysis of 1,766 near-complete Firmicutes genomes inferred from the concatenation of phylogenetically informative proteins. The inner layer surrounding the genomic tree designates taxonomic annotations, while the remaining layers depict the log normalized number of BGCs in the genome with the ascribed function. Bacterial class and order labels are displayed for those lineages in which more than 50 genomes were identified. Near-complete genomes were defined as being ≥90% complete, having ≤5% contamination, and contig N50 ≥15 kbp. The phylogenetic tree is provided as Supplementary Data 10. c A relational network of NRPS gene clusters in Firmicutes, Bacteroidota, and Euryarchaeota highlights the similarity of NRPS BGCs from Euryarchaeota and Firmicutes. Edge weight represents the similarity of two BGCs, as determined by BiG-SCAPE (i.e. darker edges demonstrate more similarity between two BGCs). Edges are only shown for BGCs with ≥0.3 BiG-SCAPE similarity. Nodes from each phylum are duplicated to illustrate intra-phylum relationships and nodes along a given axis are ordered alphabetically by taxonomic family. d The association between genome phylogeny and the similarity of NRPS gene clusters coded by near-complete Euryarchaeota genomes. BGCs designated as NRPS were clustered with BiG-SCAPE. The relationship between NRPS clusters was portrayed through the hierarchical clustering of pairwise inter-cluster similarities. The number of NRPS clusters coded by each genome (range of 0–3) is presented alongside the assigned genus. A group of Methanobrevibacter genomes, likely of the same species (≥95% ANI), possessed very similar NRPS clusters (highlighted in blue). Yet, phylogenetically closely related genomes, belonging to two different 95% ANI clusters, did not code for any identified NRPS gene clusters (highlighted in red). The phylogenetic tree is based on the concatenation of 122 phylogenetically informative archaeal proteins and is available as Supplementary Data 11.Full size imageNRPS exhibit high molecular and structural diversity resulting in a wide array of biological activities. The diversity of NRPS, combined with their proteolytic stability and selective bioactivity, has resulted in the development of many NRPS as antimicrobials and other therapeutic agents28. Given the prevalence of NRPS among the recovered MAGs (Fig. 4a), the peptides appear to be important bioactive metabolites in the rumen. To gain fundamental insight into the phylogenetic diversity of rumen NRPS, we built a network based on BGC similarity using BiG-SCAPE29. BiG-SCAPE uses protein domain content, order, copy number, and sequence identity to calculate a distance metric. We assessed the similarity of NRPS gene clusters identified in Firmicutes, Bacteroidota, and Euryarchaeota, as these three phyla coded for 96.4% of assembled NRPS gene clusters from rumen genomes. With a BiG-SCAPE similarity threshold of 0.3, the resulting network consisted of 3,436 nodes (NRPS BGCs on contigs ≥10 kbp) and 79,112 edges (Fig. 4c and Supplementary Data 4). As expected, the network analysis depicted high inter- and intra-phylum genetic diversity among the NRPS gene clusters. The median intra-phylum, -family, and -genus similarity was 0.40, 0.44, and 0.46, respectively, while the median inter-phylum, -family, and -genus similarity was 0.32, 0.34, and 0.34, respectively. Further, only 2.6% of edges were inter-phylum and 69.0% were intra-family. Of the 6,594 Euryarchaeota edges, 8.1% were Euryarchaeota-Firmicutes (median similarity of 0.32) and 2.0% of edges were Euryarchaeota-Bacteroidota (median similarity of 0.31). To further examine the phylogenetic relationships of rumen Euryarchaeota NRPS, we clustered 265 NRPS gene clusters (≥10 kbp) from 85 near-complete Euryarchaeota genomes at a higher similarity threshold of 0.75, yielding 57 NRPS clusters (Fig. 4d). The distribution of NRPS clusters amongst the genomes suggests there exists a strong relationship between methanogen phylogeny and NRPS similarity. Only Methanobrevibacter genomes contain NRPS gene clusters, and genomes of the same species often possessed many of the same NRPS clusters (see genomes highlighted in blue in Fig. 4d). However, there are instances in which closely related methanogens code for a contrasting pattern of NRPS clusters or no NRPS clusters at all (see genomes highlighted in red in Fig. 4d).Bacteriocins likely serve as regulatory elements in complex microbial communities such as the rumen. Consequently, bacteriocins have been studied and characterized for their bactericidal activity and as agents that modulate the microbiota structure and function30. In particular, lanthipeptides, a class of ribosomally synthesized and post-translationally modified peptides (RiPPs) with thioether cross-linked amino acids31, are of pharmaceutical, preservative, and agricultural interest due to their strong antimicrobial properties against gram-positive pathogens31,32,33, low levels of antimicrobial resistance34, and stability35. We identified 195 rumen lanthipeptide BGCs from the Hungate1000 genomes and MAGs from Stewart et al. and the current study. Rumen lanthipeptide BGCs were clustered with 22,870 lanthipeptide BGCs from RefSeq genomes36,37 into gene cluster families (GCFs; groups of BGCs that may generate highly similar products). Clustering with BiG-SCAPE29 yielded 4,565 GCFs, 120 of which contained a rumen lanthipeptide. The 120 GCFs were composed of 519 lanthipeptide BGCs, where 324 were from RefSeq isolates and 195 from rumen genomes (Fig. 5a). The 324 RefSeq BGCs fell into only 18 GCFs. Lanthipeptides from the Hungate1000 isolates clustered into 36 GCFs, while rumen MAG lanthipeptides belonged to 92 GCFs, 82 of which were exclusively composed of MAG lanthipeptides. Together, this evidence suggests rumen MAGs code for diverse and novel lanthipeptides not represented in cultured isolates, including the Hungate Collection.Fig. 5: Phylogenetic diversity of 195 lanthipeptide BGCs coded by rumen genomes.a Network depicting the similarity between lanthipeptide BGCs identified from complete and draft isolate genomes in RefSeq and rumen genomes of the Hungate1000 collection, Stewart et al. MAGs, and MAGs from the current study. The BGCs were clustered into gene cluster families (GCFs) with BiG-SCAPE29. Only the GCFs containing a rumen genome and at least two BGCs were visualized. Nodes in the network represent BGCs and edges connect BGCs with BiG-SCAPE defined similarity ≥0.3. b Phylogenetic relationships of 120 near-complete rumen bacterial genomes coding for lanthipeptide BGCs. Near-complete genomes were defined as being ≥90% complete, having ≤5% contamination, and contig N50 ≥15 kbp. Layers surrounding the genomic trees indicate the bacterial phyla and family, if the genome is a MAG or Hungate Collection isolate, and the class of lanthipeptide, as predicted by antiSMASH27. Genomes without an indicated lanthipeptide class were not classified by antiSMASH. The phylogenetic tree is based on the concatenation of 120 phylogenetically informative bacterial proteins and is available as Supplementary Data 12.Full size imageWe sought to further examine the differences in rumen MAG lanthipeptides relative to isolates and the taxonomic diversity of rumen microbes coding for lanthipeptides. The 195 rumen lanthipeptides were mainly found in Firmicutes genomes, with a subset from Bacteroidota and Actinobacteriota (Fig. 5b). Fifty-two of the 55 lanthipeptides from the Hungate Collection isolates were from Firmicutes (94.5%). At the family-level, these 52 Firmicutes BGCs were distributed evenly between Lachnospiraceae and Streptococcaceae. In contrast, 19.2% and 8.6% of lanthipeptides from rumen MAGs belonged to Bacteroidota and Actinobacteriota, respectively. Lanthipeptides from MAGs were also found in Muribaculaceae and Oscillospiraceae. Moreover, 26.4% of rumen MAG lanthipeptides, compared to 3.6% of Hungate Collection isolates, were found in Eubacterium genomes. The majority of Eubacterium MAG lanthipeptides (62.1%) belonged to a single GCF, suggesting they code for very similar products. Lastly, antiSMASH predicted the bulk of the rumen lanthipeptides were Class II lanthipeptides, with fewer Class I and Class III types (Fig. 5b). Nearly all of the Class I lanthipeptides were from Hungate isolates. The above analysis of lanthipeptide diversity further supports that rumen MAGs code for novel secondary metabolites not represented in cultured isolates.We aligned previously published rumen metatranscriptome data from steers characterized as having high and low feed efficiency to the BGCs to demonstrate if the identified BGCs are active and to explore potential ecological roles of secondary metabolites. Despite data from the metatranscriptome study not being applied to reconstruct genomes in the current study, we identified the expression of 554 gene clusters from rumen-specific genomes in the 20 metatranscriptomes (≥100 aligned reads). Metatranscriptome read count data were normalized independently for each genome to better account for the variation in taxonomic composition across samples38. Genome-specific normalization resulted in the identification of 17 differentially expressed gene clusters between steers with high and low feed efficiency (DESeq239 false discovery rate adjusted P  More

1.Bálint M, Domisch S, Engelhardt CHM, Haase P, Lehrian S, Sauer J, et al. Cryptic biodiversity loss linked to global climate change. Nat Clim Chang. 2011;1:313–8.Article 

Google Scholar 
2.Parmesan C, Yohe G. A globally coherent fingerprint of climate change impacts across natural systems. Nature. 2003;421:37–42.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Blois JL, Zarnetske PL, Fitzpatrick MC, Finnegan S. Climate change and the past, present, and future of biotic interactions. Science. 2013;341:499–504.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Haines A, Ebi K. The imperative for climate action to protect health. N. Engl J Med. 2019;380:263–73.PubMed 
Article 
PubMed Central 

Google Scholar 
5.Deutsch CA, Tewksbury JJ, Tigchelaar M, Battisti DS, Merrill SC, Huey RB, et al. Increase in crop losses to insect pests in a warming climate. Science. 2018;361:916–9.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Kattwinkel M, Jan-Valentin K, Foit K, Liess M. Climate change, agricultural insecticide exposure, and risk for freshwater communities. Ecol Appl. 2011;21:2068–81.PubMed 
Article 
PubMed Central 

Google Scholar 
7.Moe SJ, De Schamphelaere K, Clements WH, Sorensen MT, Van den Brink PJ, Liess M. Combined and interactive effects of global climate change and toxicants on populations and communities. Environ Toxicol Chem. 2013;32:49–61.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Potts SG, Biesmeijer JC, Kremen C, Neumann P, Schweiger O, Kunin WE. Global pollinator declines: trends, impacts and drivers. Trends Ecol Evol. 2010;25:345–53.PubMed 
Article 
PubMed Central 

Google Scholar 
9.Moran EV, Alexander JM. Evolutionary responses to global change: Lessons from invasive species. Ecol Lett. 2014;17:637–49.PubMed 
Article 
PubMed Central 

Google Scholar 
10.Harwood AD, You J, Lydy MJ. Temperature as a toxicity identification evaluation tool for pyrethroid insecticides: toxicokinetic confirmation. Environ Toxicol Chem. 2009;28:1051–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Guo L, Su M, Liang P, Li S, Chu D. Effects of high temperature on insecticide tolerance in whitefly Bemisia tabaci (Gennadius) Q biotype. Pestic Biochem Physiol. 2018;150:97–104.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Mao K, Jin R, Li W, Ren Z, Qin X, He S, et al. The influence of temperature on the toxicity of insecticides to Nilaparvata lugens (Stål). Pestic Biochem Physiol. 2019;156:80–86.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Verheyen J, Delnat V, Stoks R. Increased daily temperature fluctuations overrule the ability of gradual thermal evolution to offset the increased pesticide toxicity under global warming. Environ Sci Technol. 2019;53:4600–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Moran NA. Symbiosis as an adaptive process and source of phenotypic complexity. Proc Natl Acad Sci USA. 2007;104:8627–3863.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Kikuchi Y, Hayatsu M, Hosokawa T, Nagayama A, Tago K, Fukatsu T. Symbiont-mediated insecticide resistance. Proc Natl Acad Sci USA. 2012;109:8618–22.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Jones RM, Desai C, Darby TM, Luo L, Wolfarth AA, Scharer CD, et al. Lactobacilli modulate epithelial cytoprotection through the Nrf2 pathway. Cell Rep. 2015;12:1217–25.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Cheng D, Guo Z, Riegler M, Xi Z, Liang G, Xu Y. Gut symbiont enhances insecticide resistance in a significant pest, the oriental fruit fly Bactrocera dorsalis (Hendel). Microbiome. 2017;5:13.PubMed 
PubMed Central 
Article 

Google Scholar 
18.Pang R, Chen M, Yue L, Xing K, Li T, Kang K, et al. A distinct strain of Arsenophonus symbiont decreases insecticide resistance in its insect host. PLoS Genet. 2018;14:e1007725.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
19.Kikuchi Y, Tada A, Musolin DL, Hari N, Hosokawa T, Fujisaki K, et al. Collapse of insect gut symbiosis under simulated climate change. mBio. 2016;7:e01578–16.PubMed 
PubMed Central 
Article 

Google Scholar 
20.Corbin C, Heyworth ER, Ferrari J, Hurst GDD. Heritable symbionts in a world of varying temperature. Heredity. 2017;118:10–20.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Jia FX, Yang MS, Yang WJ, Wang JJ. Influence of continuous high temperature conditions on Wolbachia infection frequency and the fitness of Liposcelis tricolor (Psocoptera: Liposcelididae). Environ Entomol. 2009;38:1365–72.PubMed 
Article 
PubMed Central 

Google Scholar 
22.Burke G, Fiehn O, Moran N. Effects of facultative symbionts and heat stress on the metabolome of pea aphids. ISME J. 2010;4:242–52.PubMed 
Article 
PubMed Central 

Google Scholar 
23.Fan Y, Wernegreen JJ. Can’t take the heat: high temperature depletes bacterial endosymbionts of ants. Micro Ecol. 2013;66:727–33.Article 

Google Scholar 
24.Hussain M, Akutse KS, Ravindran K, Lin Y, Bamisile BS, Qasim M, et al. Effects of different temperature regimes on survival of Diaphorina citri and its endosymbiotic bacterial communities. Environ Microbiol. 2017;19:3439–49.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Engl T, Eberl N, Gorse C, Krüger T, Schmidt THP, Plarre R, et al. Ancient symbiosis confers desiccation resistance to stored grain pest beetles. Mol Ecol. 2018;27:2095–108.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Zhang XJ, Yu XP, Chen JM. High Temperature effects on yeast-like endosymbiotes and pesticide resistance of the small brown planthopper, Laodelphax striatellus. Rice Sci. 2008;15:326–30.CAS 
Article 

Google Scholar 
27.Zhang B, Zuo TQ, Li HG, Sun LJ, Wang SF, Zhang CY, et al. Effect of heat shock on the susceptibility of Frankliniella occidentalis (Thysanoptera: Thripidae) to insecticides. J Integr Agric. 2016;15:2309–18.CAS 
Article 

Google Scholar 
28.Karimzadeh R, Javanshir M, Hejazi MJ. Individual and combined effects of insecticides, inert dusts and high temperatures on Callosobruchus maculatus (Coleoptera: Chrysomelidae). J Stored Prod Res. 2020;89:10693.Article 

Google Scholar 
29.Michigan State University. Arthropod Pesticide Resistance Database (APRD). East Lansing: Michigan State University; 2020. http://www.pesticideresistance.com/.30.Ju JF, Bing XL, Zhao DS, Guo Y, Xi Z, Hoffmann AA, et al. Wolbachia supplement biotin and riboflavin to enhance reproduction in planthoppers. ISME J. 2019;14:676–87.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
31.Zhang Y, Tang T, Li W, Cai T, Li J, Wan H. Functional profiling of the gut microbiomes in two different populations of the brown planthopper. Nilaparvata lugens J Asia Pac Entomol. 2018;21:1309–14.Article 

Google Scholar 
32.Ye YH, Seleznev A, Flores HA, Woolfit M, McGraw EA. Gut microbiota in Drosophila melanogaster interacts with Wolbachia but does not contribute to Wolbachia-mediated antiviral protection. J Invertebr Pathol. 2017;143:18–25.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Yamada R, Floate KD, Riegler M, O’Neill SL. Male development time influences the strength of Wolbachia-induced cytoplasmic incompatibility expression in Drosophila melanogaster. Genetics. 2007;177:801–8.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Wari D, Kabir MA, Mujiono K, Hojo Y, Shinya T, Tani A, et al. Honeydew-associated microbes elicit defense responses against brown planthopper in rice. J Exp Bot. 2019;70:1683–96.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Miller ALE, Tindall K, Leonard BR. Bioassays for monitoring insecticide resistance. J Vis Exp. 2010;46:2129.
Google Scholar 
36.Zhang J, Zhang Y, Wang Y, Yang Y, Cang X, Liu Z. Expression induction of P450 genes by imidacloprid in Nilaparvata lugens: a genome-scale analysis. Pestic Biochem Physiol. 2016;132:59–64.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
37.Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods. 2001;25:402–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Noda H, Koizumi Y, Zhang Q, Deng K. Infection density of Wolbachia and incompatibility level in two planthopper species, Laodelphax striatellus and Sogatella furcifera. Insect Biochem Mol Biol. 2001;31:727–37.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011. https://doi.org/10.14806/ej.17.1.20040.Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from illumina amplicon data. Nat Methods. 2016;13:581–3.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Katoh K, Misawa K, Kuma KI, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome. 2018;6:90.PubMed 
PubMed Central 
Article 

Google Scholar 
43.Liu S, Ding Z, Zhang C, Yang B, Liu Z. Gene knockdown by intro-thoracic injection of double-stranded RNA in the brown planthopper, Nilaparvata lugens. Insect Biochem Mol Biol. 2010;40:666–71.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
44.Tai V, James ER, Nalep CA, Scheffrahn RH, Perlman SJ, Keelinga PJ. The role of host phylogeny varies in shaping microbial diversity in the hindguts of lower termites. Appl Environ Microbiol. 2015;81:1059–70.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
45.Bale JS, Hayward SAL. Insect overwintering in a changing climate. J Exp Biol. 2010;213:980–94.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Rahmstorf S, Cazenave A, Church JA, Hansen JE, Keeling RF, Parker DE, et al. Recent climate observations compared to projections. Science. 2007;316:709.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Radchuk V, Reed T, Teplitsky C, van de Pol M, Charmantier A, Hassall C, et al. Adaptive responses of animals to climate change are most likely insufficient. Nat Commun. 2019;10:3019.Article 
CAS 

Google Scholar 
48.Iwamura T, Guzman-Holst A, Murray KA. Accelerating invasion potential of disease vector Aedes aegypti under climate change. Nat Commun. 2020;11:2130.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Li J, Mao T, Wang H, Lu Z, Qu J, Fang Y, et al. The CncC/keap1 pathway is activated in high temperature-induced metamorphosis and mediates the expression of Cyp450 genes in silkworm, Bombyx mori. Biochem Biophys Res Commun. 2019;541:1045–50.Article 
CAS 

Google Scholar 
50.Kalsi M, Palli SR. Transcription factor cap n collar C regulates multiple cytochrome P450 genes conferring adaptation to potato plant allelochemicals and resistance to imidacloprid in Leptinotarsa decemlineata (Say). Insect Biochem Mol Biol. 2017;83:1–12.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Kalsi M, Palli SR. Transcription factors, CncC and Maf, regulate expression of CYP6BQ genes responsible for deltamethrin resistance in Tribolium castaneum. Insect Biochem Mol Biol. 2015;65:47–56.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Misra JR, Lam G, Thummel CS. Constitutive activation of the Nrf2/Keap1 pathway in insecticide-resistant strains of Drosophila. Insect Biochem Mol Biol. 2013;43:1116–24.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
53.Tang B, Cheng Y, Li Y, Li W, Ma Y, Zhou Q, et al. Adipokinetic hormone regulates cytochrome P450-mediated imidacloprid resistance in the brown planthopper, Nilaparvata lugens. Chemosphere. 2020;259:127490.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Cheng Y, Li Y, Li W, Song Y, Zeng R, Lu K. Inhibition of hepatocyte nuclear factor 4 confers imidacloprid resistance in Nilaparvata lugens via the activation of cytochrome P450 and UDP-glycosyltransferase genes. Chemosphere. 2021;263:128269.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
55.Li Y, Liu X, Wang N, Zhang Y, Hoffmann AA, Guo H. Background-dependent Wolbachia-mediated insecticide resistance in Laodelphax striatellus. Environ Microbiol. 2020;22:2653–63.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
56.Berticat C, Rousset F, Raymond M, Berthomieu A, Weill M. High Wolbachia density in insecticide-resistant mosquitoes. Proc R Soc B Biol Sci. 2002;269:1413–6.Article 

Google Scholar 
57.Zhang G, Hussain M, O’Neill SL, Asgari S. Wolbachia uses a host microRNA to regulate transcripts of a methyltransferase, contributing to dengue virus inhibition in Aedes aegypti. Proc Natl Acad Sci USA. 2013;110:10276–81.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
58.Bi J, Sehgal A, Williams JA, Wang YF. Wolbachia affects sleep behavior in Drosophila melanogaster. J Insect Physiol. 2018;107:81–88.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Roughgarden J, Gilbert SF, Rosenberg E, Zilber-Rosenberg I, Lloyd EA. Holobionts as units of selection and a model of their population dynamics and evolution. Biol Theory. 2018;13:44–65.Article 

Google Scholar 
60.Pan X, Zhou G, Wu J, Bian G, Lu P, Raikhel AS, et al. Wolbachia induces reactive oxygen species (ROS)-dependent activation of the Toll pathway to control dengue virus in the mosquito Aedes aegypti. Proc Natl Acad Sci USA. 2012;109:E23–31.PubMed 
Article 
PubMed Central 

Google Scholar 
61.Gong JT, Li Y, Li TP, Liang Y, Hu L, Zhang D, et al. Stable introduction of plant-virus-inhibiting Wolbachia into planthoppers for rice protection. Curr Biol. 2020;30:4837–45.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
62.Elzaki MEA, Li ZF, Wang J, Xu L, Liu N, Zeng RS, et al. Activiation of the nitric oxide cycle by citrulline and arginine restores susceptibility of resistant brown planthoppers to the insecticide imidacloprid. J Hazard Mater. 2020;396:122755.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
63.Werren JH. Biology of Wolbachia. Annu Rev Entomol. 1997;42:587–609.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
64.Kokou F, Sasson G, Nitzan T, Doron-Faigenboim A, Harpaz S, Cnaani A, et al. Host genetic selection for cold tolerance shapes microbiome composition and modulates its response to temperature. Elife. 2018;77:e36398.Article 

Google Scholar  More

1.Arim, M. & Marquet, P. A. Intraguild predation: A widespread interaction related to species biology. Ecol. Lett. 7, 557–564. https://doi.org/10.1111/j.1461-0248.2004.00613.x (2004).Article 

Google Scholar 
2.Polis, G. A., Myers, C. A. & Holt, R. D. The ecology and evolution of intraguild predation: Potential competitors that eat each other. Annu. Rev. Ecol. Syst. 20, 297–330. https://doi.org/10.1146/annurev.es.20.110189.001501 (1989).Article 

Google Scholar 
3.Rosenheim, J. A., Kaya, H. K., Ehler, L. E., Marois, J. J. & Jaffee, B. A. Intraguild predation among biological-control agents: Theory and evidence. Biol. Control 5, 303–335. https://doi.org/10.1006/bcon.1995.1038 (1995).Article 

Google Scholar 
4.Rosenheim, J. A. & Harmon, J. P. The influence of intraguild predation on the suppression of a shared prey population: An empirical reassessment. In Trophic and Guild in Biological Interactions Control (eds Brodeur, J. & Boivin, G.) 1–20 (Springer, 2006) https://doi.org/10.1007/1-4020-4767-3_1.Chapter 

Google Scholar 
5.Fonseca, M. M. et al. How to evaluate the potential occurrence of intraguild predation. Exp. Appl. Acarol. 72, 103–114. https://doi.org/10.1007/s10493-017-0142-x (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
6.Ferguson, K. I. & Stiling, P. Non-additive effects of multiple natural enemies on aphid populations. Oecologia 108, 375–379 (1996).ADS 
Article 

Google Scholar 
7.Hindayana, D., Meyhöfer, R., Scholz, D. & Poehling, H.-M. Intraguild predation among the hoverfly Episyrphus balteatus de Geer (Diptera: Syrphidae) and other aphidophagous predators. Biol. Control 20, 236–246 (2001).Article 

Google Scholar 
8.Denoth, M., Frid, L. & Myers, J. H. Multiple agents in biological control: Improving the odds?. Biol. Control 24, 20–30. https://doi.org/10.1016/S1049-9644(02)00002-6 (2002).Article 

Google Scholar 
9.Muştu, M., Kilinçer, N., Ülgentürk, S. & Kaydan, M. B. Feeding behavior of Cryptolaemus montrouzieri on mealybugs parasitized by Anagyrus pseudococci. Phytoparasitica 36, 360–367 (2008).Article 

Google Scholar 
10.Lucas, E. Intraguild predation among aphidophagous predators. Eur. J. Entomol. 102, 351–364 (2005).
Google Scholar 
11.Muştu, M. & Kilinçer, N. Intraguild predation of Planococcus ficus parasitoids Anagyrus pseudococci and Leptomastix dactylopii by Nephus kreissli. Biocontrol Sci. Technol. 24, 257–269. https://doi.org/10.1080/09583157.2013.856866 (2014).Article 

Google Scholar 
12.Diehl, S. & Feißel, M. Effects of enrichment on three-level food chains with omnivory. Am. Nat. 155, 200–218 (2000).Article 

Google Scholar 
13.Holt, R. D. & Polis, G. A. A theoretical framework for intraguild predation. Am. Nat. 149, 745–764 (1997).Article 

Google Scholar 
14.Kuijper, L. D. J., Kooi, B. W., Zonneveld, C. & Kooijman, S. A. L. M. Omnivory and food web dynamics. Ecol. Modell. 163, 19–32 (2003).Article 

Google Scholar 
15.Morin, P. Productivity, intraguild predation, and population dynamics in experimental food webs. Ecology 80, 752–760 (1999).Article 

Google Scholar 
16.Mylius, S. D., Klumpers, K., de Roos, A. M. & Persson, L. Impact of intraguild predation and stage structure on simple communities along a productivity gradient. Am. Nat. 158, 259–276 (2001).CAS 
Article 

Google Scholar 
17.Polis, G. A. & Holt, R. D. Intraguild predation: The dynamics of complex trophic interactions. Trends Ecol. Evol. 7, 151–154 (1992).CAS 
Article 

Google Scholar 
18.Janssen, A. et al. Intraguild predation usually does not disrupt biological control. In Trophic and Guild in Biological Interactions Control (eds Brodeur, J. & Boivin, G.) 21–44 (Springer, 2006) https://doi.org/10.1007/1-4020-4767-3_2.Chapter 

Google Scholar 
19.Skalski, G. T. & Gilliam, J. F. Functional responses with predator interference: Viable alternatives to the Holling type II model. Ecology 82, 3083–3092 (2001).Article 

Google Scholar 
20.de Villemereuil, P. B. & López-Sepulcre, A. Consumer functional responses under intra- and inter-specific interference competition. Ecol. Modell. 222, 419–426. https://doi.org/10.1016/j.ecolmodel.2010.10.011 (2011).Article 

Google Scholar 
21.Sutherland, W. J. From Individual Behaviour to Population Ecology (Oxford Series in Ecology and Evolution, 1996).
Google Scholar 
22.Pedersen, B. S. & Mills, N. J. Single vs. multiple introduction in biological control: The roles of parasitoid efficiency, antagonism and niche overlap. J. Appl. Ecol. 41, 973–984 (2004).Article 

Google Scholar 
23.Godfray, H. C. J. & Godfray, H. C. J. Parasitoids: Behavioral and Evolutionary Ecology Vol. 67 (Princeton University Press, 1994).Book 

Google Scholar 
24.Harvey, J. A., Poelman, E. H. & Tanaka, T. Intrinsic inter-and intraspecific competition in parasitoid wasps. Annu. Rev. Entomol. 58, 333–351 (2013).CAS 
Article 

Google Scholar 
25.Peri, E., Cusumano, A., Amodeo, V., Wajnberg, E. & Colazza, S. Intraguild interactions between two egg parasitoids of a true bug in semi-field and field conditions. PLoS ONE 9(6), e99876. https://doi.org/10.1371/journal.pone.0099876 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
26.Bruzzone, O. A., Logarzo, G. A., Aguirre, M. B. & Virla, E. G. Intra-host interspecific larval parasitoid competition solved using modelling and bayesian statistics. Ecol. Modell. 385, 114–123 (2018).Article 

Google Scholar 
27.Triapitsyn, S. V. et al. Complex of primary and secondary parasitoids (Hymenoptera: Encyrtidae and Signiphoridae) of Hypogeococcus spp. mealybugs (Hemiptera: Pseudococcidae) in the New World. Florida Entomol. 101, 411–434. https://doi.org/10.1653/024.101.0320 (2018).Article 

Google Scholar 
28.Aguirre, M. B. et al. Analysis of biological traits of Anagyrus cachamai and Anagyrus lapachosus to assess their potential as biological control candidate agents against Harrisia cactus mealybug pest in Puerto Rico. Biocontrol 64, 539–551. https://doi.org/10.1007/s10526-019-09956-y (2019).CAS 
Article 

Google Scholar 
29.Poveda-Martínez, D. et al. Species complex diversification by host plant use in an herbivorous insect: The source of Puerto Rican cactus mealybug pest and implications for biological control. Ecol. Evol. 10, 10463–10480. https://doi.org/10.1002/ece3.6702 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
30.Poveda-Martínez, D. et al. Untangling the Hypogeococcus pungens species complex (Hemiptera: Pseudococcidae) for Argentina, Australia, and Puerto Rico based on host plant associations and genetic evidence. PLoS ONE 14(7), e0220366. https://doi.org/10.1371/journal.pone.0220366 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
31.Thurstone, L. L. A law of comparative judgment. Psychol. Rev. 34, 273 (1927).Article 

Google Scholar 
32.Bradley, R. A. & Terry, M. E. Rank analysis of incomplete block designs: I. The method of paired comparisons. Biometrika 39, 324–345. https://doi.org/10.2307/2334029 (1952).MathSciNet 
Article 
MATH 

Google Scholar 
33.Gelman, A., Carlin, J. B., Stern, H. S. & Rubin, D. B. Bayesian data analysis. In Texts Stat. Sci. 2nd ed, 661 (CRC Press, 2003).34.Stevens, S. S. On the Theory of Scales of Measurement, vol. 103, 677–680 (1946).35.Patil, A., Huard, D. & Fonnesbeck, C. J. PyMC: Bayesian stochastic modelling in Python. J. Stat. Softw. 35, 1 (2010).Article 

Google Scholar 
36.Zwölfer, H. The structure and effect of parasite complexes attacking phytophagous host insects. In Proc. Adv. Study Inst. Dyn. Numbers Popul. 405–418 (1971).37.Zwölfer, H. Strategies and counterstrategies in insect population systems competing for space and food in flower headsand plant galls. Fortschr. Zool. 25(2/3), 331–353 (1979).
Google Scholar 
38.Vance, R. R. The stable coexistence of two competitors for one resource. Am. Nat. 126, 72–86 (1985).Article 

Google Scholar 
39.Fellers, J. H. Interference and exploitation in a guild of woodland ants. Ecology 68, 1466–1478 (1987).Article 

Google Scholar 
40.Cusumano, A., Peri, E., Vinson, S. B. & Colazza, S. Intraguild interactions between two egg parasitoids exploring host patches. Biocontrol 56, 173–184 (2011).Article 

Google Scholar 
41.Mizutani, N. Interspecific larval competition among three egg parasitoid species on the host, Riptortus clavatus (Thunberg) (Heteroptera: Alydidae). Proc. Assoc. Plant Prot. Kyushu 40, 106–110 (1994).Article 

Google Scholar 
42.Weber, C. A., Smilanick, J. M., Ehler, L. E. & Zalom, F. G. Ovipositional behavior and host discrimination in three scelionid egg parasitoids of stink bugs. Biol. Control 6, 245–252 (1996).Article 

Google Scholar 
43.Alim, M. A. & Lim, U. T. Interspecific larval competition between two egg parasitoids in refrigerated host eggs of Riptortus pedestris (Hemiptera: Alydidae). Biocontrol Sci. Technol. 21, 395–407 (2011).Article 

Google Scholar 
44.De Moraes, C. M. & Lewis, W. J. Analyses of two parasitoids with convergent foraging strategies. J. Insect Behav. 12, 571–583 (1999).Article 

Google Scholar 
45.Mackauer, M. Host discrimination and larval competition in solitary endoparasitoids. Crit. Issues Biol. Control, Intercept, Andover, Hants, UK. xvii + 330 pp (1990).46.Strand, M. R. & Godfray, H. C. J. Superparasitism and ovicide in parasitic Hymenoptera: Theory and a case study of the ectoparasitoid Bracon hebetor. Behav. Ecol. Sociobiol. 24, 421–432 (1989).Article 

Google Scholar 
47.Abram, P. K., Brodeur, J., Urbaneja, A. & Tena, A. Nonreproductive effects of insect parasitoids on their hosts. Annu. Rev. Entomol. 64, 259–276 (2019).CAS 
Article 

Google Scholar 
48.Abram, P. K., Brodeur, J., Burte, V. & Boivin, G. Parasitoid-induced host egg abortion: An underappreciated component of biological control services provided by egg parasitoids. Biol. Control 98, 52–60 (2016).Article 

Google Scholar 
49.Abram, P. K., Gariepy, T. D., Boivin, G. & Brodeur, J. An invasive stink bug as an evolutionary trap for an indigenous egg parasitoid. Biol. Invasions 16, 1387–1395 (2014).Article 

Google Scholar 
50.Steiner, A. L. Stinging behaviour of solitary wasps. In Venoms of the Hymenoptera. Biochemical, Pharmacological and Behavioural Aspects (ed Piek, T.) 63–148 (Academic Press, London, 1986) https://doi.org/10.1016/b978-0-12-554770-3.50008-5.51.Feng, Y., Wratten, S., Sandhu, H. & Keller, M. Interspecific competition between two generalist parasitoids that attack the leafroller Epiphyas postvittana (Lepidoptera: Tortricidae). Bull. Entomol. Res. 105, 426–433 (2015).CAS 
Article 

Google Scholar 
52.De Moraes, C. M. & Mescher, M. C. Intrinsic competition between larval parasitoids with different degrees of host specificity. Ecol. Entomol. 30, 564–570 (2005).Article 

Google Scholar 
53.Desneux, N., Barta, R. J., Hoelmer, K. A., Hopper, K. R. & Heimpel, G. E. Multifaceted determinants of host specificity in an aphid parasitoid. Oecologia 160, 387–398 (2009).ADS 
Article 

Google Scholar 
54.Brodeur, J. & Boivin, G. Functional ecology of immature parasitoids. Annu. Rev. Entomol. 49, 27–49 (2004).CAS 
Article 

Google Scholar 
55.Rogers, D. Random search and insect population models. J. Anim. Ecol. 41, 369–383 (1972).Article 

Google Scholar 
56.Holling, C. S. Some characteristics of simple types of predation and parasitism. Can. Entomol. 91, 385–398. https://doi.org/10.4039/Ent91385-7 (1959).Article 

Google Scholar  More

1.Zhao, F. et al. Diversity and species abundance patterns of the early Cambrian (Series 2, Stage 3) Chengjiang Biota from China. Paleobiology 40, 50–69 (2014).Article 

Google Scholar 
2.Zhu, M.-Y., Zhang, J.-M. & Li, G.-X. Sedimentary environments of the early Cambrian Chengjiang biota: sedimentology of the Yu’anshan Formation in Chengjiang County, eastern Yunnan. Acta Palaeontol. Sin. 40, 80–105 (2001).
Google Scholar 
3.Hu, S.-X. Taphonomy and palaeoecology of the early Cambrian Chengjiang Biota from eastern Yunnan, China. Berl. Palobiologische Abhandlungen 7 (2005).4.Hou, X. et al. The Cambrian Fossils of Chengjiang, China. The Flowering of Early Animal Life 2nd edn (John Wiley & Sons, 2017).5.Zhang, W.-T. & Hou, X.-G. Preliminary notes on the occurrence of the unusual trilobite Naraoia in Asia. Acta Palaeontol. Sin. 24, 591–595 (1985).
Google Scholar 
6.Luo, H.-L, Hu, S.-X, Chen, L.-Z, Zhang, S.-S & Tao, Y.-H. Early Cambrian Chengjiang Fauna from Kunming Region, China (Yunnan Science and Technology Press, 1999).7.Chen, J.-Y The Dawn of Animal World (Jiangsu Science and Technology Press, China, 2004).8.Duan, Y. et al. Reproductive strategy of the bradoriid arthropod Kunmingella douvillei from the lower Cambrian Chengjiang Lagerstätte, South China. Gondwana Res. 25, 983–990 (2014).Article 

Google Scholar 
9.Zhao, F.-C., Zhu, M.-Y. & Hu, S.-X. Community structure and composition of the Cambrian Chengjiang biota. Sci. China Earth Sci. 53, 1784–1799 (2010).Article 

Google Scholar 
10.Liu, Y. et al. Three-dimensionally preserved minute larva of a great-appendage arthropod from the early Cambrian Chengjiang biota. Proc. Natl Acad. Sci. USA 113, 5542–5546 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
11.Ou, Q. et al. Evolutionary trade-off in reproduction of Cambrian arthropods. Sci. Adv. 6, 33–76 (2020).
Google Scholar 
12.Dornbos, S. Q. & Chen, J.-Y. Community palaeoecology of the Early Cambrian Maotianshan Shale biota: ecological dominance of priapulid worms. Palaeogeogr. Palaeoclimatol. Palaeoecol. 258, 200–212 (2008).Article 

Google Scholar 
13.Fu, D. et al. The Qingjiang biota—a Burgess Shale-type fossil Lagerstätte from the early Cambrian of South China. Science 363, 1338–1342 (2019).CAS 
PubMed 
Article 

Google Scholar 
14.Caron, J.-B. & Jackson, D. A. Paleoecology of the Greater Phyllopod Bed community, Burgess Shale. Palaeogeogr. Palaeoclimatol. Palaeoecol. 258, 222–256 (2008).15.Nanglu, K., Caron, J.-B. & Gaines, R. R. The Burgess Shale paleocommunity with new insights from Marble Canyon, British Columbia. Paleobiology 46, 58–81 (2020).Article 

Google Scholar 
16.Gaines, R. R. in Reading and Writing of the Fossil Record: Preservational Pathways to Exceptional Fossilization Vol. 20 (eds Laflamme, M. et al.) 123–146 (Paleontological Research Institution, 2014).17.Zhai, D. et al. Variation in appendages in early Cambrian bradoriids reveals a wide range of body plans in stem-euarthropods. Commun. Biol. 2, 329 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
18.Isaevaa, V. V., Ozernyukc, N. D. & Rozhnov, S. V. Evidence for evolutionary changes in ontogeny: paleontological, comparative morphological, and molecular aspects. Biol. Bull. 40, 243–252 (2013).Article 

Google Scholar 
19.Liu, Y., Haug, J. T., Haug, C., Briggs, D. E. G. & Hou, X.-G. A 520 million-year-old chelicerate larva. Nat. Commun. 5, 4440 (2014).CAS 
PubMed 
Article 

Google Scholar 
20.Chipman, A. D. An embryological perspective on the early arthropod fossil record. BMC Evol. Biol. 15, 285 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
21.Wolfe, J. M. Metamorphosis is ancestral for crown euarthropods, and evolved in the Cambrian or earlier. Integr. Comp. Biol. 57, 499–509 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
22.Haug, T. J. Why the term “larva” is ambiguous, or what makes a larva? Acta Zool. 101, 167–188 (2018).Article 

Google Scholar 
23.Fu, D., Zhang, X., Budd, G. E., Liu, W. & Pan, X. Ontogeny and dimorphism of Isoxys auritus (Arthropoda) from the Early Cambrian Chengjiang biota, South China. Gondwana Res. 25, 975–982 (2014).Article 

Google Scholar 
24.Yang, X.-F., Kimmig, J., Lieberman, B. S. & Peng, S.-C. A new species of the deuterostome Herpetogaster from the early Cambrian Chengjiang biota of South China. Sci. Nat. 107, 37 (2020).CAS 
Article 

Google Scholar 
25.Zhai, D. Y. et al. Fine-scale appendage structure of the Cambrian trilobitomorph Naraoia spinosa and its ontogenetic and ecological implications. Proc. R. Soc. B 286, 20192371 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
26.Hughes, N. C. et al. Articulated trilobite ontogeny: suggestions for a methodological standard. J. Paleont. 95, 298–304 (2021).Article 

Google Scholar 
27.Chen, J.-Y. & Zhou, G.-Q. Biology of the Chengjiang fauna. Bull. Natl Mus. Nat. Sci. 10, 11–106 (1997).
Google Scholar 
28.Haug, J. T., Caron, J.-B. & Haug, C. Demecology in the Cambrian: synchronized moulting in arthropods from the Burgess Shale. BMC Biol. 11, 64 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
29.Robison, R. A., Babcock, L. E. & Gunther, V. G. Exceptional Cambrian fossils from Utah: A Window into the Age of Trilobites (Utah Geological Survey, 2015).30.Kimmig, J., Strotz, L. C., Kimmig, S. R., Egenhoff, S. O. & Lieberman, B. S. The Spence Shale Lagerstätte: an important window into Cambrian biodiversity. J. Geol. Soc. Lond. 176, 609–619 (2019).Article 

Google Scholar 
31.Paterson, J. R. et al. The Emu Bay Shale Konservat-Lagerstätte: a view of Cambrian life from East Gondwana. J. Geol. Soc. Lond. 173, 3107 (2016).32.Du, K. et al. A new early Cambrian Konservat-Lagerstätte expands the occurrence of Burgess Shale-type deposits on the Yangtze Platform. Earth Sci. Rev. 211, 103409 (2020).Article 

Google Scholar 
33.Harper, D. A. T. et al. The Sirius Passet Lagerstätte of North Greenland: a remote window on the Cambrian explosion. J. Geol. Soc. Lond. 176, 1023–1037 (2019).Article 

Google Scholar 
34.Chen, L. Z et al. Early Cambrian Chengjiang Fauna in Eastern Yunnan, China (Yunnan Science and Technology Press, 2002).35.Zhao, F. C., Caron, J.-B., Hu, S. X. & Zhu, M. Y. Quantitative analysis of taphofacies and paleocommunities in the Early Cambrian Chengjiang Lagerstätte. PALAIOS 24, 826–839 (2009).CAS 
Article 

Google Scholar 
36.Beck, M. K. et al. The identification, conservation, and management of estuarine and marine nurseries for fish and invertebrates. BioScience 51, 633–641 (2001).Article 

Google Scholar 
37.Botton, M. L. & Loveland, R. E. Abundance and dispersal potential of horseshoe crab (Limulus polyphemus) larvae in the Delaware estuary. Estuar. Coasts 26, 1472–1479 (2003).Article 

Google Scholar 
38.Watson, W. H. & Chabot, C. C. High resolution tracking of adult horseshoe crabs Limulus polyphemus in a New Hampshire estuary using a fixed array ultrasonic telemetry. Curr. Zool. 56, 599–610 (2010).Article 

Google Scholar 
39.Perry, F. A. et al. Habitat partitioning in Antarctic krill: spawning hotspots and nursery areas. PLoS ONE 14, e0219325 (2019).40.Nagelkerken, I. in Ecological Connectivity among Tropical Coastal Ecosystems (ed. Nagelkerken, I.) 357–399 (Springer, 2009).41.Kanciruk, P. in The Biology and Management of Lobsters Vol. 2 (eds Cobb, J. S. & Phillips, B. F.) 59–96 (Academic Press, 1980).42.Sandt, V. J. & Stoner, A. W. Ontogenetic shift in habitat by early juvenile queen conch, Strombus gigas: patterns and potential mechanisms. Fish. Bull. 91, 516–525 (1993).
Google Scholar 
43.Pedrotti, M. L. & Fenaux, L. Dispersal of echinoderm larvae in a geographical area marked by upwelling (Ligurian Sea, NW Mediterranean). Mar. Ecol. Prog. Ser. 87, 217–227 (1992).Article 

Google Scholar 
44.Zhai, D. et al. Spatial heterogeneity of the population age structure of the ostracode Limnocythere inopinata in Hulun Lake, Inner Mongolia and its implications. Hydrobiologia 716, 29–46 (2013).CAS 
Article 

Google Scholar 
45.Baillon, S., Hamel, J. F., Wareham, V. E. & Mercier, A. Deep cold-water corals as nurseries for fish larvae. Front. Ecol. Environ. 10, 351–356 (2012).Article 

Google Scholar 
46.Treude, T., Kiel, S., Linke, P., Peckmann, J. & Goedert, J. Elasmobranch egg capsules associated with modern and ancient cold seeps: a nursery for marine deep-water predators. Mar. Ecol. Prog. Ser. 437, 175–181 (2011).Article 

Google Scholar 
47.Rooper, C. N., Boldt, J. L. & Zimmermann, M. An assessment of juvenile Pacific Ocean perch (Sebastes alutus) habitat use in a deepwater nursery. Estuar. Coast. Shelf Sci. 75, 371–380 (2007).Article 

Google Scholar 
48.Pimiento, C., Ehret, D. J., MacFadden, B. J. & Hubbell, G. Ancient nursery area for the extinct giant shark Megalodon from the Miocene of Panama. PLoS ONE 5, e10552 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
49.Villafaña, J. A. et al. First evidence of a palaeo-nursery area of the great white shark. Sci. Rep. 10, 8502 (2020).50.Paterson, J. R., Jago, J. B., Brock, G. A. & Gehling, J. G. Taphonomy and palaeoecology of the emuellid trilobite Balcoracania dailyi (early Cambrian, South Australia). Palaeogeogr. Palaeoclimatol. Palaeoecol. 249, 302–321 (2007).Article 

Google Scholar 
51.Hartnoll, R. G. in Physiology and Behaviour of Marine Organisms (eds McLusky, D. S. & Berry, A. J.) 349–358 (Pergamon Press, 1978).52.Hartnoll, R. G. & Bryant, A. D. Size-frequency distributions in decapod Crustacea—the quick, the dead and the cast-offs. J. Crust. Biol. 10, 14–19 (1990).Article 

Google Scholar 
53.Sheldon, P. R. Trilobite size-frequency distributions, recognition of instars, and phyletic size changes. Lethaia 21, 293–306 (1988).Article 

Google Scholar 
54.Herrnkind, W. F. in The Biology and Management of Lobsters Vol. 1 (eds Cobb, J. S. & Phillips B. F.) 349–407 (Academic Press, 1980)55.Linnane, A., Dimmlich, W. & Ward, T. Movement patterns of the southern rock lobster, Jasus edwardsii, of South Australia. NZ J. Mar. Freshw. Res. 39, 335–346 (2005).Article 

Google Scholar 
56.Blazejowski, B. et al. Ancient animal migration: a case study of eyeless, dimorphic Devonian trilobites from Poland. Palaeontology 59, 743–751 (2016).Article 

Google Scholar 
57.Hughes, N. C., Kříž, J., Macquaker, J. H. S. & Huff, W. D. The depositional environment and taphonomy of the Homerian “Aulacopleura shales” fossil assemblage near Loděnice, Czech Republic (Prague Basin, Perunican microcontinent). Bull. Geosci. 89, 219–238 (2014).Article 

Google Scholar 
58.Whitaker, A. F. & Kimmig, J. Anthropologically introduced biases in natural history collections, with a case study on the invertebrate paleontology collections from the middle Cambrian Spence Shale Lagerstätte. Palaeontol. Electron. 23, a58 (2020).
Google Scholar 
59.Conway Morris, S. The community structure of the Middle Cambrian phyllopod bed (Burgess Shale). Palaeontology 29, 423–467 (1986).
Google Scholar 
60.Caron, J.-B., Gaines, R. R., Aria, C., Mángano, M. G. & Streng, M. A new phyllopod bed-like assemblage on from the Burgess Shale of the Canadian Rockies. Nat. Commun. 5, 3210 (2014).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
61.Ihaka, R. R. & Gentleman, R. A language for data analysis and graphics. J. Comput. Graph. Stat. 5, 299–314 (1996).
Google Scholar  More

1.Boudreau BP, Huettel M, Forster S, Jahnke RA, McLachlan A, Middelburg JJ, et al. Permeable marine sediments: overturning an old paradigm. Eos Trans AGU. 2001;82:133–6.
Google Scholar 
2.Huettel M, Berg P, Kostka JE. Benthic exchange and biogeochemical cycling in permeable sediments. Annu Rev Mar Sci. 2014;6:23–51.Article 

Google Scholar 
3.Huettel M, Ziebis W, Forster S. Flow-induced uptake of particulate matter in permeable sediments. Limnol Oceanogr. 1996;41:309–22.Article 

Google Scholar 
4.Huettel M, Rusch A. Transport and degradation of phytoplankton in permeable sediment. Limnol Oceanogr. 2000;45:534–49.CAS 
Article 

Google Scholar 
5.Rusch A, Forster S, Huettel M. Bacteria, diatoms and detritus in an intertidal sandflat subject to advective transport across the water-sediment interface. Biogeochemistry. 2001;55:1–27.CAS 
Article 

Google Scholar 
6.Ahmerkamp S, Winter C, Krämer K, de Beer D, Janssen F, Friedrich J, et al. Regulation of benthic oxygen fluxes in permeable sediments of the coastal ocean. Limnol Oceanogr. 2017;62:1935–54.CAS 
Article 

Google Scholar 
7.Jahnke RA Global Synthesis. In: Liu KK, Atkinson L, Quinones R, Talaue-McManus L, editors. Carbon and nutrient fluxes in continental margins. Ch. 16 Berlin: Springer; 2010.8.Joiris C, Billen G, Lancelot C, Daro MH, Mommaerts JP, Bertels A, et al. A budget of carbon cycling in the Belgian coastal zone: relative roles of zooplankton, bacterioplankton and benthos in the utilization of primary production. Neth. J. Sea Res. 1982;16:260–75.CAS 
Article 

Google Scholar 
9.Jørgensen BB, Bang M, Blackburn TH. Anaerobic mineralization in marine-sediments from the Baltic-Sea-North Sea transition. Mar Ecol Prog Ser. 1990;59:39–54.Article 

Google Scholar 
10.Middelburg JJ, Barranguet C, Boschker HTS, Herman PMJ, Moens T, Heip CHR. The fate of intertidal microphytobenthos carbon: an in situ 13C-labeling study. Limnol Oceanogr. 2000;45:1224–34.CAS 
Article 

Google Scholar 
11.Böer SI, Arnosti C, van Beusekom JEE, Boetius A. Temporal variations in microbial activities and carbon turnover in subtidal sandy sediments. Biogeosciences. 2009;6:1149–65.Article 

Google Scholar 
12.Goto N, Mitamura O, Terai H. Biodegradation of photosynthetically produced extracellular organic carbon from intertidal benthic algae. J Exp Mar Biol Ecol. 2001;257:73–86.CAS 
PubMed 
Article 

Google Scholar 
13.Rusch A, Huettel M, Reimers CE, Taghon GL, Fuller CM. Activity and distribution of bacterial populations in Middle Atlantic Bight shelf sands. FEMS Microb Ecol. 2003;44:89–100.CAS 
Article 

Google Scholar 
14.Hewson I, Vargo GA, Fuhrman JA. Bacterial diversity in shallow oligotrophic marine benthos and overlying waters: effects of virus infection, containment, and nutrient enrichment. Microb Ecol. 2003;46:322–36.CAS 
PubMed 
Article 

Google Scholar 
15.Teske A, Durbin A, Ziervogel K, Cox C, Arnosti C. Microbial community composition and function in permanently cold seawater and sediments from an Arctic fjord of Svalbard. Appl Environ Microbiol. 2011;77:2008–18.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Zinger L, Amaral-Zettler LA, Fuhrman JA, Horner-Devine MC, Huse SM, Welch DBM, et al. Global patterns of bacterial beta-diversity in seafloor and seawater ecosystems. PLoS ONE. 2011;6:e24570.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Cardman Z, Arnosti C, Durbin A, Ziervogel K, Cox C, Steen AD, et al. Verrucomicrobia are candidates for polysaccharide-degrading bacterioplankton in an Arctic fjord of Svalbard. Appl Environ Microbiol. 2014;80:3749–56.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Teeling H, Fuchs BM, Becher D, Klockow C, Gardebrecht A, Bennke CM, et al. Substrate-controlled succession of marine bacterioplankton populations induced by a phytoplankton bloom. Science. 2012;336:608–11.CAS 
PubMed 
Article 

Google Scholar 
19.Teeling H, Fuchs BM, Bennke CM, Kruger K, Chafee M, Kappelmann L, et al. Recurring patterns in bacterioplankton dynamics during coastal spring algae blooms. eLife. 2016;5:e11888.PubMed 
PubMed Central 
Article 

Google Scholar 
20.Fuhrman JA, Hewson I, Schwalbach MS, Steele JA, Brown MV, Naeem S. Annually reoccurring bacterial communities are predictable from ocean conditions. Proc Natl Acad Sci USA. 2006;103:13104–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
21.Chafee M, Fernàndez-Guerra A, Buttigieg PL, Gerdts G, Eren AM, Teeling H, et al. Recurrent patterns of microdiversity in a temperate coastal marine environment. ISME J. 2018;12:237–52.PubMed 
Article 

Google Scholar 
22.Mayer LM. Extracellular proteolytic enzyme activity in sediments of an intertidal mudflat. Limnol Oceanogr. 1989;34:973–81.CAS 
Article 

Google Scholar 
23.Middelburg J, Klaver G, Nieuwenhuize J, Wielemaker A, Haas W, Vlug T, et al. Organic matter mineralization in intertidal sediment along an estuarine gradient. Mar Ecol Prog Ser. 1996;132:157–68.24.Tabuchi K, Kojima H, Fukui M. Seasonal changes in organic matter mineralization in a sublittoral sediment and temperature-driven decoupling of key processes. Microb Ecol. 2010;60:551–60.PubMed 
Article 

Google Scholar 
25.Hoffmann K, Hassenrück C, Salman-Carvalho V, Holtappels M, Bienhold C. Response of bacterial communities to different detritus compositions in Arctic deep-sea sediments. Front Microbiol. 2017;8:266.PubMed 
PubMed Central 

Google Scholar 
26.Gobet A, Boer SI, Huse SM, van Beusekom JEE, Quince C, Sogin ML, et al. Diversity and dynamics of rare and of resident bacterial populations in coastal sands. ISME J. 2012;6:542–53.PubMed 
Article 

Google Scholar 
27.Mills HJ, Hunter E, Humphrys M, Kerkhof L, McGuinness L, Huettel M, et al. Characterization of nitrifying, denitrifying, and overall bacterial communities in permeable marine sediments of the northeastern Gulf of Mexico. Appl Environ Microbiol. 2008;74:4440–53.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
28.Probandt D, Knittel K, Tegetmeyer HE, Ahmerkamp S, Holtappels M, Amann R. Permeability shapes bacterial communities in sublittoral surface sediments. Environ Microbiol. 2017;19:1584–99.CAS 
PubMed 
Article 

Google Scholar 
29.Tait K, Airs RL, Widdicombe CE, Tarran GA, Jones MR, Widdicombe S. Dynamic responses of the benthic bacterial community at the Western English Channel observatory site L4 are driven by deposition of fresh phytodetritus. Prog Oceanogr. 2015;137:546–58.Article 

Google Scholar 
30.Wiltshire K, Kraberg A, Bartsch I, Boersma M, Franke H-D, Freund J, et al. Helgoland Roads, North Sea: 45 years of change. Estuaries and Coasts. 2010;33:295–310.CAS 
Article 

Google Scholar 
31.Probandt D. Microbial ecology of subtidal sandy sediments [PhD thesis]. Bremen: University of Bremen; 2017.32.Berge J, Renaud PE, Darnis G, Cottier F, Last K, Gabrielsen TM, et al. In the dark: a review of ecosystem processes during the Arctic polar night. Prog Oceanogr. 2015;139:258–71.Article 

Google Scholar 
33.Boehnert S, Ruiz Soto S, Fox BRS, Yokoyama Y, Hebbeln D. Historic development of heavy metal contamination into the Firth of Thames, New Zealand. Geo-Mar Lett. 2020;40:149–65.CAS 
Article 

Google Scholar 
34.Lorenzen CJ. Determination of chlorophyll and pheo-pigments: spectrophotometric eqations. Limnol Oceanogr. 1967;12:343–6.CAS 
Article 

Google Scholar 
35.Zhou J, Bruns MA, Tiedje JM. DNA recovery from soils of diverse composition. Appl Environ Microbiol. 1996;62:316–22.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Herlemann DPR, Labrenz M, Jürgens K, Bertilsson S, Waniek JJ, Andersson AF. Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. ISME J. 2011;5:1571–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Bushnell B, Rood J, Singer E. BBMerge—accurate paired shotgun read merging via overlap. PLoS ONE. 2017;12:e0185056.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
38.Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537–41.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–D596.CAS 
PubMed 
Article 

Google Scholar 
40.Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 2017;11:2639–43.PubMed 
PubMed Central 
Article 

Google Scholar 
41.Oksanen J, Blanchet F, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: Community Ecology Package. R package version. 2019;2:5–6.
Google Scholar 
42.Team R.C. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.r-project.org/; 2019.43.Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, et al. Welcome to the Tidyverse. J Open Source Softw. 2019;4:1686.Article 

Google Scholar 
44.Chapman MG, Underwood AJ. Ecological patterns in multivariate assemblages: information and interpretation of negative values in ANOSIM tests. Mar Ecol Prog Ser. 1999;180:257–65.Article 

Google Scholar 
45.Pernthaler A, Pernthaler J, Amann R. Fluorescence in situ hybridization and catalyzed reporter deposition for the identification of marine bacteria. Appl Environ Microbiol. 2002;68:3094–101.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Pernthaler J, Pernthaler A, Amann R. Automated enumeration of groups of marine picoplankton after fluorescence in situ hybridization. Appl Environ Microbiol. 2003;69:2631–7.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Bennke CM, Reintjes G, Schattenhofer M, Ellrott A, Wulf J, Zeder M, et al. Modification of a high-throughput automatic microbial cell enumeration system for shipboard analyses. Appl Environ Microbiol. 2016;82:3289–96.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
48.Ludwig W, Strunk O, Westram R, Richter L, Meier H, Yadhukumar, et al. ARB: a software environment for sequence data. Nucleic Acids Res. 2004;32:1363–71.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Snaidr J, Amann R, Huber I, Ludwig W, Schleifer K, Snaidr J, et al. Phylogenetic analysis and in situ identification of bacteria in activated sludge. Appl Environ Microbiol. 1997;63:2884–96.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Bockelmann F-D, Puls W, Kleeberg U, Müller D, Emeis K-C. Mapping mud content and median grain-size of North Sea sediments—a geostatistical approach. Mar Geol. 2018;397:60–71.Article 

Google Scholar 
51.Hoshino T, Doi H, Uramoto G-I, Wörmer L, Adhikari RR, Xiao N, et al. Global diversity of microbial communities in marine sediment. Proc Natl Acad Sci USA. 2020;117:27587–97.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Probandt D, Eickhorst T, Ellrott A, Amann R, Knittel K. Microbial life on a sand grain: from bulk sediment to single grains. ISME J. 2017;12:623.PubMed 
PubMed Central 
Article 

Google Scholar 
53.Acosta-González A, Rosselló-Móra R, Marqués S. Characterization of the anaerobic microbial community in oil-polluted subtidal sediments: aromatic biodegradation potential after the Prestige oil spill. Environ Microbiol. 2013;15:77–92.PubMed 
Article 
CAS 

Google Scholar 
54.Tian F, Yu Y, Chen B, Li H, Yao Y-F, Guo X-K. Bacterial, archaeal and eukaryotic diversity in Arctic sediment as revealed by 16S rRNA and 18S rRNA gene clone libraries analysis. Polar Biol. 2009;32:93–103.Article 

Google Scholar 
55.Zeng Y, Zou Y, Grebmeier JM, He J, Zheng T. Culture-independent and culture-dependent methods to investigate the diversity of planktonic bacteria in the northern Bering Sea. Polar Biol. 2012;35:117–29.Article 

Google Scholar 
56.Santelli CM, Orcutt BN, Banning E, Bach W, Moyer CL, Sogin ML, et al. Abundance and diversity of microbial life in ocean crust. Nature. 2008;453:653–6.CAS 
PubMed 
Article 

Google Scholar 
57.Ravenschlag K, Sahm K, Pernthaler J, Amann R. High bacterial diversity in permanently cold marine sediments. Appl Environ Microbiol. 1999;65:3982–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
58.Hunter EM, Mills HJ, Kostka JE. Microbial community diversity associated with carbon and nitrogen cycling in permeable shelf sediments. Appl Environ Microbiol. 2006;72:5689–701.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
59.Dyksma S, Bischof K, Fuchs BM, Hoffmann K, Meier D, Meyerdierks A, et al. Ubiquitous Gammaproteobacteria dominate dark carbon fixation in coastal sediments. ISME J. 2016;10:1939–53.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Allers E, Wright JJ, Konwar KM, Howes CG, Beneze E, Hallam SJ, et al. Diversity and population structure of Marine Group A bacteria in the Northeast subarctic Pacific Ocean. ISME J. 2013;7:256–68.CAS 
PubMed 
Article 

Google Scholar 
61.Hodal H, Falk-Petersen S, Hop H, Kristiansen S, Reigstad M. Spring bloom dynamics in Kongsfjorden, Svalbard: nutrients, phytoplankton, protozoans and primary production. Polar Biol. 2012;35:191–203.Article 

Google Scholar 
62.Jönsson BF, Salisbury JE, Mahadevan A. Large variability in continental shelf production of phytoplankton carbon revealed by satellite. Biogeosciences. 2011;8:1213–23.Article 
CAS 

Google Scholar 
63.Kuliński K, Kędra M, Legeżyńska J, Gluchowska M, Zaborska A. Particulate organic matter sinks and sources in high Arctic fjord. J Mar Syst. 2014;139:27–37.Article 

Google Scholar 
64.Bourgeois S, Kerhervé P, Calleja ML, Many G, Morata N. Glacier inputs influence organic matter composition and prokaryotic distribution in a high Arctic fjord (Kongsfjorden, Svalbard). J Mar Syst. 2016;164:112–27.Article 

Google Scholar 
65.Zaborska A, Włodarska-Kowalczuk M, Legeżyńska J, Jankowska E, Winogradow A, Deja K. Sedimentary organic matter sources, benthic consumption and burial in west Spitsbergen fjords—signs of maturing of Arctic fjordic systems? J Mar Syst. 2018;180:112–23.Article 

Google Scholar 
66.McGovern M, Pavlov AK, Deininger A, Granskog MA, Leu E, Søreide JE, et al. Terrestrial inputs drive seasonality in organic matter and nutrient biogeochemistry in a high Arctic fjord system (Isfjorden, Svalbard). Front Mar Sci. 2020;7:747.Article 

Google Scholar 
67.Avci B, Krüger K, Fuchs BM, Teeling H, Amann RI. Polysaccharide niche partitioning of distinct Polaribacter clades during North Sea spring algal blooms. ISME J. 2020;14:1369–83.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Braeckman U, Janssen F, Lavik G, Elvert M, Marchant H, Buckner C, et al. Carbon and nitrogen turnover in the Arctic deep sea: in situ benthic community response to diatom and coccolithophorid phytodetritus. Biogeosciences. 2018;15:6537–57.CAS 
Article 

Google Scholar 
69.Guilini K, Oevelen DV, Soetaert K, Middelburg JJ, Vanreusela A. Nutritional importance of benthic bacteria for deep-sea nematodes from the Arctic ice margin: results of an isotope tracer experi5ment. Limnol Oceanogr. 2010;55:1977–89.CAS 
Article 

Google Scholar 
70.van Oevelen D, Soetaert K, Middelburg J, Herman P, Moodley L, Hamels I, et al. Carbon flows through a benthic food web: Integrating biomass, isotope and tracer data. J Mar Res. 2006;64:453–82.Article 

Google Scholar 
71.Danovaro R, Dell’Anno A, Corinaldesi C, Magagnini M, Noble R, Tamburini C. et al. Major viral impact on the functioning of benthic deep-sea ecosystems. Nature. 2008;454:1084–7.CAS 
PubMed 
Article 

Google Scholar 
72.Miller DC. Abrasion effects on microbes in sandy sediments. Mar Ecol Prog Ser. 1989;55:73–82.Article 

Google Scholar 
73.Ahmerkamp S, Marchant HK, Peng C, Probandt D, Littmann S, Kuypers MM. et al. The effect of sediment grain properties and porewater flow on microbial abundance and respiration in permeable sediments. Sci. Rep. 2020;10:3573CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
74.Barka EA, Vatsa P, Sanchez L, Gaveau-Vaillant N, Jacquard C, Klenk HP. et al. Taxonomy, physiology, and natural products of Actinobacteria. Microbiol Mol Biol Rev. 2016;80:1–43.PubMed 
Article 

Google Scholar 
75.Schrempf H. Actinobacteria within soils: capacities for mutualism, symbiosis and pathogenesis. FEMS Microbiol Lett. 2013;342:77–78.CAS 
PubMed 
Article 

Google Scholar 
76.Giovannoni SJ, Stingl U. Molecular diversity and ecology of microbial plankton. Nature. 2005;437:343–8.CAS 
PubMed 
Article 

Google Scholar 
77.Yilmaz P, Iversen MH, Hankeln W, Kottmann R, Quast C, Glöckner FO. Ecological structuring of bacterial and archaeal taxa in surface ocean waters. FEMS Microbiol Ecol. 2012;81:373–85.CAS 
PubMed 
Article 

Google Scholar 
78.Bienhold C, Zinger L, Boetius A, Ramette A. Diversity and biogeography of bathyal and abyssal seafloor bacteria. PLoS ONE. 2016;11:e0148016.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
79.Rappé MS, Kemp PF, Giovannoni SJ. Phylogenetic diversity of marine coastal picoplankton 16S rRNA genes cloned from the continental shelf off Cape Hatteras, North Carolina. Limnol Oceanogr. 1997;42:811–26.Article 

Google Scholar 
80.Zeng Y-X, Yu Y, Li H-R, Luo W. Prokaryotic community composition in Arctic Kongsfjorden and sub-arctic northern Bering Sea sediments as revealed by 454 pyrosequencing. Front Microbiol. 2017;8:2498.PubMed 
PubMed Central 
Article 

Google Scholar 
81.Fang X-M, Zhang T, Li J, Wang NF, Wang Z, Yu LY. Bacterial community pattern along the sediment seafloor of the Arctic fjorden (Kongsfjorden, Svalbard). Antonie Van Leeuwenhoek. 2019;112:1121–36.PubMed 
Article 

Google Scholar 
82.Ziemert N, Lechner A, Wietz M, Millán-Aguiñaga N, Chavarria KL, Jensen PR. et al. Diversity and evolution of secondary metabolism in the marine actinomycete genus salinispora. Proc Natl Acad Sci USA. 2014;111:e1130–1139.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
83.Manivasagan P, Venkatesan J, Sivakumar K, Kim SK. Pharmaceutically active secondary metabolites of marine actinobacteria. Microbiol Res. 2014;169:262–78.CAS 
PubMed 
Article 

Google Scholar 
84.Kamjam M, Sivalingam P, Deng Z, Hong K. Deep sea Actinomycetes and their secondary metabolites. Front Microbiol. 2017;8:760.PubMed 
PubMed Central 
Article 

Google Scholar 
85.Lewin GR, Carlos C, Chevrette MG, Horn HA, McDonald BR, Stankey RJ. et al. Evolution and ecology of Actinobacteria and their bioenergy applications. Annu Rev Microbiol. 2016;70:235–54.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
86.Matsumoto A, Kasai H, Matsuo Y, Ōmura S, Shizuri Y, Takahashi Y. Ilumatobacter fluminis gen. nov., sp. nov., a novel actinobacterium isolated from the sediment of an estuary. J Gen Appl Microbiol. 2009;55:201–5.CAS 
PubMed 
Article 

Google Scholar 
87.Ghai R, Mizuno CM, Picazo A, Camacho A, Rodriguez-Valera F. Metagenomics uncovers a new group of low GC and ultra-small marine Actinobacteria. Sci Rep. 2013;3:2471.PubMed 
PubMed Central 
Article 

Google Scholar 
88.El Kaoutari A, Armougom F, Gordon J, Raoult D, Henrissat B. The abundance and variety of carbohydrate-active enzymes in the human gut microbiota. Nat Rev Microbiol. 2013;11:497–504.89.Berlemont R, Martiny AC. Glycoside hydrolases across environmental microbial communities. PLoS Comp. Biol. 2016;12:e1005300.Article 
CAS 

Google Scholar 
90.Becker S, Tebben J, Coffinet S, Wiltshire K, Iversen MH, Harder T, et al. Laminarin is a major molecule in the marine carbon cycle. Proc Natl Acad Sci USA. 2020;117:6599–607.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
91.Coutinho MCL, Teixeira VL, Santos CSG. A review of “Polychaeta” chemicals and their possible ecological role. J Chem Ecol. 2018;44:72–94.CAS 
PubMed 
Article 

Google Scholar 
92.Arnosti C. Functional differences between Arctic seawater and sedimentary microbial communities: contrasts in microbial hydrolysis of complex substrates. FEMS Microbiol Ecol. 2008;66:343–51.CAS 
PubMed 
Article 

Google Scholar 
93.Krüger K, Chafee M, Francis TB, Del Rio TG, Becher D, Schweder T, et al. In marine Bacteroidetes the bulk of glycan degradation during algae blooms is mediated by few clades using a restricted set of genes. ISME J. 2019;13:2800–16.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
94.Reintjes G, Arnosti C, Fuchs BM, Amann R. An alternative polysaccharide uptake mechanism of marine bacteria. ISME J. 2017;11:1640–50.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
95.Arnosti C, Jørgensen BB. High activity and low temperature optima of extracellular enzymes in Arctic sediments: implications for carbon cycling by heterotrophic microbial communities. Mar Ecol Prog Ser. 2003;249:15–24.CAS 
Article 

Google Scholar 
96.Arnosti C, Jørgensen BB. Organic carbon degradation in Arctic marine sediments, Svalbard: a comparison of initial and terminal steps. Geomicrobiol J. 2006;23:551–63.CAS 
Article 

Google Scholar  More

Data sourcesNo single comprehensive dataset of planktonic foraminiferal distributional records currently exists. Instead, these data are available from a wide range of sources in many different structures. Some of these sources are compilations of existing data (e.g., Neptune14,15,16, ForCenS21), and others derive from individual sampling sites (e.g. ocean drilling expeditions). Triton combines these disparate sources (Fig. 1) to produce a single spatio-temporal dataset of Cenozoic planktonic foraminifera with updated and consistent taxonomy, age models, and paleo-coordinates.Neptune is currently the most comprehensive database of fossil plankton data, with records exclusively from the DSDP, ODP and IODP representing planktonic foraminifera, calcareous nannofossils, diatoms, radiolaria and dinoflagellates14,15,16. A subset of these sites is included in Neptune, representing those with the most continuous sampling through time. The raw data from Neptune form the core of our dataset. All foraminiferal occurrences for the Cenozoic (i.e. last 66 Ma) were downloaded using the GTS 2012 timescale. In the download options, all questionable identifications and invalid taxa were removed, as were records that had been identified as reworked.In addition to Neptune, three other compilation datasets were included in Triton: ForCenS21, which consists of global core-top samples; the Eocene data from Fenton, et al.8 created based on literature searches for planktonic foraminiferal data in the Eocene; and the land-based records from Lloyd, et al.22 that were created from literature searches. The marine records in Lloyd, et al.22 were not included, as they were obtained from Neptune.Following preliminary compilation of existing datasets, we identified all legacy DSDP, ODP and IODP cores missing from Triton. The online DESCLogik (http://web.iodp.tamu.edu/DESCReport/) and Pangaea17 databases were then mined for .csv files containing planktonic foraminiferal species count data for the missing cores, supplemented with data from AWI_Paleo (URI: http://www.awi.de/en/science/geosciences/marine-geology.html), GIK/IFG (URI: http://www.ifg.uni-kiel.de/), MARUM (URI: https://www.marum.de/index.html), and QUEEN (URI: http://ipt.vliz.be/eurobis/resource?r=pangaea_2747). All additional cores were assessed individually by inspecting the scientific drilling proceedings to determine whether sites were suitable to contribute to our dataset. The primary assessment criterion was identification of continuous sedimentary sections, wherein two or more confidently assigned consecutive chronostratigraphic tie points existed to allow for construction of age models.In addition to these longer cores, many sediment sampling projects have produced planktonic foraminiferal distribution data from shorter cores that tend to correspond in age to the last few million years. The website PANGAEA17 (www.pangaea.de) has been used as a repository for most of these occurrence data. This website was searched using the terms “plank* AND foram”, with resulting datasets downloaded using the R package ‘pangaear’23. These datasets were filtered to exclude records collected using multinets, sediment traps or box cores, as these methods produce samples not easily correlated to sediment cores. Column names allowed for further filtering to exclude records with no species-level data, records that had only isotopic data (rather than abundance data), or records with no age controls.Data processingThe data sources underpinning Triton serve their records in different formats. Therefore, processing was necessary to convert records into a unified framework, with one species per row for each sample and associated metadata (see below for details). Some metadata could be used without modification when available (e.g. water depth, data source), whereas other data needed processing to ensure consistency (e.g. abundance, paleo-coordinates, age). Without this processing, samples from different sources were not directly comparable. Where data were not available, they were set to NA. Those records with missing data in crucial columns (species name, abundance, age, and paleo-coordinates) were removed from the final dataset. All data processing was performed using R v. 3.6.124.Taxonomic consistency is essential to enable comparison of datasets created at different times. The species and synonymy lists used in Triton are based on the Paleogene Atlases20,25,26, with additional information from mikrotax27 (http://www.mikrotax.org/pforams/). These sources were supplemented, when necessary with more up to date literature including Poole and Wade28 and Lam and Leckie29. (A full list of the taxonomic sources can be found in the PFdata.xlsx file18.) A synonymy list was generated to convert species names to the senior synonym. At the same time, typographic errors were corrected. For example, Globototalia flexuosa should be Globorotalia flexuosa. Exclusively Mesozoic taxa were omitted, as were all instances when species names were unclear or imprecise (i.e. not at the species level). Junior synonyms were merged with their senior synonyms and their abundances summed, although the original names and abundances are also retained in the processed dataset. For presence/absence samples, these numerical merged abundances were set to one (i.e. present). The full species list and list of synonyms can be found in the accompanying data.Abundance data for planktonic foraminifera are provided in different formats: presence/absence, binned abundance, relative abundance, species counts, and number of specimens per gram. These metrics were converted into numeric relative abundance to make comparisons easier, although both the original abundance value and its numeric version are retained, as is a record of the abundance type. Presence/absence data were converted to a binary format (one for present; zero for absent). Species counts were converted to relative percent abundances based on the total number of specimens in the sample (this was calculated where it was not already recorded). When full counts were not performed, binned abundances were frequently used. These binned abundances were converted into numeric abundances based on the sequence. So, for example, the categorical labels of N, P, R, F, C, A, D (indicating none, present, rare, few, common, abundant, dominant) were converted to a numerical sequence of 0 to 6. As the meaning of letters can depend on the context (e.g. ‘A’ could be absent or abundant), conversion was done in a semi-automated fashion on a sample-by-sample basis. A value of 0.01 was assigned to records where an inconsistent abundance was recorded (e.g. samples with mostly numeric counts but a few species were designated ‘P’, indicating presence). Samples with zero abundance were retained in the full dataset to provide an indication of sampling.The age of samples were recorded in multiple ways. For some samples, age models provide precise numerical estimates of the age (e.g., those in Neptune). Other samples are dated relative to stratigraphic events such as biostratigraphic zones (including benthic and planktonic foraminifera, diatoms, radiolarians and nannofossils) or magnetic reversals. In this case, ages sometimes needed to be converted to reflect revised age estimates. The start and end dates of biostratigraphic zones are defined in relation to events in marker species, e.g. their speciation, extinction or acme events. All such marker events were updated to their most recent estimates and tuned to the GTS 2020 timescale19. The process of updating included correction of synonymies. Additional care was taken to ensure the correct interpretation of abbreviations (e.g. determining whether LO meant lowest occurrence or last occurrence) based on the entire list of events for a study. Where up-to-date ages were not available or events were ambiguous, they were removed from the age models.The marker events defining a zone can depend on the zonal scheme used. For example, Berggren30 defined the base of the planktonic foraminifera zone M8 as the first occurrence of Fohsella fohsi. Wade et al.31 used this same event to define the base of M9. Therefore, the zonal scheme was recorded when collecting age models, to accurately convert ages to the GTS 2020 time scale. Some marker events have different ages depending on the ocean basin or latitude, and these differences are not necessarily well studied31,32. Where these differences in marker events have been recorded, the coordinates of a site were used to determine whether sites were in the Atlantic or Indo-Pacific Ocean, and whether they were tropical or temperate (with the division at 23.5° latitude). However, this is an area where more research is needed to improve the accuracy of higher-latitude dating32. Magnetostratigraphic ages were also tuned to the GTS 2020 timescale.We constructed new age models for samples not already assigned a numeric age. Where the depths of biostratigraphic events were already recorded, these were converted directly to GTS 2020. Where samples were not given any ages, often the case for the cores collected in the early days of ocean drilling, ages were reconstructed from the shipboard and post-cruise biostratigraphic data available in DESCLogik, Pangaea, and drilling publications. For holes where no tie point data were retrievable, biostratigraphic count data were extracted directly from drilling publications, and biostratigraphic events were assigned via GTS 2020. The first and last occurrences in raw shipboard biostratigraphic data often do not represent true datums, and careful assessment of the shipboard, and post-cruise literature was a prerequisite to confidently assigning chronostratigraphic datums. Tie point depths were assigned as the midpoint depth between the core sample before and after an event. For example, for an extinction event, the recorded depth was the midway point between the last recorded occurrence of a species and the first sample from which the species is absent. All sites were assessed individually to determine the age of the seafloor. Where IODP reports or sample-based publications strictly stated that the sediment surface (i.e. 0.00 meters below seafloor (mbsf)) was deemed to be “Holocene”, “Recent” or “Modern” in age, an additional 0 Ma tie point was assigned appropriately. All samples present outside the maximum/minimum age tie points for that site were removed, as they could not be confidently assigned an age. During assessment, individual drilling reports were investigated for geological structures. Where features such as unconformities, reverse faults, stratigraphic inversions, décollements, and major slips and slumps were identified, separate age models were generated for individual intact stratal intervals to account for potential externally emplaced or repeated strata (see “Age models” and “Triton working” in the figshare data repository18). Similarly, age gaps of greater than 10% of the age range of the core were classified as hiatuses, leading to separate age models (see Fig. 4). Cores of denser sediments that have been sampled using rotary drilling will often have only ~50–60% recovery in a core (9.5 m)33. As it is not possible to determine where the recovered core material came from within this length, all intact core pieces are grouped together as a continuous section from the section top, regardless of where the pieces were sourced (e.g. 4.5 m of recovered material will be recorded as 0–4.5 m of cored interval even if some came from 9–9.5 m). Consequently, age estimates within cores where recovery was low, typically the samples collected longer ago, will necessarily be less certain.Fig. 4Different age model estimates applied to core material from IODP Site U1499A in the South China Sea. Mag – mean age based only on the magnetostratigraphic marker events. Zones – mean age based on all the marker events. Int Mag – interpolation of the points between the magnetostratigraphic marker events. Interp – interpolation between the full set of marker events. Model – the model of age as a function of depth. Note the hiatus between 50 and 100 m. For the shallower section of the dataset, with only three data points, a simple linear model was used. For the deeper section, a GAM smooth was fitted. For this site, the model predictions were chosen as the best fit.Full size imageUsing the updated marker event ages, we created age-depth plots and modelled the best fit to the data. There are different ways of creating these models, and multiple methods were applied to each core. The one that provided the best fit to the original data was chosen (the different age models are available in “Age models” in the figshare data repository18). These choices were confirmed manually (see Fig. 4). The simplest age model used interpolation of the marker events to create ‘zones’ and assign estimated ages assuming a continuous sedimentation rate between the start and end of each of these zones. Where the events do not provide a continuous sequence (e.g. gaps in the zonal markers), age estimates were assigned as the mean of that zone with error estimates of the width of the zone. Where magnetostratigraphic events were present they were given preference. This method leads to different estimates of sedimentation rate for each zone. The more complex age model estimates a smoother sedimentation rate. When there were fewer than 5 marker events, a linear model of age as a function of depth was fitted for the entire core. For larger datasets, generalised additive models (GAMs) for the same variables were used, to allow for variation in sedimentation rates through time. GAMs were run using the mgcv R library, with a gamma value of 1.134. The type of age model used in the analysis was recorded. Where appropriate, the number of points and the r2 of the model are recorded to give an indication of the accuracy of the age model.The latitude and longitude coordinates of samples were recorded in decimal degrees. For all samples except modern ones, plate tectonic reconstructions were necessary to determine the coordinates at which the sample was originally deposited. Reconstructions were performed using the Matthews, et al.35 plate motion model, which is an updated version of the Seton, et al.36 model used by Neptune. Comparisons of age models35,36,37,38,39 suggest this model is most appropriate for the deep sea environment where most of the samples occur, and is able to assign coordinates to significantly more sites than the Scotese39 GPlates model. This test was performed with a subset of the data (10633 unique sites); the Matthews, et al.35 model provided paleocoordinates for 95% of the data, whilst the GPlates model only provided coordinates for 17% of the data. The calculation of paleocoordinates was automated using an adaptation of https://github.com/macroecology/mapast.When sediment samples are derived from multiple sources, duplication will inevitably occur. All such duplicated records, identified based on the combination of species, abundance, sample depth, and coordinate values, were removed. Additionally, working on an individual record level, species that occurred significantly outside their known ranges were flagged (following updated age models) on the assumption these records were misidentifications, contamination or re-working. Records were classified as falling significantly outside their known range if they were more than 5 Ma outside the species’ range in the Palaeogene (66-23 Ma) and more than 2 Ma in the Neogene (23-0 Ma). These values were chosen based on the tradeoff between removing reworked specimens and allowing for some errors in the age estimates. Age estimates for older samples tend to be less precise. Ages were obtained from Lamyman et al. (in prep) and are available in “PFdata” in the figshare data repository18. In total, 10,990 suspect records were flagged (~2% of all records). More

Based on the present-day distribution of photosynthetic bacteria31, we assume a competitive advantage for anoxygenic photosynthetic bacteria in early environments where electron donors such as Fe2+, H2S, or H2 were present. We also assume the contemporaneous existence of environments where cyanobacterial populations could thrive, providing a seedbed for migration. Non-marine waters provide an example of the latter, supported by the branching of non-marine taxa from basal nodes in cyanobacterial phylogenies44,45 and also by the presence of stromatolites in Archean lacustrine successions46, despite the likelihood that many Archean lakes and rivers had low levels of potential electron donors such as Fe2+ and H2S47.Following Jones et al.40 and Ozaki et al.42, we use Fe (iron) and P (phosphorus) to represent the environment, which is similar to the H2 and P employed in other studies48,49. The logic of this choice is that in Archean oceans, Fe2+ is thought to have been the principal electron donor for anoxygenic photosynthesis50,51, whereas P governed total rates of photosynthesis. (Kasting14 argued that H2 was key to photosynthesis on the early Earth, a view supported by low iron concentrations in some early Archean stromatolites52.). In any event, under the conditions of low P availability thought to have characterized early oceans25,40,49,53,54,55, anoxygenic photosynthesis would have depleted limiting nutrients before alternative electron donors were exhausted. In consequence, rates of photosynthetic oxygen production would be low. As iron availability declined and/or P availability increased, the biosphere would inevitably reach a point where P would remain after Fe2+ had been depleted, expanding the range of environments where cyanobacteria are favored by natural selection42.Our model keeps track of the abundances of anoxygenic photosynthetic bacteria (APB), x1, cyanobacteria, x2, and three crucial chemicals: iron(II) (Fe2+), y1, phosphate (PO43−), y2, and dioxygen (O2), z. Both types of bacteria require phosphate for reproduction. APB needs iron(II) (or some other suitable reductant) as an electron donor in photosynthesis. The following five equations describe the reproduction and death of APB and of cyanobacteria as well as the dynamics of iron(II), phosphate, and dioxygen:$${rm{APB}}: {dot{x}}_{1} ={x}_{1}{y}_{1}{y}_{2}-{x}_{1}+{u}_{1}\ {rm{Cyano}}: {dot{x}}_{2} =c{x}_{2}{y}_{2}-{x}_{2}+{u}_{2}\ {{rm{Fe}}}^{2+}: {dot{y}}_{1} ={f}_{1}-{y}_{1}-{x}_{1}{y}_{1}{y}_{2}-{y}_{1}z\ {{rm{PO}}_{4}}^{3-}: {dot{y}}_{2} ={f}_{2}-{y}_{2}-{x}_{1}{y}_{1}{y}_{2}-{x}_{2}{y}_{2}\ {{rm{O}}}_{2}: dot{z} =a{x}_{2}{y}_{2}-bz-{y}_{1}z$$
(1)
Here, we have omitted to write symbols for those rate constants that, for understanding the GOE, can be set to one without loss of generality (Supplementary Note 1). Each remaining rate constant is a free parameter. Equations (1) thus satisfy redox balance by construction. We are left with a system that has five main parameters: c specifies the rate of reproduction of cyanobacteria; f1 and f2 denote the rates of supply of iron(II) and phosphate, respectively; a denotes biogenic production of oxygen; b denotes geochemical consumption of oxygen. Note that iron(II) and phosphate are also removed by geochemical processes at a rate proportional to their abundance. In addition, iron(II) is used up during anoxygenic photosynthesis, and iron(II) reacts with oxygen and is thereby removed from the system. Phosphate is used up during the growth of APB and cyanobacteria. (We investigate extensions of the model that incorporate bounded bacterial growth rates and organic carbon in Supplementary Note 2 and Supplementary Note 3, respectively.)We posit iron(II) as the primary electron donor for anoxygenic photosynthesis, and for simplicity of presentation, we refer to y1 and f1 in this context. However, as noted above, y1 and f1 can similarly represent the abundances and influxes of other alternative electron donors, especially dihydrogen (H2)56,57 and hydrogen sulfide (H2S)58. Our model, its analytical solution, and the conclusions that follow hold equally well by considering any of these electron donors or all together.We also include small migration rates, u1 and u2, which allow for the possibility that APB and cyanobacteria persist in privileged sites from which they can migrate into the main arena of competition. On the Archean Earth, these parameters could have been affected by the flow of water and by surface winds. For the mathematical analysis presented in the main text, we assume that these rates are negligibly small.The GOE represents the transition from a world dominated by APB (Equilibrium E1) to one that is dominated by cyanobacteria (Equilibrium E2) (Figs. S1, S2). On a slowly changing planet, the abundances of APB and cyanobacteria and of the three chemicals are approximately in steady state. Therefore, we consider the fixed points of Eqs. (1).Pure equilibriaIn the absence of APB and cyanobacteria, the abiotic equilibrium abundances of iron(II) and of phosphate are given by f1 and f2, respectively, and there is no oxygen in the system. If f1f2  > 1, then APB can emerge. Subsequently, the system settles to Equilibrium E1, where only APB are present and there is still no oxygen. E1 is stable against invasion of cyanobacteria if$${f}_{1}-{f}_{2}, > ,frac{(c+1)(c-1)}{c}.$$
(2)
This condition can be fulfilled if the influx of iron, f1, is large enough, or if the influx of phosphate, f2, is small enough. The term on the right-hand side of the inequality is an increasing function of the reproductive rate, c, of cyanobacteria.If cf2  > 1, then the system admits another equilibrium, E2, where only cyanobacteria are present and oxygen is abundant. Equilibrium E2 is stable against invasion of APB if$$a(c{f}_{2}-1), > ,(b+c)({f}_{1}-c).$$
(3)
The left-hand side of the inequality is positive. If the right-hand side is negative (that is, if f1  ,c(a-1).$$
(4)
Condition (4) is understood as follows. If b is sufficiently large, then there is not enough atmospheric oxygen for rusting to render E2 stable against invasion of APB before E1 loses stability; the result is stable coexistence. But if b is sufficiently small, then rusting causes E2 to become stable before E1 becomes unstable. The critical value of b therefore depends on the input of atmospheric oxygen for Equilibrium E2; it is an increasing function of the reproductive rate of cyanobacteria and of their rate of production of oxygen.If a  c(a − 1). Figure 3 shows gradual oxygenation due to decreasing f1. In this case, the transition occurs via the mixed equilibrium, (hat{E}), where both types of bacteria coexist (Fig. 4). A subsequent increase in f1 can cause APB to regain dominance (Fig. S3a).Fig. 3: The GOE can be triggered by a decline in the influx of iron(II) and is gradual if b  > c(a − 1).Equilibrium E1 (APB dominate) loses stability and Equilibrium E2 (cyanobacteria dominate) gains stability when f1 drops below ({f}_{1}^{* }) and (f_1^{prime}), respectively. We set f2 = 80, c = 10, a = 10, b = 100, and u1 = u2 = 10−3. a We simulate Eqs. (8) from Supplementary Note 1 with α1 = α2 = β1 = β2 = 1, and we set f1 = 100 − 40(t/105). t* denotes the time at which Equilibrium E1 loses stability. b There is stable coexistence of both types of bacteria for (f_1^{prime} , More