Increased nutrient diversity can induce loss of microbial diversity through enhanced resource uptake

Media and bacterial culture

Bacteria were precultured from frozen stocks in nutrient medium (NM) composed of 1% yeast extract and 1% soytone (both Sigma-Aldrich) for 24 h at 30 °C in shaking (225 rpm). Following preculturing, cultures were diluted 1:5 in PBS medium and inoculated into M9 medium supplemented with the relevant carbon sources (Supplementary Table 1 gives a full list and vendor details) and 1:100 diluted NM. The NM supplement was included on the basis of preliminary tests (Supplementary Fig. 35), which showed that this small addition substantially increased the number of strains reaching measurable optical density, whereas only a few strains exhibited limited growth on the NM supplement alone in the absence of added carbon sources. This indicates a synergistic effect with defined carbon sources, probably through the alleviation of micronutrient or other abiotic limitations that would otherwise restrict growth. Inoculation was performed using an acoustic liquid handler (Echo 525, Beckman Coulter) by transferring 25 nl of 1:5 diluted culture into 40 µl of media per well (corresponding to a dilution of 1:8,000). For both community and pairwise experiments (Figs. 1 and 2), cultures were passaged daily by diluting 1:27 into fresh media (1.5 µl into 39 µl). Communities were sampled on days 9 and 10 (after eight and nine passages, respectively) and pairwise interactions were sampled on day 5 (after four passages). All experiments were conducted in 384-well plates. Community and pairwise experiments were carried out in transparent plates with lids (VWR), while sentinel strain experiments were performed in black plates with lids (Sarstedt, Lumox model) optimized for luminescence readouts.

Resources and their mixtures

A total of 22 individual resources were used (Supplementary Table 1). These resources were selected to encompass a broad spectrum of microbial growth substrates, including sugars, amino acids, organic acids and alcohol sugars, thereby capturing the main metabolic classes that distinguish ‘sugar specialists’ and ‘acid specialists’ commonly observed among heterotrophic bacteria36. Resources were dispensed into 384-well plates containing M9 media without carbon sources using an acoustic liquid handler (Echo 525, Beckman) before microbial inoculation. For all resource conditions (individual or mixtures), a final concentration of 0.2% (mass/volume) was maintained. We used concentration (mass/volume) rather than molarity to keep the total available energy comparable across resources of different molecular weights, consistent with previous resource-diversity studies17,18,37. Mixture conditions ranged from 2 to 16 resources and were assigned randomly.

Strain library and whole-genome sequencingCommunity and pairwise experiments

For all community and pairwise experiments, we used 14 bacterial strains with distinct 16S rRNA gene sequences, enabling unambiguous taxonomic identification (strain names and 16S sequences listed in Supplementary Data 2). These strains are part of a Caenorhabditis elegans-associated isolate collection of the Ratzke laboratory and represent a subset of the 16 strains previously used in a similar context38.

Strain collection in sentinel strain experiments

For the sentinel strain experiments, we used 298 bacterial strains isolated from the gut of C. elegans collected from diverse natural environments, including rotting fruits, soil and compost (Supplementary Data 5). Strains were provided by C. Ratzke and H. Schulenburg. This collection spans a broad range of bacterial taxa (Supplementary Data 5 and Supplementary Fig. 19) and represents a highly diverse and ecologically relevant set of C. elegans-associated microbes suitable for controlled experimental studies.

Bacterial whole-genome sequencing

