Bacterial phylogeny structures soil resistomes across habitats

Abstract

Ancient and diverse antibiotic resistance genes (ARGs) have previously been identified from soil, including genes identical to those in human pathogens. Despite the apparent overlap between soil and clinical resistomes, factors influencing ARG composition in soil and their movement between genomes and habitats remain largely unknown. General metagenome functions often correlate with the underlying structure of bacterial communities. However, ARGs are proposed to be highly mobile, prompting speculation that resistomes may not correlate with phylogenetic signatures or ecological divisions. To investigate these relationships, we performed functional metagenomic selections for resistance to 18 antibiotics from 18 agricultural and grassland soils. The 2,895 ARGs we discovered were mostly new, and represent all major resistance mechanisms. We demonstrate that distinct soil types harbour distinct resistomes, and that the addition of nitrogen fertilizer strongly influenced soil ARG content. Resistome composition also correlated with microbial phylogenetic and taxonomic structure, both across and within soil types. Consistent with this strong correlation, mobility elements (genes responsible for horizontal gene transfer between bacteria such as transposases and integrases) syntenic with ARGs were rare in soil by comparison with sequenced pathogens, suggesting that ARGs may not transfer between soil bacteria as readily as is observed between human pathogens. Together, our results indicate that bacterial community composition is the primary determinant of soil ARG content, challenging previous hypotheses that horizontal gene transfer effectively decouples resistomes from phylogeny.

Figure 1

Functional selections of 18 soil libraries yields diverse ARGs. (a) Bar chart depicting contigs >500bp across all antibiotic selections from CC (red) and KBS (blue) libraries. (b) Amino acid identity between all ORFs (black) or ARGs only (red) and their top hit in NCBI protein. (c) Total ARGs by antibiotic class; y-axis shows number of ORFs, ARG types are in the boxed legend.

Figure 2

Resistance is encoded by diverse soil phyla. (a) Network of predicted bacterial phyla for each ARG. Edge thickness indicates number of ARGs within an ARG family (diamonds) from a predicted phylum (rounded squares). Phyla containing >15 ARGs are labeled, and are shaded dark grey at >3% 16S rRNA abundance. (b) Simplified network of general ARG mechanisms; edge thickness represents significance of phylum and ARG mechanism co-occurrence (Fisher’s exact test, line width increases with ranked significance). Node size indicates number of ARGs (diamonds) or contig count (rounded squares).

Figure 3

Resistomes correlate with phylogeny across soil type and nitrogen amendment. (a to b) PCoA plots depict Bray-Curtis distances between soils, using unique ARG counts. (a) Resistomes from CC (red) and KBS (blue) soils cluster separately (p<10⁻⁵, ANOSIM). Asterisk denotes two soils with near-identical coordinates. (b) CC soils amended with high N-levels cluster separately from other CC soils (p=0.01, ANOSIM). (c to d) Procrustes analyses depict significant correlation between ARG content (Bray-Curtis) and bacterial composition (Bray-Curtis) for (c) CC (red) and KBS (blue) soils and (d) only CC soils.

Figure 4

Pathogen ARGs show higher HGT potential than soil ARGs. (a) Mobility elements syntenic with ARGs are proportionally higher in pathogens than soil genomes or soil functional selections. (*) indicates significance determined from 1000 Monte Carlo simulations; (**) indicates significance determined by Student’s T test. Error bars depict two standard deviations from mean. (b) Pathogens show significantly increased HGT potential relative to soil genomes and soil selections at all distances (20bp intervals) greater than 580bp (dashed line, p < 0.05, Fisher’s exact test). Inset depicts mobility elements encountered within 1.5Kb of ARGs, demonstrating that data from soil selections resembles soil genomes.

Extended Data Figure 1

(a) Results of selections of 18 soil metagenomic libraries for antibiotic resistance to 18 compounds. A dark gray cell means a resistance phenotype was observed whereas white cells indicate the absence of any drug-tolerant transformants. Grassland soils from Cedar Creek (CC) are labeled in red and agricultural soils from Kellogg Biological Station (KBS) in blue. (b to c) Alpha diversity representations. On the left is depicted the number of distinct ARG annotations observed as increasing numbers of ARGs are sampled from each soil. On the right, Shannon diversity scores are shown at each rarefaction step.

Extended Data Figure 2

