Horizontal Gene Transfer Is the Main Driver of Antimicrobial Resistance in Broiler Chicks Infected with Salmonella enterica Serovar Heidelberg

ABSTRACT The overuse and misuse of antibiotics in clinical settings and in food production have been linked to the increased prevalence and spread of antimicrobial resistance (AR). Consequently, public health and consumer concerns have resulted in a remarkable reduction in antibiotics used for food animal production. However, there are no data on the effectiveness of antibiotic removal in reducing AR shared through horizontal gene transfer (HGT). In this study, we used neonatal broiler chicks and Salmonella enterica serovar Heidelberg, a model food pathogen, to test if chicks raised antibiotic free harbor transferable AR. We challenged chicks with an antibiotic-susceptible S. Heidelberg strain using various routes of inoculation and determined if S. Heidelberg isolates recovered carried plasmids conferring AR. We used antimicrobial susceptibility testing and whole-genome sequencing (WGS) to show that chicks grown without antibiotics harbored an antimicrobial resistant S. Heidelberg population at 14 days after challenge and chicks challenged orally acquired AR at a higher rate than chicks inoculated via the cloaca. Using 16S rRNA gene sequencing, we found that S. Heidelberg infection perturbed the microbiota of broiler chicks, and we used metagenomics and WGS to confirm that a commensal Escherichia coli population was the main reservoir of an IncI1 plasmid acquired by S. Heidelberg. The carriage of this IncI1 plasmid posed no fitness cost to S. Heidelberg but increased its fitness when exposed to acidic pH in vitro. These results suggest that HGT of plasmids carrying AR shaped the evolution of S. Heidelberg and that antibiotic use reduction alone is insufficient to limit antibiotic resistance transfer from commensal bacteria to Salmonella enterica. IMPORTANCE The reported increase in antibiotic-resistant bacteria in humans has resulted in a major shift away from antibiotic use in food animal production. This shift has been driven by the assumption that removing antibiotics will select for antibiotic susceptible bacterial taxa, which in turn will allow the currently available antibiotic arsenal to be more effective. This change in practice has highlighted new questions that need to be answered to assess the effectiveness of antibiotic removal in reducing the spread of antibiotic resistance bacteria. This research demonstrates that antibiotic-susceptible Salmonella enterica serovar Heidelberg strains can acquire multidrug resistance from commensal bacteria present in the gut of neonatal broiler chicks, even in the absence of antibiotic selection. We demonstrate that exposure to acidic pH drove the horizontal transfer of antimicrobial resistance plasmids and suggest that simply removing antibiotics from food animal production might not be sufficient to limit the spread of antimicrobial resistance.

