This study was licensed by Teagasc Animal Ethics Committee and was carried out on commercial pig farms in the Republic of Ireland. Ten farms were selected, five farms with routine use of in-feed prophylactic antibiotics and therapeutic ZnO (3,000 ppm) during first two weeks post-weaning (Treated, n = 5) and 5 farms that had not used either strategy during the previous three years (ZnO-free, n = 5). The farms used in this study ranged in size a size between 200 and 3000 sows, were all farrow-to-finish operations, and weaned piglets between 28 and 32 days of age. The antibiotics used on the Treated farms were amoxicillin (Stabox, Virbac, France) or trimethoprim and sulfadiazine in combination (Sulfoprim, Univet Limited, Ireland) at the dose indicated by the manufacturer. On each farm, two pens from two different rooms were sampled. Pen size varied depending on the farm and ranged between 12 and 72 pigs per pen with similar stocking densities, slightly above legal requirement, throughout. Environmental sampling was performed in empty clean pens immediately before the pigs were moved into the rooms on weaning day (0 day post-weaning, dpw) using sponge swabs (3M™ Sponge-Stick in sample bag with 10 mL D/E neutralizing broth, 3M Deutschland GmbH, Neuss, Germany). One swab was used to sample the feeders and the drinkers in the pen (FD) and another swab was used to sample two sections of 50 cm2 of the walls and the floor of the pen (WF). Pig fecal samples were collected at 0, 7 and 14dpw. In addition, a diarrheic fecal sample was collected at 7dpw if available. For the fecal sampling, one random freshly voided fecal sample from one pig per pen was collected and transferred to 1.5mL microcentrifuge tube. Each sample was collected using a sterile 140x7mm conical steel spatula avoiding the part in direct contact with the floor. Samples were transported to the laboratory under cooling conditions (in less than 2 hours) where swabs were processed extracting the sampled material from the swabs using 5 mL of sterile Phosphate Saline Buffer (PBS) 1X. Approximately 8 mL and 16 mL were recovered from each WF and FD swab, respectively, and transferred to a 20mL centrifuge tube. These tubes were centrifuged at 3000 x g for 15 min at 4° C, the supernatant was discarded, and the pellet was suspended in 1mL of PBS and transferred to a 1.5 mL tube. The processed environmental samples and the fecal samples were stored at -80° C until DNA extraction.

The DNA was extracted using the QIAamp PowerFecal Pro DNA Kit (Qiagen, Crawley, West Sussex, UK) following the manufacturer’s instructions, using 200 ± 50 mg of fecal content from samples. Environmental samples were previously thawed on ice, centrifuged at 15000 rpm for 1 min at 4° C, the supernatant was discarded and pellet was used for DNA extraction. A Qubit fluorometer (Qubit 3, BioSciences, Dublin, Ireland) was used to determine the total DNA concentration. The 2 samples from the different rooms for each type and time point were pooled by adding 5µL of each sample at a concentration of 1ng/µL. Paired-end sequencing libraries were prepared from the extracted DNA using the Illumina Nextera XT Library Preparation Kit (Illumina Inc., San Diego, CA) followed by sequencing on the Illumina NextSeq 500 platform using high-output chemistry (2 × 150 bp) according to the manufacturer’s instructions. Library size from each sample was assessed on an Agilent Technology 21000 Bioanalyzer using a High Sensitivity DNA chip.

