Standardized multi-omics of Earth’s microbiomes reveals microbial and metabolite diversity

Dataset description

Sample collection

Our research complies with all relevant ethical regulations following policies at the University of California, San Diego (UCSD). Animal samples that were sequenced were not collected at UCSD and are not for vertebrate animals research at UCSD following the UCSD Institutional Animal Care and Use Committee (IACUC). Samples were contributed by 34 principal investigators of the Earth Microbiome Project 500 (EMP500) Consortium and are samples from studies at their respective institutions (Supplementary Table 1). Relevant permits and ethics information for each parent study are described in the ‘Permits for sample collection’ section below. Samples were contributed as distinct sets referred to here as studies, where each study represented a single environment (for example, terrestrial plant detritus). To achieve more even coverage across microbial environments, we devised an ontology of sample types (microbial environments), the EMP Ontology (EMPO) (http://earthmicrobiome.org/protocols-and-standards/empo/)¹, and selected samples to fill out EMPO categories as broadly as possible. EMPO recognizes strong gradients structuring microbial communities globally, and thus classifies microbial environments (level 4) on the basis of host association (level 1), salinity (level 2), host kingdom (if host-associated) or phase (if free-living) (level 3) (Fig. 1a). As we anticipated previously¹, we have updated the number of levels as well as states therein for EMPO (Fig. 1b) on the basis of an important additional salinity gradient observed among host-associated samples when considering the previously unreported shotgun metagenomic and metabolomic data generated here (Fig. 3c,d). We note that although we were able to acquire samples for all EMPO categories, some categories are represented by a single study.

Samples were collected following the Earth Microbiome Project sample submission guide⁵⁰. Briefly, samples were collected fresh, split into 10 aliquots and then frozen, or alternatively collected and frozen, and subsequently split into 10 aliquots with minimal perturbation. Aliquot size was sufficient to yield 10–100 ng genomic DNA (approximately 10⁷–10⁸ cells). To leave samples amenable to chemical characterization (metabolomics), buffers or solutions for sample preservation (for example, RNAlater) were avoided. Ethanol (50–95%) was allowed as it is compatible with LC–MS/MS although it should also be avoided if possible.

Sampling guidance was tailored for four general sample types: bulk unaltered (for example, soil, sediment, faeces), bulk fractionated (for example, sponges, corals, turbid water), swabs (for example, biofilms) and filters. Bulk unaltered samples were split fresh (or frozen), sampled into 10 pre-labelled 2 ml screw-cap bead beater tubes (Sarstedt, 72.694.005 or similar), ideally with at least 200 mg biomass, and flash frozen in liquid nitrogen (if possible). Bulk fractionated samples were fractionated as appropriate for the sample type, split into 10 pre-labelled 2 ml screw-cap bead beater tubes, ideally with at least 200 mg biomass, and flash frozen in liquid nitrogen (if possible). Swabs were collected as 10 replicate swabs using 5 BD SWUBE dual cotton swabs with wooden stick and screw cap (281130). Filters were collected as 10 replicate filters (47 mm diameter, 0.2 um pore size, polyethersulfone (preferred) or hydrophilic PTFE filters), placed in pre-labelled 2 ml screw-cap bead beater tubes, and flash frozen in liquid nitrogen (if possible). All sample types were stored at –80 °C if possible, otherwise –20 °C.

To track the provenance of sample aliquots, we employed a QR coding scheme. Labels were affixed to aliquot tubes before shipping when possible. QR codes had the format ‘name.99.s003.a05’, where ‘name’ is the PI name, ‘99’ is the study ID, ‘s003’ is the sample number and ‘a05’ is the aliquot number. QR codes (version 2, 25 pixels × 25 pixels) were printed on 1.125’ × 0.75’ rectangular and 0.437’ circular cap Cryogenic Direct Thermal labels (GA International, DFP-70) using a Zebra model GK420d printer and ZebraDesigner Pro 3 software for Windows. After receipt but before aliquots were stored in freezers, QR codes were scanned into a sample inventory spreadsheet using a QR scanner.

Sample metadata

Environmental metadata were collected for all samples on the basis of the EMP Metadata Guide, which combines guidance from the Genomics Standards Consortium MIxS (Minimum Information about any Sequence) standard⁷⁴ and the Qiita Database (https://qiita.ucsd.edu)⁵¹. The metadata guide provides templates and instructions for each MIxS environmental package (that is, sample type). Relevant information describing each PI submission, or study, was organized into a separate study metadata file (Supplementary Table 1).

Metabolomics

LC–MS/MS sample extraction and preparation

To profile metabolites among all samples, we used LC–MS/MS, a versatile method that detects tens of thousands of metabolites in biological samples. All solvents and reactants used were LC–MS grade. To maximize the biomass extracted from each sample, the samples were prepared depending on their sampling method (for example, bulk, swabs, filter and controls). The bulk samples were transferred into a microcentrifuge tube (polypropylene, PP) and dissolved in 7:3 MeOH:H₂O using a volume varying from 600 µl to 1.5 ml, depending on the amounts of sample available, and homogenized in a tissue lyser (QIAGEN) at 25 Hz for 5 min. Then, the tubes were centrifuged at 2,000 × g for 15 min, and the supernatant was collected in a 96-well plate (PP). For swabs, the swabs were transferred into a 96-well plate (PP) and dissolved in 1.0 ml of 9:1 ethanol:H₂O. The prepared plates were sonicated for 30 min, and after 12 h at 4 °C, the swabs were removed from the wells. The filter samples were dissolved in 1.5 ml of 7:3 MeOH:H₂O in microcentrifuge tubes (PP) and sonicated for 30 min. After 12 h at 4 °C, the filters were removed from the tubes. The tubes were centrifuged at 2,000 × g for 15 min, and the supernatants were transferred to 96-well plates (PP). The process control samples (bags, filters and tubes) were prepared by adding 3.0 ml of 2:8 MeOH:H₂O and recovering 1.5 ml after 2 min. After the extraction process, all sample plates were dried with a vacuum concentrator and subjected to solid phase extraction (SPE). SPE was used to remove salts that could reduce ionization efficiency during mass spectrometry analysis, as well as the most polar and non-polar compounds (for example, waxes) that cannot be analysed efficiently by reversed-phase chromatography. The protocol was as follows: the samples (in plates) were dissolved in 300 µl of 7:3 MeOH:H₂O and put in an ultrasound bath for 20 min. SPE was performed with SPE plates (Oasis HLB, hydrophilic-lipophilic-balance, 30 mg with particle sizes of 30 µm). The SPE beds were activated by priming them with 100% MeOH, and equilibrated with 100% H₂O. The samples were loaded on the SPE beds, and 100% H₂O was used as wash solvent (600 µl). The eluted washing solution was discarded, as it contains salts and very polar metabolites that subsequent metabolomics analysis is not designed for. The sample elution was carried out sequentially with 7:3 MeOH:H₂O (600 µl) and 100% MeOH (600 µl). The obtained plates were dried with a vacuum concentrator. For mass spectrometry analysis, the samples were resuspended in 130 µl of 7:3 MeOH:H₂O containing 0.2 µM of amitriptyline as an internal standard. The plates were centrifuged at 30 × g for 15 min at 4 °C. Samples (100 µl) were transferred into new 96-well plates (PP) for mass spectrometry analysis.

LC–MS/MS sample analysis

The extracted samples were analysed by ultra-high performance liquid chromatography (UHPLC, Vanquish, Thermo Fisher) coupled to a quadrupole-Orbitrap mass spectrometer (Q Exactive, Thermo Fisher) operated in data-dependent acquisition mode (LC–MS/MS in DDA mode). Chromatographic separation was performed using a Kinetex C₁₈ 1.7 µm (Phenomenex), 100 Å pore size, 2.1 mm (internal diameter) × 50 mm (length) column with a C₁₈ guard cartridge (Phenomenex). The column was maintained at 40 °C. The mobile phase was composed of a mixture of (A) water with 0.1% formic acid (v/v) and (B) acetonitrile with 0.1% formic acid. Chromatographic elution method was set as follows: 0.00–1.00 min, isocratic 5% B; 1.00–9.00 min, gradient from 5% to 100% B; 9.00–11.00 min, isocratic 100% B; followed by equilibration 11.00–11.50 min, gradient from 100% to 5% B; 11.50–12.50 min, isocratic 5% B. The flow rate was set to 0.5 ml min⁻¹.

The UHPLC was interfaced to the orbitrap using a heated electrospray ionization source with the following parameters: ionization mode, positive; spray voltage, +3,496.2 V; heater temperature, 363.90 °C; capillary temperature, 377.50 °C; S-lens RF, 60 arbitrary units (a.u.); sheath gas flow rate, 60.19 a.u.; and auxiliary gas flow rate, 20.00 a.u. The MS¹ scans were acquired at a resolution (at m/z 200) of 35,000 in the m/z 100–1500 range, and the fragmentation spectra (MS²) scans at a resolution of 17,500 from 0 to 12.5 min. The automatic gain control target and maximum injection time were set at 1.0 × 10⁶ and 160 ms for MS¹ scans, and set at 5.0 × 10⁵ and 220 ms for MS² scans, respectively. Up to three MS² scans in data-dependent mode (Top 3) were acquired for the most abundant ions per MS¹ scans using the apex trigger mode (4–15 s), dynamic exclusion (11 s) and automatic isotope exclusion. The starting value for MS² was m/z 50. Higher-energy collision induced dissociation (HCD) was performed with a normalized collision energy of 20, 30 and 40 eV in stepped mode. The major background ions originating from the SPE were excluded manually from the MS² acquisition. Analyses were randomized within plate and blank samples analysed every 20 injections. A quality control mix sample assembled from 20 random samples across the sample types was injected at the beginning, the middle and the end of each plate sequence. The chromatographic shift observed throughout the batch was estimated as less than 2 s, and the relative standard deviation of ion intensity was 15% per replicate.

LC–MS/MS data processing

The mass spectrometry data were centroided and converted from the proprietary format (.raw) to the m/z extensible markup language format (.mzML) using ProteoWizard (ver. 3.0.19, MSConvert tool)⁷⁵. The mzML files were then processed with MZmine 2 toolbox⁷⁶ using the ion-identity networking modules⁷⁷ that allow advanced detection for adduct/isotopologue annotations. The MZmine processing was performed on Ubuntu 18.04 LTS 64-bits workstation (Intel Xeon E5-2637, 3.5 GHz, 8 cores, 64 Gb of RAM) and took ~3 d. The MZmine project, the MZmine batch file (.XML format) and results files (.MGF and .CSV) are available in the MassIVE dataset MSV000083475. The MZmine batch file contains all the parameters used during the processing. In brief, feature detection and deconvolution was performed with the ADAP chromatogram builder⁷⁸ and local minimum search algorithm. The isotopologues were regrouped and the features (peaks) were aligned across samples. The aligned peak list was gap filled and only peaks with an associated fragmentation spectrum and occurring in a minimum of three files were conserved. Peak shape correlation analysis grouped peaks originating from the same molecule and annotated adduct/isotopologue with ion-identity networking⁷⁷. Finally, the feature quantification table results (.CSV) and spectral information (.MGF) were exported with the GNPS module for feature-based molecular networking analysis on GNPS⁷⁹ and with SIRIUS export modules.

LC–MS/MS data annotation

The results files of MZmine (.MGF and .CSV files) were uploaded to GNPS (http://gnps.ucsd.edu)⁵² and analysed with the feature-based molecular networking workflow⁷⁹. Spectral library matching was performed against public fragmentation spectra (MS²) spectral libraries on GNPS and the NIST17 library.

For the additional annotation of small peptides, we used the DEREPLICATOR tools available on GNPS^80,81. We then used SIRIUS⁸² (v. 4.4.25, headless, Linux) to systematically annotate the MS² spectra. Molecular formulae were computed with the SIRIUS module by matching the experimental and predicted isotopic patterns⁸³, and from fragmentation trees analysis⁸⁴ of MS². Molecular formula prediction was refined with the ZODIAC module using Gibbs sampling⁸⁵ on the fragmentation spectra (chimeric spectra or those with poor fragmentation were excluded). In silico structure annotation using structures from biodatabase was done with CSI:FingerID⁸⁶. Systematic class annotations were obtained with CANOPUS⁴¹ and used the NPClassifier ontology⁸⁷.

The parameters for SIRIUS tools were set as follows, for SIRIUS: molecular formula candidates retained, 80; molecular formula database, ALL; maximum precursor ion m/z computed, 750; profile, orbitrap; m/z maximum deviation, 10 ppm; ions annotated with MZmine were prioritized and other ions were considered (that is, [M+H3N+H]+, [M+H]+, [M+K]+, [M+Na]+, [M+H-H2O]+, [M+H-H4O2]+, [M+NH4]+); for ZODIAC: the features were split into 10 random subsets for lower computational burden and computed separately with the following parameters: threshold filter, 0.9; minimum local connections, 0; for CSI:FingerID: m/z maximum deviation, 10 ppm; and biological database, BIO.

To establish putative microbially related secondary metabolites, we collected annotations from spectral library matching and the DEREPLICATOR+ tools and queried them against the largest microbial metabolite reference databases (Natural Products Atlas⁸⁸ and MIBiG⁸⁹). Molecular networking⁷⁹ was then used to propagate the annotation of microbially related secondary metabolites throughout all molecular families (that is, the network component).

LC–MS/MS data analysis

We combined the annotation results from the different tools described above to create a comprehensive metadata file describing each metabolite feature observed. Using that information, we generated a feature-table including only secondary metabolite features determined to be microbially related. We then excluded very low-intensity features introduced to certain samples during the gap-filling step described above. These features were identified on the basis of presence in negative controls that were universal to all sample types (that is, bulk, filter and swab) and by their relatively low per-sample intensity values. Finally, we excluded features present in positive controls for sampling devices specific to each sample type (that is, bulk, filter or swab). The final feature-table included 618 samples and 6,588 putative microbially related secondary metabolite features that were used for subsequent analysis.

We used QIIME 2’s⁹⁰ (v2020.6) ‘diversity’ plugin to quantify alpha-diversity (that is, feature richness) for each sample and ‘deicode’⁹¹ to quantify beta-diversity (that is, robust Aitchison distances, which are robust to both sparsity and compositionality in the data) between each pair of samples. We parameterized our robust Aitchison principal components analysis (RPCA)⁹¹ to exclude samples with fewer than 500 features and features present in fewer than 10% of samples. We used the ‘taxa’ plugin to quantify the relative abundance of microbially related secondary metabolite pathways and superclasses (that is, on the basis of NPClassifier) within each environment (that is, for each level of EMPO 4), and ‘songbird’ v1.0.4⁹² to identify sets of microbially related secondary metabolites whose abundances were associated with certain environments. We parameterized our ‘songbird’ model as follows: epochs, 1,000,000; differential prior, 0.5; learning rate, 1.0 × 10⁻⁵; summary interval, 2; batch size, 400; minimum sample count, 0; and training on 80% of samples at each level of EMPO 4 using ‘Animal distal gut (non-saline)’ as the reference environment. Environments with fewer than 10 samples were excluded to optimize model training (that is, ‘Animal corpus (non-saline)’, ‘Animal proximal gut (non-saline)’, ‘Surface (saline)’). The output from ‘songbird’ includes a rank value for each metabolite in every environment, which represents the log fold change for a given metabolite in a given environment⁹². We compared log fold changes for each metabolite from this run to those from (1) a replicate run using the same reference environment and (2) a run using a distinct reference environment: ‘Water (saline)’. We found strong Spearman correlations in both cases (Supplementary Table 8), and therefore focused on results from the original run using ‘Animal distal gut (non-saline)’ as the reference environment, as it has previously been shown to be relatively unique among other habitats. In addition to summarizing the top 10 metabolites for each environment (Supplementary Table 3), we used the log fold change values in our multi-omics analyses described below.

We used the RPCA biplot and QIIME 2’s⁹⁰ EMPeror⁹³ to visualize differences in composition among samples, as well as the association with samples of the 25 most influential microbially related secondary metabolite features (that is, those with the largest magnitude across the first three principal component loadings). We tested for significant differences in metabolite composition across all levels of EMPO using PERMANOVA implemented with QIIME 2’s ‘diversity’ plugin⁹⁰ and using our robust Aitchison distance matrix as input. In parallel, we used the differential abundance results from ‘songbird’ described above to identify specific microbially related secondary metabolite pathways and superclasses that varied strongly across environments. We then went back to our metabolite feature-table to visualize differences in the relative abundances of those pathways and superclasses within each environment by first selecting features and calculating log-ratios using ‘qurro’⁹⁴, and then plotting using the ‘ggplot2’ package⁹⁵ in R⁹⁶ v4.0.0. We tested for significant differences in relative abundances across environments using Kruskal–Wallis tests implemented with the base ‘stats’ package in R⁹⁶.

GC–MS sample extraction and preparation

To profile volatile small molecules among all samples in addition to what was captured with LC–MS/MS, we used gas chromatography coupled with mass spectrometry (GC–MS). All solvents and reactants were GC–MS grade. Two protocols were used for sample extraction, one for the 105 soil samples and a second for the 356 faecal and sediment samples that were treated as biosafety level 2. The 105 soil samples were received at the Pacific Northwest National Laboratory and processed as follows. Each soil sample (1 g) was weighed into microcentrifuge tubes (Biopur Safe-Lock, 2.0 ml, Eppendorf). H₂O (1 ml) and one scoop (~0.5 g) of a 1:1 (v/v) mixture of garnet (0.15 mm, Omni International) and stainless steel (0.9–2.0 mm blend, Next Advance) beads and one 3 mm stainless steel bead (Qiagen) were added to each tube. Samples were homogenized in a tissue lyser (Qiagen) for 3 min at 30 Hz and transferred into 15 ml polypropylene tubes (Olympus, Genesee Scientific). Ice-cold water (1 ml) was used to rinse the smaller tube and combined into the 15 ml tube. Chloroform:methanol (10 ml, 2:1 v/v) was added and samples were rotated at 4 °C for 10 min, followed by cooling at −70 °C for 10 min and centrifuging at 150 × g for 10 min to separate phases. The top and bottom layers were combined into 40 ml glass vials and dried using a vacuum concentrator. Chloroform:methanol (1 ml, 2:1) was added to each large glass vial and the sample was transferred into 1.5 ml tubes and centrifuged at 1,300 × g. The supernatant was transferred into glass vials and dried for derivatization.

The remaining 356 samples received from UCSD that included faecal and sediment samples were processed as follows: 100 µl of each sample was transferred to a 2 ml microcentrifuge tube using a scoop (MSP01, Next Advance). The final volume of the sample was brought to 1.5 ml, ensuring that the solvent ratio is 3:8:4 H₂O:CHCl₃:MeOH by adding the appropriate volumes of H₂O, MeOH and CHCl₃. After transfer, one 3 mm stainless steel bead (QIAGEN), 400 µl methanol and 300 µl H₂O were added to each tube and the samples were vortexed for 30 s. Then, 800 µl chloroform was added and samples were vortexed for 30 s. After centrifuging at 150 × g for 10 min to separate phases, the top and bottom layers were combined in a vial and dried for derivatization.

The samples were derivatized for GC–MS analysis as follows: 20 µl of a methoxyamine solution in pyridine (30 mg ml⁻¹) was added to the sample vial and vortexed for 30 s. A bath sonicator was used to ensure that the sample was completely dissolved. Samples were incubated at 37 °C for 1.5 h while shaking at 1,000 r.p.m. N-methyl-N-trimethylsilyltrifluoroacetamide (80 µl) and 1% trimethylchlorosilane solution was added and samples were vortexed for 10 s, followed by incubation at 37 °C for 30 min, with 1,000 r.p.m. shaking. The samples were then transferred into a vial with an insert.

An Agilent 7890A gas chromatograph coupled with a single quadrupole 5975C mass spectrometer (Agilent) and an HP-5MS column (30 m × 0.25 mm × 0.25 μm; Agilent) was used for untargeted analysis. Samples (1 μl) were injected in splitless mode, and the helium gas flow rate was determined by the Agilent Retention Time Locking function on the basis of analysis of deuterated myristic acid (Agilent). The injection port temperature was held at 250 °C throughout the analysis. The GC oven was held at 60 °C for 1 min after injection, and the temperature was then increased to 325 °C at a rate of 10 °C min⁻¹, followed by a 10 min hold at 325 °C. Data were collected over the mass range of m/z 50–600. A mixture of FAMEs (C8–C28) was analysed each day with the samples for retention index alignment purposes during subsequent data analysis.

GC–MS data processing and annotation

The data were converted from vendor’s format to the .mzML format and processed using GNPS GC–MS data analysis workflow (https://gnps.ucsd.edu)⁹⁷. The compounds were identified by matching experimental spectra to the public libraries available at GNPS, as well as NIST 17 and Wiley libraries. The data are publicly available at the MassIVE depository (https://massive.ucsd.edu); dataset ID: MSV000083743. The GNPS deconvolution is available in GNPS (https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=d5c5135a59eb48779216615e8d5cb3ac), as is the library search (https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=59b20fc8381f4ee6b79d35034de81d86).

GC–MS data analysis

For multi-omics analyses including GC–MS data, we first removed noisy (that is, suspected background contaminants and artifacts) features by excluding those with balance scores <50%. Balance scores describe compositional consistency of deconvoluted spectra across the dataset, where high values indicate reproducible spectral patterns and thus high-quality spectra. We then used QIIME 2’s ‘deicode’⁹¹ plugin to estimate beta-diversity for each dataset using robust Aitchison distances. The final feature-table for GC–MS beta-diversity analysis included 460 samples and 216 features.

Metagenomics

DNA extraction

For each round of DNA extractions described below for both amplicon and shotgun metagenomic sequencing, a single aliquot of each sample was processed for DNA extraction. DNA was extracted following the EMP 96-sample, magnetic bead-based DNA extraction protocol⁹⁸ following refs. ^99,100,101 and using the QIAGEN MagAttract PowerSoil DNA KF kit (384-sample) (that is, optimized for KingFisher, 27100-EP). Importantly, material from each sample was added to a unique bead tube (containing garnet beads) for single-tube lysis, which has been shown to reduce sample-to-sample contamination common in plate-based extractions¹⁰¹. For bulk samples, 0.1–0.25 g of material was added to each well; for filtered samples, one entire filter was added to each well; for swabbed samples, one swab head was added to each well. The lysis solution was dissolved at 60 °C before addition to each tube, then capped tubes were incubated at 65 °C for 10 min before mechanical lysis at 6,000 r.p.m. for 20 min using a MagNA lyser (Roche). Lysate from each bead tube was then randomly assigned and added to wells of a 96-well plate, and then cleaned-up using the KingFisher Flex system (Thermo Fisher). Resulting DNA was stored at –20 °C for sequencing. We note that whereas QIAGEN does not offer a ‘hybrid’ extraction kit allowing for single-tube lysis and plate-based clean-up, the Thermo MagMAX Microbiome Ultra kit does, and was recently shown to be comparable to the EMP protocol used here¹⁰².

Amplicon sequencing

We generated amplicon sequence data for variable region four (V4) of the bacterial and archaeal 16S rRNA gene, variable region nine (V9) of the eukaryotic 18S rRNA gene, and the fungal internal transcribed spacer one (ITS1). For amplifying and sequencing all targets, we used a low-cost, miniaturized (that is, 5 µl volume), high-throughput (384-sample) amplicon library preparation method implementing the Echo 550 acoustic liquid handler (Beckman Coulter)¹⁰³. The same protocol was modified with different primer sets and PCR cycling parameters depending on the target. Two rounds of DNA extraction and sequencing were performed for each target to obtain greater coverage per sample. For a subset of 500 samples, we also generated high-quality sequence data for full-length bacterial rRNA operons following a previously published protocol¹⁰⁴, which is briefly outlined below.

The protocol for 16S is outlined fully in ref. ¹⁰⁵. To target the V4 region, we used the primers 515F (Parada) (5’-GTGYCAGCMGCCGCGGTAA-3’) and 806R (Apprill) (5’-GGACTACNVGGGTWTCTAAT-3’). These primers are updated from the original EMP 16S-V4 primer sequences^106,107 to (1) remove bias against Crenarchaeota/Thaumarchaeota¹⁰⁸ and the marine freshwater clade SAR11 (Alphaproteobacteria)¹⁰⁹, and (2) enable the use of various reverse primer constructs (for example, the V4-V5 region using the reverse primer 926R¹¹⁰) by moving the barcode/index to the forward primer¹⁰⁸. We note that while we previously named these updated primers ‘515FB’ and ‘806RB’ to distinguish them from the original primers, the ‘B’ may be misinterpreted to indicate ‘Barcode’. To avoid ambiguity, we now use the original names suffixed with the lead author name (that is, ‘515F (Parada)’, ‘806R (Apprill)’ and ‘926R (Quince)’). We highly recommend to always check the primer sequence in addition to the primer name. For Qiita users, studies with ‘library_construction_protocol’ as ‘515f/806rbc’ used the original primers, whereas ‘515fbc/806r’ indicates use of updated primers, where ‘bc’ refers to the location of the barcode.

To facilitate sequencing on Illumina platforms, the following primer constructs were used to integrate adapter sequences during amplification^106,107,111. For the barcoded forward primer, constructs included (5’ to 3’): the 5’ Illumina adapter (AATGATACGGCGACCACCGAGATCTACACGCT), a Golay barcode (12 bp variable sequence), a forward primer pad (TATGGTAATT), a forward primer linker (GT) and the forward primer (515F (Parada)) (GTGYCAGCMGCCGCGGTAA). For the reverse primer, constructs included (5’ to 3’): the reverse complement of 3’ Illumina adapter (CAAGCAGAAGACGGCATACGAGAT), a reverse primer pad (AGTCAGCCAG), a reverse primer linker (CC) and the reverse primer (806R (Apprill)) (GGACTACNVGGGTWTCTAAT).

For each 25 µl reaction, we combined 13 µl PCR-grade water (Sigma W3500, or QIAGEN 17000-10), 10 µl Platinum Hot Start PCR master mix (2X) (Thermo Fisher, 13000014), 0.5 µl of each primer (10 µM) and 1 µl of template DNA. The final concentration of the master mix in each 1X reaction was 0.8X and that of each primer was 0.2 µM. Cycling parameters for a 384-well thermal cycler were as follows: 94 °C for 3 min; 35 cycles of 94 °C for 1 min, 50 °C for 1 min and 72 °C for 105 s; and 72 °C for 10 min. For a 96-well thermal cycler, we recommend the following: 94 °C for 3 min; 35 cycles of 94 °C for 45 s, 50 °C for 1 min and 72 °C for 90 s; and 72 °C for 10 min.

We amplified each sample in triplicate (that is, each sample was amplified in three replicate 25 µl reactions) and pooled products from replicate reactions for each sample into a single volume (75 µl). We visualized expected products between 300–350 bp on agarose gels, and note that while low-biomass samples may yield no visible bands, instruments such as a Bioanalyzer or TapeStation (Agilent) can be used to confirm amplification. We quantified amplicons using the Quant-iT PicoGreen dsDNA Assay kit (Thermo Fisher, P11496) following the manufacturer’s instructions. To pool samples, we combined an equal amount of product from each sample (240 ng) into a single tube and cleaned the pool using the UltraClean PCR Clean-Up kit (QIAGEN, 12596-4) following the manufacturer’s instructions. We checked DNA quality using a Nanodrop (Thermo Fisher), confirming that A260/A280 ratios were between 1.8–2.0.

For sequencing, the following primer constructs were used. Read 1 constructs included (5’ to 3’): a forward primer pad (TATGGTAATT), a forward primer linker (GT) and the forward primer (515F (Parada)) (GTGYCAGCMGCCGCGGTAA). Read 2 constructs included (5’ to 3’): a reverse primer pad (AGTCAGCCAG), a reverse primer linker (CC) and the reverse primer (806R (Apprill)) (GGACTACNVGGGTWTCTAAT). The index primer sequence was AATGATACGGCGACCACCGAGATCTACACGCT, which we highlight as having an extra GCT at the 3’ end compared to Illumina’s index primer sequence, to increase the melting temperature for read 1 during sequencing.

The protocol for 18S is outlined fully in ref. ¹¹². To target variable region nine (V9), we used the primers 1391f (5’-GTACACACCGCCCGTC-3’) and EukBr (5’-TGATCCTTCTGCAGGTTCACCTAC-3’). These primers are based on those of ref. ^113,114 and are designed for use with Illumina platforms. The forward primer is a universal small-subunit primer, whereas the reverse primer favours eukaryotes but with mismatches can bind and amplify Bacteria and Archaea. In addition to deviations from the 16S protocol above with respect to primer construct sequences and PCR cycling parameters, we included a blocking primer that reduces amplification of vertebrate host DNA for host-associated samples, on the basis of the strategy outlined in ref. ¹¹⁵. We note that the blocking primer is particularly useful for host-associated samples with a low biomass of non-host eukaryotic DNA.

The following primer constructs were used to integrate adapter sequences during amplification. For the barcoded forward primer, constructs included (5’ to 3’): the 5’ Illumina adapter (AATGATACGGCGACCACCGAGATCTACAC), a forward primer pad (TATCGCCGTT), a forward primer linker (CG) and the forward primer (Illumina_Euk_1391f) (GTACACACCGCCCGTC). For the reverse primer, constructs included (5’ to 3’): The reverse complement of 3’ Illumina adapter (CAAGCAGAAGACGGCATACGAGAT), a Golay barcode (12 bp variable sequence), a reverse primer pad (AGTCAGTCAG), a reverse primer linker (CA) and the reverse primer (806R (Apprill)) (TGATCCTTCTGCAGGTTCACCTAC). The construct for the blocking primer is as such and is formatted for ordering from IDT: ‘GCCCGTCGCTACTACCGATTGG/ideoxyI//ideoxyI//ideoxyI//ideoxyI//ideoxyI/TTAGTGAGGCCCT/3SpC3/’.

Reaction mixtures without the blocking primer (that is, those for non-vertebrate hosts or free-living sample types as defined by EMPO) were prepared as described for 16S. For reactions including the blocking primer, we combined 9 µl PCR-grade water, 10 µl master mix, 0.5 µl of each primer (10 µM), 4 µl of blocking primer (10 µM) and 1 µl of template DNA. The final concentration of the master mix in each 1X reaction was 0.8X, that of each primer was 0.2 µM and that of the blocking primer was 1.6 µM. Without blocking primers, cycling parameters for a 384-well thermal cycler were as follows: 94 °C for 3 min; 35 cycles of 94 °C for 45 s, 57 °C for 1 min and 72 °C for 90 s; and 72 °C for 10 min. With blocking primers, cycling parameters for a 384-well thermal cycler were as follows: 94 °C for 3 min; 35 cycles of 94 °C for 45 s, 65 °C for 15 s, 57 °C for 30 s and 72 °C for 90 s; and 72 °C for 10 min. Expected bands ranged between 210–310 bp.

For sequencing, the following primer constructs were used. Read 1 constructs (Euk_illumina_read1_seq_primer) included (5’ to 3’): a forward primer pad (TATCGCCGTT), a forward primer linker (CG) and the forward primer (1391f) (GTACACACCGCCCGTC). Read 2 constructs (Euk_illumina_read2_seq_primer) included (5’ to 3’): a reverse primer pad (AGTCAGTCAG), a reverse primer linker (CA) and the reverse primer (EukBr) (TGATCCTTCTGCAGGTTCACCTAC). The index primer construct (Euk_illumina_index_seq_primer) included (5’ to 3’): the reverse complement of the reverse primer (EukBr) (GTAGGTGAACCTGCAGAAGGATCA), the reverse complement of the reverse primer linker (TG) and the reverse complement of the reverse primer pad (CTGACTGACT).

The protocol for ITS is outlined fully in ref. ¹¹⁶. To target the fungal internal transcribed spacer (ITS1), we used the primers ITS1f (5’-CTTGGTCATTTAGAGGAAGTAA-3’) and ITS2 (5’-GCTGCGTTCTTCATCGATGC-3’). These primers are based on those of ref. ¹¹⁷, and we note that primer ITS1f used here binds 38 bp upstream of ITS1 reported in that study.

The following primer constructs were used to integrate adapter sequences during amplification. For the barcoded forward primer, constructs included (5’ to 3’): the 5’ Illumina adapter (AATGATACGGCGACCACCGAGATCTACAC), a forward primer linker (GG) and the forward primer (ITS1f) (CTTGGTCATTTAGAGGAAGTAA). For the reverse primer, constructs included (5’ to 3’): the reverse complement of 3’ Illumina adapter (CAAGCAGAAGACGGCATACGAGAT), a Golay barcode (12 bp variable sequence), a reverse primer linker (CG) and the reverse primer (ITS2) (GCTGCGTTCTTCATCGATGC).

Reaction mixtures were prepared as described for 16S. Cycling parameters for a 384-well thermal cycler were as follows: 94 °C for 1 min; 35 cycles of 94 °C for 30 s, 52 °C for 30 s and 68 °C for 30 ; and 68 °C for 10 min. Expected bands ranged between 250–600 bp^118,119.

For sequencing, the following primer constructs were used. Read 1 sequencing primer constructs included (5’ to 3’): a forward primer segment (TTGGTCATTTAGAGGAAGTAA) and a region extending into the amplicon (AAGTCGTAACAAGGTTTCC). Read 2 sequencing primer constructs included (5’ to 3’): a reverse primer segment (CGTTCTTCATCGATGC) and a region extending into the amplicon (VAGARCCAAGAGATC). The index sequencing primer construct included (5’ to 3’): the reverse complement of the region extending into the amplicon (TCTC), the reverse complement of the reverse primer (GCATCGATGAAGAACGCAGC) and the reverse complement of the linker (CG).

The protocol for generating bacterial full-length rRNA operon data is described in ref. ¹⁰⁴. The method uses a unique molecular identifier (UMI) strategy to remove PCR errors and chimeras, resulting in a mean error rate of 0.0007% and a chimera rate of 0.02% of the final amplicon data. Briefly, the bacterial rRNA operons were targeted with an initial PCR using tailed versions of 27f (AGRGTTYGATYMTGGCTCAG)¹²⁰ and 2490r (GACGGGCGGTGWGTRCA)¹²¹. The primer tails contained synthetic priming sites and 18-bp-long patterned UMIs (NNNYRNNNYRNNNYRNNN). The PCR reaction (50 µl) contained 1–2 ng DNA template, 1 U Platinum SuperFi DNA Polymerase High Fidelity (Thermo Fisher) and a final concentration of 1× SuperFi buffer, 0.2 mM of each deoxynucleotide triphosphate, and 500 nM of each tailed 27f and tailed 2490r. The PCR cycling parameters consisted of an initial denaturation (3 min at 95 °C) and two cycles of denaturation (30 s at 95 °C), annealing (30 s at 55 °C) and extension (6 min at 72 °C). The PCR product was purified using a custom bead purification protocol ‘SPRI size selection protocol for >1.5–2 kb DNA fragments’ (Oxford Nanopore Technologies). The resulting product consists of uniquely tagged rRNA operon amplicons. The uniquely tagged rRNA operons were amplified in a second PCR, where the reaction (100 µl) contained 2 U Platinum SuperFi DNA Polymerase High Fidelity (Thermo Fisher) and a final concentration of 1X SuperFi buffer, 0.2 mM of each dNTP, and 500 nM of each forward and reverse synthetic primer targeting the tailed primers from above. The PCR cycling parameters consisted of an initial denaturation (3 min at 95 °C) and then 25–35 cycles of denaturation (15 s at 95 °C), annealing (30 s at 60 °C) and extension (6 min at 72 °C), followed by final extension (5 min at 72 °C). The PCR product was purified using the custom bead purification protocol above. Batches of 25 amplicon libraries were barcoded and sent for PacBio Sequel II library preparation and sequencing (Sequel II SMRT Cell 8M and 30 h collection time) at the DNA Sequencing Center at Brigham Young University. Circular consensus sequencing (CCS) reads were generated using CCS v.3.4.1 (https://github.com/PacificBiosciences/ccs) using default settings. UMI consensus sequences were generated using the longread_umi pipeline (https://github.com/SorenKarst/longread_umi) with the following command: longread_umi pacbio_pipeline -d ccs_reads.fq -o out_dir -m 3500 -M 6000 -s 60 -e 60 -f CAAGCAGAAGACGGCATACGAGAT -F AGRGTTYGATYMTGGCTCAG -r AATGATACGGCGACCACCGAGATC -R CGACATCGAGGTGCCAAAC -U ‘0.75;1.5;2;0’ -c 2.

Amplicon data analysis

For multi-omics analyses including amplicon sequence data, we processed each dataset for comparison of beta-diversity. For all amplicon data except that for bacterial full-length rRNA amplicons, raw sequence data were converted from bcl to fastq, and then multiplexed files for each sequencing run uploaded as separate preparations to Qiita (study: 13114).

For each 16S sequencing run, in Qiita, data were demultiplexed, trimmed to 150 bp and denoised using Deblur¹²² to generate a feature-table of sub-operational taxonomic units (sOTUs) per sample, using default parameters. We then exported feature-tables and denoised sequences from each sequencing run, used QIIME 2’s ‘feature-table’ plugin to merge feature-tables and denoised reads across sequencing runs, and placed all denoised reads into the GreenGenes 13_8 phylogeny¹²³ via fragment insertion using QIIME 2’s⁹⁰ SATé-Enabled Phylogenetic Placement (SEPP)¹²⁴ plugin to produce a phylogeny for diversity analyses. To allow for phylogenetically informed diversity analyses, reads not placed during SEPP (that is, 513 sOTUs, 0.1% of all sOTUs) were removed from the merged feature-table. We then used QIIME 2’s ‘feature-table’ plugin to exclude singleton sOTUs and rarefy the data to 5,000 reads per sample. Rarefaction depths for all amplicon analyses were chosen to best normalize sampling effort per sample while maintaining ≥75% of samples representative of Earth’s environments, and also to maintain consistency with the analyses from EMP release 1. We then used QIIME 2’s⁹⁰ ‘diversity’ plugin to estimate alpha-diversity (that is, sOTU richness) and beta-diversity (that is, unweighted UniFrac distances). The final feature-table for 16S beta-diversity analysis included 681 samples and 93,260 features. We performed a comparative analysis of the data including and excluding the reads not placed during SEPP, and note that both alpha-diversity (that is, sOTU richness) and beta-diversity (that is, sample–sample RPCA distances) were highly correlated between datasets (Spearman r = 1.0) (Supplementary Fig. 5). We thus proceeded with the SEPP-filtered dataset and used phylogenetically informed diversity metrics where applicable.

For 18S data, we used QIIME 2’s⁹⁰ ‘demux’ plugin’s ‘emp-paired’ method^125,126 to first demultiplex each sequencing run, and then the ‘cutadapt’ plugin’s¹²⁷ ‘trim-paired’ method to trim sequencing primers from reads. We then exported trimmed reads, concatenated R1 and R2 read files per sample, and denoised reads using Deblur’s^122,128 ‘workflow’ with default settings, trimming reads to 90 bp, and taking the ‘all.biom’ and ‘all.seqs’ output, for each sequencing run. We then used QIIME 2’s ‘feature-table’ plugin to merge feature-tables and denoised sequences across sequencing runs, and then the ‘feature-classifier’ plugin’s ‘classify-sklearn’ method to classify taxonomy for each sOTU via pre-fitted machine-learning classifiers¹²⁹ and the SILVA 138 reference database¹³⁰. We then used QIIME 2’s⁹⁰ ‘feature-table’ plugin to exclude reads assigned to bacteria and archaea, singleton sOTUs and samples with a total frequency of <5,500 reads, and the ‘deicode’⁹¹ plugin to estimate beta-diversity for each dataset using robust Aitchison distances⁹¹. The final feature-table for 18S beta-diversity analysis included 461 samples and 14,839 features.

For fungal ITS data, we used QIIME 2⁹⁰ to generate and merge feature-tables and denoised sequences across sequencing runs, as for 18S data but trimming reads to 150 bp. We then classified taxonomy for each sOTU as for 18S data, but using the UNITE 9 reference database¹³¹. We then used QIIME 2’s ‘feature-table’ plugin to exclude singleton sOTUs and samples with a total frequency of <500 reads, and the ‘deicode’⁹¹ plugin to estimate beta-diversity for each dataset using robust Aitchison distances⁹¹. The final feature-table for fungal ITS beta-diversity analysis included 500 samples and 10,966 features.

For full-length rRNA operon data, per-sample fasta files were reformatted for importing to QIIME 2 as ‘SampleData[Sequences]’ (that is, with each header as ‘>{sample_identifier}_{sequence_identifier}’), concatenated into a single fasta file and imported. We then used QIIME 2’s ‘vsearch’ plugin¹³² to dereplicate sequences and then cluster them at 65% similarity (that is, due to rapid evolution at bacterial ITS regions). The 65% OTU feature-table had 365 samples and 285 features. The concatenated fasta file and 65% OTU feature-table were uploaded to Qiita as distinct preparations (study: 13114). We then used QIIME 2’s⁹⁰ ‘feature-table’ plugin to exclude singleton OTUs and samples with a total frequency of <500 reads, and the ‘deicode’⁹¹ plugin to estimate beta-diversity for each dataset using robust Aitchison distances⁹¹. The final feature-table for full-length rRNA operon beta-diversity analysis included 242 samples and 196 features.

Shotgun metagenomic sequencing

One round of DNA extraction was performed as above for shotgun metagenomic sequencing. Sequencing libraries were prepared using a high-throughput version of the HyperPlus library chemistry (Kapa Biosystems) miniaturized to approximately 1:10 reagent volume and optimized for nanolitre-scale liquid-handling robotics¹³³. An exhaustive step-by-step protocol and accompanying software can be found in ref. ¹³³. Briefly, DNA from each sample was transferred to a 384-well plate and quantified using the Quant-iT PicoGreen dsDNA Assay kit (P7589, Thermo Fisher), and then normalized to 5 ng in 3.5 µl of molecular-grade water using an Echo 550 acoustic liquid-handling robot (Labcyte). For library preparation, reagents for each step (that is, fragmentation, end repair and A-tailing, ligation and PCR) were added at 1:10 the recommended volumes using a Mosquito HTS micropipetting robot (SPT Labtech). Fragmentation was performed at 37 °C for 20 min and A-tailing at 65 °C for 30 min.

Sequencing adapters and barcode indices were added in two steps¹³⁴. First, the Mosquito HTS robot was used to add universal adapter ‘stub’ adapters and ligase mix to the end-repaired DNA, and the ligation reaction performed for 20 °C for 1 h. Adapter-ligated DNA was then cleaned-up using AMPure XP magnetic beads and a BlueCat purification robot (BlueCat Bio) by adding 7.5 µl magnetic bead solution to the total sample volume, washing twice with 70% ethanol and resuspending in 7 µl molecular-grade water. Then, the Echo 550 robot was used to add individual i7 and i5 indices to adapter-ligated samples without repeating any barcodes, and iterate the assignment of i7 to i5 indices to minimize repeating unique i7:i5 pairs. Cleaned adapter-ligated DNA was then amplified by adding 4.5 µl of each sample to 5.5 µl PCR master mix and running for 15 cycles, and then purified again using magnetic beads and the BlueCat robot. Each sample was eluted into 10 µl water, and then transferred to a 384-well plate using the Mosquito HTS robot. Each library was quantified using qPCR and then pooled to equal molar fractions using the Echo 550 robot. The final pool was sequenced at Illumina on a NovaSeq6000 using S2 flow cells and 2 × 150 bp chemistry (Illumina). To increase sequence coverage for certain samples, libraries were re-pooled and a second sequencing run performed as above.

Shotgun data analysis

Raw sequence data were converted from bcl to fastq and demultiplexed to produce per-sample fastq files. The mean sequencing depth was 7,580,347 ± 7.82 × 10¹³ reads per sample. We processed raw reads with Atropos (v1.1.24)¹³⁵ to trim universal adapter sequences, poly-G tails introduced by the NovaSeq instrument (that is, from use of two-colour chemistry) and low-quality bases from reads. Atropos parameters included poly-G trimming (nextseq-trim=30), inclusion of ambiguous bases (match-read-wildcards), a maximum error rate for adapter matching (error-rate=0.1, default), removal of low-quality bases at 3’ and 5’ ends before adapter removal (quality-cutoff=15), a maximum error rate for adapter matching (insert-match-error-rate=0.2, default), discarding of short trimmed reads (minimum-length=100) and discarding of paired reads if even one fails filtering (pair-filter=any). Trimmed reads were then mapped to the Web of Life database of microbial genomes (release 1)¹³⁶ using bowtie2 v2.3.2¹³⁷ in very-sensitive mode to produce alignments that were used for taxonomic and exploratory functional analysis of microbial communities. Bowtie2 settings included maximum and minimum mismatch penalties (mp=[1,1]), a penalty for ambiguities (np=1; default), read and reference gap open- and extend penalties (rdg=[0,1], rfg=[0,1]), a minimum alignment score for an alignment to be considered valid (score-min=[L,0,−0.05]), a defined number of distinct valid alignments (k=16), and the suppression of SAM records for unaligned reads, as well as SAM headers (no-unal, no-hd). The Web of Life database is particularly attractive as it includes a phylogeny that can be used for diversity analyses, and was curated to represent phylogenetic breadth of Bacteria and Archaea¹³⁶, ideal for analyses across diverse environments. We compared mapping to the Web of Life to Rep200, a curated database of NCBI representative and reference microbial genomes (that is, corresponding to RefSeq release 200, released 14 May 2020) and found little difference across environments (Supplementary Fig. 6). We therefore chose the Web of Life as it allows for phylogenetically informed analyses.

For taxonomic analysis, we generated a feature-table of counts of operational genomic units (OGUs) for each sample using a reference-based approach. We chose this method over the de novo or reference-free approach, as the latter uses assembly/clustering to deconvolute short reads into larger sequence units; the reference-free approach allows for the direct observation of the actual organisms in the community, but alone does not allow their meaningful characterization⁶. Reference-based approaches use reference sequences from described organisms, allowing us to find the closest matches and use them to describe the taxa in a community⁶. This strategy is advantageous as results are not dependent on the samples included and it is less difficult because sequences can more easily be aligned to a reference vs assembled into MAGS^138,139. Most importantly, it allows for comparisons of results across samples and studies, therefore representing a standardized method. Specifically, we used Woltka’s v0.1.4¹⁴⁰ ‘classify’ function, with per-genome alignments and default parameters. Woltka’s default normal mode is such that for one query sequence mapped to k genomes, each genome receives a count of 1/k. To permit examination of rare taxa across environments, no genomes were excluded. For diversity analyses, to best normalize sampling effort per sample while maintaining ≥75% of samples representative of Earth’s environments, we rarefied the OGU feature-table to 6,550 reads per sample. The final feature-table for analyses of shotgun metagenomic taxonomic diversity included 612 samples and 8,692 OGUs.

For alpha-diversity, we quantified three metrics, in part to see which had the strongest correlations with microbially related metabolite richness. We used the R package ‘geiger’¹⁴¹ to quantify weighted Faith’s PD for each sample following the method of Swensen¹⁴². We used QIIME 2’s ‘diversity’ plugin⁹⁰ to quantify richness and Faith’s PD (that is, unweighted), as well as beta-diversity (that is, using weighted UniFrac distance) between each pair of samples. We performed PERMANOVA on that distance matrix to test for significant differences in microbial community composition across the various levels of EMPO, and verified that differences were robust across sampling depths spanning three orders of magnitude (Supplementary Table 9). We then used principal coordinates analysis (PCoA) and EMPeror⁹³ to visualize differences in microbial community composition among samples. We used ‘songbird’⁹² to identify sets of microbial taxa whose abundances were associated with certain environments, and parameterized our songbird model as above for our LC–MS/MS data. We then mapped the differential abundance results from songbird onto a phylogeny representing all microbial taxa using ‘empress’¹⁴³ to visualize phylogenetic relationships related to log fold changes in abundance relative to specific environments.

For the functional analysis, we initially generated two sets of annotations for comparison of read mapping across environments. First, we generated a feature-table of counts of Gene Ontology (GO) terms (that is, for biological process, molecular function and cellular compartment) for each sample using Woltka’s ‘collapse’ function, inputting per-gene alignments and with default parameters for mapping to GO terms through MetaCyc. For subsequent analysis, we used QIIME 2’s⁹⁰ ‘feature-table’ plugin to exclude singleton features and rarefy the data to 5,000 sequences per sample. The final feature-table included 517 samples and 3,776 features (that is, GO terms). We also generated a feature-table of counts of KEGG^144,145,146 EC features (that is, enzymes) for each sample using PRROMenade¹⁴⁷. Trimmed, quality-controlled reads were mapped to the PRROMenade index of bacterial and viral protein domains via the IBM Functional Genomics Platform¹⁴⁸ following ref. ¹⁴⁹, searching for maximal exact matches with a length ≥11 amino acids and retaining samples with ≥10,000 annotated reads (that is, summed across R1 and R2 read files). Annotated read counts were pushed to leaf level nodes in the four-level EC hierarchy (for example, EC 1.2.3.4). For diversity analysis, we used QIIME 2’s⁹⁰ ‘feature-table’ plugin to exclude singleton features and samples with fewer than 150,000 reads. The final feature-table included 616 samples (representing 18 environments) and 1,250 enzymes (that is, KEGG ECs). We performed a comparative analysis comparing the Woltka GO-term analysis and the PRROMenade KEGG EC analysis, and found PRROMenade to more efficiently map reads across the majority of environments (Supplementary Fig. 4). We therefore proceeded with our analysis of microbial functions using PRROMenade. With that table, we used QIIME 2’s ‘deicode’⁹¹ plugin to estimate beta-diversity for each dataset using robust Aitchison distances⁹¹ and EMPeror⁹³ to visualize differences in microbial community composition among samples. We then performed PERMANOVA as above to test for significant differences in microbial functional composition across the various levels of EMPO.

Nestedness analysis of metabolites and microbial taxa

As our analysis of turnover (replacement) of microbial taxa suggested a degree of nestedness (gain or loss of taxa promoting differences in richness) among environments in line with previous observations based on EMP 16S release 1, we tested for nestedness in our shotgun metagenomics data for microbial taxa. We used the NODF statistic¹⁵⁰ to quantify nestedness on the basis of the degree to which less diverse communities are subsets of more diverse communities, which we quantified at each major taxonomic level from phylum to species. We used the rarefied feature-table described above and a null model (that is, equiprobable rows, fixed columns) for assessing observed values of NODF, which we considered at each taxonomic level, and for all of the samples and each subset of the samples at EMPO 2. To compute standardized effect sizes and P values for significance, we used simulated results (n = 10,000 iterations) to find the expectation and variance of the NODF statistic under the null model. Standardized effect sizes were large (>90).

Multi-omics

Alpha-diversity correlations

Using the alpha-diversity metrics for LC–MS/MS (that is, richness) and shotgun metagenomic taxonomic data (that is, richness, unweighted Faith’s PD and weighted Faith’s PD), we performed correlation analysis to better understand relationships therein. We used the function ‘multilevel’ available in the R package ‘correlation’¹⁵¹ to perform Spearman correlations for each environment (that is, based on EMPO 4), treating study (that is, the variable representing distinct PI submissions of samples) as a random effect and adjusting for multiple comparisons using the Benjamini-Hochberg correction. We performed additional correlations with our shotgun metagenomics data rarefied to sampling depths across three orders of magnitude, and confirmed that patterns observed with our focal sampling depth of 6,550 are robust, although loss of samples at higher sampling depths results in reduced effects (Supplementary Table 10).

Machine-learning analyses

To better understand community composition of microbes and metabolites across environments and specifically which features are predictive of certain habitats, we performed machine learning. For analyses of LC–MS/MS and shotgun metagenomic taxonomic and functional data, additional samples were filtered from the feature-tables noted previously to exclude environments with relatively low sample representation (that is, <9 samples). For the LC–MS/MS feature-table, we excluded samples in the four EMPO environments (that is, ‘Animal corpus (non-saline)’, ‘Animal proximal gut (non-saline)’, ‘Soil (saline)’ and ‘Surface (saline)’). The final feature-table included 605 samples (representing 15 environments) and 6,588 microbially related metabolites. For the shotgun metagenomic feature-table for taxonomic analysis, we excluded samples in four EMPO environments (that is, ‘Animal corpus (non-saline)’, ‘Fungus corpus (non-saline)’, ‘Surface (saline)’ and ‘Subsurface (non-saline)’). The final feature-table included 598 samples (representing 15 environments) and 8,587 microbial taxa (that is, Woltka OGUs). For the shotgun metagenomic feature-table for functional analysis, we used QIIME 2’s⁹⁰ ‘feature-table’ plugin to exclude samples in three EMPO environments (that is, ‘Animal corpus (non-saline)’, ‘Surface (saline)’ and ‘Subsurface (non-saline)’), exclude singleton features and normalize the total count per sample to 10,000 sequences. The final feature-table included 706 samples (representing 16 environments) and 1,133 enzymes (that is, KEGG ECs).

For each feature-table, we trained an auto-AI classifier¹⁵² with SHAP explanations¹⁵³ and the hyper-tuned XGBoost method¹⁵⁴ for predicting environments (on the basis of EMPO 4). Each dataset was split into a training set (80%) and a testing set (20%), with similar environmental distributions in each iteration for the classification of samples. We evaluated the predictive performance of each classifier by quantifying accuracy statistics across 20 randomized iterations, and specifically by using resulting confusion matrices to quantify the overall and per-environment precision, recall and F1 score. To identify the most important features contributing to the classification, we examined SHAP explanations, which we used to describe the impact of each feature for prediction. For features with an impact in at least one of 20 iterations examined, we assigned absolute ranks for each feature per iteration, and then assigned final ranks on the basis of the mean of absolute ranks across iterations. For the top 20 ranked features per feature-table, we visualized the environment for which each feature was impactful, as well as the direction of impact. Direction was determined by assessing differences in the mean relative abundances of the focal environment vs all other environments combined. Positive impact indicates that a feature was predictive of the focal environment when it was more abundant there vs the other environments.

Metabolite–microbe co-occurrence analysis

To begin to explore co-occurrences between microbes and metabolites across environments, we implemented an approach that generates co-occurrence probabilities between all metabolite and microbial features, clusters metabolites on the basis of their co-occurrence with the microbial community and highlights individual microbial features driving global patterns in metabolite distribution in this space. For co-occurrence analyses of LC–MS/MS metabolites and genomes profiled from shotgun metagenomic data, feature-tables were further filtered to retain only the 434 samples found in both datasets. For the LC–MS/MS feature-table of microbially related secondary metabolites, we excluded 172 samples lacking shotgun metagenomics data, resulting in a final set of 6,501 microbially related metabolites. For the shotgun metagenomics feature-table for taxonomy, we excluded 150 samples lacking LC–MS/MS data, resulting in a final set of 4,120 OGUs.

Specifically, we obtained co-occurrence probabilities and ordinated metabolites in microbial taxon space using ‘mmvec’ v1.0.6, which uses the probabilities (that is, log conditional probabilities, or co-occurrence strength) to predict metabolites on the basis of microbial taxa from neural-network, compositionally robust modelling⁵⁶. The model was trained on 80% of the 434 samples, which were selected to balance environments (that is, EMPO 4), and used the following parameters: epochs, 200; batch size, 165; learning rate, 1.0 × 10⁻⁵; summary interval, 1; and with ‘equalize-biplot’. For training and testing, we filtered to retain only those features present in at least 10 samples (that is, min-feature-count, 10), and restricted decomposition of the co-occurrence matrix to 10 principal components (PCs) (that is, latent-dim, 10). The model predicting metabolite–microbe co-occurrences was more accurate than one representing a random baseline, with a pseudo-Q² value of 0.18, indicating much reduced error during cross-validation.

To relate these metabolite–microbe co-occurrences to the distribution of metabolites across environments, we calculated the Spearman correlation between the loadings of metabolites on each co-occurrence PC vs (1) log fold changes in metabolite abundances for each environment (that is, from ‘songbird’), (2) loadings for metabolites on the first three axes from the ordination corresponding to clustering of samples by environment (that is, from RPCA) and (3) a vector representing the global magnitude of metabolite importance across all three axes from that same ordination. To explicitly highlight metabolite–microbe co-occurrences specific to particular environments, we visualized the relationships between metabolite–microbe co-occurrences and (1) by considering the first three PCs of the co-occurrence ordination (that is, from mmvec) and colouring metabolites by their log fold change values for a focal environment (for example, Fig. 4b, Extended Data Fig. 10 and Supplementary Table 5). Then, focusing on the co-occurrence PC exhibiting the strongest correlation with log fold changes in metabolite abundances with respect to the focal environment, we manually selected one subset of metabolites highly abundant with respect to the focal environment but similar with respect to co-occurrences with microbes (that is, high values on both axes, the focal group of metabolites) and one subset of metabolites lowly abundant with respect to the focal environment but similar with respect to co-occurrences with microbes (that is, low values on both axes, the reference group of metabolites)¹⁵⁵. Each select group of metabolites was chosen to represent a single pathway. Then, depending on the focal environment, we chose either the top 10 or top 10% of co-occurring microbes (that is, on the basis of co-occurrence strength) for each of the focal and reference groups of metabolites¹⁵⁴. Finally, we visualized differences in the log-ratio of the focal group to the reference group between the focal environment and all other environments, separately for metabolites and microbes¹⁵⁴.

Mantel correlations between datasets

To explore the relationships between sample–sample distances for any two datasets (for example, LC–MS/MS vs shotgun metagenomics for taxonomy), we used QIIME 2’s ‘diversity’ plugin⁹⁰ to perform Mantel tests on all pairings of the datasets using Spearman correlations. Input distance matrices are those described above for each dataset.

Statistics and reproducibility

Samples and studies were crowd-sourced to span microbial environments described by EMPO version 1. Before acceptance as an EMP500 study, scientific justification was required (Supplementary Table 1). Sample sizes for each study were determined by each EMP500 PI on the basis of sample availability (that is, no statistical method was used to predetermine sample size, but our sample sizes are similar to those reported in previous studies^{59,60,61,62,156}). No samples were excluded from analysis, except when inclusion violated assumptions or best practices of statistical tests, which we detail for each method used above. As each sample was split into 10 aliquots, samples from several studies are available for future use. Similarly, as we used standardized protocols and methods throughout from sample collection to data analysis, the results are reproducible. No experiments requiring randomization or blindness were carried out. For each analysis, we used non-parametric statistical tests unless tests for normality and equal variances showed that these assumptions were met.

Permits for sample collection

For all animal, geological and international sample collection, the proper procedures for sampling, exporting and importing were followed. In accordance with the genetic resource sharing component of the Nagoya Protocol, we have made all sequence data publicly available at NCBI. Here we provide specific statements for sample collection where relevant.

Animal specimens and geological samples used in this study were collected for a range of different parent studies and were contributed to the project at UCSD (Supplementary Table 1). Therefore, based on IACUC policies, this project was not considered vertebrate animal research at UCSD. Here we provide relevant ethical information for samples from parent studies that included animals: studies 9, 18, 63 and 72 include only lower-level invertebrates and are thus exempt from animal use protocol based on IACUC guidelines; study 50 did not require handling of animals; collection for studies 51–53 and 54 was approved by the University of Colorado, Boulder (UCB); collection for study 54 was approved by UCSD (protocols S12219, S09392); study 76 did not require handling of animals; collection for study 81 was approved by UCB (protocol 08-04-AK-01); collection for study 88 was approved by the Animal Experiment Board in Finland (ESAVI/7256/04.10.07/2014).

For samples from Costa Rica, permits were granted by the Institutional Biodiversity Commission of the University of Costa Rica (UCR, resolution number 055-2016) and authorized by the Organization for Tropical Studies (OTS) and the Central Pacific Conservation Area (ACOPAC) of the Ministry of Energy and the Environment (MINAE), Costa Rican government, under UCR project B6-656.

For samples from Ukraine, all procedures were performed in accordance with legal requirements and regulations from Ukrainian authorities (957-i/16/05/2016), and the Animal Experiment Board in Finland (ESAVI/7256/04.10.07/2014). The samples were transported to Finland for research purposes on the basis of the import permission from the Evira (3679/0460/2016).

Samples from Namibia were collected under the Republic of Namibia – Ministry of Mines and Energy permit number ES30246 and transported to South Africa for research purpose with import permit P0067933 from the Department of Agriculture, Forestry and Fisheries of the Republic of South Africa.

For samples from Singapore, permit (No:NP/RP18-086) was granted by National Parks Board (NParks) and sampling was conducting according to stipulations of the permit.

All coral samples were collected by AAUS-certified scientific divers, in accordance with local regulations. Relevant permit numbers are: CITES (PWS2014-AU-002155, 12US784243/9), Great Barrier Reef Marine Park Authority (G12/35236.1, G14/36788.1), Lord Howe Island Marine Park (LHIMP/R/2015/005), New South Wales Department of Primary Industries (P15/0072–1.0, OUT 15/11450), US Fish and Wildlife Service (2015LA1632527, 2015LA1703560), and Western Australia Department of Parks and Wildlife (SF010348, CE004874, ES002315).

How to access and contribute to the EMP500

All methods and protocols can be accessed at www.earthmicrobiome.org and GitHub (https://github.com/biocore/emp/). All data are available as indicated below. We note that in parallel to future sample collection efforts directed by the EMP Consortium, all projects adhering to the EMP standardized protocols for sample collection and sample processing can be analysed using meta-analyses with the data provided here, and all other data generated by following those protocols. Announcements for future sample collection directed by the EMP500 Consortium will be made via https://earthmicrobiome.org.

Reporting summary

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

Source: Ecology - nature.com