Mapping Microbiome Signatures Associated With Antibiotic Usage On A Swine Production Farm Using Amplicon And Shotgun Data

Background: High antimicrobial usage in swine production has the potential to create reservoirs of antimicrobial resistance (AMR) genes which are transferable to human pathogens via mobile genetic elements. Understanding microbial community responses to antibiotic use is central to unravelling transfer of such resistance genes. Our previous investigation revealed a scenario of optimal antibiotic activity associated with saturation of AMR genes on this farm. Here, we use amplicon and shotgun sequence data to investigate the microbiome signatures that underwrite such a phenomenon. Results: We generated 1.24 and 576 million high quality 16S rRNA gene amplicon and shotgun sequences from 24 porcine faecal samples, respectively. The ratio of taxa detection at genus level between the two methods was 1:24. Using shotgun sequence data, 235 unique AMR genes, 122 modes of action and 17 antibiotic classes were identied using the MEGARes AMR database. Antibiotic usage in growing pigs was signicantly associated with microbial and AMR resistome structural and compositional changes detectable two weeks after antibiotic initiation. These were characterised by a down regulation of MDR eux pumps and an up regulation of macrolide-specic eux pumps in the growing pigs (treated-group) linked to lower abundance of Verrucomicrobiaeceae. In the sows (non-treated group), a potentially undetected infection, was characterised by a high abundance of pathogenic viral sequences, microbial structural changes i.e. family Alcaligenaceae, and an up regulation of beta-lactamases, including MDR eux pumps. We assembled 682 near complete bacterial genomes revealing that a large proportion of the resistome is carried by Firmicutes and Proteobacteria, specically multi-class gene carriage by Clostridium species and Escherichia coli, which occurred exclusively in the treatment group. Conclusion: Microbiome signatures i.e. microbial structure, composition and resistome carriage associated with antibiotic-use can be cost effectively screened with amplicon sequencing but their granularity unravelled using shotgun metagenomic data.


Background
Antibiotic usage is integral to livestock production globally [1,2]. Until 2006, drugs such as tylosin and tetracycline were used to improve feed conversion e ciency in swine and poultry production systems in Europe [3]. Today, their use is restricted to therapeutic and where justi able, prophylactic use, with the aims of limiting welfare and production impacts of infectious diseases [4]. However, global projections indicate that twothirds of the growth in antimicrobial usage is expected to be within the veterinary sector, with a doubling in use for pigs [5]. This raises concerns about the role pig farms can play as reservoirs of antimicrobial resistance (AMR) genes that could be transferred to human and pig pathogens via mobile genetic elements [6]. Unravelling the microbial characteristics that underwrite these genetic transfer mechanisms is important to inform the development of control strategies [7], however, recent evidence suggest that the vast majority of unexplored AMR is harboured by non-pathogenic and non-culturable bacteria present in human and animal gut microbiomes [8,9]. Advances in sequencing technologies [8,10] and improved e ciency in data processing have provided access and ability to examine this non-culturable fraction in epidemiologically relevant settings. However, there are few longitudinal studies that have interrogated AMR gene generation and dissemination in livestock production systems to date. Such approaches are essential to disentangle the contribution of temporal and contextual microbiome variation in high antibiotic use systems such as swine production [11]. Our previous work [11] indicated a state of saturation of AMR genes in a swine production system, with no apparent impact on the utility of antibiotics to limit infections and mortality. Here, we seek to unravel the microbiome characteristics and the association with antimicrobial resistance genes, referred to as the resistome, in the pig faecal samples on this farm. To do so, the resolution of shotgun sequence data was exploited and this allowed comparison with 16S rRNA gene sequence data from the same samples. We targeted four critical time points as part of the pig production cycle, i.e. at farrowing, weaning, growing and nishing where batch antibiotic treatment (metaphylaxis) is considered necessary to limit the impact of infectious diseases. It is noteworthy that welfare concerns precluded the use of within group controls, therefore sows on the same farm were sampled at the same timepoints as the production animals as a comparative control for changes which were also analysed over time within groups.

