Community composition and development of the post-weaning piglet gut microbiome

We report on the largest metagenomic analysis of the pig gut microbiome to date. By processing over 800 faecal time-series samples from 126 piglets and 42 sows, we generated over 8Tbp of metagenomic shotgun sequence data. Here we describe the animal trial procedures, the generation of our metagenomic dataset and the analysis of the microbial community composition using a phylogenetic framework. We assess the effects of intramuscular antibiotic treatment and probiotic oral treatment on the diversity of gut microbial communities. We found differences between individual hosts such as breed, litter, and age, to be important contributors to variation in the community composition. Treatment effects of the antibiotic and probiotic treatments were found but were subtle, while host age was the dominant factor in shaping the gut microbiota of piglets after weaning. The post-weaning pig gut microbiome appears to follow a highly structured developmental program with characteristic post-weaning changes that can distinguish hosts that were born as little as two days apart in the second month of life.

while host age was the dominant factor in shaping the gut microbiota of piglets after weaning. The current study shows, by means of a phylogenetic diversity framework, that the post-weaning pig gut microbiome appears to follow a highly structured developmental program with characteristic postweaning changes that can distinguish hosts that were born as little as two days apart in the second month of life.

Background
As the world population grows, there is an accompanying demand for animal derived products. In a seminatural environment pig weaning occurs between the 12 th and the 17 th weeks from birth, whereas in the farm this typically occurs at 4 weeks of age 1 . Intensive animal husbandry and early weaning practices are commonly used to maximise production rates while minimizing costs. Both practices increase the risk of infections with pathogenic organisms, and thereby the need for antimicrobial strategies, which has included the common use of antibiotics. [2][3][4][5][6][7][8][9][10][11][12] Although antibiotics may kill some pathogens, the surviving bacteria can develop antimicrobial resistance (AMR) against the class of antibiotic used [13][14][15][16][17][18][19][20] , as well as PhyloSift 86 was employed as a means to study microbial community diversity among the samples, and to test for associations with treatment, time of sampling, and differences among hosts at the start of the trial. To this end, analysis of the unrooted phylogenetic diversity (PD) 87 , the balance weighted phylogenetic diversity (BWPD) 88 and principal component analysis (PCA) of the Kantorovich-Rubinstein distances 89 (beta diversity analysis) were performed.
Taxonomic analysis of the technical replicates of the probiotic ColiGuard® also showed a species pro le consistent with the expected pro le, with Lactobacillus salivarius and Lactobacillus plantarum in a 9:1 ratio (Lactobacillus salivarius: mean ± SD: 93.52 ± 1.617; Lactobacillus plantarum: mean ± SD: 6.10 ± 1.134) across the replicates (S8). ColiGuard® contained a total of 20 contaminants, of which 16 and 4 were identi ed at the species and the genus level, respectively. Contaminants were present at a higher level in two technical replicates (R5, R7), with R7 displaying the most diverse and highest contamination rate (R7: 14 taxa; total contaminating reads: 2.67%; R5: 9 taxa; total contaminating reads: 0.30%). (S12).
Phylogenetic diversity of positive controls and how it compares to the taxonomic pro le Alpha diversity of the positive controls re ects the expected alpha diversity, with the mock community, ColiGuard® and D-Scour™ positive controls, displaying a progressively higher unrooted PD (mean ± SD: Mock community: 31.53 ± 29.50; ColiGuard®: 58.52 ± 21.70; D-Scour™: 64.84 ± 21.30). BWPD for the positive control ColiGuard® appears on the far left of the plot (mean ± SD: 0.82 ± 0.15), distant from the mock community (mean ± SD: 1.58 ± 0.12) and the positive control D-Scour™ (mean ± SD: 1.89 ± 0.26), which contains 2, 7 and 8 main species, respectively (S13; Supplementary le 1). The unrooted PD of the positive control ColiGuard® (mean ± SD: 58.52 ± 21.70) is closer to the unrooted PD of the Mock Community (mean ± SD: 31.53 ± 29.50) and to that of the positive control D-Scour™ (mean ± SD: 64.84 ± 21.30), than the BWPD.
Based on Kruskal-Wallis one-way analysis of variance (Hommel adjusted p values to correct for multiple testing), phylogenetic diversity of the piglet samples did not cluster signi cantly by breed cross type (p > 0.05) or by pig line (p > 0.05) at the start of the trial, but a signi cant difference was found with breed and line in the fourth week of the trial (breed: p = 0.04; pig line: p = 0.02) (Fig. 2). Additionally a correlation of breed cross type was found with beta diversity principal components (PC) in the second week of the trial (PC3: p = 0.0357) and in the fth week of the trial (PC1: p = 0.027). Alpha diversity signi cantly correlated with pig line in the fourth week of the trial (BWPD: p = 0.017). (Fig. 2; Supplementary le 5) Notably, we found a signi cant correlation between alpha diversity and the age of the piglets at the rst sampling time point (unrooted PD: p = 0.0016; BWPD: p = 0.0355) ( Fig. 2; Supplementary le 5). While a correlation with age was found for alpha diversity only at the start of the trial, age of the piglet was observed to be weakly associated with differences in community composition in week two (PC3: p = 0.0507) (S14; Supplementary le 5).
As age groups were confounded with cross-breed types (i.e. not all age groups are represented by each of the four cross-breed types), we compared the phylogenetic diversity of age groups within each breed. As cross-breed types "Landrace × cross bred (LW × D)" and "Large White × Duroc" had only a small number of piglets in each age group, we tested for an association between phylogenetic diversity and age in crossbreeds "Duroc × Landrace" and "Duroc × Large White". In these cross-breeds, alpha and beta diversity are correlated with age in the "Duroc × Landrace" piglets (unrooted PD: p = 0.0059; BWPD: p = 0.0226; PC2: p = 0.0263; PC5: p = 0.0063) and to a lesser extent in the "Duroc × Large White" piglets (unrooted PD: p = 0.0529; PC5: p = 0.0310) during the rst week of the trial. In the "Duroc × Landrace" piglets, a correlation between age and beta diversity was detected in week two (PC2: p = 0.0347) ( A strong effect of aging on phylogenetic diversity Beta diversity analysis revealed a distinct and consistent change of the microbial community over time in all piglets, regardless of the treatment. Samples clustered in PC2 (accounting for 21.68% of the variation), showing a higher representation of Bacteroidetes, Gammaproteobacteria and Prevotellaceae from day 0-20, and a higher representation of Firmicutes, Mollicutes and Ruminococcaceae from day 20-40 ( Fig. 3; S20). In PC1 (accounting for 46.99% of the variation) samples shifted towards a higher representation of Lactobacilli from day 0 to day 20, and towards a higher representation of Actinobacteria, Clostridiales and Mollicutes from day 20 to day 40 ( Fig. 3; Supplementary le 5). A temporary shift towards a higher representation of Bi dobacteriaceae (Bi dobacterium being a component of the probiotic D-Scour™) is seen in PC5 (3.04% of the variation explained) one and two weeks from the start of the trial (S21).
Beta diversity analysis was performed separately for samples within each time point in order to nd taxa associated with variation within each time point. Extent of variation was derived from the product of branch width by the variation explained by the principal component (S22). The following taxa were responsible for variation only at the start of the trial: Gammaproteobacteria (t0 = 0.12), Enterobacteriaceae (t0 = 0.08) and Archaea (t0 = 0.03). The following taxa were responsible for variation throughout the trial: Clostridiales (min = 0.11; max = 0. In comparing four timepoints at one week intervals from the start of the trial (intervals labeled as A, B, C, D), changes in alpha diversity among all the piglets were tested for and signi cance was determined using the Bonferroni correction. Unrooted phylogenetic diversity increased after the start of the trial (

Effect of antibiotic and probiotic treatment on alpha diversity
We hypothesized that the probiotic treatments, whether alone (D-Scour™ and ColiGuard®) or administered after neomycin (neomycin + D-Scour™ and neomycin + ColiGuard®) would cause a change in the microbial community composition that would be measurable via phylogenetic diversity. We tested whether the treatments correlated with a change in phylogenetic diversity independently of the changes occurring with time. Given the differences in alpha and beta diversity detected among the subjects at the start of the trial, we analyzed the deltas of phylogenetic diversity instead of relying on the absolute means, similar to the procedure applied by Kembel et al (2012) 91 . Time-point measurements of alpha diversity were taken and deltas were computed for each piglet. Delta means were compared between cohorts, where the control cohort would serve as a control group for neomycin, D-Scour™ and ColiGuard® cohorts, whereas the neomycin cohort would serve as a control group for the neomycin + D-Scour™ and neomycin + ColiGuard® cohorts.
After the rst week of the trial, the majority of the piglets displayed an increase in unrooted PD (85%) and The only signi cant difference between the neomycin and the control cohort was detected in BWPD in the second week of the trial (control mean: 7.10; neomycin mean: 0.28) p = 0.033751) (S26 B-C).
No signi cant differences were detected between the control and the D-Scour™ cohort (p < 0.05) (Supplementary le 5). Instead, a signi cant decrease in BWPD was detected in the ColiGuard® cohort between the start of the trial and the end of the ColiGuard® treatment, where nearly all piglets (88%) in the ColiGuard® cohort displayed a decrease in BWPD, whereas the control cohort was split between piglets increasing (53%) and decreasing (47%) in BWPD (control mean: 0.85; ColiGuard® mean: -5.49; p = 0.040041) (S26 A-C).

Effects of antibiotic and probiotic treatment on beta diversity
To investigate the treatment effect on beta diversity, principal component analysis (PCA) of the Kantorovich-Rubenstein distances (beta diversity analysis) was performed on all samples and, additionally, on samples within individual time points. This analysis is conceptually similar to the weighted Unifrac approach for beta diversity analysis, but is designed to work with phylogenetic placement data 89

Effect of treatments on weight gain
Overall weight gain from initial to nal weight (S34 A-E; Supplementary le 5) was not signi cantly affected by any treatment. However, the probiotic ColiGuard® was found to have a partial effect on piglet weight gain (S34; Supplementary le 5). Weight was measured weekly for a total of six measurements. Based on Tukey adjusted p values, a lower weight gain was detected in the ColiGuard® cohort in week 3 compared to the control cohort (p = 0.0008) (S34 C-D; Supplementary le 5). Similarly, a lower weight gain was detected in the neomycin + ColiGuard® cohort in C-E compared to the neomycin cohort (p = 0.0393) (S34 C-E; Supplementary le 5). Breed and age differences among piglets were not associated with weight gain (Supplementary le 5).

Discussion
The consistent trend in community composition over time, across all the cohorts, indicates that an agerelated process of ecological succession is the largest factor shaping the microbial community of postweaning piglets, as found in this study where animals aged 20-63 days and were fed the same diet.  115 . In this study, piglets reached a comparable alpha diversity to the sows after the rst week of the trial, at which time the piglets were aged between 3.8 and 4.6 weeks. Unrooted PD did not reach higher levels at later sampling time points. The highest BWPD accompanied by a high unrooted PD was reached after the second week of the trial when piglets were aged between 4.9 and 5.6 weeks. Age-dependent physiological changes could explain i) the major shifts we detected in alpha diversity during the rst two weeks of the trial and, ii) the distinct differences in community composition with age, even with a narrow age difference between piglets (1-6 days). We were able to appreciate a signi cant trend of increasing unrooted PD and decreasing BWPD with age in piglets that separated a 6 days maximum by day of birth from each other. Since age groups were confounded with breeds in our study, we attempted to determine the correlation within single breeds. Unfortunately, although the correlations with age could still be detected, we could not determine the association at later time points due to the introduction of treatment effects.
Animal trials are often conducted in controlled environments so as to minimize environmental effects. However, individual variations such as breed and age are often unavoidable in large animal trials.
Previously reported confounding factors include: individual variation 44 44 , behavioural differences between breeds (e.g. coprophagy, mouth to mouth contact) and extent of long-term behavioural adaptation, which can differ between breeds for reasons not attributable to genetics 44,92,93 . A litter effect was found in piglets at the start of the trial and was lost at later time points during the trial. This could be due to either of the aforementioned factors. Cohousing, aging and the splitting of the piglets in separate rooms to receive a different treatment, are possible causes for loss of the litter effect. In this study we con rm the importance of these factors in the contribution to inter-individual variability of gut microbial composition. Motta et al (2019) report a correlation of beta diversity with age and no correlation of genotype and litter effect with either alpha or beta diversity 9 . On the contrary, we found the piglet samples to signi cantly cluster by litter, breed and by age up to the second and the fourth week post weaning, in alpha diversity and beta diversity, respectively. We conclude that even small age differences among post-weaning piglets, down to the day, must be accounted for in an experimental set up.
Three groups of piglets (cohorts neomycin, neomycin + D-Scour™ and neomycin + ColiGuard®) underwent 5 days of treatment with the broad-spectrum antibiotic neomycin, via intramuscular administration.
Intramuscular neomycin poorly diffuses (< 10%) into a healthy gastrointestinal tract 96 , therefore a direct effect of neomycin on the gut microbiome may not be expected. However, neomycin showed a different trend in unrooted PD between the second and the third week of the trial, corresponding to the week following the neomycin treatment period for the neomycin cohort. Taking this time frame into consideration, the neomycin cohort did not increase in BWPD to the extent of the Control cohort. Although statistically signi cant differences between neomycin and Control in alpha diversity were not reached, the BWPD of the neomycin cohort appears to follow a different trend to the Control from the rst week shifts of the Firmicutes/Bacteroidetes ratio following treatment (length of the treatment not reported) 98 .
The effects of intramuscular administration of neomycin (aminoglycoside class) on the gut microbiota have to our knowledge not been investigated. Based on our results we conclude that a mild effect on phylogenetic diversity is appreciable post IM neomycin treatment, up to two weeks after termination of the treatment. Additional compositional and functional analysis is necessary to determine the source of this mild variation. Differences were not detected at later time points, based on our phylogenetic diversity analyses, suggesting a full recovery of the microbial communities after two week from the end of the treatment.
It is possible that the large shifts in phylogenetic diversity taking place in the rst two weeks irrespective of the treatment (an increase, then decrease of unrooted PD, and an opposite trend of BWPD) have masked the milder effects of the treatment, despite our efforts to control for the effects of aging. This could be the reason why a signi cantly distinct alpha diversity trend was found in the neomycin + D-Scour™ cohort compared to the neomycin cohort, but not in the D-Scour™ cohort compared to the Control cohort. The neomycin + D-Scour™ cohort underwent 5 days of neomycin treatment followed by 2 weeks of D-Scour™ treatment. A signi cant increase of BWPD was detected in the two-week period of D-Scour™ treatment, indicating a possible enhancement of microbiome evenness following neomycin treatment. To our knowledge there are no studies reporting an increased evenness in piglet gut community composition following a speci c probiotic treatment. There are instead multiple studies reporting bene cial effects of probiotic treatment in sucker and weaner piglets in terms of improved gut mucosal integrity 66,80 , growth rate 80− 83 , digestibility of proteins and water absorption 80,83 , reduction of pathogen invasion e ciency 76,79,80 , and decreased mortality 80,82 . Although the assessment of physiologic changes from probiotic treatments was outside the scope of this study, we found signi cant separation of neomycin + D-Scour™ cohort samples to neomycin cohort samples in beta diversity 3 and 10 days after D-Scour™ treatment, where neomycin + D-Scour™ samples showed a higher representation of Lactobacillales compared to neomycin samples, suggesting a transient establishment of the probiotic strains in the piglet guts.
The second probiotic in this study, ColiGuard®, did not have an effect on alpha diversity, but clustering was detected in beta diversity, where ColiGuard® samples separated from Control cohort samples in the third principal component (explaining 7.88% of the variation) two weeks post probiotic treatment. Additionally, the ColiGuard® treatment correlated with a lower weight gain, whether or not it was preceded by the antibiotic treatment. However, when comparing the overall weight gain (from the start to the end of the trial) the weight gain in the cohorts receiving ColiGuard® did not differ from the other cohorts.
We extracted the 16S rRNA gene hypervariable regions from our dataset, obtained the counts, and ran a correlation analysis to discover taxa that correlated with the weight of the piglets. As a consequence of the library size normalization step, the use of correlation with compositional data can in ate the false discovery rate 99,100 . For this reason it can be expected that some of the taxa we found to correlate with the weight of the piglets (eighty-three distinct species) could be spurious while other correlations may have been missed.

