Clusters of Antibiotic Resistance Genes Enriched Together Stay Together in Swine Agriculture.

UNLABELLED: Antibiotic resistance is a worldwide health risk, but the influence of animal agriculture on the genetic context and enrichment of individual antibiotic resistance alleles remains unclear. Using quantitative PCR followed by amplicon sequencing, we quantified and sequenced 44 genes related to antibiotic resistance, mobile genetic elements, and bacterial phylogeny in microbiomes from U.S. laboratory swine and from swine farms from three Chinese regions. We identified highly abundant resistance clusters: groups of resistance and mobile genetic element alleles that cooccur. For example, the abundance of genes conferring resistance to six classes of antibiotics together with class 1 integrase and the abundance of IS6100-type transposons in three Chinese regions are directly correlated. These resistance cluster genes likely colocalize in microbial genomes in the farms. Resistance cluster alleles were dramatically enriched (up to 1 to 10% as abundant as 16S rRNA) and indicate that multidrug-resistant bacteria are likely the norm rather than an exception in these communities. This enrichment largely occurred independently of phylogenetic composition; thus, resistance clusters are likely present in many bacterial taxa. Furthermore, resistance clusters contain resistance genes that confer resistance to antibiotics independently of their particular use on the farms. Selection for these clusters is likely due to the use of only a subset of the broad range of chemicals to which the clusters confer resistance. The scale of animal agriculture and its wastes, the enrichment and horizontal gene transfer potential of the clusters, and the vicinity of large human populations suggest that managing this resistance reservoir is important for minimizing human risk. IMPORTANCE: Agricultural antibiotic use results in clusters of cooccurring resistance genes that together confer resistance to multiple antibiotics. The use of a single antibiotic could select for an entire suite of resistance genes if they are genetically linked. No links to bacterial membership were observed for these clusters of resistance genes. These findings urge deeper understanding of colocalization of resistance genes and mobile genetic elements in resistance islands and their distribution throughout antibiotic-exposed microbiomes. As governments seek to combat the rise in antibiotic resistance, a balance is sought between ensuring proper animal health and welfare and preserving medically important antibiotics for therapeutic use. Metagenomic and genomic monitoring will be critical to determine if resistance genes can be reduced in animal microbiomes, or if these gene clusters will continue to be coselected by antibiotics not deemed medically important for human health but used for growth promotion or by medically important antibiotics used therapeutically.

INTRODUCTION

Antibiotic resistance is recognized as a worldwide health threat (1–5), especially from Gram-negative pathogens (6). Antibiotic resistance is an ancient trait (7) that has coevolved with natural antibiotic production (8) to result in widespread resistance in nature (9, 10). Humanity’s use of thousands of tons of antibiotics annually (11) has provided a selective advantage for resistant bacteria to flourish in the clinic (2, 12) and the environment (13). Subinhibitory antibiotic concentrations are sufficient to enrich for resistant bacteria (14), and mixtures of subinhibitory doses of antibiotics, such as are common in agriculture, select for the highest level of resistance (15). The fight against infectious disease is approaching a crisis situation (16), especially considering that while antibiotic use (12) and resistance increase, the process of antibiotic discovery has nearly halted (6).

Governments are increasingly regulating the use of antibiotics in agriculture. The European Union (17), led by Denmark, has banned the use of antibiotics as growth promoters (18). In the United States, the use of antibiotics important for human medicine will cease for animal growth promotion after 2016, and veterinary oversight will be required for therapeutic uses (5). Observations that agricultural antibiotics increase the abundance of resistance genes in manure (19, 20) and manure-amended soil (21, 22) and potentially in the general public (23) provide a clear scientific basis for the regulation of agricultural uses of antibiotics. Compliance with new regulations will likely require alternative management practices (24) to overcome the losses of antibiotic growth-promoting and disease-controlling benefits (25, 26) and to maintain agricultural productivity (18).

Horizontal gene transfer (HGT) can uncouple gene content and phylogeny (27) and allows for rapid (28) acquisition of resistance or multidrug resistance (29), in some cases from environmental to pathogenic bacteria (30, 31). Genetic elements mobilizing genes within a cell (integrons and transposons) and between cells (plasmids and phage) are known to be enriched in agriculture due to the use of antibiotics (19, 32–35). We will refer to both types of elements as mobile genetic elements (MGEs). Insertion sequences confer a DNA cut-and-paste functionality at inverted repeat sequences (36). Integrons allow for the recruitment, genetic integration, and promotion of exogenous genes as integron cassettes (34). Class 1 integrons in clinical strains are thought to have all arisen from a single source due to widespread occurrence of identical intI1 sequences (34, 37). Both of these within-cell mobile genetic elements can lead to the evolution of multidrug resistance plasmids (28, 38), which can be transferred horizontally by conjugation (27).

In this study, we identify high correlations in the cooccurrence of clusters of identical antibiotic resistance genes (ARGs) and mobile genetic element sequences in Chinese swine farms and farm-impacted soils as well as U.S. laboratory swine. Despite consistency of detection and abundance of ARGs across sample types, the phylogenetic membership of these communities is distinct between farms and sample types. These results provide a strong impetus to further study the effects that swine agriculture has in possibly developing multidrug resistance islands and plasmids and their dispersal throughout these communities. This is the first study to perform high-throughput quantitative PCR (qPCR) followed by sequencing of all amplicon products to describe the effects of antibiotic use in animal agriculture.

RESULTS

Samples and DNA extraction.