S. Heidelberg detected in their ceca was 80% for the oral and seeder methods and 100% for the cloacal method during trial 1. For trial 2, 100% of the chicks in the oral and cloacal treatment were positive compared with 40% for the seeder group (n = 10, 10, and 15 chicks for oral, cloacal, and seeder treatments, respectively, in both trials). The ceca (n = 10 chicks) and litter of control chicks were negative for S. Heidelberg by culture; however, 1 litter sample was positive after a 24-h enrichment in buffered peptone water. Litter pH and moisture were higher for trial 1 (pH, 6.86 6 0.27; moisture, 25% 6 0.03%) than those of trial 2 (pH, 6.54 6 0.17; moisture, 21% 6 0.03%) (pH, W = 14, P = 0.11; moisture, W = 13, P = 0.18) (see Fig. S1 in the supplemental material).
S. Heidelberg perturbed the ceca and litter microbiome in a challenge routespecific manner. Next, we determined if S. Heidelberg challenge perturbed the microbiome of broiler chicks. We profiled the microbial community present in the ceca and litter using 16S rRNA gene amplicon sequencing. The bacterial alpha diversity was higher in the litter than in the ceca for the oral treatment for all alpha diversity measures (observed, Chao1, Shannon, and InvSimpson) examined (P , 0.01) (see Fig. S2a in the supplemental material). In addition, observed and Chao1 alpha diversity indices were higher for litter than ceca for cloacal and control treatments (P = 0.007). Contrastingly, the ceca from orally treated chicks had the lowest alpha diversity for all measures used compared with the control (P , 0.05).
Beta diversity differed between the litter and ceca, and the route of S. Heidelberg challenge affected the beta diversity of the litter (F = 43.7, df = 11, P = 0.001) more than the ceca (F = 2.21, df = 29, P = 0.007) (Fig. 1d). Furthermore, cecal samples were more dispersed than litter samples. Firmicutes sp. dominated the cecal microbiome of broiler chicks (Fig. 1e) and Ruminococcaceae and Lachnospiraceae (class Clostridia) were the major families. Amplicon sequence variants (ASVs) matching Subdolinagrum, Lactobacillus, and Family_Lachnospiraceae were the most abundant ASVs in oral and cloacal chicks compared with controls in the ceca (Fig. S2b).
Firmicutes, Proteobacteria, Actinobacteria, and Bacteriodetes were the dominant phyla present in litter (Fig. 1e), and Escherichia/Shigella and Klebsiella were the most abundant ASVs in litter (Fig. S2b). Escherichia/Shigella was higher in cloacal and control chicks than in oral and seeder chicks, and Klebsiella was higher in control than in oral, cloacal, and seeder chicks (Escherichia/Shigella, x 2 = 9.9744, P = 0.018; Klebsiella, x 2 = 9.6667, P = 0.021). This result suggested that S. Heidelberg perturbed the bacterial community present in the litter and ceca of broiler chicks, especially members of the family Enterobacteriaceae.
Determining if S. Heidelberg acquired AR plasmid after broiler chicks were challenged. To determine if S. Heidelberg acquired AR after 14 days of challenge, we used the National Antimicrobial Resistance Monitoring System Sensititre panel for Gram-negative bacteria to conduct antimicrobial susceptibility testing (AST). We performed AST on 250 S. Heidelberg isolates recovered from the ceca and litter of broiler chicks. SH-2813 nal R carries a gyrA mutation for nal resistance; therefore, all isolates recovered were resistant to nal. S. Heidelberg also carries a chromosome-carried fosA7 gene that confers resistance to fosfomycin (18).
Next, we found the genetic determinants responsible for the acquired AR phenotype using whole-genome sequencing (WGS). We performed WGS on 69 S. Heidelberg isolates from trial 2 and 2 isolates from trial 1. We focused our resources on trial 2, for which AR acquisition was higher than trial 1. The isolates sequenced were either susceptible (n = 33) or carried resistance to at least one antibiotic in two antimicrobial drug classes (n = 38). The ancestral SH-2813 nal R (isolate used for challenge) harbors an IncX1 plasmid, and carriage of this plasmid does not confer AR in S. Heidelberg. The IncX1 incompatibility (inc) region was not detected by PlasmidFinder (19) in 42 of the 71 isolates sequenced.
Thirty-nine percent of S. Heidelberg isolates acquired Col plasmids, but none carried a known antibiotic resistance gene (ARG). ARGs known to confer resistance to aminoglycosides [aadA1 and aac(3)-Via] and tetracyclines (tetA) were detected in all AR isolates from trial 2. The AR isolates from trial 1 harbored a sulfonamide resistance gene (sul1) in addition to aminoglycoside genes but no tet gene. The ARGs were found on a plasmid belonging to inc group IncI1 (Fig. 2c) and plasmid multilocus sequence type 26 (pST26). A change in AR profile was seen in four isolates from the seeder group; for these isolates, AST confirmed that they acquired AR but after WGS the AR determinants were not detected. Upon AST retest, these isolates were found to be susceptible. Acquired ARG was absent in four isolates with resistance to either tetracycline or streptomycin.
pST26 differed in its AR genetic context and genome architecture. We questioned if the pST26s present in the sequenced S. Heidelberg isolates were similar or different in gene content or structure. To answer this question, we used hybrid assemblies (Illumina short reads and MinION or PacBio long reads) to achieve complete circular pST26 plasmids for two AR S. Heidelberg isolates recovered from trial 1 and trial 2. The pST26 plasmids from trial 1 (p1ST26) and trial 2 (p2ST26) were ;112 kbp long and ;87% identical ( Fig. 2c). Both carried the atypical IncI1 backbone, including regions encoding replication, stability, leading, and conjugative transfer (20). To determine the IncI1-complex group, we used the classical traY and excA protein sequences of plasmids R64 (IncI1) and R621a (IncI1-gamma). A reconstructed tree confirmed that both plasmids belong to the IncI1 group (Fig. 2d). A 19.8-kb variable region encoding transposons and AR was found in both plasmids. For p1ST26, this region encoded aadA1, aac(3)-Via, sul1 and quaternary (quats) ammonium compound resistance (qacED1) (Fig. 2c, Fig. S3b). In contrast, this region carried tetA and mer operons in addition to aadA1 and aac(3)-Via in p2ST26 (Fig. 2c).
To investigate whether gene flux influenced IncI1 genome architecture, we used p1ST26 and p2ST26 as reference genomes for the rest of the AR isolates carrying IncI1. By aligning raw reads to complete pST26 genomes, we were able generate a consensus IncI1 plasmid contig for each isolate. A tree built with these plasmids resulted in two clades represented by p1ST26 and p2ST26 plasmids, respectively (Fig. 3a). Multiple alignment of protein sequences revealed that pST26 plasmids differed in the number of genes carried and gene alleles. This finding was pronounced for regions encoding AR, pilus/shufflon assembly proteins (pil), and IS66 family of transposases (see Data Set S1 in the supplemental material; Data Set S1 to S5are available in the Zenodo repository online at https://doi.org/10.5281/zenodo.4976002 and in the supplemental material). A pangenome analysis revealed that the pST26 from this study carried 70 core genes (genes present in $95% of the plasmids) and 90 accessory genes (genes present in ,95% of the plasmids).
A tree reconstructed with the core genes and accessory genes divided the pST26 into 5 major groups ( Fig. 3b and c). For both trees, the p1ST26 formed one clade and the p2ST26 made up the other clades. No single nucleotide polymorphisms (SNPs) were found between p2ST26s, but p1ST26 differed from p2ST26 with 21 SNPs and 2 insertions (Data Set S1, see Table S1 in the supplemental material). The insertions were present in only one p1ST26 plasmid. One insertion was found on tniA (a DDE-type integrase/transposase/recombinase protein) and the other was between tn3-like transposon and yadA. Inserted DNA showed significant homology to tniA present in E. coli plasmid EcPF5 (GenBank accession no. CP054237) and IS5057 present in E. coli plasmid pIOMTU792 (GenBank accession no. LC542972.1), respectively.
Determining if the acquisition of new genes affected the genome of S. Heidelberg. We hypothesized that successful S. Heidelberg colonization will be dependent on gene flux (gain and loss of genes) and mutations. To test this hypothesis, we reconstructed a maximum likelihood (ML) tree based on the pangenome and mutations of S. Heidelberg strains recovered. The core genome (genes present in $95% of the strains) and accessory genome (genes present in ,95% of the strains) was composed of 3,729 and 1,531 genes, respectively. A total of 91 new mutations (SNPs and indels) were on the chromosome of S. Heidelberg strains, and no single mutation was present in every strain. Rather, mutations were unique to individual strains or shared between 2 and 26 isolates (Data Set S2). Consequently, the core genome and SNP-based trees did not provide a clear division of strains based on plasmid acquired, and clades had low bootstrap support values (see Fig. S4a and b in the supplemental material).
Contrastingly, the accessory genome tree grouped the strains based on AR phenotype, plasmid, and ARG carried (Fig. 4a). Acquisition of IncI (IncI1 and IncI2) and Col plasmids and the presence/absence of IncX1 plasmids defined the clades seen on the accessory genome tree. Susceptible strains carrying only IncX1 plasmids represented the ancestral clade and AR strains dominated nested clade II. Clade II strains carried only either IncI1 or IncI1 plus a Col plasmid or both were present with an IncX1 plasmid (Fig. 4a). A few susceptible and strains with no plasmid (n = 12) were nested within subclade IIa, and the genome of these strains showed signs of substantial gene loss (Fig. 4a, see Fig. S5a in the supplemental material). One S. Heidelberg strain (og8-05a) harbored an untypeable circular episome encoding uncharacterized proteins and two genes encoding NADH-quinone oxidoreductase (nuo) (;22 kbp) (Fig. S5b) and an IncX1 plasmid with divergent protein sequences from the IncX1 present in the ancestor (Fig. S5c). FIG 2 Legend (Continued) antibiotic resistance grouped by the route the broiler chicks and were challenged (n = 54, 62, and 42 for oral, cloaca, and seeder treatment, respectively). (c) BLASTn alignment of the IncI1 plasmid acquired during trial 1 (p1ST26) and trial 2 (p2ST26); dashed black lines highlight mobile region carrying antimicrobial resistance genes, while blue dashed lines show other regions that are different between p1ST26 and p2ST26. (d) Maximum likelihood tree constructed using TraY and ExcA protein sequences from a representative IncI1 plasmid from this study, R64 (IncI1), and R621a (IncI1-gamma). The GTR model of nucleotide substitution and the GAMMA model of rate heterogeneity were used for sequence evolution prediction. Tree was rooted with the traY sequence of R621a. Subclade IIa strains also had a higher number of contigs, higher misassemblies, and lower genome size than the rest of the strains sequenced for this study (test statistic for Kruskal-Wallis test [H] = 18.82, df =1, P , 0.001) (Fig. S5d to f).
To determine if erroneous assembly affected the clustering of this clade, we performed long read sequencing (PacBio) on one strain from this clade and used it in a FIG 3 Gene flux contributed to the diversity of IncI1 plasmids. (a to c) Maximum likelihood tree of pST26 IncI1 plasmids from this study (n = 36) constructed using complete plasmid DNA sequences (a), core genes (b), and accessory genes (c). The GTR1I, GTR, and JC1I models of nucleotide substitution and the GAMMA model of rate heterogeneity were used for sequence evolution prediction for a, b, and c, respectively (tree was rooted with the closest relative; GenBank CP016585), found through NCBI BLASTn search. Clade numbering was assigned arbitrarily to show the number of clades found. Numbers shown next to the branches represent the percentage of replicate trees where associated isolates cluster together based on ;100 bootstrap replicates. Heidelberg isolates (n = 72) recovered from the ceca and litter of chicks colonized with S. Heidelberg. The GTR model of nucleotide substitution and the GAMMA model of rate heterogeneity were used for sequence evolution prediction. Numbers shown next to the branches represent the percentage of replicate trees where associated isolates cluster together based on ;100 bootstrap replicates. All S. Heidelberg strains were assembled using Illumina short reads except OG8-05A, IC9B, and IC4-3A and SH2813-ancestral that were assembled using both Illumina short reads and PacBio or MinION long reads. The tree was rooted with the ancestral susceptible S. Heidelberg strain. Clade numbers were assigned arbitrarily to ease discussion on isolates that acquired antibiotic resistance. A blue rectangular box highlights the subclade with susceptible strains nested within antibiotic resistance strains due to misassembly bias. (Nal, nalidixic acid; Gen, gentamicin; Str, streptomycin; Tet, tetracycline; Fis, sulfisoxazole; NA, not applicable). (b) Mauve visualization of the inverted genome of S. Heidelberg strain og8-05a. The chromosomal contigs of og8-05a were aligned and ordered to the complete chromosome of the S. Heidelberg ancestor. A colored similarity plot is shown for each genome, of which the height is proportional to the level of sequence identity in that region. When the similarity plot points downward, it shows an alignment to the reverse strand of the S. Heidelberg ancestor genome, i.e., inversion. A segment highlighted with solid black rectangular box denotes the inverted region in og8-05a, while dashed blue and red rectangular boxes denote segments with mutations. hybrid approach with Illumina short reads. This procedure resulted in a partially complete chromosome (longest contig, 4.7 Mbp), but 30 noncircular contigs were not scaffolded. An alignment of the ordered contigs with the ancestor showed that the strain did not suffer significant gene loss as suggested through short or long read only assembly (see Table S2 in the supplemental material); however, a 3.2-Mbp chromosomal region was inverted (Fig. 4b). To confirm if chromosomal inversion influenced the assembly of these genomes, we aligned protein sequences of two complete circular S. Heidelberg chromosomes from trial 1 and 2 with the ancestor. This analysis showed that the isolate from trial 1 harbored an ;2-Mbp chromosomal inversion (see Fig. S6a in the supplemental material).
Gene inversion changes the leading or lagging strand sequence to its reverse complement, thus altering the GC skew of the affected gene or genome from a positive value to a negative value or negative to a positive (21). Accordingly, we confirmed that there was a GC skew inversion in these regions ( Fig. S6b and c). Gene inversions occur after rejoining DNA breaks and inverted genes can introduce sequencing bias and misassembly (22). Here, deletions affecting the cytochrome c maturation gene cluster (ccm) contributed to the inversion of genes in S. Heidelberg strain og8-05a ( Fig. 4b, see Fig. S6b in the supplemental material). Furthermore, the oriC of og8-05a was reoriented compared with the ancestor. On the other hand, deletion of a putative YebC/ PmpR transcriptional regulator and a deletion between cysK and cysZ led to the reorientation of oriC in strain ic9b ( Fig. S6a and c).
The inverted region in these genomes encoded 34% to 60% of the 168 virulence genes present in the ancestor, including Salmonella pathogenicity islands 1 and 2, fimbrial and adherence genes, and type 1 and 2 secretion system proteins and prophages (Data Set S3; Fig. S6b and c). Gene inversions can increase the diversity and virulence of pathogens (21), and our bioinformatic analysis confirmed that these S. Heidelberg strains inverted their genomes. Nevertheless, in vivo testing will be required to confirm their virulence potential. Our results suggested that gene flux and homologous recombination shaped the genome of S. Heidelberg strains after they gained entry into the gut of broiler chicks.
Determining the commensal reservoir of IncI1 plasmids in broiler chicks. The application of the proximity-ligation method (Hi-C) has improved the assembly of metagenomes and made it possible to detect plasmid-host associations (23). Therefore, we used Hi-C metagenomics to find the bacterial species carrying IncI1 plasmids in the microbiome of broiler chicks. First, we checked if Hi-C correctly assembled S. Heidelberg genome present in the ceca of challenged chicks. We selected two cecal samples from cloacal chicks from trial 2 (here referred to as Hi-IC-FL1 and Hi-IC-FL2). We chose cloacal chicks because all chicks in this group were colonized with S. Heidelberg compared to oral or seeder. The two chicks carried S. Heidelberg at ;100 CFU/g of ceca.
A total of ;125 million shotgun reads and ;206 million Hi-C reads were obtained per sample. Taxonomic assignment of Hi-C metagenome-assembled genomes (MAGs) using Mash (24,25) assigned 226 MAGs to a taxon for Hi-IC-FL1, while 87 MAGs were assigned a taxon for Hi-IC-FL2. One MAG from Hi-IC-FL1 and one from Hi-IC-FL2 were classified as Salmonella enterica serovar Cubana and S. Heidelberg, respectively. However, when we used a Salmonella serovar prediction tool (26) to validate the result from Mash, the MAGs (;1.4-and 1.2-Mbp total contig size for Hi-IC-FL1 and Hi-IC-FL2, respectively) were not Salmonella sp. Furthermore, aligning these MAGs to a complete S. Heidelberg genome yielded poor alignment results (data not shown). CheckM also revealed that the Salmonella MAGs were incomplete for both Hi-IC-FL1 (42% completeness) and Hi-IC-FL2 (12.4% completeness). Similarly, taxonomic assignments performed with the Genome Taxonomy Database Toolkit (GTDB-Tk) (27) did not assign any MAG as Salmonella sp. Contrastingly, 23,610 6 7,156 raw shotgun reads from each sample were classified as Salmonella enterica using the Kraken 2 database (28) (Fig. 5a). This (Continued on next page) Oladeinde et al. result suggested that our metagenomic approach did not have the discriminatory power to identify the S. Heidelberg strain in the cecal samples.
Next, we investigated if the MAGs harbored IncI1 plasmids. The IncI1 plasmid contigs detected in both samples were linked to E. coli MAGs ( Fig. 5b and c; see Table S3 in the supplemental material) and multiple MAGs were classified as E. coli by Mash for both Hi-IC-FL1 (n = 3) and Hi-IC-FL2 (n = 12), and one MAG was ;95% complete. Some bacterial MAGs assigned to the genus Firmicutes showed a low level of linkage to IncI1 (;100Â less compared than E. coli in Hi-IC-FL1). In addition, we found contigs (n = 151; mean contig size, 17,607 bp; smallest contig, 1,027 bp; largest contig, 282,389 bp) with protein sequences identical to pST26 in both Hi-IC-FL1 and Hi-IC-FL2 (Data Set S4). This result suggested that E. coli was the main reservoir of IncI1 plasmids.
Broiler chicks carried E. coli strains harboring pST26 IncI1 plasmids. To confirm if broiler chicks from this study harbored E. coli populations that could serve as a reservoir of pST26 plasmids, we retrospectively screened the cecal contents of the broiler chicks used for Hi-C on CHROMagar supplemented with gentamicin and tetracycline and selected five colonies for whole-genome sequencing. As expected, we found pST26 in two strains (phylogroup F, multilocus sequence type [MLST] 6858), while the other three strains (phylogroup D, MLST 69) carried an IncI1 that harbored no ARG and an unknown MLST (Table 1). Next, we compared the pST26 in S. Heidelberg to that found in E. coli using complete plasmids. The p1ST26 and p2ST6 IncI1 plasmids of S. Heidelberg were 86.8% and 99.4% identical to pST26 present in E. coli. The main difference was in the region encoding AR and transposases in p1ST26 (see Fig. S7 in the supplemental material; Fig. S7 to S10 are available in Zenodo repository online at https://doi.org/10.5281/ zenodo.4976002 and in the supplemental material).
pST26 carriage poses a variable fitness cost and benefit under selection pressure. pST26 carried accessory genes for AR, metal resistance, or disinfectants (i.e., quaternary ammonium compounds). Furthermore, S. Heidelberg strains carrying pST26 differed by route of exposure, e.g., higher proportion in oral versus cloacal exposure. To this end, we questioned if any of these factors could exert a selective pressure for pST26. We did not administer antibiotics to broiler chicks throughout study; therefore, we did not consider antibiotic selective pressure. We detected metals, including nickel, chromium, copper, and zinc, in the starter feed in parts per million; and we found contaminants, FIG 5 Legend (Continued) of the metagenomic reads) were used and zoomed in to show 0% to 1% of total reads within phylum Proteobacteria (96.7% of reads within phylum Proteobacteria were assigned to Komagataeibacter). (b and c) IncI1 contig hosts found by Hi-C contacts. IncI1-derived contigs (horizontal axis of heatmap) show specific Hi-C associations with metagenome-assembled genomes (MAGs) present in samples Hi-IC4 (a) and Hi-IC6 (b) (vertical axis of heatmap). MAGs are derived from Hi-C deconvolution of the metagenome assembly and placed into a bacterial phylogeny using Mash and CheckM. Cluster.10 in each sample is an Escherichia coli genome. In Hi-IC4, one MAG (cluster.82) representing an extremely fragmented and partially contaminated Escherichia coli genome is omitted for clarity (due to its contamination, this cluster was placed in Clostridia by CheckM). Heatmap values indicate transformed counts of Hi-C read contacts (indicating intracellular physical proximity of IncI1 contigs to those genomes). Heatmap values were pseudocounted to facilitate plotting of log-transformed data, including zeroes. such as arsenic, lead, and cadmium, in trace amounts (see Table S4 in the supplemental material). Disinfectants, including quaternary ammonium compounds, were used for cleaning broiler houses and equipment; thus, we considered disinfectant as another source of selective pressure in this study.
To determine their effect on S. Heidelberg fitness, we used phenotype MicroArray (PM) 96-well plates supplemented with selected metal chlorides and disinfectants at various concentrations to compare a p2ST26-carrying strain (here referred to as gen R ) to a susceptible evolved strain. There was no significant difference in metabolic activity for S. Heidelberg strain carrying gen R compared with its susceptible counterpart for metals and disinfectant, but gen R exhibited a higher metabolism for benzethonium chloride, chromate, tellurite, and zinc, whereas the susceptible strain showed higher metabolism for dequalinum chloride, chromium, and cesium (x 2 = 1.83, df = 1, P = 0.1762) ( Fig. S8a  and b).
Cox and colleagues have suggested that the acidic pH of the upper gastrointestinal (GI) tract was associated with the lower number of chicks positive for Salmonella and Campylobacter sp. after oral gavage than that after cloacal inoculation (29)(30)(31). Based on this premise, we hypothesized that S. Heidelberg strains carrying gen R will survive better under exposure to acidic pH. To address this hypothesis, we first compared the survival of the strains under different pH levels (3.5 to 10) with or without nitrogen sources and amino acids using PM plates. An S. Heidelberg strain carrying gen R showed lower metabolic activity than the susceptible strain at pH 3.5 to 4 (V-statistic for Wilcoxon signed-rank test [V] = 3; P = 0.5), at pH 4.5 with and without nitrogen sources (V = 599; P = 3.426e-06), and at pH 9.5 to 10 with and without sources of nitrogen (V = 714; P = 6.547e-07) (Fig. 6a).
Next, we compared the survival of the susceptible ancestor, susceptible evolved strains, and gen R evolved strains (n = 3 strains for each population) in pine shaving extract (PSE) adjusted to a pH of 2.5. We acclimatized the bacterial population to PSE (pH 6.5) for 2 h before exposing it to acidified PSE (Fig. S8c). We screened for gen R -carrying populations using the gentamicin resistance marker on p2ST26. We found that p2ST26 carriage produced higher variable fitness effects at all time points than the susceptible strains (Fig. 6b). For instance, two gen R strains showed a decrease in fitness (values of ,1 show a reduction in fitness [w], while a value of .1 indicates an increase in fitness) after 30 minutes in PSE, while one strain showed an increase (w = 2.7). After 2 h of exposure, there was no significant difference in fitness between susceptible and gen R populations, but cells carrying gen R exhibited higher fitness than their susceptible counterparts (x 2 = 4.35, df = 2, P = 0.113) (Fig. 6b). Similarly, there was no significant difference in the final population of susceptible and evolved strains; however, the gen R strain had a lower population size (x 2 = 6.37, df = 3, P = 0.095) (Fig. 6c). This result suggested that p2ST26 carriage poses a variable fitness cost and benefit on the S. Heidelberg host when exposed to acidic pH.

DISCUSSION
In this study, neonatal Cobb 500 broiler chicks challenged with a susceptible S. Heidelberg strain carried susceptible and antimicrobial resistant S. Heidelberg populations 2 weeks after inoculation in their ceca and litter. This resistance phenotype was conferred by the acquisition of IncI1-pST26 plasmids. These plasmids are present in S. Heidelberg, S. Typhimurium, Salmonella enterica serovar Anatum, Salmonella enterica serovar Derby, Salmonella enterica serovar Schwarzengrund, Salmonella enterica serovar Saintpaul, and E. coli strains isolated from animal and human sources (Fig. S9). The closest IncI1 plasmids to pST26 from this study are pST26 present in E. coli (CP018625) and S. Typhimurium (CP027409). The plasmids carry aminoglycoside resistance genes with or without tetracycline, sulfonamide, and qac and mer resistance genes, making it an antibiotic/disinfectant/ metal resistance plasmid.
The acquisition of pST26 was higher in chicks challenged orally than in chicks inoculated via the cloaca in trial 2. This route-specific rate of conjugation suggested that the selection for pST26 in our study was associated with the differential pH in the upper GI versus the lower GI tract. This finding was corroborated by in vitro experiments, where acidic pH reduced the metabolism of S. Heidelberg cells carrying p2ST26 compared with S. Heidelberg cells with no p2ST26. Furthermore, carriage of p2ST26 imposed a variable fitness effect on S. Heidelberg exposed to acidic pH. This type of acid-imposed selection for p2ST26 alludes to the "negative frequency-dependent selection hypothesis"-where the fitness of the population is higher when a plasmid is present at a low frequency in the population (32) and also supports the tenet that variability of fitness effects contributes to plasmid maintenance and persistence in bacterial communities (33).
The length of the GI tract and the route traveled could make it challenging for S. Heidelberg to colonize the ceca of chicks or acquire AR via HGT. When we measured the distance of the cloaca and esophagus to the ceca of a 49-day-old Cobb 500 broiler chicken (photo not shown), we estimated that S. Heidelberg will have to travel ;20Â farther to get to the ceca of chicks gavaged orally than chicks inoculated cloacally. It is plausible that this longer "residence" time in the upper GI tract allows more opportunities for S. Heidelberg to make contact and compete with other members of the chicken microbiome. In this study, such an event could also explain the higher transfer rate of p2ST26 to S. Heidelberg in chicks that were gavaged.
It is noteworthy that AR acquisition occurred at a lower rate during trial 1 than during trial 2. On one hand, the lower pH of the litter in trial 2 than that in trial 1 supports our acid-imposed selection hypothesis; alternatively, several other factors can be responsible for the observed difference, including the season we conducted the experiments. Trial 1 was done in the fall, and trial 2 was performed in the spring. The temperature outside for trial 1 was higher (average, 70°F; high, 90°F; low, 53°F) than trial 2 (average, 57°F; high, 81°F; low, 34°F) (www.timeanddate.com/weather/usa/ athens/). The temperature and relative humidity inside a broiler house are correlated with the weather outside, and these environmental parameters are controlled by manipulating the ventilation system installed in the house. Broiler house ventilation rates have a direct effect on the moisture, pH, ammonia levels, and microbiome of the litter (34). The higher levels of S. Heidelberg in the litter from trial 1 than that from trial 2 suggest that high house temperature, litter moisture, and litter pH selected for an S. Heidelberg population that had no plasmid-borne AR. Lastly, we did not administer antibiotics to the 1-day-old chicks used for both trials and do not have the information to speculate on the antimicrobial practices of the hatchery or parent/breeder farms of the chicks. One or all of these factors could explain the differences between trials and indicate that further studies are needed before the results can be generalized.
The IncI1 complex is widespread in Enterobacteriaceae and has been reported to be the major driver of AR in Salmonella sp. and E. coli. Yet, our understanding of its mode of transfer and AR evolution has been limited to studies on its prevalence, mating experiments, and in vivo challenge with mice (20,35,36). In this study, IncI1 plasmid contigs found using metagenomics were linked to E. coli MAGs, and this result was corroborated by WGS of E. coli isolates. The E. coli isolates sequenced carried IncI1 plasmids that were similar to the pST26 acquired by S. Heidelberg. Fischer et al. (37) showed that E. coli populations in broiler chickens acquired b-lactam resistance after 4-days-old broiler chicks were gavaged with an E. coli strain carrying bla CTX-M-1 on an IncI1 plasmid. In a twin pair experiment, Hagbo et al. (38) revealed that IncI1 carrying bla CMY-2 , ColE1, and P1 bacteriophage are often transferred between microbiota present in the preterm infant gut. Salmonella resistance to b-lactams in retail meat was shown to correlate with AR found in E. coli, and horizontal transfer was the hypothesis put forward (39). This study and others have demonstrated that IncI1 is the most prevalent plasmid carrying b-lactam in Enterobacteriaceae and that they are transferable from Salmonella to E. coli and vice versa.
The nutrient-niche hypothesis suggests that S. Heidelberg will be able to invade a new niche only if it can metabolize a growth-limiting resource or if it can outcompete a resident species and a vacant niche arises due to a new component in diet or by eliminating the competitor (40). Metagenomics revealed that the two most abundant bacterial phyla in the ceca and litter from this study were Firmicutes and Proteobacteria. Therefore, bacterial species belonging to these phyla are the resident bacteria with which S. Heidelberg will have to compete. The most abundant genera of Firmicutes found by aligning shotgun reads to the Kraken database were Staphylococcus, Lachnoclostridium, Streptococcus, Enterococcus, Faecalibacterium, and Lactobacillus. Firmicutes have been shown to be the dominant phylum in chicken ceca in earlier studies, with relative abundances usually well above .50% as found by 16S rRNA gene sequencing (41,42). Comparable results have also been demonstrated with shotgun metagenomics (43)(44)(45). Interestingly, S. Heidelberg inoculation did not have an impact on the abundance of Firmicutes in this study. Contrastingly, we saw an effect of S. Heidelberg inoculation on Proteobacteria and was significantly associated with the genera Escherichia/Shigella and Klebsiella. The most abundant species of these genera found in ceca using shotgun metagenomics were E. coli, Shigella dysenteriae and Klebsiella pneumoniae.
Competition for oxygen was shown to be the limiting factor that contributed to the virulence of S. Enteritidis toward E. coli in neonatal chickens (46). Similarly, competition between S. Heidelberg strains and members of the family Enterobacteriaceae played a key role in selecting for antibiotic-resistant S. Heidelberg in our study. This conclusion was drawn from results showing that S. Heidelberg challenge perturbed the bacterial community of broiler chicks and S. Heidelberg acquired plasmids found only in E. coli populations. For S. Heidelberg and E. coli to coevolve and persist in the same new niche, a recombination of genes occurred and led to the emergence of new S.
Heidelberg strains. In our study, S. Heidelberg strains carrying pST26 appeared after conjugating with E. coli populations. As expected, S. Heidelberg genomes were grouped into nine core genome MLST clusters, with eight clusters representing new strains (Fig. S10). The S. Heidelberg strains that acquired pST26 exhibited variable fitness under acid exposure experiments and are expected to be persist longer in the new niche. Therefore, we conclude that simply removing antibiotics from food animal production might not be sufficient to reverse the spread of antimicrobial resistance.

MATERIALS AND METHODS
S. Heidelberg inoculum preparation. We have previously described the characteristics of the S. Heidelberg strain (SH-2813nal R ) used for this study and how the inoculum was prepared (17). Briefly, the strain belongs to the multilocus sequence type 15; carries a 37-kb conjugative plasmid; and is resistant to erythromycin, tylosin, and fosfomycin. The strain was recovered from a broiler chicken carcass in 2013 and made resistant to 200 ppm of nalidixic acid (nal R ) for selective enumeration. The nal R phenotype is conferred by a serine-to-tyrosine substitution at position 83 of DNA gyrase subunit A protein (GyrA). The S. Heidelberg inoculum was grown overnight in poultry litter extract, centrifuged, and resuspended in 1Â phosphate-buffered saline (PBS). The resuspended cells were used as inocula. The genome of the ancestor to SH-2813nal R was sequenced using Illumina and PacBio sequencing to achieve a complete circular chromosome and plasmid (GenBank accession number CP066851). SH-2813nal R carries six mutations (Data Set S2), including the mutation affecting gyrA, that are not present in the ancestor.
Challenging broiler chicks with S. Heidelberg. One-day-old Cobb 500 broiler chicks were purchased from a commercial hatchery in Cleveland, Georgia. Upon purchase, chicks were placed in plastic crates lined with brown paper and transported to the University of Georgia, Poultry Research Center (33.90693362620281, 283.37871698746522). One hundred chicks were either uninoculated (n = 25), gavaged orally (n = 25), or inoculated cloacally (n = 25) with a 100ml-volume of S. Heidelberg inoculum (Fig. 1a). We included a seeder colonization method, whereby five gavaged chicks were mingled with 20 uninoculated chicks. Each inoculated chick received ;10 6 colony forming units (CFU) of SH-2813nal R . Afterward, chicks were placed in floor pens at a stocking density of 0.65 m 2 /chick on fresh pine shaving litter. Broiler chicks were given water and feed ad libitum and raised antibiotic-free on starter diet for 2 weeks (starter feed was synthesized by the University of Georgia poultry research center feed mill). The concentration of metals and priority pollutants in feed was determined by the University of Georgia's feed and environmental water laboratory (Athens, GA, USA) (see Table S4 in the supplemental material). Husbandry and management followed commercial broiler chicken industry guidelines. After 2 weeks, 10 chicks from oral, cloacal, and control groups and 15 chickens from seeder (5 seeder chicks and 10 uninoculated pen mate) were sacrificed to determine the extent of S. Heidelberg colonization in the ceca. The experiments were done in two trials conducted in September 2017 (trial 1) and April 2018 (trial 2). The study was approved by the University of Georgia Office of Animal Care and Use under Animal Use Protocol A2017 04-028-A2.
Cecal and litter bacteriological analyses. Ceca were aseptically removed from the eviscera of 2week-old broiler chicks, placed in a stomacher bag, and transported on ice to the U.S. National Poultry Research Center for analysis. Ceca were weighed, and buffered peptone water (BPW) (BD Difco, MD, USA) was added 3Â volume to the weight (v/w) and stomached for 60 s. Serial dilutions were made and plated onto Brilliant green sulfa agar (BGS) (BD Difco, MD, USA) containing 200 ppm nal. In addition, a 10-ml inoculating loop was used to streak cecal slurry (cecal contents in BPW) onto xylose lysine tergitol-4 agar (BD Difco, Sparks, MD) supplemented with 4 ppm tetracycline for trial 1 and BGS supplemented with 32 ppm ampicillin (amp) or 32 ppm streptomycin for trial 2. All antibiotics were purchased from Sigma-Aldrich (St. Louis, MO), unless otherwise noted. Plates were incubated along with the cecal slurry for 24 h. All bacterial incubations were carried out at 37°C, unless otherwise noted. After incubation, colonies were manually counted, and serial dilution plates with 2 to 100 colonies were used for CFU per gram calculation. If no colonies appeared on serial dilution plates, a cecal slurry was streaked onto a new BGS plus nal plate and incubated overnight. After incubation, plates were examined for the presence/absence of Salmonella colonies.
Broiler chicken litter was collected as grab samples from 7 locations (4 corners of the pen and under the waterer) in each pen after chicks were removed. The litter samples were pooled, and 30 g was processed in duplicates from each pen as previously described (47). Serial dilutions of the litter slurry were made and plated onto BGS agar containing 200 ppm nal. Plates were incubated overnight, and colonies were manually counted and reported per gram litter dry weight. Litter pH and moisture were determined as described previously (47). We selected randomly 2 to 6 single colonies from BGS plates supplemented with nal from each ceca and litter sample and archived them in 30% Luria Bertani broth (LB) glycerol at 280°C. In addition, cecal slurry was saved at a 4:1 ratio in Luria Bertani broth (BD Difco, MD, USA) containing 30% glycerol at 280°C, whereas litter samples were stored in vacuum-sealed whirl pak bags at 220°C.
Antimicrobial susceptibility testing. We performed antimicrobial susceptibility testing (AST) on 2 50 S. Heidelberg isolates recovered from the ceca and litter of broiler chicks following the National Antimicrobial Resistance Monitoring System (NARMS) protocol for Gram-negative bacteria. MICs for isolates were determined by broth microdilution using the Sensititre semiautomated antimicrobial susceptibility system (Thermo Fisher Scientific, Inc., MA, USA) and interpreted according to clinical and laboratory standards institute guidelines when available; otherwise, breakpoints established by NARMS were used. AST was also done on E. coli isolates (n = 5) recovered from the ceca of two broiler chicks challenged intracloacally with S. Heidelberg. The Kirby-Bauer disk diffusion assay for tobramycin, kanamycin, neomycin, and netilmicin was done on four gentamicin-resistant S. Heidelberg isolates as previously described (47).
DNA extraction. DNA was extracted and purified from bacterial colonies using a FastDNA spin kit (MP Biomedicals, LLC, CA, USA), whereas 250 mg of cecal and litter was extracted with the Qiagen DNeasy power soil DNA kit (Hilden, Germany). The modifications we made to the manufacturer's protocol for DNA extraction have been reported (47,48).
16S rRNA gene sequence processing and analysis. Cecal and litter DNA were used for bacterial community analysis through the sequencing of the V4 hypervariable region of the 16S rRNA gene of bacterial genomes. Sequencing was done using the paired-end (2 Â 250 bp) method on the Illumina MiSeq platform. Statistical analysis of microbial communities was performed in the R environment using the packages "phyloseq," "Ampvis2," and "vegan." Alpha diversity indices were calculated using the function "estimate_richness" from phyloseq with a data set rarefied to the minimum sample size (8,740 sequences) (49). After assessing normal distribution by qqplots, histograms, and the Shapiro-Wilk normality test, the not normally distributed groups were compared using the Wilcoxon signed-rank test. Beta diversity was calculated with "amp_ordinate" in the ampvis2 package (50). Data have been Hellinger transformed, and ASVs that were not present in more than 0.1% relative abundance in any sample have been removed. The principal-coordinate analysis was based on Bray-Curtis distances. The permutational multivariate analysis of variance (PERMANOVA) was calculated with the "adonis" function and 5,000 permutations in vegan (51). All raw FASTQ reads for 16S rRNA gene sequences have been deposited under NCBI accession number PRJNA669215.
Whole-genome sequencing and processing. Illumina short read sequencing was performed on DNA extracted from S. Heidelberg isolates recovered from ceca and litter. In addition, five E. coli isolates recovered from two cecal samples were sequenced. Libraries were prepared using either Nextera XT or Nextera DNA flex library preparation kits (Illumina, Inc., San Diego, CA) following the manufacturer's protocol. Libraries were sequenced on the Illumina MiSeq platform with 150-or 250-bp paired-end reads. Additionally, five S. Heidelberg and two E. coli isolates were selected for long-read sequencing using the Sequel II system (PacBio Biosciences Inc.) or MinION device (Oxford Nanopore Technology).
Preparation and sequencing of long read libraries were done by next-generation sequencing core centers of University of Georgia and Colorado State University. MinION sequencing of samples was done using R9.5 chemistry on a 1D flowcell. Lima (https://github.com/PacificBiosciences/pbbioconda) was used to demultiplex barcoded PacBio samples and to convert the split BAM files into FASTQ format. Genome assembly, resistome characterization, and quality assessment of long reads were done using Reads2Resistome pipeline v1.1.1 (53). Reads2Resistome performed hybrid assemblies using either Illumina reads and PacBio or Illumina and MinION reads, using both Unicycler (54) and SPAdes (55) (see Data Set S5 for assembly statistics). In addition, hierarchical genome assembly process (HGAP) assembly was performed using the PacBio single-molecule real-time (SMRT) Link v9.0 analysis software suite using default settings. Assembly quality was assessed by QUAST v5.0.2 (56), using default settings, and genome annotation was done using Prokka (57), Rapid Annotation using Subsystem Technology (RAST) (58), and BlastKOALA (59). We confirmed that all S. Heidelberg isolates were Salmonella enterica serovar Heidelberg using Salmonella In Silico Typing Resource (SISTR) (26) and used ClermonTyping (60) and multilocus sequence typing (MLST) (61) for E. coli strains.
For resistome characterization, Reads2Resistome used ARG-ANNOT (62), the Comprehensive Antibiotic Resistance Database (63), MEGARes (64), AMRFinderPlus (65), PlasmidFinder (19), ResFinder (66), and VirulenceFinder database (67). For plasmid typing and IncI1 clonal complex determination, we used plasmid MLST (19). PHAST (68) was used to identify prophages present in chromosomal contigs, and the predicted prophage DNA sequence was annotated with RAST and BlastKOALA. progressiveMAUVE v1.1.1 (69) and MAFFT v1.4.0 (70) implemented in Geneious Prime v2020.0.1 were used for aligning and comparing sequences. In addition, a pangenome analysis of annotated assemblies was conducted with Roary (71). Phylogenetic trees based on the core genome and accessory genome were reconstructed using the maximum likelihood (ML) method implemented in RAxML-NG 2.0.0 (72). When computationally possible, the best model of sequence evolution predicted by jModelTest (73) was used for tree reconstruction; otherwise, the GTR1GAMMA model was implemented. Lastly, we used the bacterial genomes sequenced in the study to create a BLAST database that can be searched for sequences of interest in Geneious Prime.
Illumina short reads were assembled de novo into contigs using Unicycler v0.4.7 and characterized with bioinformatic tools described for long reads. Illumina reads were used to determine single nucleotide polymorphisms (SNPs) and indels present in S. Heidelberg isolates. Alignment of raw FASTQ reads to the genome of the S. Heidelberg ancestor was done using Burrows-Wheeler aligner (BWA) (74), and SNPs/ indels were called using the Genome Analysis Toolkit (75) as described previously (48). Variant call format (VCF) files of identified SNPs/indels and the Linux/Unix shell script used have been deposited in Dryad Digital Repository online at https://doi.org/10.5061/dryad.4tmpg4f8d. The SNP-based ML tree was reconstructed by converting merged VCF files to the PHYLogeny Inference Package format using PGDSpider (76). Afterward, isolates with duplicated SNPs/indels were removed, and an ML tree was drawn using the Jukes-Cantor model of nucleotide substitution and GAMMA model of rate heterogeneity.
To determine the IncI1-pST26 consensus, first we used BWA (74) or Bowtie v 7.2.1 (77) (implemented in Geneious Prime) to align raw FASTQ files against a complete circular p1ST26 or p2ST26 and used the binary alignment map file generated for consensus determination in Geneious Prime. To identify identical IncI1-pST26 plasmids, we performed a BLAST search against the NCBI nonredundant database and selected the top 50 plasmids matching p2ST26 from this study. To construct a whole plasmid-based maximum likelihood tree, we downloaded the FASTA files for the top 50 plasmids from NCBI and used them for whole-genome alignment. All raw FASTQ reads for sequenced bacterial genomes are publicly available under NCBI accession numbers PRJNA683658, PRJNA684578, and PRJNA684580.
Hi-C and cecal metagenome library preparation and analysis. A shotgun and Hi-C DNA library of two cecal samples was created using the Nextera XT library preparation kit and Phase Genomics (Seattle, WA) ProxiMeta Hi-C Microbiome kit following the manufacturer's instructions. Library sequencing was performed by Novogene corporation (Sacramento, CA, USA) on the Illumina HiSeq platform using 150-bp paired-end reads. Two libraries were sequenced per HiSeq flow cell lane. Metagenomic FASTQ files were uploaded to the Phase Genomics cloud-based bioinformatics portal for subsequent analysis.
Shotgun reads were filtered and trimmed for quality using bbduk (78), normalized using bbtools (78), and then assembled with metaSPAdes (79) using default options. Hi-C reads were then aligned to the assembly following the Hi-C kit manufacturer's recommendations (https://phasegenomics.github.io/ 2019/09/19/hic-alignment-and-qc.html). Briefly, reads were aligned using BWA-MEM with the -5SP and -t 8 options specified and all other options default. SAMBLASTER (80) was used to flag PCR duplicates. Alignments were then filtered with SAMtools (81) using the -F 2304 filtering flag to remove nonprimary and secondary alignments. Lastly, Hi-C read alignments to the assembly were filtered using Matlock (https://github.com/phasegenomics/matlock) with default options (removing all alignments with MAPQ of ,20, edit distance greater than 5, and read duplicates). Metagenome deconvolution was performed with ProxiMeta (82,83), resulting in the creation of putative genome and genome fragment clusters.
Metagenome-assembled genomes/clusters were assessed for quality using CheckM (84) and assigned preliminary taxonomic classifications with Mash (24) and GTDB-Tk (27) (database release 202). Quality-controlled shotgun reads were additionally classified using Kraken2 (v2.0.8 beta) (28). To identify the plasmids and ARGs present in metagenome-assembled genomes, we used Minimap2 (85) to align assembled metagenome contigs to the PlasmidFinder and ResFinder database. We used a 75% identity threshold for the alignment match. To search the metagenome for ARG and plasmids of interest, we used Geneious Prime to create a BLAST metagenome database. Metagenome contigs with significant homology (E value of ,0.05 and at least 1,000 bp of aligned sequence) were considered contigs of the respective plasmid or ARG. Hi-C data were then used to link identified sequences to host genomes and genome fragments within ProxiMeta (23). Phylogenetic visualizations of Hi-C linkages used the placement of clusters in a large prokaryotic phylogeny as estimated by CheckM (84) and Mash (24). Shotgun and Hi-C reads are publicly available under NCBI accession number PRJNA688069.
Isolation of E. coli from ceca of broiler chicks. To isolate E. coli, we retrospectively screened two cecal slurry samples from chicks challenged intracloacally on CHROMagar plates supplemented with gentamicin (8 ppm) and tetracycline (8 ppm). Frozen vials of cecal contents were thawed on ice and vortexed. For each sample, a 10-ml aliquot was struck for isolation, and a 100-ml aliquot was spread plated onto CHROMagar with relevant antibiotics. Plates were incubated overnight at 37°C. Blue-green and blue-cream colonies were counted as presumptive E. coli, and colonies were streaked again for isolation on CHROMagar supplemented with gentamicin and tetracycline. After incubation, 2 to 3 colonies from each sample were selected to represent different forms of blue-green/blue-cream colonies, streaked in parallel onto sheep blood agar (SBA) and eosin methylene Bbue (EMB) (Remel Inc., KS, USA), and incubated overnight. Growth on EMB was the characteristic green metallic sheen after incubation. Colonies from SBA were stored in 30% LB glycerol at 280°C and were used for AST and WGS.
Bioinformatic tools used for visualization of high-throughput sequence data. The gene presence/absence heatmap, with gene presence/absence data obtained from Roary, was generated using the pheatmap v1.0.12, tidyverse v1.3.0, and viridis v0.5.1 packages in R v4.0.2. BLAST Ring Image Generator (86) was used for genome comparison visualization, including GC skew change, and Phandango (87) was used for visualizing phylogenetic trees with their associated metadata. PHYLOViZ 2.0 (88) was used to generate a minimum spanning tree of S. Heidelberg isolates from core genome information obtained through SISTR (i.e., wgMLST_330:complete-alleles, wgMLST_330:missing-alleles, wgMLST_330:partial-alleles, closest-public-genome-alleles-matching, and cgMLST_cluster_level). SnapGene was used for drawing linear maps of the plasmid, prophages, and other regions of interest in bacterial genomes. Geneious Prime and SnapGene were used to view MAFFT alignments of DNA and amino acid sequences.
Determining S. Heidelberg metabolism under selective pressure. To determine if exposure to metals, disinfectants, or acidic pH poses a selection on S. Heidelberg strains carrying IncI1 plasmids, we used 96-well phenotype microarray (PM) MicroPlates (PM10, PM12B, PM13B, PM15B, and PM16A) (Biolog, Inc., Hayward CA) to compare the metabolic profile of one gentamicin-resistant strain harboring IncI1 and one susceptible strain. PM plates use cell respiration via NADH production to determine cell metabolic activity. If the phenotype is positive in a well, the cells respire actively, reducing a tetrazolium dye and forming a strong color (Biolog, Inc.). Briefly, S. Heidelberg strains were subcultured twice in universal growth agar (Biolog Universal Growth Medium 1 5% sheep blood) for 24 h at 37°C. Cells were removed with a sterile swab and transferred to 16 ml inoculating fluid 0 (IF-0) to achieve a final turbidimetric transmission of 42% transmittance (42%T). Afterward, 15 ml of the 42%T cell suspension was transferred to 75 ml of IF-01 dye to achieve a final transmission of 85%T before 600 ml of the cell suspension was transferred to 120 ml IF-101 dye. The suspension was mixed, and each well in a microplate was inoculated with 100 ml of the suspension. After inoculation, microplates were covered with a sterile plastic film and monitored automatically for color development every 15 min for 48 h at 37°C using and an OmniLog reader for 48 h. To identify phenotypes, the kinetic curves of the gentamicin-resistant and -susceptible strains were compared using OmniLog PM software. For each PM plate, the respiratory unit for each well at 22 h was extracted and log 10 transformed. All supplies used were purchased from Biolog, Inc.
Exposing S. Heidelberg to acidic pH. To determine the fitness of evolved S. Heidelberg strains that acquired antibiotic resistance compared with that of evolved susceptible strains after exposure to acidic pH, we exposed gentamicin-resistant (n = 3) and -susceptible (n = 3) strains recovered from the litter of chicks gavaged with S. Heidelberg to acidified filter-sterilized pine shaving extract (PSE). PSE was prepared as described for poultry litter extract (47), using fresh pine shavings collected from the University of Georgia, Poultry Research Center, Athens, GA, USA. Three single colonies of each strain, including the ancestral SH-2813nal R , were selected from overnight cultures grown on sheep blood agar and transferred to a microcentrifuge tube containing 900 ml of PSE (pH 6.52), i.e., one tube per strain.
After transfer, tubes were vortexed, covered with a gas-permeable paper strip, and incubated at 41°C under microaerophilic conditions (5% O 2 , 10% CO 2 , and 85% N 2 ) for 2 h. After incubation, tubes were vortexed, and 100 ml of the suspension was transferred to a microcentrifuge tube containing 900 ml of PSE (pH of PSE was adjusted to 2.54 using 1 M HCl [Spectrum Chemical Mfg. Corp., CA, USA] and 1 M NaOH [Fisher Chemical, NJ, USA]) Afterward, tubes were vortexed and incubation was continued for another 2 h as was done for PSE at pH 6.5. To determine the S. Heidelberg population in PSE, one replicate tube per strain was removed and serially diluted in 1Â PBS at time points 0 and 2 h for pH 6.5 and 0.5 and 1 and 2 h for pH 2.5. Serial dilutions were plated onto BGS agar plates supplemented with or without 16 ppm gentamicin. S. Heidelberg colonies were counted 18 to 24 h after incubation using a calibrated automated colony counter.
The fitness of each evolved S. Heidelberg population relative to the SH-2813nal R population was determined after 2 h at pH 2.54 as described by San Millan et al. (89) where W strain is the fitness of the evolved susceptible (total susceptible) or gentamicin-resistant populations (total resistant and gen R ); N initial , evol and N final , evol are the numbers of cells (in CFU) of the evolved susceptible or gentamicin-resistant population at time point 0 and 2 h after exposure to PSE; and N initial , anc and N final , anc are the numbers of cells of SH-2813nal R population at time point 0 and 2 h. Statistical analyses. Continuous variables did not meet the assumption of a normal distribution; therefore, nonparametric testing for direct comparisons was performed using Wilcoxon rank-sum and signed-rank tests, and the Kruskal-Wallis rank-sum test was used for one-way analysis of variance tests. Furthermore, continuous variables were log-transformed before any statistical tests were performed. Statistical analyses were performed using R (v3.4.1).
Ethics statement. All animal experiments were approved by the University of Georgia Office of Animal Care and Use under Animal Use Protocol A2017 04-028-A2.
Data availability. All raw FASTQ reads, including short and long reads for sequenced bacterial genomes, are publicly available under NCBI accession numbers PRJNA683658, PRJNA684578, and PRJNA684580. Shotgun and Hi-C reads are publicly available under NCBI accession number PRJNA688069 and 16S rRNA gene sequences under NCBI accession number PRJNA669215. The whole-genome assemblies for S. Heidelberg strains SH-2813-ancestor and ic9b harboring p1ST26 have been made available under NCBI GenBank accession number CP066851 and DDBJ/ENA/GenBank JAEMHU000000000, respectively. The whole-genome assemblies for E. coli strain Ec-FL1-1X and Ec-FL1-2X carrying IncI1 are available under GenBank accession numbers CP066836 and JAFCXR000000000, respectively. Variant call format (VCF) files of identified SNPs/indels and the Linux/Unix shell script used have been deposited in Dryad Digital Repository online at https://doi.org/10.5061/dryad.4tmpg4f8d.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.

ACKNOWLEDGMENTS
We are grateful to Marlo Sommers, Jasmine Johnson, Carolina Hall, Jeromey Jackson, and Latoya Wiggins for their logistical and technical assistance. This work was supported by USDA Agricultural Research Service (project number 6040-32000-010-00-D), nonassistance cooperative agreement (58-6040-6-030) between USDA Agricultural Research Service and University of Georgia, Research Foundation, and research service agreement (58-6040-8-035) between USDA Agricultural Research Service and Colorado State University. M.O.P. was supported in part by grant R44AI150008 from NIAID to Phase Genomics. B.Z. was supported by the competence center. The competence center (FFoQSI) is funded by the Austrian ministries BMVIT and BMDW and the Austrian provinces Niederoesterreich, Upper Austria, and Vienna within the scope of COMET-Competence Centers for Excellent Technologies. The Austrian Research Promotion Agency FFG handles the program COMET. This study was supported in part by resources and technical expertise from the Georgia Advanced Computing Resource Center, a partnership between the University of Georgia's Office of the Vice President for Research and Office of the Vice President for Information Technology.