Technical controls in metagenomic studies and methodological limitations
Taxonomic assignment of the raw reads from the positive controls was performed with MetaPhlAn2 101 which relies on a ca. 1M unique clade-speci c markers derived from 17,000 reference genomes. Such a database to map against the positive controls su ces as these organisms are cultivable, and for this reason they are widely studied hence the sequences are known. This is not the case for real-world samples where mapping against a database (which completeness relies on studied and often cultivable organisms) would narrow the view on the true diversity within the sample.
Positive controls with well-studied members and known ratios within the samples, has proven to be a valuable tool to assess consistency among technical replicates across batches and to detect possible biases derived from the DNA extraction method.
Systematic taxonomic bias in microbiome studies, resulting from differences in cell wall structures between Gram positive and Gram negative bacteria, have previously been reported; sample treatment with enzymatic cocktails can modestly reduce this bias [102][103][104] . Although we implemented this step in our work ow, it seems that, from the read abundance of our mock community, which contained three Gram negative and four Gram positive strains, a bias towards Gram negative taxa may still be present.
In terms of contamination we concluded that: a) contamination in our study was not batch speci c; b) a problem of sample cross-contamination existed at the DNA extraction step between neighbouring wells.
During the bead beating step of DNA extraction, the deep-well plate is sealed with a sealing mat, rotated and placed in a plate shaker for the bead beating to take place. We consider that sample crosscontamination is most likely to occur during this step.