Farm description
This experiment was approved by the Royal (Dick) School of Veterinary Studies Veterinary Ethics Research Committee. A 600 sow (Landrace x Large White) commercial farrowing to nishing unit in the United Kingdom was recruited for this study. On this unit, piglets (antibiotic treated group) were batch farrowed every four weeks and the sample collection commenced one week before the October 2016 batch farrowing process. Dry sows (Nontreated group) were housed in straw yards (referred to as the "sow barn") and nursing sows and piglets were housed on slats. Detailed herd and antimicrobial usage records have been published in our previous work [11].
Study design and sampling The main study design and sampling procedures were reported in our previous publication [11]. Here, we focus on a subset of faecal samples collected at four time points − 9th November 2016, 30th November 2016, 14th December 2016 and 5th April 2017 -labelled as T1, T2, T3 and T4, respectively. These time points were selected on the basis of differing growth stages and antimicrobial treatment phases during the production cycle of the growing pigs. In parallel, faecal samples from the dry sow barn (where routine antimicrobial administration did not occur) were taken at each of these time points (n = 3) as shown in Fig. 1a. Faecal samples were collected in spooned universal tubes prior to being stored immediately at -20 o C onsite. Samples were batch transported to the laboratory on dry ice for DNA extraction and dry matter calculations, as described previously [11]. Here we refer to this analysis as the selected amplicon resistome.

Sequencing
The extracted DNA was submitted to Edinburgh Genomics (United Kingdom) for shotgun metagenomic sequencing. Brie y, Illumina TruSeq DNA libraries were prepared from the DNA extracts on the HiSeq 4000 platform to generate 150bp paired end reads (Illumina, United States). The raw sequences for this study are publicly available on EMBL-EBI under accession number PRJEB34736. Reads generated from host DNA and Phix (Phix174) were ltered out by mapping the reads to the Sus Scrofa reference genome GCA_000003025 version 11.1 and Phix174 respectively using the run_contaminant_ lter.pl script in Microbiome Helper suite [12]. Adapters and low-quality reads were then trimmed and removed respectively using Trimmomatic [13], the ltered reads were then taken forward for analysis. Detailed description of the 16S rRNA gene sequencing is reported here [11] but in brief; the paired end reads from the sequencing step above were processed using the Quantitative Insights Into Microbial Ecology-2™ (QIIME-2) software. The reads were then de-replicated and chimeric sequences removed before the denoising was done in DADA2 with output here as a feature table used for the downstream analysis. OTU taxonomic classi cation was done using a naïve Bayes classi er trained on the SILVA database (version March 2020) at 97% similarity. Visualisation of taxonomic abundances and alpha diversity was done using ggplot [14] and plotly in R.
We also track the copy number of four selected AMR genes i.e. tetQ, tetB, ermB and dfra1 as using qPCR across time as earlier reported [11].

Mapping contigs to AMR, plasmid and taxonomic databases
In order to maximize the amount of sequences used, rather than annotating full genomes, we mapped the Megahit-assembled contigs to databases and then lter output using stringent thresholds. The aim here is to detect contigs carrying AMR genes (full resistome), plasmid and identi able taxonomic sequences. Whilst this approach would not provide details of taxa carrying genes lower than Phylum, it allows us to see a more comprehensive picture for the composition of these genetic elements on a sample-wide basis. To identify contig carrying AMR genes ,we used Abricate[16] to map the contigs against the MEGARes[17] -https://megares.meglab.org. Here the output was ltered by contig length (n > = 800bp) and percentage identity ( > = 75%).
Contig taxonomic assignment -Further taxonomic identi cation of contigs was performed using the CAT and BAT pipeline[18] which identi es assembled sequences up to genus level. Here, contigs were ltered according to length and retained taxonomic assignment based on at least three known markers [18]. This taxonomic assignment was validated using the Plas ow pipeline [19] which also predicts plasmid sequences and taxonomy at phylum level.
Merging output for analysis-The AMR gene and taxonomic mapping and assignment process generated a tab separated value (tsv) le with eight standard BLAST output columns for the former http://www.metagenomics.wiki/tools/blast/blastn-output-format-6 Using the UCI, AMR gene and taxonomic assignment data was merged into one le for downstream data analysis.
Analysis and Visualization-To visualise the entire resistome we use geom-tile in ggplot [14] while the taxonomic hierarchical clustering of 16S rRNA gene and shotgun data was done using the complexheatmaps [20] package in R. To detect differential abundance in the resistome and microbial populations, the negative binomial additive model in Deseq2 was applied [21]. Note that analysis of taxonomic abundance was done for bacterial, viral and eukaryotic sequences by ltering at kingdom level.