We continued the investigation of two classes of swine manure samples that were collected and have been described and studied previously by our laboratories (19, 20). Briefly, we sampled individual swine feces from the USDA National Animal Disease Center (NADC), Ames, IA, from swine living in highly controlled environments, sampled immediately preceding (0 days) and 14 days following continued antibiotic-free feeding (not medicated [NM], n = 3) or introduction to in-feed antibiotics (medicated [M], n = 3), which contained a mixture of chlortetracycline, sulfamethazine, and penicillin (ASP250). These sample names were abbreviated with the treatment group (NM or M) and the day sampled (0D or 14D). The other class of samples comes from industrial-size swine farms in north (Beijing [B]), central (Jiaxing [J]), and southern (Putian [P]) China as well as a pristine Chinese soil (CS). At each of these farms, the manure piles (M), composted manure (C), and soil amended with composted manure (S) were each sampled in triplicate. Samples were given a two-letter abbreviation: the first is the city (B, J, or P) and the second is the sample type (M, C, or S). In addition, we obtained feral pig feces (n = 4) as a reference sample (39) for antibiotic-free swine without contact with modern agricultural systems. These were abbreviated F (feral) with a pig identification (ID) number.

Cooccurrence of ARG and MGE sequences.

Among the Chinese farm sites, despite their distant geographic locations or sample types, the abundance and occurrence of the most abundant allele of multiple genes were significantly correlated. For example, the abundance of the most prevalent allele of intI1 compared with that of qacEΔ1 and IS6100 resulted in r2 values of 0.92 and 0.87, respectively. When we perform Spearman correlations with all the complete pairs of genes from the Chinese farms, we observe extensive clusters of statistically significant abundance-correlated alleles (Fig. 1A). In general, there are two major clusters of genes that cooccur in the Chinese farms, with little cross-correlation between the two clusters. The first cluster is composed of intI1, qacEΔ1, IS6100, aminoglycoside phosphotransferases and nucleotidyltransferases, tetracycline efflux, sul2, dfrA1, cmlA1, and incompatibility group IncW plasmids. We designate this group of genes the intI1-IS6100 cluster (Fig. 1A, left cluster). The average correlation coefficient (ρ) within the cluster of significant cooccurrences is 0.86, showing very high between-gene correlations.

FIG 1 

Cooccurrence network with primer names as node labels of the ΔΔAC values of the most abundant allele of each antibiotic resistance gene and mobile genetic element together with genus-level classification of 16S rRNA sequences within all Chinese farm samples (n = 25; JS1 and JS2 were excluded) (A) and all NADC pigs (n = 12) (B). Nodes connected by a line have a statistically significant Spearman correlation and are cooccurring. Various requirements were required to call two alleles cooccurring: codetection in at least half the samples (for cooccurrence to a genus, this requirement was relaxed to n > 4), false-discovery correction q value of <0.05, and ρ of >0.75. Node size is dependent on the number of connections to other nodes (degree). Shaded circles were added post hoc to clusters of alleles that have high degrees of abundance among all members of that cluster and limited abundance outside the cluster. (C to E) Representative correlations from the intI1-IS6100, IS1216-tet, and blaTEM-sul2 resistance clusters, respectively. The intI1-IS6100 cluster is enriched and diluted differently in all farms, while the IS1216-tet cluster steadily declines from manure to soil in all farms. Note that the ΔΔAC detection limit for ARGs was −15. Axis labels are followed by _1 to indicate that they are the most abundant allele of the amplicons obtained.

Tetracycline ribosomal protection protein resistance genes and two transposase alleles (IS614 and IS1216) dominate a second abundance-correlation cluster. We designate this group of genes the IS1216-tet cluster (Fig. 1A, right cluster). Also included in this cluster is mefA (macrolide efflux pump) and aphA3, an aminoglycoside phosphotransferase gene. Throughout the manure management process, i.e., from fresh manure to composted manure to soils receiving composted manure, the IS1216-tet cluster steadily declines, while the abundance of the intI1-IS6100 cluster is much more dynamic between sample types (Fig. 1C and D; see also Fig. S1 in the supplemental material). The sequences from the IS1216-tet cluster are most commonly associated with Gram-positive bacteria, while the intI1-IS6100 alleles are associated with Gram-negative bacteria as determined using the Basic Local Alignment Search Tool (BLASTn) with the nr database from the National Center for Biotechnology Information (NCBI).

The NADC research animals exhibited a different set of cooccurring ARGs and MGEs than the Chinese pigs: the blaTEM-sul2 cluster (Fig. 1B). The gene cooccurring with the greatest number of other genes within the cluster was blaTEM, and it was shown to cooccur with nearly every other gene of the cluster, including aminoglycoside phosphorylases (strA and strB), nucleotidyltransferases (aadA), sulfonamide resistance (sul2), and multiple mobile genetic elements (intI2, IS6100, IS4, and incN plasmids). The average ρ within the blaTEM cluster was 0.79. Tetracycline resistance (tetO and tetW) clustered separately. Figure 1A and B show only the most abundant allele of each gene. We also observed high intragene cooccurrence (cooccurrence between alleles of the same gene), especially within the Chinese samples (see Fig. 4A and B), when more alleles were included in the analysis.

FIG 4 

(A and B) Average percent abundance of the five most abundant alleles (100% identity OTUs) from all primer sets in the Chinese farms (A) and the NADC pigs (B). Allele abundance was determined for all samples combined (NADC and China). For this reason, the same allele number in both the Chinese and NADC samples represents the same allele. For a detailed heat map of the abundance of the 40 most abundant alleles in all samples individually, see Fig. S4 in the supplemental material. The genes are organized based on clustering structure of the data in the cooccurrence network (Fig. 1). Numbers on the bars indicate the number of intergene edges (number of edges between the top five alleles of each gene; number of maximum intergene edges possible = 20) when the network analysis is performed with the top five alleles rather than what is shown with only the most abundant allele. (C) Shannon diversity indices for both China and NADC samples.