Conclusion
Our study provides a publicly available databank of the pig gut metagenome. Our ndings further stress the importance of confounding factors such as breed, age and maternal effects when assessing the effect of treatment on the gut microbiome. We found that age, even within a narrow age span (1-6 days) can have an impact on microbial shifts and should be accounted for in microbiome studies. Intramuscular neomycin treatment correlated with a clustering in alpha diversity and a higher representation of Mollicutes compared to control. D-Scour™ treated piglets displayed a mild shift in alpha diversity compared to control, and a transient establishment of Lactobacillales. ColiGuard® treated piglets displayed a clustering in beta diversity and a transient lower weight gain compared to control.
Weight correlated with the abundance of a number of taxa. Age was the strongest factor shaping phylogenetic diversity of the piglets.
As previously mentioned, phylogenetic diversity is based on distinct taxa (richness) and their collective structure (proportions re ected by BWPD) and not on a direct assessment of composition and function.
These types of analyses will be necessary to further describe the effects of the treatments.

Pig trial
Animal studies were conducted at the Elizabeth Macarthur Agricultural Institute (EMAI) NSW, Australia and were approved by the EMAI Ethics Committee (Approval M16/04). The trial animals comprised 4- week old male weaner pigs (n = 126) derived from a commercial swine farm and transferred to the study facility in January 2017. These were cross-bred animals of Landrace, Duroc and Large White breeds and had been weaned at approximately 3 weeks of age.
The pig facility consisted of four environmentally controlled rooms (Rooms 1-4) with air conditioning, concrete slatted block ooring with underground drainage and open rung steel pens (S1). Each room had nine pens, consisting of a set of six and a set of three pens, designated a-f and g-i respectively, with the two sets of pens being physically separate, i.e. animals could come in contact with each other through the pen's bars within each set of pens, but not between sets. The rooms were physically separated by concrete walls and contamination between rooms was minimized by using separate equipment (boots, gloves, coveralls) for each room. In addition, under-oor drainage was ushed twice weekly and the ushed faeces/urine was retained in under-oor channels that ran the length of the facility, so that Rooms 1, 2 were separate from Rooms 3, 4 and ushing was in the direction 1 to 2 and 3 to 4.
The pigs were fed ad libitum a commercial pig grower mix of 17.95% protein free of antibiotics, via selffeeders. On the day of arrival (day 1) 30, 18, 18, and 60 pigs were allocated randomly to Rooms 1, 2, 3 and 4 respectively in groups of 6, 6, 6 and 6-7 pigs per pen respectively (S1A). Pigs were initially weighed on day 2, and some pigs were moved between pens to achieve an initial mean pig weight per treatment of approximately 6.5 kg (range: 6.48-6.70; mean ± SD: 6.53 ± 0.08). Pigs were weighed weekly throughout the trial (Supplementary le 1). Behaviour and faecal consistency scores were taken daily over the 6-week period of the trial (Supplementary le 2). Developmental and commercial probiotic paste preparations ColiGuard® and D-Scour™ from International Animal Health, Australasia, were used in some treatment groups.
The animals were acclimatised for 2 days before the following treatments were administered: Room 1oral 1 g/pig of placebo paste daily for 14 d; Room 2 -oral 1 g/pig of D-Scour™ paste daily for 14 d; Room 3 -oral 1 g/pig of ColiGuard® paste daily for 14 d; Room 4 -intramuscular (IM) injection of antibiotic administered at 0.1 mL per pig daily from a 200 mg/mL solution for a total treatment duration of 5 d.
On the day following the nal neomycin treatment (day 8), 36 pigs were moved from Room 4 to Room 2 (n = 18, 6 in each pen, pens g-i), and to Room 3 (n = 18, 6 in each pen, pens g-i) (S1). The following day (day 9), oral administration of D-Scour™ (1 g/pig) and of ColiGuard® (1 g/pig) commenced for pigs in Room 2 pens g-i and in Room 3 pens g-i, respectively, and continued for a period of 14 days. Assignment of the 36 neomycin-treated pigs to the treatment groups neomycin + D-Scour™ (n = 18; Room 2 pens g-i) and neomycin + ColiGuard® (n = 18; Room 3 pens g-i), was carried out by distributing them so that the mean weight of the animals distributed across pens and rooms was similar (Supplementary le 1). By this time point, each occupied pen in the trial housed six pigs. (S1B) From that time, twelve piglets from the original 126 were no longer present, as they had been euthanised as pre-treatment controls at the start of the trial.
Faecal samples were collected from all piglets once per week and from a subset (n = 48 pigs; 8 from each of the six cohorts) twice per week over the 6-week study period (Fig. 1).