AMR gene content inference
We developed a mixed effects Poisson regression models [22] implemented using the lme4 package [23] in R to analyse within and between groups variation in abundance of AMR genes at the four time points. To do so we make the following assumptions as showing in (Fig. 1b); a) The piglet's microbiome gut will exhibit more variation over time than the sows, b) sows will have a higher absolute count/abundance and richness of microbiome units (taxa or genes), c) Therefore, the gut microbiome changes (Δm = Δx + Δy) will be a function of the change in temporal signal (Δx), developmental changes and routine activities such as antibiotic use (Δy). We then de ned the outcome as the read-normalised gene abundance. To compare between groups, we created subsets of the data i.e. Piglets and Sows (Farrowing + Dry), then developed a mixed effects Poisson model for each data set, with the xed component as antibiotic class (AC) and time (T) point while the random component using the sample identi cation (SI). On its own, the random component represents within sample variation, i.e. the total variance in gene count for a sample including developmental effects, antibiotic use and seasonal variability) (Fig. 1b). The difference between the total variance and that which is explained by the xed component is the unexplained variance. Comparing this between groups allows us to partition the variance and account for the antibiotic usage. For example; Here we assume that most of the unexplained variance is likely due to antibiotic use. To compare the outputs from the models from the subsets above, we use jtools [24]. To assess the validity of the model, we assess if the data supports our assumptions. In addition, we use the modes of action [17] as the label in the same model to infer the functional characteristics of the gene abundance estimates. Model interpretation; a) Statistically signi cant variables in xed component of the model for a dataset represented explanatory factors for within group variations, b) variables that were statistically signi cant between models (between datasets) represent factors with an overlapping effect between groups and c) distinct lack of overlap of variable estimates between group (between datasets) represents factors with a differential affect at group level.

Results
Sample sequence characteristics on amplicon and shotgun 16S rRNA gene sequencing of 24 faecal samples generated 1.55 million 150bp paired-end reads, 95% of these yielding at least Q30 values per sample. Assembly and ltering by sequencing depth (4000 sequences) retained 1.24 million sequences, which clustered into 1284 operational taxonomic units, 12 phyla and 141 genera. The mean alpha diversity based on the Shannon index was 2.2 across all the samples, 2.6(SD = 1.8-3.1) for piglets, 2.3(SD = 1.28-3.0) for dry sows and 1.8(SD = 1.3-1.8) for farrowing sows. With the exception of T3, alpha diversity was stable among sows but a signi cantly higher diversity for piglets at T3 and T4 (Fig S1).
Shotgun sequencing generated 456 gigabytes of data, i.e. 27.8 million sequences of 150bp paired end reads per sample. After quality control, each sample on average retained 25 million sequences, which were assembled into 9.2 million contigs. 47% of these contigs were ≥ 800bp of which only 9434 mapped to the MEGARes AMR database with a percentage identity ≥ 75%. On the other hand, 650,368 unique contigs of ≥ 800bp were taxonomically classi ed to at least genus level.