Of the 298 strains used in this study, 92 had been previously sequenced and assembled (Supplementary Data 5; courtesy of H. Schulenburg, collected by Johnke et al.39). The remaining 206 strains were sequenced as part of this work. Genomic DNA was extracted using the peqGOLD Bacteria DNA Kit (VWR). Whole-genome libraries were prepared with the NEBNext Ultra II FS DNA Library Prep Kit for Illumina (New England Biolabs), using dual index oligos for barcoding. To reduce reagent use and cost, the protocol was miniaturized to 1:7 of the standard volume. This scaled-down protocol was benchmarked against the standard 1:1 protocol and produced comparable results. Libraries were sequenced on an Illumina NovaSeq platform to a mean depth of ~120×. Raw reads were trimmed using fastp (v.0.23.2)40 and assembled de novo with SPAdes (v.3.15.5)41 using the –isolate flag and default error correction. Assemblies were generated using 24 threads per sample. Assembly quality was assessed with QUAST (v.5.2)42 and metrics were summarized across strains. Contamination was evaluated at the contig level using DIAMOND43 followed by MEGAN44 with GTDB taxonomy. Genomes were excluded if ≥10% of contigs were assigned to taxa inconsistent with the expected lineage, based on contradictions at or above the genus level. Sequencing statistics are in Supplementary Data 6.

Strain phylogeny calculation

Core-genome phylogeny was constructed using panX45, which clusters orthologous genes, aligns core-gene families and infers a maximum-likelihood tree from the concatenated alignment. We used default bacterial parameters and defined core genes as those present in ≥95% of strains. Taxonomy classification was done as described next.

Taxonomic classification of bacterial strains

Taxonomic classification of all bacterial strains was performed on the basis of their assembled genomes using a similarity-based approach against the AnnoTree reference database (v.2021). Each genome was queried with DIAMOND v.2.1.6 (blastx mode) against the AnnoTree protein database to identify the closest reference taxa, using a maximum of 25 target sequences per query and default-sensitive alignment parameters. Resulting alignment files (.daa) were filtered to retain only the top 90% of high-confidence matches. For each strain, the corresponding RefSeq or GenBank assembly accession (GCF_/GCA_) most frequently detected across all hits was selected as the best-supported taxonomic representative. These assembly identifiers were then mapped to National Center for Biotechnology Information (NCBI) taxonomy identifiers (TaxIDs) and organism names using official NCBI assembly summary files from both RefSeq and GenBank. To ensure consistency and completeness of the taxonomic hierarchy, TaxIDs were further resolved into full seven-rank lineages (kingdom, phylum, class, order, family, genus and species) using TaxonKit v.0.16.0 with the most recent NCBI taxdump. The resulting consensus table provides, for each strain, its genome accession, matched assembly reference, NCBI TaxID, organism name and the complete hierarchical taxonomic lineage.

Luciferase tagging of the sentinel strain

P. aeruginosa PAO1 was genomically transformed with the luxCDABE operon using the mini-Tn7 insertion system as previously described38. Briefly, the plasmids pUC18T-mini-Tn7T-Gm-lux (Addgene no. 64953) and pTNS2 (Addgene no. 64968) were curated from Escherichia coli DH5a using a plasmid extraction kit (peqGOLD Plasmid MiniPrep Kit II; VWR). PAO1 cells were made competent for transformation by treatment with a sucrose 300 mM solution and were then electroporated with 50 ng of the curated pUC18T-mini-Tn7T-Gm-lux and pTNS2 plasmids. After recovery in 1 ml of LB medium at 30 °C, the electroporated cells were plated on an LB-agar+Gm30 plate. Plates were incubated for 24 h at 30 °C and colonies were tested for their bioluminescence activity using a plate reader. Positive colonies were then further grown in a selective medium (LB+Gm30) and were stocked in 25% glycerol at −80 °C.

Community composition based on 16S rRNA gene sequencing

Amplicon-based 16S rRNA gene sequencing was used to profile bacterial community composition in pairwise and multispecies experiments.

DNA extraction