Raw reads were filtered using trimmomatic v0.38 (Bolger et al., 2014). An average quality threshold score of 15 in a sliding window of 4 base pairs was used to trim reads below the threshold. A minimum length of 40 base pairs was ensured for all reads. Bowtie2 v2.4.4 (Langmead and Salzberg, 2012) was used to map the reads against host and human reference genomes, keeping the unmapped reads for the downstream analysis. Reference genomes were downloaded from Illumina iGenomes (https://support.illumina.com/sequencing/sequencing_software/igenome.html). Read duplicates were removed using the clumpify.sh tool in bbmap 38.22 (Bushnell, 2014). Analysis of microbial composition was carried out using Metaphlan v3.0 (Beghini et al., 2021). Functional profiles were assigned using HUMAnN v3.0 (Beghini et al., 2021). Gene families identified by HUMAnN were regrouped in MetaCyc metabolic reactions and Pfam protein domains functional categories using the utility script ‘humann_regroup_table’. MetaCyc is a database of metabolic pathways from all domains of life, containing pathways involved in primary and secondary metabolism, associated metabolites, reactions, enzymes, and genes (Caspi et al., 2020). Pfam is a comprehensive Database of Protein Domain Families (Finn et al., 2010). Abundance of MetaCyc metabolic reactions and Pfam protein domains was obtained from UniRef90 gene families using HUMAnN v3.0. Processed reads were assembled into contigs using the metaSPAdes pipeline from SPAdes v3.15.3 (Nurk et al., 2017). Mass screening of assembled contigs for E. coli virulence factors was performed with ABRicate (v1.0.1; https://github.com/tseemann/abricate), using the Ecoli_VF database (https://github.com/phac-nml/ecoli_vf).

All analyses were carried out in R v4.0.2 (R Core Team, 2022) with alpha level for significance of 0.05 and trend between 0.05 and 0.10 unless otherwise indicated. The fixed factors to be studied were the type of sample (FD, WF, feces or diarrhea), treatment or not with in-feed ZnO and antibiotics (Treated or ZnO-free) and day post weaning (0dpw, 7dpw or 14dpw). The type of sample and day post weaning were merged into a unique factor named “type-dpw”, having 6 different levels for this factor: feces 0dpw, feces 7dpw, feces 14dpw, diarrhea 7dpw, FD and WF. The farm was included in all the clustering analyses.

The effect of the treatment on the microbiota was studied between and within each type-dpw level. Alpha and beta diversities were computed at both the species and functional level using the R package Vegan v2.5-7 (Jari Oksanen et al., 2020). For alpha diversity estimation, Species richness, Inverse Simpson, and Shannon and Pielou evenness indices of diversity were calculated. Statistical differences in alpha diversity indexes were tested, after checking their normal distribution, with ANOVA and pairwise compared with Tukey (car v3.0.10) (Broom et al., 2006); multcompView v0.1.8, (Graves et al., 2019) and lsmeans v2.30.0 (Lenth, 2016) R packages or otherwise by using the Kruskal-Wallis test followed by pairwise Wilcoxon tests (R Core Team, 2022) R package. Beta diversity and ordination of samples were performed by non-metric multidimensional scaling (NMDS) of previously calculated Weighted Unifrac and Aitchison distances between samples by Species and functional abundance data, respectively. Weighted Unifrac distances on the species abundance table were calculated using the utility R script calculate_unifrac.R from Metaphlan. Aitchison distances were computed calculating the Euclidean distances of the CLR transformed pathways. Ordination of virulence factors data was performed by PCoA of distances calculated using Simple matching coefficient (Sokal and Michener, 1958). Separation between groups was tested with PERMANOVA (adonis2 and pairwise adonis) (Jari Oksanen et al., 2020; Martinez Arbizu, 2020). Factors and species influencing the ordination were assessed by linear models fitting on the ordination results (envfit function in Vegan R package). All P-values were adjusted by Benjamini-Hochberj (BH) procedure (Benjamini and Hochberg, 1995). For fitting species in ordination space, taxa and pathways were filtered, keeping the top 15 species and functions with the highest mean abundance across samples.

MetaCyc pathways obtained using HUMAnN were regrouped into MetaCyc superclasses using ‘humann2meco’ function from microeco package in R (Liu et al., 2021). Bacteria and function abundance analyses among type, dpw and treatment were performed using Linear Discriminant Analysis Effect Size (LDA LEfSE) (Segata et al., 2011). Data grouped in variables type_dpw, treatment, or type_dpw_treatment were used as classes selecting an alpha cut-off of 0.05 and a LDA threshold of 4 for type_dpw species composition analysis and type_dpw_treatment metacyc grouped superclass2 analysis, and 2 for species composition comparison between treatments. Species and functional genes explaining differences between classes were determined by LEfSE using Kruskal-Wallis test (P < 0.05) followed by linear discriminant analysis. Plots were built using ggplot2 v3.3.3 and pheatmap v1.0.12 in R (Wickham, 2016; Raivo, 2019) and the figures produced were subsequently arranged using Inkscape software v1.0.2 (Inkscape Project, 2020).

Results obtained from computing richness (Species richness) and diversity (Shannon, Inverse Simpson, Pielou) did not reveal any difference by the treatment (ZnO-free vs Treated) or farm (farms 1 to 10) variables, while there were differences when data was analysed by sample-type and days post-weaning (dpw) factors. Figure 1A summarises the α-diversity taxonomic results while functional analyses are shown in Figure 1B . Fecal diversity evenness measured by Pielou index decreased from weaning (0dpw) to the last sampling performed at 14dpw (P < 0.05). Interestingly, the lowest diversity value across the three indexes analysed was obtained in Diarrhea 7dpw samples, which was lower than Feces 0dpw in both Pielou evenness and Shannon indexes (P = 0.022 and P = 0.009, respectively). Diversity values in environmental samples collected at weaning revealed a large species richness in FD and to a lower extent in WF categories. While taxonomic α-diversity was not impacted by treatment, the analysis of diversity in metabolic pathways revealed differences between Treated and ZnO-free samples that were more evident in Diarrhea 7dpw samples. These functional differences were also significant for sample-type and days post-weaning variables ( Figure 1B ).

Similar to the results observed in α-diversity, neither treatment nor farm factors impacted the β-diversity ordination of the microbiota ( Supplementary Table S1 ) whereas both sample-type and days post-weaning variables contributed to microbiota spatial distribution ( Supplementary Table S2 , P < 0.05), at both species composition and metabolic pathways profiles, showing marked differences in diversity between fecal and environmental samples ( Figures 1C, D ). The ordination analyses revealed a clear separation by treatment group (Treated vs ZnO-free) in Diarrhea 7dpw samples, which were more apparent in functional than taxonomic data. The most abundant species and pathways that influenced the microbiota ordinations are shown in Figures 1C, D .

We further explored the sample metagenome structure for their similarity in microbial abundance using Ward clustering of species weighted Unifrac phylogenetic distances ( Figure 2 ; Supplementary Figure S1 ). The analysis clearly split samples in two branches, i.e. environmental samples (branch A) and fecal samples (branch B). The branch A composition was characterised by increased relative abundance of the Aerococcus and Corynebacterium genera, similar to the results observed in the ordination analyses ( Figure 1C ). This branch was split into two sub-branches; i.e. sub-branch A1, dominated by species of the genera Lactobacillus, Limosilactobacillus, Corynebacterium, Aerococcus and Staphylococcus, and sub-branch A2, with Aerococcus, Corynebacterium and Acinetobacter as the main representatives, with lower Lactobacillus abundance. Sub-branch A1 included six out of the ten samples from WF and sub-branch A2 included mainly FD samples (six out of ten). Most of samples in branch B were fecal samples and further sub-clustering in this branch was influenced by the species present rather than differences in genera and reflected the time-point factor ( Figure 2 ). Thus, the Feces 0dpw samples showed higher similarity with some Feces 7dpw samples (branch B2), and with Diarrhea 7dpw samples from ZnO-free group (branch B4). The remaining Feces 7dpw and Feces 14dpw samples were allocated to branch B3.

Figure 3A summarises the mean relative abundance of the most representative species in each sample type collected in the study. Environmental samples and Feces 0dpw showed similar or homogeneous abundance of different species while there was a hierarchy in microbiota relative abundance values in Feces 7dpw (including Diarrhea 7dpw) and Feces 14dpw samples. Focusing on fecal results, during the first 14 days after weaning, we observed an increase in abundance of the dominant species, shifts in species from the same genus and the bloom of anaerobes. Lactobacillus amylovorus and Limosilactobacillus reuteri, the two dominant species in fecal samples, increased in abundance (P < 0.05) during the post-weaning period ( Supplementary Table S3 ). Feces 0dpw showed an even species composition, with the exclusive presence of Prevotella sp. CAG:873 and Anaeromassilibacilus sp. An172 (both with relative abundance values over 2%) and higher abundance of Prevotella sp. CAG:520 ( Figure 3B ). Other species associated with 0dpw feces were Escherichia coli, Phascolarctobacterium succinatutens, Collinsella aerofaciens, Lactobacillus johnsonii and Phocaeicola vulgatus. In subsequent samplings, i.e., 7dpw and 14dpw, the abundance of Prevotella sp. CAG:873 decreased, while the abundance of two other species of Prevotella increased, Prevotella sp. P3-122 at 7dpw and Prevotella copri at 14dpw ( Supplementary Table S3 ). Two species exhibited higher abundance by linear discriminant analyses (LDA) at 7dpw; Butyricicoccus porcorum and the virus Lactobacillus phage phiAQ133. Megasphaera elsdenii increased in abundance towards the end of the of the study period with higher abundance at 14dpw ( Figure 3B ), sampling time at which we also observed higher abundance of Catenibacterium mitsuokai. While the relative abundance M. elsdenii in the piglets fecal microbiota was 3.08% at 0dpw, it had increased to 7.89% by 14dpw. In contrast, we observed the opposite trend for E. coli, from 11% at 0dpw to <2% at 14dpw. Diarrhea 7dpw samples showed a pattern which resembled those of 7dpw feces samples but with increased abundance of E. coli, and higher abundance of already dominant species such as L. amylovorus, L. reuteri and Ruminococcus torques. Both environmental samples showed similar compositions, with a high relative abundance of aerobic species such as Aerococcus spp. (A. viridans, A. urinaeequi) and Corynebacterium spp (C. stationis, C. xerosis, C. glutamicum).

Figure 3C shows the comparison of the taxonomic profiles at species level between the samples from the Treated and ZnO-free farms for each sample type. The LEfSe LDA analyses of differential abundance run globally revealed different species associated with treatment ( Figure 4A ). Overall, dominance of Lactobacillus spp. was not affected by the treatment and global analysis of species abundance by LDA confirmed these differences among treatments for M. elsdenii and P. succinatutens abundance. Analysis by sample type and dpw variables ( Figures 4B–D ) showed that P. vulgatus was the only species found to be more abundant in feces from farms using in-feed antibiotics and ZnO at weaning, i.e. before the treatment began ( Figure 4B ). We observed a notable change in dominance of environmental species by treatment variable. Species associated with treatment in the FD microbiota were A. viridans and A. urinaeequii and Enterococcus casseliflavus, whereas Sanguibacter sp. Leaf3 and Corynebacterium efficiens were associated with FD samples from ZnO-free farms ( Figure 4C ).

The species P. succinatutens, Roseburia inulinivorans, Ruminococcaceae bacterium D16 and Clostridium sp. CAG:242 were also higher in abundance in feces from Treated farms both at 7dpw and 14dpw ( Figures 4D, E ), whereas different species of Prevotella spp. were associated with the fecal microbiota in Treated farms at different time points. In contrast, species more abundant in ZnO-free farms were: Eubacterium rectale along with the virus Lactobacillus phage Lj711 in both fecal samples at 7dpw and 14dpw ( Figures 4D, E ); L. johnsonii and Prevotella sp. CAG:873 in both Feces 14dpw and Diarrhea 7dpw ( Figures 4E, F ); and Ligilactobacillus salivarius and Ligilactobacillus agilis in Feces 14dpw and Diarrhea 7dpw samples, respectively. Other species enriched in feces from ZnO-free farms at 14dpw were Eubacterium eligens, Clostridium sp. CAG:590, Streptococcus hyointestinalis, Acidaminococcus fermentans, Roseburia faecis and M. elsdenii.

Analysis of E. coli abundance revealed a trend and a high within-group variability ( Figure 4G ). Analysis of variance stabilizing AST transformed E. coli abundance revealed greater E. coli abundance in Diarrhea 7dpw samples from ZnO-free farms ( Figure 4G ). Further analyses of virulence factors linked to E. coli identified adhesins (fimbriae factors) and toxins associated to ETEC as shown in Figure 5 . Enterotoxigenic E. coli associated genes were present in all diarrhea samples, regardless the treatment group. Hierarchical cluster analysis based on the virulence factor profile revealed three different clusters. Cluster 1 included ZnO-free Diarrhea 7dpw and 0dpw fecal samples, and was characterised by the presence of F18 fimbrial, heat-stable (ST) enterotoxin, and EAST1 toxin genes. Cluster 3 defined by F4 (K88), heat-labile (LT) and ST toxins genes, consisted of Treated Diarrhea 7dpw samples ( Figure 5A ). The ordination of virulence factor profiles grouped the samples according to the observed clusters, with cluster 2 composed of environmental samples and some 7 and 14dpw fecal samples ( Figure 5B ). Analysis of abundance of sequenced assigned to Pfam protein families revealed higher abundance of heat stable enterotoxin ST and heat-stable E. coli enterotoxin 1 in ZnO-free diarrhea samples ( Supplementary Figure S2 ).

Analysis of the relative abundance of functional genes, grouped by superclass 2 level of MetaCyc metabolic pathways database, revealed the dominance of three functional superclasses, which accounted for more than 50% of the relative abundance of functional genes within the fecal microbiota: (i) nucleoside and nucleotide biosynthesis, (ii) amino acid biosynthesis, and iii) cofactor, carrier and vitamin biosynthesis ( Figure 6A ). The microbiota functional profiles, shown in Figure 6 , are ordered according to Ward clustering using Aitchison distance. Most of the samples analysed were clustered by their environmental or fecal origin (branches A and B, respectively). We observed that in Feces 0dpw or ZnO-free Diarrhea 7dpw samples ( Figure 6B , sub-branch B1), metabolic categories were ascribed to a lower number of species; whereas microbiota functions in Feces 7dpw and Feces 14dpw were built by a higher number of species ( Figure 6B , sub-branch B2). Species contribution to metabolic categories varied across clusters, particularly in ZnO-free Diarrhea 7dpw samples, but also in species contribution of less dominant categories ( Figure 6B ).

Supplementary Figure S3 shows which species significantly contributed to each superclass2 grouped metacyc pathway associated with the three factors under study, i.e., sample type, dpw and treatment. As already mentioned, the dominant functions were associated with different species to depending on the sample type, dpw and treatment. Notably, the Diarrhea 7dpw samples from ZnO-free farms were enriched in aerobic processes encoded by E. coli.