Enrichment of ARG and MGE allele sequences.

Individual resistance gene alleles (sequences) were enriched with a relative change in the adjusted cycle threshold (ΔΔAC) up to 15, which is equivalent to about 32,000-fold enrichment compared to the pristine soil. The level of enrichment of the resistance cluster gene alleles is unified across sample types and farms (see Fig. S1 in the supplemental material). Interestingly, in the Beijing farm, enrichment in manure is highest, followed by compost and soil, while in the Putian farm, some gene alleles follow this same trend (blue background in the figure), while for others (green background) enrichment in the manure is lowest and compost is the highest (see Fig. S1). These genes make up a large portion of the community; for example, the most abundant intI1 allele was 5.7% as abundant as 16S rRNA in the BM, PC, and JC samples and 0.16% as abundant in the soil across all the farms but 0.93% as abundant in the PS sample. The aphA3 allele was the most abundant gene, up to 34% as abundant as 16S rRNA. This allele-specific analysis aided in the detection of enrichment of the most abundant blaTEM allele. When only the results of qPCR are used, it appears that the abundance of blaTEM is decreased (see Fig. S2). However, this is due to the presence of a different blaTEM allele in the control soil, and a different allele was enriched in both the NADC animals and the Chinese farms.

ARG-MGE clusters in sequenced genomes.

We searched known sequenced genomes to ascertain if there was a precedent for resistance cluster alleles that we observed (Fig. 1) to cooccur in individual microbial genomes, and if the genes colocalize into resistance islands. We found supporting evidence for both cases. Genes from the intI1-IS6100 gene cluster were found to colocalize most often in plasmids or resistance islands from Gram-negative organisms at 100% identity to the most abundant alleles that we obtained (Table 1). Acinetobacter baumannii strain AYE contains an extensive 86-kb resistance island, likely obtained through horizontal gene transfer (40). Amplicons from our study map with 100% identity along the entire amplicon length for 18 of the genes in this resistance island. Acinetobacter was found to correlate with the abundance of the intI1-IS6100 resistance cluster in manure (Fig. 2C). The colocalization of intI1 and IS6100, together with resistance genes from our survey, was identified in multiple Escherichia coli plasmids as well as a Klebsiella pneumoniae plasmid (see Fig. S3 in the supplemental material) that also encodes the BlaNDM-1 protein in another of its plasmids. These strains also carry up to 11 additional genes from the intI1-IS6100 cluster mapping at least 50 kb away from the identified intI1 and IS6100 genes.

TABLE 1 

Obtained amplicons align with 100% identity in colocalized groups within known sequences

Example no.	Species and strain (accession no.)	Island or plasmid	Amplicon	Location (kb)
1	Acinetobacter baumannii strain AYE (CU459141.1)	Resistance island	qacEΔ1	3621.2
			aadA	3622.0
			intI1	3624.9
			IS26	3626.6
			aphA1	3627.6
			IS26	3628.9
			IS26	3633.2
			tetR	3643.7
			IS6100	3649.0
			qacEΔ1	3651.7
			qacEΔ1	3655.8
			aadA	3656.6
			strA	3658.1
			tetG	3672.3
			qacEΔ1	3676.1
			dfrA1	3677.3
			intI1	3677.9
			strB	3679.4
2	Proteus mirabilis strain PmGUE (JX121641.1)	Genomic island	intI1	26.3
			qacEΔ1	28.4
			IS6100	30.4
			tetR	35.9
			strB	42.1
			IS26	45.3
			blaTEM	46.7
			IS26	48.4
			aphA1	49.4
			IS26	50.3
			intI1	51.7
			qacEΔ1	53.0
			tetG	56.9
			qacEΔ1	62.5
			IS6100	65.2
3	Salmonella enterica strain ST06-53 (KT334335.1)	Plasmid pHK0653	IS26	72.1
			IS26	78.2
			IS26	80.6
			IS26	84.3
			aadA2	85.0
			IS26	87.0
			aphA1	87.9
			IS26	88.8
			aadA	94.1
			cmlA1	94.6
			aadA2	96.4
			intI1	98.9
			Tn21	102.8
			sul2	106.9
			IS4	110.6
			IS26	114.5
			qacEΔ1	118.2
			IS26	120.1
4	Klebsiella pneumoniae strain JM45 (CP006657.1)	Plasmid p1	Tn21	41.9
			incN	51.0b
			IS26	63.9b
			IS6100	67.9
			intI1	69.8
			dfrA2	70.5
			aadA2	71.6
			qacEΔ1	72.3
			tetRb	76.9
			Tn23	129.4
			IS26	140.5
			blaTEM	141.9
			tetG	148.0
			qacEΔ1	151.8
			IS6100	154.0
			IS26	154.8
			aphA1	155.8
			IS26	156.7
			IS26	158.6
5	Escherichia coli (NC_010558.1)	Plasmid pIP1206	IS26	90.9
			tetB	96.1
			Tn21	101.5
			IS26	104.8
			qacEΔ1	107.9
			aadA5	108.4
			intI1	110.1
			IS26	111.5
			IS26	112.9
			blaTEM	115.8
			intI1	123.0
			IS26	123.7
			IS26b	135.0