Three prominent ARG classes are present in nearly all bacterial genomes and can provide antibiotic resistance when overexpressed. (a) Generalized as red circles are dihydrofolate reductases (DHFRs), D-alanine—D-alanine (Dala-Dala) ligases, which are the molecular targets of the drugs trimethoprim (TR) and D-cycloserine (CY) respectively (black stars), and thymidylate synthases (TSs), which can provide trimethoprim resistance by circumventing the need for an active DHFR. When overexpressed in functional selections, these genes can provide antibiotic resistance. We found substantial diversity in these genes (average pairwise amino acid identity 39.3 ± 12.2%), suggesting that variants were captured from many bacterial lineages. (b) Relative to other ARG mechanisms, large numbers of DHFRs, TSs, and Dala-Dala ligases were found in all soils, with these ARGs representing 92.5% of resistance genes identified from selections containing trimethoprim- or D-cycloserine antibiotics. Therefore, these selections encompass large genetic diversity, but constrained functional diversity, with a broad range of genes encoding limited functional traits. (c) When considered in isolation, these functions were not different between the Kellogg Biological Station (KBS) and Cedar Creek (CC) soils (p>0.05, ANOSIM), indicating that trimethoprim and D-cycloserine resistance function is similarly distributed across the surveyed soil types.

Extended Data Figure 3

Total counts of β-lactamases recovered from antibiotic selections using all soils (black), CC soils (red), and KBS soils (blue).

Extended Data Figure 4

Total counts of ARGs categorized by their predicted phylogenetic origin. The number of ORFs are indicated on the y-axis and the ARG types in the boxed legend.

Extended Data Figure 5

Principal coordinate analysis (PCoA) plots of Bray-Curtis distances between soil resistomes. The PCoA was calculated using all ORFs captured from functional selections without trimethoprim- and D-cycloserine, and shows significant separation between CC (red) and KBS (blue) resistomes (p<10⁻⁵, ANOSIM).

Extended Data Figure 6

PCoA across CC (red, grassland) and KBS (blue, agricultural) soils. (a to c) PCoA generated from all 16S data available from ref. , using (a) Bray-Curtis, (b) weighted Unifrac, and (c) unweighted Unifrac dissimilarity metrics. Samples cluster by soil location and N-level, as previously demonstrated. (d to f) The same PCoA plots generated using only samples with sufficient 16S and resistome data (i.e. those used in Procrustes and Mantel analyses). Excluding the two high-N KBS soils with insufficient resistome data eliminates the clustering pattern observed for KBS soils in (a to c). The asterisk denotes the high-N KBS soil common to both sets of analyses.

Extended Data Figure 7

Phylum level relative abundance of combined Cedar Creek (CC) and Kellogg Biological Station (KBS) datasets for major soil bacteria. (a) 16s rRNA data is depicted in black. Phylogenetic inferences based on the sequence composition of the assembled, resistance-conferring DNA fragments are depicted in red. Actinobacteria and Acidobacteria relative abundance represent the largest discrepancies between datasets. (b) Actinobacteria are most dramatically enriched in resistance-conferring DNA fragments, in accord with their role in producing antibiotics, but despite their high GC-content and predicted transcriptional incompatibilities with E. coli. Levels of Proteobacteria, the phylum to which E. coli belongs, are largely unchanged following functional selection, suggesting that any potential bias introduced to the selections by heterologous expression in E. coli is minimal compared to the effect of ARG-content of the source organisms.

Extended Data Figure 8

Procrustes analysis demonstrates that when soils cluster by bacterial composition, resistomes aggregate with phylogenetic groupings. (a to c) Procrustes analysis of the ARG content (Bray-Curtis) of CC (red) and KBS (blue) soils compared to community composition calculated by (a) Bray-Curtis, (b) weighted Unifrac, and (c) unweighted Unifrac dissimilarity metrics. (d to f) The same Procrustes transformations for CC soils only. For a given soil, black lines connect to functional resistome data while the green lines connect to points generated from 16S gene sequence data. The M² fit reported is from a Procrustes transformation over the first two principal coordinates while the p-value is calculated from a distribution of empirically determined M² values over 10,000 Monte Carlo label permutations. For M²/p-values calculated using all principal coordinates, refer to table S8.

Extended Data Figure 9

Procrustes analysis demonstrates that when soils do not form distinct phylogenetic clusters, we are unable to detect significant correlation between ARG content and phylogenetic architecture. See Extended Data figure 6 for the phylogenetic relationships between these soils. (a to c) Procrustes analysis of the ARG content (Bray-Curtis) of KBS (agricultural, blue) soils compared to 16s rRNA gene sequence using (a) unweighted Unifrac, (b) weighted Unifrac, and (c) Bray-Curtis similarity metrics. (d to f) The same Procrustes transformations for the CC soils (grassland, red) without high-N amendment, showing that soil groupings must be distinguishable by bacterial composition to detect correlations with resistome content, regardless of soil type. For a given soil, black lines connect to functional resistome data while the green lines connect to points generated from 16S rRNA gene sequence data. The M² fit reported is from a Procrustes transformation over the first two principal coordinates while the p-value is calculated from a distribution of empirically determined M² values over 10,000 Monte Carlo label permutations.

Extended Data Figure 10

Histogram of nucleotide percent identity from pairwise alignments of all predicted mobility elements, suggesting assembly does not inappropriately condense mobile DNA elements into too few sequences. The blue trace depicts a normal distribution with the same mean and standard deviation empirically observed across all pairwise comparisons.