Genomic DNA was extracted using a miniaturized protocol in 384-well plates. Pellets were resuspended in TES buffer and 10 µl was transferred to fresh PCR plates. Lysis was initiated by adding metapolyzyme (in PBS) and ReadyLyse (Lucigen, 1:5 in TES), dispensed via an ECHO acoustic liquid handler, followed by overnight incubation at 30 °C with shaking (1,350 rpm). The next day, Proteinase K (New England Biolabs, 1:2 dilution of 20 mg ml−1 stock) and 20% SDS (1% final) were added, followed by incubation at 55 °C for 30 min. RNase A (Omega Bio-Tek, 1:10 dilution) was added and incubated for 5 min at room temperature. Proteinase K was inactivated at 95 °C for 5 min. Samples were centrifuged (2,800g, 10 min) and supernatants were collected. Lysates were diluted 1:6 and then 1:3 in nuclease-free water. DNA concentrations were measured using the PicoGreen assay and 0.1–10 ng was used as template for 16S amplification.

16S library preparation and sequencing

Full-length 16S rRNA genes were amplified using a barcoded PCR protocol with primers 27 F and 1492 R (barcode sequences listed in Supplementary Data 7), followed by sequencing on an Oxford Nanopore platform. PCR-grade water and 2× Hot-Start PCR Master Mix (Biotechrabbit) were combined in ECHO-qualified plates and 4.4 µl of reaction mix was dispensed into 384-well PCR plates using an ECHO acoustic liquid handler (Beckman Coulter). After centrifugation, 50 nl of barcode-specific primers and 500 nl of genomic DNA were added to each well using the same system. Thermocycling was performed as follows: 94 °C for 5 min; 29 cycles of 94 °C for 60 s, 55 °C for 60 s and 72 °C for 2.5 min; followed by 72 °C for 10 min and a final hold at 4 °C. Barcoded amplicons were pooled per tube and libraries were prepared using Oxford Nanopore short-read amplicon workflow according to the manufacturer’s protocol, with halved reaction and bead-cleanup volumes. Sequencing was performed on a MinION device (Oxford Nanopore Technologies).

Read mapping and relative abundance estimation