6	Escherichia coli strain 2009C-3133 (CP013027.1)	Plasmid unnamed3	sul2	123.1
			IS26	125.5
			strAb	137.1
			aadA	138.2
			IS4	142.4
			IS26	144.7
			blaTEM	148.1
			IS26	148.9
			IS26b	172.7
			blaTEM	173.4

The most abundant of each allele of each primer was included in the query.

Indicates a single mismatch between the most abundant amplicon and its local alignment with the NCBI sequence.

Six examples from the NCBI database of the obtained amplicons aligning in clusters within genomes. Examples 1 to 5 include genes from the intI1-IS6100 cluster and some intercluster genes, while example 6 includes genes from the blaTEM-sul2 cluster (IS26, while found in the sequence, was not cooccurring with this cluster in the swine samples). Co-localization of genes within 10 kb from another alignment position is common as indicated by the differences in genomic locations. Sequences from the NCBI database with the highest numbers of total alignments are shown. Examples from the IS1216-tet cluster are discussed in the text.

FIG 2 

Correlations between the ΔΔAC values of the most abundant ARG alleles and representative cooccurring phylogenetic groups. (A) Lactobacillus cooccurs with six of the eight genes in the IS1216-tet cluster (Fig. 1); however, Lactobacillus was present in all nine manure samples but was detected in only one of the 16 soil and compost samples. (B) Unclassified Xanthomonadaceae sequences cooccurred with two genes (Fig. 1) within the intI1-IS6100 cluster and were detected in 10 of 25 samples. (C) Acinetobacter was detected in 5 manure samples and showed high cooccurrence with the intI1-IS6100 cluster but only within manure samples. Note that the ΔΔAC detection limit for taxa was 7.8 and that the one for ARGs was −15. Gene names in the axis label are followed by _1 to indicate that they are the most abundant allele.

Genes from the IS1216-tet cluster were not observed together in known genomes as frequently as those of the intI1-IS6100 cluster. Only pairs of these genes were frequently observed in Gram-positive organisms. In one instance, colocalization of four of the genes (mefA, aphA3, tetQ, and IS614) was observed in a transposon of Bacteroides fragilis (accession number KJ816753), a Gram-negative organism, but these genes all came from a region of the transposon which was horizontally transferred, likely from Gram-positive organisms (41).

Community analysis using V4 16S rRNA gene amplicons.

Bacterial communities from all the samples clustered by sample type but with a moderate level of heterogeneity between replicates (Fig. 3A). The NM (a composite of NM0D, M0D, and NM14D), manure, compost, and soil sample cluster centers were statistically distinct based on an analysis of molecular variance (AMOVA) Bonferroni-corrected P value of <0.01 for all comparisons between these samples. The stress calculated in this ordination was 0.24, which indicates only a fair representation of the data based on the two axes. The feral samples did not cluster distinctly from any of the other sample types.

FIG 3 

Taxonomic analysis based on sequences of the V4 region of the 16S rRNA gene. (A) Nonmetric multidimensional scaling ordination plot of the 16S rRNA data. The stress calculated in this ordination was 0.24, which indicated only a fair representation of the data based on the two axes. Colored ellipses were added to indicate sample type separation. (B) Average percent abundance of phylum-level-classified 16S rRNA sequences. In all panels, borders are included on the six most abundant phyla only for clarity in matching the legend. (C and D) Average number of sequences of genus-level sequences within Gammaproteobacteria and Bacilli, respectively. The data were subsampled to 967 sequences for all samples to normalize the number of sequences. Two Jiaxing soil replicates (JS1 and JS2) and a feral pig (F6) did not have sufficient numbers of sequences and thus were excluded from this analysis. Nonitalicized taxonomic names are unclassified sequences within that group. Only the 12 most abundant genera are shown.

All manure samples, independent of their treatment, were dominated by the Firmicutes and Proteobacteria phyla, and soil samples were dominated by Proteobacteria and Acidobacteria (Fig. 3B). The Access Array platform required primer annealing at 60°C for all reactions. This annealing temperature limited amplification of members of the Bacteroidetes phylum (see Text S1 in the supplemental material), which was expected to be a major phylum in manure samples.

While there are taxonomic differences between individual samples across sample types, noise within sample types resulted in only the few following statistical differences (q value, <0.05). Soil samples were enriched in Acidobacteria and decreased in Firmicutes compared with the NM, manure, and compost samples. The NM samples were enriched in Spirochaetes (Fig. 3B) compared with the Chinese samples. At the genus classification level, there is even greater variation within and between sample types, but some statistically significant shifts were observed. Within Gammaproteobacteria (Fig. 3C), Succinivibrio was enriched in the NM compared to the Chinese samples. While there were elevated levels of Acinetobacter in manure and compost compared to NM, this relationship was not statistically significant (P value, <0.05; q value, >0.5). Within Bacilli (Fig. 3D), compared to all other samples, Lactobacillus and Aerococcus were enriched in manure, while Gracilibacillus and unclassified Bacillaceae 2 were enriched in compost. Streptococcus enrichment in manure was nearly significant compared to compost (q value, 0.07). Within Clostridia, unclassified Ruminococcaceae and unclassified Lachnospiraceae were enriched in manure and NM compared with compost and soil, while Clostridium sensu stricto and unclassified Clostridiales were enriched in manure compared to soil.

ARG-MGE cluster cooccurrence with phylogenetic groups.