Taxonomic pro le
In general, the taxonomic pro les were comparable between the two methods but the difference in granularity was evident at family level (Fig 2A & C), i.e. The ratio of taxa detection at genus level between the two methods was 1:24, none the less, there were some discernible differences between pig groups across time (Fig 2B). Using 16S rRNA at phylum level, Firmicutes, Bacteriodetes and Proteobacteria constituted nearly 94% of gut microbiota with a marked increase in Bacteriodetes for dry sows at T3. The granularity of these is revealed by shotgun sequence data (Fig 2D), however, here the change in Bacteriodetes abundance was less dramatic and showed that same phylum was differentially higher among piglets at T3. Moreover, Shotgun sequence data revealed distinct clustering by time and group, for example; cluster 4 (Fig 2 A & C) represented samples with the highest richness and these exclusively belonged to dry sows at T3. Cluster 3 ( Fig 2C) exclusively belonged to piglets at T4, therefore, richness peaks at different time points in piglets and sows. This observation supports our model assumptions. Cluster 1 exclusively belonged to piglets distinctly separated into T3 and T2 while cluster 2 mostly belonged to sows, this too reveals distinct clusters by time. Overall, this clustering indicates a strongtaxonomic signal associated with group and time.
At family level, Shotgun data shows that Bacteriodaceae, Verrumicrobiaceae and Bradymonadaceae were differentially abundant between piglets and sows ( Fig S2a), with the Verrumicrobiaceae and Bradymonadaceae more abundant among sows (Fig S2b). Analysis of viral sequences revealed a steady increase in viral sequences in piglets over time, while that of sows was generally stable (Fig S3a and b). This was with exception of the dry sows at T3 which exhibited a marked increase in viral sequences which mapped to Siphoviridae, Phycodnaviridae and Podoviridae. Similarly, the hierarchical clustering of these viral sequences also distinctly separates dry sows at T3, however these characteristics were not mirrored by the Eukaryotic sequence data (Fig S3c).
The antibiotic resistome Selected amplicon resistome presented as AMR gene copy number normalised by 16S rRNA gene copy number [11] is shown in Fig. 3A. At farrowing (T1), the sows carried signi cantly more tetQ and ermB gene copies with the former exhibiting the least variability in gene copy (P < 0.001). This trend was mirrored by dry sows at the same time point. Among piglets, the abundance of tetQ and ermB genes was generally higher than sows across the three time points but the highest level was observed at T4 (P < 0.01), speci cally, the abundance of these two genes was 2-fold higher than that for dry sows. T3 was characterised by limited variability across all the ve genes but the copy number of tetQ and ermB remained markedly higher among piglets however, this was not the case for dfra1. Analysis of the same genes using shotgun sequence data Fig. 3B showed similar levels of tetQ across time and groups. The trend of tetQ is similar to that observed with amplicon sequence data, especially at T3 and T4.
Full shotgun sequence resistome is shown in Fig. 3C and is comprised of 235 unique genes clustering into 17 antibiotic classes with 122 corresponding modes of action. In piglets, the abundance of genes encoding resistance for tetracycline generally increased over time; speci c examples include tetW, tetQ, tetO, tet44 and tet32. By contrast, all of these genes had exhibited the opposite trend among sows (Fig S4a). This too supports our model assumptions. Macrolides, lincosamides and streptogramines (MLS) resistance-encoding genes exhibited comparatively low variation and all but lnuc increased over time (Fig S4a). In general, T4 was characterised by an increase in mphB, macB and ermTmacrolide resistance genes -in piglets. Beyond the administered antibiotics, uoroquinolone and glycopeptide-associated resistance genes also exhibited a gradual increase in abundance over time. For example; among sows, the abundance of genes that encode resistance for uoroquinolones (parc gyrA and gyrB), beta-lactams (imp and carB), aminoglycosides (ant3-prime) and CCR (sme,mexW,mexK & mexB) differentially increased at T3 (Fig S4b) AMR gene count inferences Almost all attributable variance at group level belonged to piglets, and after accounting for the temporal effect, approximately 66% of the variance remained explained, which is likely the effect size due to antibiotic usage in this group. There was a signi cant difference in within and between groups at T4 ( Fig. 4A and Fig S5). The changes in gene count at T3 were statistically signi cant within the dry sows.
At antibiotic class level, salient differences in genes within-groups were detectable for all but cationic antimicrobial peptides, metronidazole in piglets and tetracyclines in sows. It is noteworthy that the difference in variance between sows and piglets also validates our model assumptions ( Fig. 1 &  Fig. 4b). In addition, Fig S3b suggests a gradual increase in abundance and diversity of viral sequences in piglets over time. Finally, although the resistome (Fig. 3C) shows an increase in diversity of glycopeptide resistance genes, it did not translate into a signi cant differential abundance between piglets and sows.
When we consider the mode of action for resistance rather than the class of antibiotic (Fig. 4C), the results also indicate differential abundance within groups at T3 and T4, but a clear difference between groups at T3. Generally, the abundance here was skewed toward sows, but multi drug e ux pumps, polymyxin B resistance regulator and class B beta-lactamases were markedly higher in sows (Fig S6). Among piglets, the abundance of dihydrofolate reductase, vancomycin G type resistance protein, penicillin binding protein, MDR regulator and aminocoumarin-resistant DNA topoisomerases signi cantly varied (Fig. 4C & S6). On the other hand, EF-Tu inhibition, uoroquinolone DNA resistant topoisomerases, and ATPbinding cassettes ABC e ux pumps differentially varied within sows over time. Finally, mechanisms such as resistant 30S ribosomal subunit protein S12, class A beta-lactamases, macrolide resistance e ux pumps and tetracycline ribosomal protection protein were differentially abundant within both groups.