Raw nanopore reads were split into 100 parts using seqkit46 and filtered by quality (Q ≥ 15) and length (1,500–1,700 base pairs) using NanoFilt (v.2.8.0)47. Demultiplexing was performed with a custom barcode set using minibar (https://github.com/calacademy-research/minibar), allowing edit distances up to six bases. Reverse complementation, where needed, was performed using a custom Python script (rev_cmplt_minibar_output.py). Strain-specific 16S consensus sequences were generated by pooling reads across replicates, filtering (Q ≥ 18) and aligning the top 1,000 reads per strain with MAFFT (v.7.508)48. Consensus sequences were computed using a custom script applying a 90% agreement threshold and IUPAC ambiguity codes and then combined into a FASTA file to serve as a custom reference database. Chimaeric reads were removed using vsearch (v.2.21.1)49 with –uchime_ref against this custom reference database. Non-chimaeric reads were mapped to the database using minimap2 (v.2.26)50 with -ax map-ont. Alignment files were sorted and indexed using samtools (v.1.17)51 and per sample read counts were generated using samtools idxstats. The final count matrix was compiled with a custom script (make_count_table_updated.py). All custom scripts and workflow code are available at https://github.com/orshalevsk/Emergent_Biodiv_loss.

OD600 and bioluminescence measurementOD600

A 384-well plate containing bacterial cultures was mixed for 15 s at 2,000 rpm using a Mixmate plate shaker (Eppendorf), then transferred to a FLUOstar Omega plate reader (BMG Labtech) to measure absorbance at 600 nm.

Bioluminescence

Bacterial cultures were transferred to a Lumox 384-well plate (Sarstedt) for bioluminescence measurement. First, 2.5 µl of NM solution was dispensed into each well. Then, 20 µl of bacterial culture was reverse-dispensed into the same wells using a VIAFLO 384 liquid handler (Integra) and mixed with NM for 30 s at 1,500 rpm using a Mixmate plate shaker (Eppendorf). The Lumox plate was incubated at 30 °C in a FLUOstar Omega plate reader (BMG Labtech) for 2 min to allow NM to stimulate bacterial metabolism (notably in P. aeruginosa), enhancing bioluminescence signal consistency. Bioluminescence was then recorded using the luminescence emission filter of the plate reader, with a gain setting of 3,600 and an exposure time of 1 s.

Experimental designCommunities, pairwise and individual strains experiments

Precultures were prepared as described above (section ‘Media and bacterial culture’), then dispensed into 384-well plates containing M9 media with the appropriate carbon sources and 1:100 NM. Metadata for all conditions and strain combinations used in the community experiment are provided in Supplementary Data 1; metadata for pairwise and single-strain experiments are provided in Supplementary Data 3. All experimental conditions were fully randomized across wells and plates. At the end of each experiment, bacterial density (OD600) was measured for all wells, including individual strains, pairwise combinations and communities. Only pairwise and community samples were further processed for 16S sequencing: samples were centrifuged at 2,800g for 5 min, supernatants were removed and pellets were stored at −80 °C. DNA extraction and 16S rRNA gene sequencing procedures are described below (DNA extraction and 16S sequencing). Individual strain growth was measured after 72 h (carrying capacity, OD600), pairwise cultures were serially passaged for 5 days, followed by a final 24-h incubation before OD600 measurement and 16S profiling, and communities were serially passaged for 9 days, followed by a final 24-h incubation before OD600 measurement. The impact factor was computed from experimental OD600 data as an empirical measure of total inhibition acting on the sentinel strain. It was defined as the relative reduction in OD600 compared with monoculture controls (Fig. 2a). In the framework of generalized Lotka–Volterra model this impact factor corresponds to the steady-state population density of species 2 in the presence of species 1 divided by the steady-state population density of species 2 grown alone, which equals \({\alpha }_{22}\left({\alpha }_{11}-{\alpha }_{21}\right)/\left({\alpha }_{11}{\alpha }_{22}-{\alpha }_{12}{\alpha }_{21}\right)\). Accordingly, the impact factor depends on all four alphas of a pairwise interaction.

Sentinel strain selection rationale

We selected P. aeruginosa PAO1 as the sentinel because it supports stable, single-copy chromosomal integration of the luxCDABE reporter via the mini-Tn7 system—ensuring consistent reporter dosage, avoiding plasmid-borne variability and yielding a strong, reliable bioluminescence signal compatible with our plate-reader workflow. PAO1 grows robustly at 30 °C—matching the growth range of our strain collection—and exhibits broad metabolic niche breadth (generalist physiology), so its growth is sensitive to a wide range of changes in metabolic space. These practical and biological properties make PAO1 a reliable readout strain for quantifying competitive balance across diverse resource conditions. We note that any single sentinel may introduce bias.

Sentinel strain experiments

The full experimental setup for sentinel strain co-cultures is described above (sections ‘Media and Bacterial Culture‘ and ‘Resources and their mixtures’). Briefly, the sentinel strain and 298 co-cultured strains were grown overnight in NM medium and diluted 1:5 in PBS (without washing). Cultures were then dispensed into 384-well plates containing the appropriate media and resource conditions. Six experimental batches were performed: three using individual resources (biological replicates) and three using distinct mixed-resource conditions (no replicates). In the individual-resource experiments, all 298 strains were co-incubated with the sentinel strain across 22 single-resource conditions. The sentinel strain was also grown alone with 16 replicates per condition. In the mixed-resource experiments, one batch included all 298 strains co-incubated with the sentinel strain in all conditions. In the other two batches, 80 strains were randomly selected per condition to maximize the number of unique resource combinations tested. In these batches, the sentinel strain alone was included with eight replicates per condition. All conditions were fully randomized across wells and plates, with each experiment distributed over 20–24 separate 384-well plates per batch. Plates were incubated at 30 °C for 24 h before luminescence was measured. The average luminescence of the sentinel strain grown alone in each condition was used to calculate the impact factor of co-cultured strains, as shown in Figs. 2a and 3a. Impact factor values for individual resources represent the average across three replicate experiments (‘5’, ‘6’ and ‘7’). Full experimental conditions and metadata are provided in Supplementary Data 4.

Supernatant experiments

The full setup for supernatant experiments is described above (sections ‘Media and Bacterial Culture‘ and ‘Resources and their mixtures’). Briefly, selected highly influenced and mildly influenced strains were grown overnight in NM medium and diluted 1:5 in PBS (without washing). Cultures were dispensed into 384-well plates containing M9 medium supplemented with either individual- or mixed-carbon sources, specifically: sucrose, glutamic acid, citric acid, maltose, mannose, xylose, lysine, ribose, proline, fructose, histidine and asparagine. Each strain inoculated alone into four adjacent wells per condition (for example, A1, A2, B1 and B2) to ensure sufficient volume for supernatant extraction. Plates were incubated at 30 °C for 24 h. To extract supernatants, technical replicates (four wells per condition) were pooled into 96-well plates using a VIAFLO 384 liquid handler (Integra). Each set of four wells (40 µl per well) yielded a combined volume of 160 µl. The pooled cultures were filtered through 0.2-µm filter plates (PALL) into a new sterile 96-well plate, removing cells and leaving only the cell-free supernatant. This design enabled combinatorial mixing of supernatants to reconstruct mixed-carbon environments. Mixing was performed using a VIAFLO 96 liquid handler by combining supernatants from wells on different plates (for example, A1 from plate 1 with carbon A and A1 from plate 2 with carbon B were combined 1:1 to generate A + B mixtures). For the sentinel growth assay (Fig. 4c), eight strains were randomly selected—four highly influenced and four mildly influenced—were tested across 66 mixed-resource conditions. For the MS experiment (Fig. 4b), a subset of four of these strains (two highly influenced and two mildly influenced) was tested across nine mixed-resource conditions. Full strain lists and conditions are provided in Supplementary Data 8.

GC–MS sample preparation and analysis

A total of 72 samples (20 µl each) were thawed on ice and mixed with 50 µl of internal standard solution containing 60 µM 13 C6-glucose and 60 µM 3-hydroxybenzoic acid (equivalent to 3,000 pmol each). Samples were dried overnight using a vacuum concentrator (Eppendorf Concentrator) to remove water, then placed in a desiccator over phosphorus pentoxide for 1 h. Derivatization was performed following a protocol adapted from ref. 52. Briefly, 50 µl of methoxylamine hydrochloride in pyridine (20 mg ml−1, freshly prepared) was added to each sample. Samples were incubated in an ultrasonic bath for 10 min, followed by shaking at 30 °C for 90 min at 1,400 rpm. After a brief centrifugation, 70 µl of N-methyl-N-trimethylsilyltrifluoroacetamide (Sigma-Aldrich) was added. Samples were then incubated at 40 °C for 60 min at 1,200 rpm and left at room temperature for an additional 2 h. Final centrifugation was performed at 15,000g for 10 min at 4 °C (Hettich tabletop centrifuge, swing-out rotor) and 60 µl of supernatant was transferred to GC vials for injection. GC–MS analysis was carried out on a Shimadzu TQ 8040 triple quadrupole mass spectrometer coupled to a high-performance liquid chromatography system. Chromatographic separation was achieved on a Restek Rxi-5Sil MS column. Detailed GC and MS parameters are provided in Supplementary Data 9. Data acquisition and peak integration were performed using LabSolutions Insight GC–MS software, with manual curation of peak boundaries.

EMP data handlingAlpha diversity versus resource complexity analysis

Data were retrieved from the EMP study by ref. 19 via Qiita (study ID 13114: https://qiita.ucsd.edu/study/description/13114). We used both 16S-based alpha diversity and predicted microbially derived metabolite richness, as previously computed in ref. 19 (Fig. 2c and Supplementary Fig. 1c,d). Alpha diversity was taken from the column alpha_16s_deblur_nosingletons_rar5k_shannon, which reports Shannon diversity after singleton removal and rarefaction to 5,000 sequences per sample. This depth was identified by the authors as optimal for cross-sample comparison. Resource diversity was taken from alpha_lcms_fbmn_microbial_richness, defined as the number of HPLC–MS peaks >0 in each sample following feature-based molecular networking and interpreted as microbially associated metabolite richness. Samples were filtered to include only those with non-missing alpha diversity values (477 of 618 total). To ensure statistical robustness, we retained only environments (as defined by EMP level empo_4) with ≥6 samples, resulting in 462 samples across 13 distinct ecosystems. We used the empo_4 level of environmental classification, which provides the most specific ecosystem labels available in the EMP dataset and avoids confounding that can arise at broader levels (for example, empo_1), where ecologically distinct environments—such as animal guts and plants—may be grouped together.

Mapping metagenomic reads to highly and mildly influenced strains and abundance estimation

To assess the environmental representation of highly and mildly influenced strains, we mapped metagenomic shotgun reads from the EMP19 to our full genome collection of 298 experimentally characterized strains. Genome assemblies were merged into a single reference database and a minimap2 index (v.2.26)50 was created. Reads were aligned using minimap2 with the -ax sr –secondary=no parameters to exclude secondary alignments (pipeline: minimap2_pipline_no_multi_match.sh). Alignment files (SAM) were processed using a custom script (make_count_table_updated_with_mistmatch_filter_db.py) to filter for high-quality matches (–min_mapq 50) and compute per strain read counts, along with unmapped reads (Supplementary Data 10). To ensure sufficient coverage and minimize noise, only samples with a total mapped read count ≥10,000 were retained. This filtering yielded 453 metagenomes spanning 18 distinct ecosystems. For each sample, raw counts were normalized to relative abundances. Strains were categorized as highly or mildly influenced based on their experimentally measured impact on biodiversity and reads were binned accordingly into highly influenced, mildly influenced or unmapped groups. Relative abundances for each category were then calculated (Supplementary Data 11). These values were integrated with matched metabolomic profiles from the EMP19 for downstream ecological analyses. All custom scripts used in this workflow are available at https://github.com/orshalevsk/Emergent_Biodiv_loss.

Consumer–resource model

For the classic consumer–resource model, we consider a community of 20 microbial species competing for 30 resources in a continuous dilution (chemostat-like) environment. Species take up available resources and either turn them into biomass or transform them into secreted metabolic byproducts (cross-feeding) and die through dilution. Resources are depleted through species uptake, replenished through secretion from species and are subject to dilution, which replenishes supplied resources and washes out others. For each of 20 independently parameterized communities, we varied resource complexity across 1, 2, 4, 8, 16 supplied resources (equal total supply split evenly among the chosen resources). Initial abundances were uniform (N0 = 1/20 per species). Simulations used scipy.integrate.solve_ivp v.1.16.3 in Python v.3.14.0 as ODE solver (rtol=1e-3, atol=1e-6, method = ”RK45”), over t = 0–800 (a.u.) with dilution of 0.1. Cross-feeding is modelled as a metabolic transformation tensor (bᵢ,β → α) which allows a fraction of consumed resource β to be secreted as resource α by species i, thereby generating new metabolites available to other species. This formulation follows standard extensions of MacArthur-type models incorporating metabolic byproduct exchange and saturating resource uptake kinetics8,17,53. For full details, see Supplementary Text 1, Section 1 and the full code (under ‘CR_model_simulations’).

For the highly influenced strains modified model version see Supplementary Text 1, Section 2. In this updated model, we introduce a species-specific sensitivity parameter γ. When γ > 0, the species increases its resource uptake affinity (reduces Kₘ) as resource diversity increases. We set γ = 1 for highly influenced strains (HIS) and γ = 0 for mildly influenced strains (MIS) by default. For each of the 20 independently parameterized communities, we vary the HIS fraction f = 0, 0.05, …, 1.0 (14 fractions) (Fig. 4e and Supplementary Figs. 28–34). To assign the species as either HIS or MIS, we first assigned a tendency to be highly influenced for each species i, which is correlated with resource generalism. We then rank the species by their tendency to be highly influenced and assign the ones with highest values as HIS according to the preassigned fraction f and the rest as MIS. We simulated the 20 communities under 14 HIS fractions, across the same five resource complexities (1–16) with equal total supply, t = 0–800, dilution of 0.1 and uniform initial abundances.

We further performed robustness tests by systematically varying key parameters and assumptions. These include: (1) cross-feeding rate; (2) correlation between highly influenced behaviour and resource generalism; (3) γ value for HISs; (4) a continuous instead of binary distribution of MIS and HIS behaviour; (5) an alternative mechanism where resource complexity influences HISs by increasing assimilation efficiency instead of resource affinity; and (6) the number of initial species for each community. For full details see Supplementary Text 1, Section 3 and Supplementary Figs. 30–34 and the full code (https://github.com/orshalevsk/Emergent_Biodiv_loss/).

Inferring interaction matrices from experimental data using the generalized Lotka–Volterra model

Data from single-species growth (Supplementary Fig. 9), pairwise interactions (Fig. 2a) and community compositions (Fig. 1d,e and Supplementary Fig. 3) at steady state were used to infer Lotka–Volterra interaction matrices of the 14 strains comprising these communities. Interaction matrices were constructed for three datasets corresponding to 1, 8 and 16 available carbon sources. Since each of these three datasets contains data for several different nutrient compositions, we treated interaction strength as a random variable by using Bayesian regression to infer their values.

The steady state of the generalized Lotka–Volterra system is obtained by setting equation (1) to zero, taking into account only those species whose population densities remain non-zero. The obtained equations were fit to the data with Bayesian regression (Monte Carlo Markov Chain sampling, NUTS solver in pyro, python). As priors we used Gaussian distributions with a mean of 1 for the off-diagonal entries and a mean equal to the inverse carrying capacity for the diagonal entries. This is based on the fact that the self-interaction coefficient (α_self) corresponds to the inverse of the carrying capacity (K), a standard notation of generalized Lotka–Volterra models32,54. Posterior sampling was performed with 30 chains, 400 warm-up steps and 200 sampling steps. For full details, see Supplementary Text 2. Note that error bars for αself in Fig. 2b (middle) may be overestimated as the interaction values are not independent (Supplementary Fig. 6a).

Estimating single-species Shapley contributions to biodiversity loss

We used the interaction matrices obtained from experimental data for 1 and 16 carbon sources (methodology described in Supplementary Text 2; results in Supplementary Fig. 4) and added Gaussian noise (mean of 0, variance of 0.7) to each entry to generate a synthetic set of 14 species resembling the 14 experimentally measured strains. We produced in total 800 such sets of 14 species and took from each set randomly six or eight species, respectively (400× six species and 400× eight species). We then simulated all possible subcommunities of these six or eight species using the generalized Lotka–Volterra model under both 1- and 16-resource conditions. Since the interaction matrices were obtained from steady-state data, we could not estimate the per capita growth rates of the single species and set them to be one instead. Since the growth rate may decide for multistable cases in which state the system ends up, this choice may result in different steady states compared with experimentally parametrized or randomized growth rates. (Supplementary Fig. 10). For each simulated subcommunity, we computed the steady-state Shannon diversity index under both conditions and calculated the change in diversity as: ΔShannon = Shannon16 − Shannon1. We then computed Shapley values, representing the contribution of each strain in a six- or eight-species community to the overall change in Shannon diversity between 1 and 16 resources. Next, using the same interaction matrices (Supplementary Fig. 4), we calculated the mean impact factor for each species on the other members of its community under both resource conditions—analogous to the analysis shown in Fig. 2d. Finally, we computed the change in mean impact factor with increasing resource diversity and plotted these values against the corresponding Shapley values for each species (Fig. 2c).

Reporting summary

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