We sought to determine if the abundance of resistance genes correlated with abundance of specific phylogenetic groups. Lactobacillus cooccurred with the Chinese IS1216-tet cluster (average ρ = 0.78) within the Chinese samples (Fig. 2A). Comparatively, the ARG-ARG correlation coefficient cluster average (ρ = 0.90) was significantly higher than that for Lactobacillus (P value, 1.0 × 10−8). There are only three cooccurrence relationships between genera and the intI1-IS6100 cluster (average ρ = 0.76) (Fig. 1A). A housekeeping gene (gapA) allele was identified as E. coli, which cooccurs with three genes of the IS1216-tet cluster, but this allele is fully integrated with the blaTEM-sul2 cluster. Another housekeeping gene, mdh, also cooccurs with the blaTEM-sul2 cluster, but this allele does not have 100% identity or match well with any distinct taxonomic group.

We explored possible relationships by considering only a single farm or a single sample type at a time. When considered alone, the Beijing farm intI1-IS6100 and IS1216-tet Beijing clusters are cocorrelated. Further, Clostridium sensu stricto and unclassified Clostridiales are well correlated with this Beijing cluster. Within manure, Acinetobacter and unclassified Clostridiales are associated with the intI1-IS6100 cluster (Fig. 2C). Within Jiaxing, Putian, compost, or soil samples, each considered individually, there were only a combined eight cooccurrence relationships (edges) between any genera and the most abundant ARG or MGE alleles, indicating nearly no ARG-taxonomic cooccurrence.

ARG diversity.

On average, the top 40 functional gene alleles account for 91% of all sequences of each sample of the genes (see Fig. S4 in the supplemental material). For many of the genes, especially within the intI1-IS6100 cluster, a single sequence represented more than 90% of all the sequences obtained (Fig. 4). In silico PCR of the reference sequences also commonly showed dominance of a single sequence, which indicated a conserved amplicon region. The Shannon diversity of all alleles of functional genes was low, often ~1 (Fig. 4C). However, the diversity of some genes (e.g., tetG, intI2, and mefA) was much higher. Despite the dominance of the sequences by a single allele, the less abundant alleles were maintained across samples. When cooccurrence analysis was performed with the top five alleles of every gene, if all five alleles cooccur with every other allele, there would be 20 lines between these intergene alleles. Often, the top five alleles of a gene were cooccurring with each other across the Chinese samples (noted on the bars in Fig. 4A and B).

DISCUSSION

Clusters of highly abundant cooccurring resistance genes and mobile genetic elements of identical sequence were found in swine farms across China and another resistance cluster in laboratory swine at the U.S. National Animal Disease Center. Cooccurrence resistance clusters could be the result of a number of actual scenarios, including simultaneous enrichment of numerous strains with single or few resistance genes from the clusters, and/or enrichment of strains with multiple resistance genes scattered or clustered within their genome. While our results are congruent with the selection for multidrug-resistant bacteria in animal agriculture, our experimental design did not test this hypothesis directly. Several lines of evidence from this study highlight a need to understand the phylogenetic and genomic context of resistance genes and mobile genetic elements in animal agriculture environments. (i) The abundance of resistance genes and mobile genetic elements with no sequence divergence are tightly correlated across many sample types in and between swine farms. (ii) The resistance genes cooccur with integrases and insertion sequences which are known for the assembly of resistance islands (28, 38). (iii) Identical sequences in our resistance clusters are found in known resistance islands (40) and multidrug-resistant plasmids (42), and most of the intI1-IS6100 cluster genes (including aadA, aadA2, cmlA1, qacEΔ1, aadA5, dfrA1, aadA9, sul2, and dfrA2) have been observed previously as integron cassettes (43, 44). The cooccurrence of intI1, qacEΔ1, and aadA is not unexpected, as these are the canonical integron cassettes (37), together with sul1 (which was not included in our set of primers). Colocalization of intI1 with IS6100 has been previously identified in multidrug-resistant plasmids in the soil bacterium Corynebacterium (42), enterobacterial plasmids (45), Aeromonas from a Spanish river (46), Acinetobacter baumannii (40), and Pseudomonas aeruginosa clinical isolates (47). The combination of intI1 and IS6100 allows for the recruitment of resistance genes by the integrase and the duplication and transfer of large chromosomal inversions by IS6100 (47). Consequently, resistance genes tend to colocalize, and the development of complex resistance gene loci has been shown in single strains (reviewed in reference 48) and recently to a small degree in metagenomic sequences (49). Less abundant resistance alleles cooccur with these clusters as well (numbers in Fig. 4A and B), indicating rare instances of different allelic composition of resistance clusters. (iv) The resistance genes confer resistance to antibiotics both fed and not fed to the animals. With certainty, the NADC swine were never fed aminoglycosides; however, the aminoglycoside phosphotransferase strA [synonymous with aph(3″)-Ib] was found to be statistically enriched with ASP250 feed (20). The aph(3″)-Ib sequences obtained in this study were amplified with the same primer as in reference 20 and are shown to cooccur with blaTEM and in the same cluster as sul2, genes that were enriched presumably because of the use of penicillin and sulfamethazine (20). This potentially indicates their coselection by colocalization of blaTEM and aph(3″)-Ib in single genomes, as has been observed in E. coli strain 2009C-3133 (Table 1). Interestingly, E. coli gapA alleles cooccur with the blaTEM-sul2 cluster, and E. coli was shown previously to be enriched with ASP250 feed (20). Now, we show the cooccurrence of E. coli, blaTEM, sul2, and aph(3″)-Ib, implicating E. coli, as the carrier of all the resistance genes, as the possible mechanism of the coenrichment. In a similar manner, coselection of resistance genes likely occurred with the intI1-IS6100 resistance cluster. While each Chinese farm used different antibiotic cocktails for therapy and growth promotion (antibiotic usage reported in reference 19), no farm reported the use of all the antibiotics for which this cluster confers resistance (e.g., no farm used chloramphenicols or trimethoprim, Jiaxing did not use sulfonamides, and Putian did not use aminoglycosides, while resistance to these antibiotics is found across all farms). Thus, independent of antibiotics used, the same genes, which may be ubiquitously distributed in China, are maintained by selection and coselection of genetically linked resistance genes.