Antibiotic resistance gene carriage by taxa
Using shotgun metagenomics, we observe that the microbiome structure based on resistome and microbial taxonomy indicate that groups and time point accounted for 63% and 60% of the variation respectively (Fig S2a, i and iii). This observation highlights a direct relationship between microbiota and its resistome. Indeed, the clusters observed here also map to groups at speci c time points. To unravel this relationship, 2,674 genomes were assembled from 24 metagenomes, and after quality control, 638 genomes (Fig. 5) were retained. Of these, 599 (93.9%) were taxonomically classi ed as bacteria and 10 (1.6%) as archaea at kingdom level. Of the bacteria, 29 (4.6%) were unclassi ed and therefore represented potentially unculturable and /or novel prokaryotic taxa. At the lower taxonomic ranks i.e. family and genus levels, there represented 163 (25.6%) and 231 (36.21%) of the microbiota respectively. In the phylogenetic tree (Fig. 5), most of these unclassi ed phyla clustered with Bacteriodetes, Firmicutes and Lentisphaerae. Similarly, Bacteriodetes, Firmicutes and Proteobacteria accounted for 75.4 % of observed taxa. We observed that AMR genes were generally carried by a wide range of bacteria across phyla but a small cluster carried genes across class of antibiotics. This cluster was exclusively observed in piglets at T2 and T3, it mostly contained phylum Proteobacteria and dominated by Escherichia coli and Clostridium species

Discussion
In this study, we assess the relationship between the resistome and microbiota, exploiting the resolution provided by shotgun sequence data. Here, we compare the granularity of two data sets -a) contigs and genomes assembled from shotgun metagenome sequencing data and b) OTUs from 16S rRNA gene sequence data as well as gene copy count using qPCR. We target four critical time points on a swine production farm, i.e. farrowing, weaning, growing and nishing, where batch antibiotic treatment (metaphylaxis) was administered. The analysis was motivated by our previous work [11], which concluded that AMR gene levels on this farm were in a state of equilibrium, or saturated, but the antibiotics remained relatively effective in reducing the frequency of infections and mortality. We therefore sought to unravel the microbiome characteristics that underwrite such a phenomenon using amplicon and shotgun sequence data. In doing so, we generate critical insights needed to design strategies that limit welfare impact of disease whilst controlling the spread of antimicrobial resistance on farms.

Characteristics of the gut resistome
The sow microbiome exhibited limited variation in abundance and diversity of antibiotic resistance genes, therefore most of the variation was attributed to piglets [35]. Partitioning this variance suggests that up to 66% may be due to chlorotetracycline and tylosin usage in these piglets. Resistome variability in piglets has previously been reported [35] and its reduction as piglets mature is in part due to natural microbial succession [1,35]. In this study, we were unable to fully delineate between changes due to the latter and antibiotic usage. Tylosin usage in piglets is known to accelerate microbial succession [35]. However, unravelling this speci c aspect was beyond the scope of this study. Fundamentally, we show the direct relationship between microbiota and its resistome, furthermore, the characteristics of the resistome, mode of action for resistance and taxonomic pro les have led us to hypothesise that the marked changes seen at T3 in both piglets and sows may have different precursors.
Changes due to chlorotetracyline and tylosin use Antibiotic-driven changes of the pig gut microbiome take up 14 days to manifest and may last for months [36,37]. Here, the piglets had been administered chlortetracycline and acidi ed water at T2, but the dramatic changes in the selected amplicon resistome, i.e. abundance were observed 28 days later at T3. This may indicate a delayed response attributable to acidi ed water. At T3, post chlortetracycline administration the full resistome was characterised by an increase in diversity and abundance of tetracycline resistance encoding genes (P > 0.01), for example, tet44, tetW and tetO. Indeed this is re ected as an upregulation of tetracycline ribosomal protection protein was observed, which works by protecting the bacterial ribosome from binding the antibiotic tetracycline [38]. All of this points towards tetracycline-driven changes as reported elsewhere [39].
Beyond this, there was a down regulation of the MDR regulator mechanism which modulates resistance to aminoglycosides and aminocumarins, suggesting a counter selection. This would be expected to impede the functionality of e ux pump mechanisms which confer multi-drug resistance via the bacterial two-component signal transduction system CpxAR [40,41]. Such an effect would likely compromise the survival of certain microbial populations [42]. Finally, the resistome of piglets at T3 was also characterised by an increase in diversity of genes that encode resistance to glycopeptides and metronidazole, although this indicates co-selection, the mechanism behind this remains unclear. We had expected that the use of antibiotics would result in a reduction of both microbial and resistome diversity, however the use of chlortetracycline was associated with a signi cant increase in Shannon index [39], suggesting an increase in unique taxa and their abundance. Indeed, this increase is exhibited by a signi cant reduction of Verrumicrobiaceae and an increase of Bacteroidaceae. Low abundance of the former in piglets under similar conditions has been reported elsewhere [43].
Earlier studies have shown that sub therapeutic administration of tylosin in pigs resulted in a much larger effect than chlortetracycline [35,43,44]. In this study, therapeutic levels were used (265.8 mg/PCU tylosin, 103.2 mg/PCU chlortetracycline), and most of the variance in piglets was at T4 after tylosin administration. This could represent a combined effect of both antibiotics, however, if we assume that most of variation attributable to tetracycline was manifested at T3, then tylosin would account for a much larger effect [35,44]. In this case, counter-selection i.e. chlortetracycline causing a suppression of e ux pump mechanism, may have resulted in microbial population that are more susceptible to antibiotics, leading to a much larger effect when tylosin was subsequently used. The response to this was an upregulation of e ux pumps speci c for macrolides and this supports export of chemicals associated with macrolide activity out of the bacterial cell [45]. Interestingly, ermB and lnuC remain abundantly stable across the three time points, which is possibly an indication of legacy use of this antibiotic on this farm. It is also possible that the full extent of tylosin's impact on this system was in the weeks beyond our selected time points.
The changes to microbiome structure and taxonomic distribution at T4 are comparable with that observed for the resistome, for example; although lower than T3 the Shannon index at T4 is signi cantly higher than in sows. It is likely that changes at T4 were more taxa-speci c, which is why they account for signi cantly more variation in the resistome [35,44]. To this effect, a reduction in the abundance of phyla Bacteriodetes, Actinobacteria and Proteobacteria but an increase in Firmicutes might explain the level of change.