Library preparation
Sample index barcode design using a previously introduced method 105 yielded a set of 96 × 8nt sequences with a 0.5 mean GC content and none of the barcodes containing 3 or more identical bases in a row. Nine hundred-sixty different combinations of i5 and i7 primers were used to create a uniquely barcoded library for each sample. The detailed sample-to-barcode assignment is given in Supplementary le 3. Library preparation was carried out using a modi cation of the Nextera Flex protocol to produce low bias, low cost shotgun libraries, as described in a previous manuscript 105 . Following the ampli cation step, samples were centrifuged at 280 x g for 1 min and stored between 1 and 5 days at 4 °C.

Size selection and puri cation
Samples from the same 96-well plates were pooled into one tube by taking 5 µL from each library. This generated 10 pooled samples, one for each plate. A master pool was created by pooling 5 µL from the pool of each plate into a single pool. Forty microliters from each of the 10 plate pools and 40 µL from the master pool underwent library size selection and puri cation using equal volumes of SPRIselect beads (Beckman Coulter) and ultrapure water (Invitrogen). Sample cleaning with SPRI-beads was performed as described previously 105 . A puri ed master pool comprising samples from all plates, and puri ed pools of individual plates to check for plate-speci c anomalies, were diluted to 4 nM and fragment size distribution was assessed using the High Sensitivity DNA kit on the Bioanalyzer (Agilent Technologies, USA).

Normalisation and sequencing
The master pool was sequenced on an Illumina MiSeq v2 300 cycle nano ow cell (Illumina, USA). Read counts were obtained and used to normalise libraries. The liquid handling robot OT-One (Opentrons) was programmed to re-pool libraries based on read counts obtained from the previous MiSeq run. The script used to achieve the normalization is available through our Github repository.
The read count distribution after normalisation is displayed in supplementary gure (S3
Additionally, SortmeRNA 126 (version 2.1) was applied to extract 16S ribosomal RNA genes from raw reads. The rRNA reference database silva-bac-16 s-id90 was used. Hits were ltered (e-value < = 1e- 30) and PCA was computed with R 108 . Sample counts were normalized for library size by proportions and were tested with the Spearman's Rank correlation coe cient method to nd taxa correlating with the weight of the piglets across the trial.

Batch effects
A randomized block design was adopted to mitigate batch effects. Because samples were distributed across ten 96-well plates during DNA extraction and library preparation, plate effects were expected.
Although samples did not visibly cluster by DNA extraction plate across the rst ve principal components (S5), a batch effect was found by multiple comparison analysis with ANOVA and by applying Tukey post-hoc correction to pairwise comparisons. Batch effects were detected (ANOVA, alpha