Persistence of resistance genes despite their discontinued use will be an obstacle in reversing the prevalence of antibiotic resistance genes now present in the environment, on animal farms, and in clinical settings, as has been observed in human populations (50). Coselection for resistance genes can occur when an antibiotic (49, 51, 52), heavy metal (19, 53, 54), or disinfectant (53) is introduced for which resistance to that agent and other resistance genes are genetically clustered. The abaR1 resistance island in Acinetobacter baumannii strain AYE is a case in point: it has 45 antibiotic resistance genes, operons encoding arsenic and mercury resistance, and the disinfectant efflux pump qacEΔ1 (40). Our data highlight two important different cases of resistance clusters. First, the IS1216-tet cluster, while very abundant in the manure (the most abundant of all resistance alleles), is well eliminated by composting and land application. This may be due to the degradation of tetracycline in these sites (55) and other IS1216-tet cluster-related selective agents. However, the intI1-IS6100 cluster remained abundant in all farms, and composting sometimes significantly increased its enrichment, potentially due to chemical persistence of sulfonamide molecules that were enriched by composting (19) or by some other nondegraded coselective agent, such as copper. As the use of many antibiotics will be phased out as growth promoters in agriculture, monitoring of all antibiotic resistance genes should continue in order to estimate the new practice’s impact on the total resistome.

Perhaps the most significant question in the field of antibiotic resistance gene ecology is the frequency of horizontal resistance gene transfer in the real world. Swine manure has been shown to contain individual strains with transferable resistance plasmids (21). Several studies have shown specific instances of horizontal gene transfer of resistance genes in the human gut (56), swine (57), and soil (30). Smillie et al. showed that the horizontal gene transfer of resistance genes is more common than that of all other genes and 25 times more likely in the gut than in other environments (58). Agricultural antibiotics have been shown to induce resistance gene transfer by prophages (35). Metagenomic sequencing projects have shown that phylogeny is a strong predictor of the antibiotic resistance profile in soil (59). In this study, Lactobacillus was an indicator organism of the IS1216-tet resistance cluster; however, there was no phylogenetic predictor of intI1-IS6100. It is possible that our 16S rRNA primers did not amplify or that the amplicon libraries were not sequenced sufficiently deeply to identify the strains that cooccur with (and, by inference, carry) these resistance genes. Another possible explanation is that a decoupling of resistance genes and phylogeny occurred in these farms via horizontal gene transfer. Resistance clusters persisted and were even enriched in multiple environments (especially in the Putian soil) with different oxygen, nutrient, and water contents where the original bacterial communities were vastly different (Fig. 3). Taken together, these results provide a strong impetus to more fully understand the genetic context of antibiotic resistance genes in high-selection-pressure environments and the distribution of resistance elements throughout the associated bacterial community.

In human-impacted environments, integrons have been shown repeatedly to be important elements in assembling arrays of antibiotic resistance cassettes, often on plasmids (37, 60, 61). Clinical class 1 integrase (intI1) sequences showed very high enrichment in the Chinese farms while the gene is hardly detected (3 to 4 orders of magnitude less abundant than Chinese manure) in the pristine soil, NADC, and feral fecal samples (see Fig. S1 in the supplemental material). The obtained intI1 amplicon sequences were dominated (89.4% of 125,847 total intI1 sequences) by a single sequence which is 100% identical to “clinical” (34) type 1 integron integrase and is the allele that overwhelmingly cooccurs with antibiotic resistance genes. When including all other operational taxonomic units (OTUs) within a single nucleotide polymorphism of the dominant sequence, these sequences account for 98.2% of the obtained sequences. These findings reinforce the potential to use intI1 as a marker for anthropogenic antibiotic use (34) and support the theory of broad dissemination of a single genotype of class 1 integrons associated with antibiotic resistance genes (37).

Antibiotic resistance genes are becoming recognized as contaminants independent of their bacterial host (28, 62). To pose a human health risk, an agricultural ARG must be in or transferred to a human pathogen, and the ARG must be transported to a location with which humans have direct physical contact (63). The probability of a resistance cluster from these farms contaminating surrounding heavily populated eastern China (~1 billion people), through routes previously observed (aerial [64], water [65], or farmer transport [66]), seems to be elevated by both the thousands-fold enrichment resulting in up to 0.1 to 5% of the bacterial community having resistance genes and the proximity to high-population centers. Furthermore, the probability for these resistance clusters to cross the commensal-pathogen barrier (due to coenrichment with MGEs) also seems to be elevated (63). Horizontal gene transfer elements like the intI1-IS6100 combination allow not only resistance gene recruitment by the integrase but genetic relocation, potentially to plasmids, by the transposases. Transposons with both intI1 and IS6100 have been identified in plasmids or chromosomes in multiple phyla (40, 42, 45–47), demonstrating their potential for future horizontal gene transfer. On the Chinese farms, one mobile genetic element, IS4, cooccurs with many genes both within and outside the intI1-IS6100 cluster. This gene might allow a future colocalization of its cooccurring genes (IS26, aphA1, dfrA2, tetM, and IS613) together with the intI1-IS6100 cluster. This genetic platform could allow for future integration and subsequent horizontal gene transfer of troublesome resistance genes, should these strains be introduced to the human microbiome or a hospital environment. This highlights the risk and importance of these microbial communities and manure disposal management.