Changes due to a potentially undetected infection in sows
The sows did not receive antibiotics during this experiment, however they exhibited marked microbiome structural, resistome and taxonomic changes indicating an undetected incursion. Certain characteristics point towards an infection. Firstly, despite comparable microbial diversity, resistome and taxonomic characteristics were markedly different between sows and piglets at T3 i.e. sows exhibited an up regulation in Class A & B betalactamases which catalyse the hydrolysis of the beta lactam ring using serine and zinc-based enzymes respectively [46]. On the other hand, multidrug e ux pump which are usually down-regulated by upstream pump operon transcription [1,47]. Such a change is likely to represent an increase in microbial multidrug resistance phenotype for beta-lactams among the sows but the opposite would be expected in piglets [48]. Furthermore, an over representation of genes that encode resistance speci cally those that modulate uoroquinolone resistance, is likely an indication of legacy of use on the farm. Finally, a gradual increase of genes in piglets suggests continuum of exposure via microbial succession hastened or otherwise by antibiotic use [35].
Secondly, viral sequences reveal a marked increase in the abundance of Siphoviridae, Phycodnaviridae, Podoviridae and Smascoviridae. Siphoviridae is associated with Bordetella bronchiseptica, the causative agent of atrophic rhinitis [49]. Indeed, we observe a marked increase (T3) of sequences of the family to which Bordetella species belong. (Fig S7). On the other hand, Phycodnaviridae primarily affects algae [50], suggesting a potential environmental incursion. Smascoviridae, speci cally of genus Porprismacovirus has previously been reported in pigs, however its potential to cause disease remains unknown. It is noteworthy that Adenoviridae sequences were also detectable in sows at T3 and these are known to cause mild gastroenteritis [51]. Crucially, the signi cantly higher abundance of Verrucomicrobiaceae in sows supports the notion of an infection as members in this family i.e. Akkermansia are associated with attenuation of gut in ammatory response [43].
Characterising AMR gene carrying taxa While it is widely accepted that the vast majority of AMR genes are carried by the uncultured fraction of a microbiome, here we have used 682 near complete genomes across major culturable taxa to show that majority of the AMR genes are carried by Firmicutes and Proteobacteria. This observation supports the notion of a direct relationship between microbiota and its resistome given that these phyla account for 58% of microbial population in this study. Furthermore, we show that a group of Proteobacteria i.e. Escherichia coli and Clostridium sp, were exclusively recovered from piglets and carried a variety of genes encoding resistance to classes of antibiotics i.e. the phenotype of MDR would be expected in this group of bacteria. This further underscores the impact of antibiotic use in this group.
Relevance to swine production Unravelling gut microbiome characteristics is an essential step in developing strategies that maximise antibiotic utility, nutrient extraction [52], mucosal immune-modulation [52] to improve welfare outcomes [37] for sustainable livestock production systems. To do so, we ought to use methodology that accurately capture these characteristics. Here, we show some differences in the utility of 16S rRNA gene metabarcoding and shotgun metagenomic sequencing and conclude that; a) the former could cost effectively be used to screening for patterns in longitudinal studies, b) and the latter for granular examination of what they represent.
While most studies have focused on bacterial shotgun sequences, here we show that viral sequences can be used to support hypothesis investigations. To this effect we used viral sequences to investigate a suspected infection among sows, and the ndings suggest that changes due to infections or otherwise could induce far larger microbiome changes than antibiotic use. It is therefore essential to account for such factors when analysing microbiome data.
Our original hypothesis here was that antibiotic use would restrict microbial diversity due to its bactericidal effect and this would result in an overrepresentation of genes encoding the corresponding resistance. The ndings show that antibiotic use in pigs is associated with changes in the microbial structure, abundance, and typically an increase in diversity. The changes in the resistome can take between 14-28 days to manifest and depending on the method of detection, remain detectable for a long time. All this highlights the integral relationship between microbial populations and the genes they carry, and the opportunities and challenges of exploiting them.

