DNA extraction, 16S rRNA sequencing and data analysis

DNA extraction, 16S rRNA sequencing and analysis pipeline

DNA extraction was performed with the QIAGEN DNeasy PowerSoil HTP 96 Kit following the provided protocol. The obtained DNA was used for 16S amplicon sequencing of the V4 region. Library preparation and sequencing, which was done on an Illumina MiSeq platform.

We used the R package DADA2 to obtain the amplicon sequence variants (ASVs). following the workflow described in Callahan et al. Taxonomic identities were assigned to ASVs using the SILVA version 132 database. 

To monitor the dynamics of the microbial communities, we measured community composition via 16S ribosomal RNA (rRNA) amplicon sequencing. DNA extraction was performed with the QIAGEN DNeasy PowerSoil HTP 96 Kit following the protocol provided by the manufacturer. The obtained DNA was used for 16S (V4 region) amplicon sequencing. Library preparation and Illumina MiSeq sequencing were performed by the Environmental Sample Preparation and Sequencing Facility at Argonne National Laboratory. We used the R package DADA2 to obtain the amplicon sequence variants (ASVs) as described by Callahan et al.(39). Taxonomic identities were assigned to the ASVs by using SILVA (version 132) as a reference database. For each sample. species richness was calculated as the number of ASVs with a relative abundance ≥0.1%, which corresponds to the 0.1% extinction threshold used in simulation (Fig. S1, S2). The phylogenetic tree (Fig. S14) was constructed using Simple Phylogeny (40) by the EMBL’s European Bioinformatics Institute. Taxonomic identities were assigned to ASVs using Randomized Axelerated Maximum Likelihood (RAxML) using default parameters. All plot of relative ASV abundances with stack bars in this paper show the results of one replicate.

Each step in the workflow to assess community composition through amplicon sequencing presents its own biases (41). This includes taxononomical biases in DNA extraction, PCR amplification (e.g., differences in 16S gene copy number), sequencing, and bioinformatics processing. In our study, such biases can significantly compromise the quantitative accuracy of the reported relative abundance of community members, although they are unlikely to significantly compromise our results qualitatively. In particular, these biases could lead to underestimations in community diversity if species fall below the extinction threshold as a result. Similarly, quantitative measures of abundance fluctuations could also be affected. However, these quantitative changes (e.g., under- or overestimation of specific community member) should be comparable across samples, making our qualitative results (e.g., transitions between dynamical modes in a specific order across the phase space) robust to these biases. 

In our sequencing dataset, sequencing depth varied from 3579 to 67354 reads, with an average of 21609 reads. This means that we could not effectively resolve any species abundance on the order of .01% or below. Our main observables, diversity and fluctuation fraction, were calculated (Methods) only from species abundances that exceed a threshold of 0.1% (the extinction threshold). On the one hand, we were able to detect abundances for all the members of each species pool in all the data points for Day 0 (Fig. 2 and S20-S22). Considering that community inoculation consisted in mixing monocultures at equal volumes at Day 0, this suggests that species-specific, sequencing-associated errors are relatively modest in our dataset. On the other hand, our communities are composed of members of a defined set of 48-species library. We did not detect any reads from any ASV that does not correspond to a member of the 48-species library, which suggests an absence of significant contamination during experimental data acquisition and processing and is an additional indicator of reliability of the sequencing data for the purposes of this study.