Future sequencing projects should use advanced technologies to sequence longer DNA amplicons to increase resolution of the gene identity. Shotgun metagenomics cannot yet rival the cost, depth, specificity, and quantitative value of this highly parallel targeted-metagenomics approach. The future direction for antibiotic resistance ecology research is clear: to directly establish the genetic context surrounding resistance genes in microbial communities. This will allow us to understand the mobility of individual genes and groups of resistance genes, the assemblage of resistance islands and plasmids, and their distribution across the microbial community. This should be central to more informed stewardship programs and contribute to the U.S. National Action Plan (5) and international plans to combat antimicrobial resistance.

MATERIALS AND METHODS

Amplicon library preparation and sequencing.

All purified DNAs (46 samples in total) were stored at −20°C until use in this study. The 46 samples (50 ng/µl) and 47 primer sets were input into a 48.48 Access Array integrated fluidic circuit (Fluidigm Corporation, CA, USA) according to the manufacturer’s instructions for 4-primer amplicon tagging in 2,209 individual 30-nl qPCR mixtures. EvaGreen dye was added to the reaction mixture to allow for real-time quantification of amplification products during cycling. Threshold cycle values were obtained following cycling, and the barcoded amplification products were pooled into one mixture. The amplicon pool was prepared using Fluidigm FL1 and FL2 workflow and sequencing primers according to the manufacturer’s protocol. Sequencing was performed using 150-bp paired-end reads plus the barcode indexing read on an Illumina MiSeq sequencer, and 6.7 million raw reads were obtained, of which 6.0 million passed the MiSeq quality filter and had a barcode detected. Each sample had an average of 128,679 sequences matched to it by the barcode.

Genes targeted.

We used a novel highly parallel qPCR and amplicon library-generating platform (Fluidigm Access Array). Using this system, we obtained quantitative information on the abundance of the genes in the original sample and harvested sequencing-ready barcoded amplicons. We selected primer sets that target antibiotic resistance genes, transposases, integrases, plasmid mobility genes, and housekeeping genes (see Table S1 in the supplemental material). We used two universal 16S rRNA gene primers: one amplifies a 60-bp region, ideal for qPCR, while the other amplifies the 16S rRNA gene V4 region to determine bacterial community membership. We selected for this study the genes that were found to be enriched in the NADC or China samples from our previous studies (19, 20), as well as primer sets to target additional mobile genetic elements (32, 54, 67) and housekeeping genes (68, 69). Among the primer sets are some that target the same genes. Amplicons of the same gene obtained with different primers are differentiated with a gene suffix (for example, aadA and aadA.1).

Resistance gene quantification and analysis.

Threshold cycle values obtained during library generation were quality checked by the Fluidigm software and compared to values obtained previously with other high-throughput qPCR platforms (19) to access the reliability of the gene quantification. Following these quality checks, the data were processed as described previously (19) to calculate the average ΔΔC for each sample type, using the Chinese pristine soil as the reference sample. Heat maps were generated for the ΔΔC using RStudio (R version 3.1.2) using the heat map.2() function.

Sequence processing and analysis.

In total, 5,490,078 forward and reverse sequences (90.8%) were assembled with the RDP paired-end read assembler (70) using default parameters and with the minimum overlap set to 10 bp and a minimum read quality score of 25. Amplicons less than 195 bp from each primer had similar assembly rates ranging from 83% to 97% of raw reads. Longer amplicons were assembled at a lower rate (see Fig. S5 in the supplemental material). Assembled sequences were analyzed to detect and remove chimeras using the Uchime software (71). RDP FrameBot was used to correct frameshifts caused by indels. Sequences were quality filtered with primers removed by RDP’s Initial Process tool (70), allowing a minimum length of 20 bp, 0 N base calls, and 1 mismatch to the forward or reverse primers. The resulting 5,289,504 quality-filtered sequences in all samples and primer sets were clustered at 100% nucleotide similarity using the RDP tool McClust. Project-wide single and double sequences were removed, and the 16S rRNA gene sequences were processed using mothur according to the Schloss standard analysis protocol (72). OTUs were binned by their taxonomic classification within mothur v.1.36.0 by comparison to the RDP database using the phylotype() and classify.otu() commands. Metastats implemented within mothur v.1.34.4 was used to find statistical differences (false discovery rate q value, <0.05) in phylogenetic groups between samples.

Sequences were compared to the NCBI database both for primer specificity (see the supplemental material) and for multiple amplicon colocalization in known genomes. Amplicon colocalization analysis was accomplished by submitting a concatenated sequence of all the most abundant amplicons obtained from each primer set.

Cooccurrence analysis and network generation.

In order to compare the abundances of individual functional gene alleles to each other, we adjusted the C value of each functional gene to calculate an allele-specific adjusted C (AC): AC = C −log2(allele % abundance), where C is the C value of the primer set (the abundance of all alleles) specific to the allele in question and the difference with log2(allele % abundance) adjusts that C value to estimate the contribution of each allele individually. This calculation assumes 100% PCR efficiency and assumes that an allele with a percent abundance of 0.5 will increase the C by 1. For genes that were not detected by qPCR or sequencing, the limit of detection was used to calculate the AC, and the percent abundance was set to 1/2,500, mimicking the library sizes that were typically obtained. After generating the AC, the traditional ΔΔC approach was used, using the AC only in place of the C to calculate a ΔΔAC rather than a ΔΔC. One deviation from this protocol was ΔΔAC = ΔCT(control sample) − ΔACT(experimental sample), where the ΔC used was that of the detection limit ΔC of the pristine Chinese soil. This use of a single value for ΔCT(control sample) was used to avoid correlation artifacts between gene data sets. This resulted in a no-detection ΔΔAC value of 7.8 for 16S OTUs and −15 for non-16S rRNA genes. Abundance relative to 16S rRNA was calculated using log2(−ΔAC). The ΔΔAC for all 16S rRNA gene OTUs and the top 5 alleles for all other genes were concatenated into a single file and used as input for cooccurrence analysis as described previously (73). Networks were additionally rendered using Cytoscape v. 3.0.2. Networks were organized by the software’s “preferred layout” option, with node size dependent on number of cooccurring genes (edges).