Conclusions
Gut microbiota and its resistome changed in response to antibiotic usage and these changes remain detectable for a long time. However, we ought to be prudent as similar change can be elicited by other factors, here we show that these are delineable using shotgun sequences including but not limited to the bacterial kingdom. We therefore recommend using 16S rRNA amplicon to screen for signatures and shotgun sequence data for their granular analysis.  a. Shows the study design, two groups i.e. piglets and sows longitudinally sampled at T1, T2, T3 and T4. After DNA extraction, sequencing and quality control, the reads are assembled to feed two analytical approaches. In green, the reads are assembled and ltered to retain >=800bp contigs, which are then annotated with AMR genes as well as taxonomic classi cation. In Blue, the reads are assembled into near complete genomes >75% and AMR genes, plasmids as well as taxonomic classi cation done. Fig 1b. Shows the hypothesis and assumption on which the analysis is done. In this study we have measured time represent as categorical increment (x), analysed the microbiome therefore changes between individuals over time (y) can be established and m is the rate of change in the microbiome characteristics. We assume that m reduces as piglets grow which is why sows in general will show limited variation (Fig 4b). We assume that total sample variation will follow the black line, antibiotic use and other routine activities can affect the m therefore shift the line to the grey and broken line positions, and thus accounting for the variation attributable to these factors. This assumption is validated in gure (Fig 4b)

Figure 2
Shows the taxonomic characteristic of the microbiome based on 16sRNA and shotgun metagenomic sequencing data. Panel A and C Shows the hierarchical clustering of family level abundance based on 16sRNA and shotgun sequencing data respectively. The latter's taxonomic assignment is based on contigs (>= 800bp). Panel B and D shows the taxonomic pro le at phyla level base on 16sRNA and shotgun sequencing data for pig groups across the four time points. Primarily, this gure demonstrates the difference in granularity of the two sequencing methods and the strength of the taxonomic signal based on how well it clusters samples group and time. Shows the resistome characteristics. Panel A shows the selected amplicon resistome, the copy-count of ve genes normalised by 16sRNA copy number. Panel B is an analysis of the same genes based on shotgun data measures as contigs per million reads carrying the selected AMR genes) based on qPCR and shotgun sequencing data. These two panels provide a comparison in characteristics as detected using both sequence methods. Panel C shows the full resistome i.e. AMR genes clustered by antibiotic class. The ll colour shows the number of contigs per million reads. CCRcross class resistance, MLS-macrolides Lincomycines and streptogramins, CAP-Cationic antimicrobial peptides. Others = Elfamycin, Rifamycins  Shows the taxonomic units are carrying the AMR genes and plasmids. The taxonomic units are Metagenomics assembled genomes (n=682) with genome completeness > 75% and contamination < 15%. Here the genomes are coloured by phyla (inner large ring). The larger two rings from the inner ring show the group and Time points. The smaller rings centripetally show hits on plasmids and antibiotic classes per genome colour intensity corresponds to abundance.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.