Sequence accession numbers.

Assembled sequences separated by sample and gene have been deposited in the NCBI Sequence Read Archive (SRA) under accession numbers SAMN04523341 to SAMN04523367 in association with BioProject PRJNA313201.

Allele-specific relative abundance compared to 16S rRNA abundance. Each gene listed represents the most abundant allele of that gene from each of the different farms, Beijing (A), Jiaxing (B), and Putian (C), and the NADC and feral animals (D). Colored highlights indicate genes belonging to the intI1-IS6100 cluster (blue), the IS1216-tet cluster (green), or the blaTEM-sul2 cluster (gray). In each panel, manure, compost, and soil samples are represented with orange, blue, and green symbols, respectively. It is clear that the abundance of alleles belonging to the same gene cluster fluctuates as a group. In instances where symbols are missing, the genes were not detected. The calculated abundance of the genes, using the adjusted ΔC method described here, at the limit of detection (threshold cycle = 34 and only 1 sequence of 2,500 total sequences for the gene) present would be about 1 × 10−11. Download

Figure S1, PDF file, 0.2 MB

Heat map of the averaged qPCR results. Both axes are clustered based on Euclidean distance. Samples cluster largely by sample type. The feral samples are the most distinct of all samples. intI1 was largely not detected in the NADC or feral animals. Download

Figure S2, PDF file, 0.03 MB

Examples of genes from the intI1-IS6100 resistance cluster in known genomes. Genes with bold black outlines are identical to alleles of the intI1-IS6100 cluster. The same integron and cassettes were found in two E. coli strains, including IS6100, as the terminal cassette. A similar intI1-IS6100 arrangement was found in the Klebsiella strain. Download

Figure S3, PDF file, 0.1 MB

Sequence diversity of all gene types across all the samples. There are 40 rows next to each gene name indicating the 40 most abundant alleles, sorted by abundance, shown as percent abundance in each sample. The shading to indicate the allele abundance is based on a log scale. Sample triplicates are separated by vertical gray lines. Horizontal black lines were not added post hoc; these are the most abundant alleles and often appear as an unbroken black horizontal line because of their dominance in all samples. The reference sequence sample was not an experimental sample but in silico PCR products from the FrameBot reference set, analyzed together with the experimental sequences. Download

Figure S4, PDF file, 0.7 MB

Success rate of assembly of the 47 amplicon sets. The percentage of the raw forward and reverse reads that were successfully assembled was nearly 100% when the entire amplicon was sequenced by both the forward and reverse reads. When the overlap between the forward and reverse reads is shorter, the percentage of the reads that can be assembled drops quickly. The major limitation to assembly of forward and reverse reads is the low-quality score of the reverse read and low-quality scores at the distal ends of both reads. Low-quality scores at the distal ends of the reads can be overcome by overlapping with the read from the other end of the amplicon, for when the base call from the two reads is in agreement, the calculated quality score increases dramatically due to the agreement of the two reads, despite individual read quality scores approaching the low cutoff. One primer set (aadA) produced a 101-bp amplicon, which had an unexpectedly low assembly rate. Download

Figure S5, PDF file, 0.03 MB

Percent identity of the assembled amplicons to the reference sequences. Assembled amplicons were mapped to reference sequences as part of the FrameBot software to correct frameshift errors. From this analysis, we observed high percent identity of the obtained amplicon reads to the reference sequences. Each of the thin lines represents the percentage of one gene amplicon pool as a percent identity to the nearest reference sequence. The bold gray line is the overall average of all the amplicons. Over 98% of the assembled reads were >90% identical to a reference sequence, and the reference sequence to which the reads were mapped was always the gene that was targeted. One outlier was the strA primer set, of which about 50% of the amplicons were only 75 to 79% similar to our FrameBot reference set. When these sequences are compared with the NCBI database (BLAST nr database), the amplicons are shown to be aph(3″)-Ib, which was the intended target. This shows that our reference set was incomplete and that the amplicons were indeed specific as intended. These data show that the amplicons are specific to the genes that were intended to be amplified and serve as a validation of the primer sets. Download

Figure S6, PDF file, 0.03 MB

Performance of the qPCR primer set and the 16S rRNA V4 degenerate primer set in amplifying Bacteroidetes. When the temperature was high, the V4 primer did not amplify Bacteroidetes efficiently. Download

Figure S7, PDF file, 0.03 MB

Primer sets used in this study. Primer name indicates the targeted gene.

Table S1, PDF file, 0.1 MB

qPCR threshold cycle results for all primer sets and samples used in our analysis.

Table S2, XLSX file, 0.05 MB

Results are shown including validation of the 4-primer amplification chemistry with our primers, validation of the qPCR results by comparison with previous results on the Open Array platform (Applied Biosystems), amplicon specificity and diversity, and detection of Bacteroidetes by the V4 16S rRNA primer set. Download

Text S1, PDF file, 0.1 MB