Differential Analysis of Longitudinal Methicillin-Resistant Staphylococcus aureus Colonization in Relation to Microbial Shifts in the Nasal Microbiome of Neonatal Piglets

ABSTRACT Methicillin-resistant Staphylococcus aureus (MRSA) is an important human pathogen and often colonizes pigs. To lower the risk of MRSA transmission to humans, a reduction of MRSA prevalence and/or load in pig farms is needed. The nasal microbiome contains commensal species that may protect against MRSA colonization and may be used to develop competitive exclusion strategies. To obtain a comprehensive understanding of the species that compete with MRSA in the developing porcine nasal microbiome, and the moment of MRSA colonization, we analyzed nasal swabs from piglets in two litters. The swabs were taken longitudinally, starting directly after birth until 6 weeks. Both 16S rRNA and tuf gene sequencing data with different phylogenetic resolutions and complementary culture-based and quantitative real-time PCR (qPCR)-based MRSA quantification data were collected. We employed a compositionally aware bioinformatics approach (CoDaSeq + rmcorr) for analysis of longitudinal measurements of the nasal microbiota. The richness and diversity in the developing nasal microbiota increased over time, albeit with a reduction of Firmicutes and Actinobacteria, and an increase of Proteobacteria. Coabundant groups (CAGs) of species showing strong positive and negative correlation with colonization of MRSA and S. aureus were identified. Combining 16S rRNA and tuf gene sequencing provided greater Staphylococcus species resolution, which is necessary to inform strategies with potential protective effects against MRSA colonization in pigs. IMPORTANCE The large reservoir of methicillin-resistant Staphylococcus aureus (MRSA) in pig farms imposes a significant zoonotic risk. An effective strategy to reduce MRSA colonization in pig farms is competitive exclusion whereby MRSA colonization can be reduced by the action of competing bacterial species. We complemented 16S rRNA gene sequencing with Staphylococcus-specific tuf gene sequencing to identify species anticorrelating with MRSA colonization. This approach allowed us to elucidate microbiome dynamics and identify species that are negatively and positively associated with MRSA, potentially suggesting a route for its competitive exclusion.

during initial time points. Piglets from litter A were more often positive for S. aureus than piglets from litter B (24/52 versus 15/52 samples). The mean number of S. aureus in positive samples was 4.60 Â 10 4 and 3.0 Â 10 4 log CFU-equivalents (CFUeq)/swab for litter A and litter B, respectively. Notably, 31 of the 39 S. aureus-positive nasal samples were also MRSA-positive. MRSA was detected in 20 S. aureus-negative samples (Fig. S1). The mean number of MRSA CFU in positive samples was 138 and 19 CFU/swab in piglets from litter A and litter B, respectively.
Microbiome analysis was carried out on 7.07 and 7.34 million error-corrected, nonchimeric amplicon sequence variant (ASV) reads with a mean count of 80,282 6 26,624 standard deviation (SD) and 71,272 6 28,735 SD reads for the 16S and tuf data sets, respectively. In addition, 4 negative-control samples were sequenced, but considerably fewer error-corrected reads were generated, with an average of 401 and 2,351 reads for the 16S and tuf data sets, respectively ( Figure S2a). Overall, 2,787 unique 16S ASVs were identified, 23 of which were detected as potential contaminants. For tuf, 39 of the 1,278 ASVs were identified as potential contaminants and were removed. However, this did not have any effect on overall sequencing depth ( Figure S2b). Samples with less than 5,000 (n = 8) reads in the 16S and 13,000 (n = 1) reads in the tuf data sets were excluded from further analysis. The numbers of reads retained after abundance-based filtering, i.e., after excluding ASVs present in less than 10% of samples with less than 0.001% abundance, were 95.24 6 8.75% and 97.17 6 6.45% for 16S and tuf sequencing, respectively ( Figure S2b). Thus, the contribution of the excluded ASVs to the overall number of reads per sample was found to be very small or negligible.
tuf gene sequencing provided species-level resolution of Staphylococcus taxa. It has previously been noted that the V4 hypervariable region of the 16S rRNA gene alone is not sufficiently discriminative for the identification of species within the Staphylococcus genus (24,25). Therefore, to increase the resolution of the Staphylococcus genus, which was one of the aims of the current study, we also carried out amplicon sequencing of the tuf gene, which better discriminates between different species of Staphylococcus, Streptococcus, and Enterococcus (26,27). In particular, tuf gene sequencing led to the identification of 22 different Staphylococcus species, while 16S rRNA gene sequencing only detected the Staphylococcus sciuri species (Fig. 1A and B). When examining the data at the even more granular ASV level, we identified 137 and 10 different sequence variants belonging to Staphylococcus taxa from the tuf and 16S data sets, respectively (Fig. S5). As expected, the abundance of ASVs assigned to the Staphylococcus genus in samples sequenced for both 16S rRNA and the tuf gene correlated significantly (rmcorr coefficient [r rm ] = 0.75, confidence interval [CI] = 0.68 to 0.86, P = 1 Â e 217 ).
Longitudinal development of piglet nasal microbiota. The composition of the piglet nasal microbiota based on 16S rRNA gene sequencing shows clear segregation of the samples based on time points (Fig. 1C). The piglets exhibited a gradual shift in collective microbiome composition from day 0 to day 42, with permutational multivariate analysis of variance (PERMANOVA) analysis showing a significant time point-associated variance on microbiota structure (R 2 = 0.55, P = 1 Â e 24 ). In alignment with the 16S rRNA gene data, significant changes in community membership and a clear structural shift from day 0 to day 42 (PERMANOVA: R 2 = 0.49, P = 1 Â e 24 ) were also apparent from tuf analysis (Fig. 1D). Litter membership explained significant but less of the variation of the microbiome community structure (PERMANOVA: 16S R 2 = 0.03, P = 1 Â e 24 ; tuf R 2 = 0.04, P = 1 Â e 24 ). High interindividual differences and separate clustering observed for 0-h and 8-h samples might be related to the effect of unstable microbiota derived from fecal or soil contamination in newborn piglets ( Fig. S4A and S6). Interestingly, we observed that as the piglets age, Proteobacteria (i.e., Moraxellaceae in tuf) increase, while taxa belonging to Firmicutes (i.e., Staphylococcaceae in tuf) decrease ( Fig. 1C and D).
We observed higher microbiota alpha diversity in 0-h and 8-h samples in both 16S and tuf data sets, but these levels then dropped dramatically at 16 h, after which the diversity gradually increased for 16S data (Fig. 1E), while it successively decreased for tuf data after a peak at day 7 (Fig. 1F). There was a high presence of the genus Clostridium, a strict anaerobe often described in the gut, at 0 h and 8 h (Fig. S6). The higher microbial richness and diversity observed in the 0-h and 8-h samples might be related to bacteria introduced from the birth canal or fecal or soil contaminants in the newborn piglets.
These observations were used as evidence to exclude the first two time points from the anticorrelative analysis against S. aureus. Next, we investigated dynamic changes in alpha diversity over time, and we found a significant increase in chao1 (richness; 16S: r rm = 0.67, CI = 0.50 to 0.80; tuf: r rm = 20.62, CI = 20.73 to 20.45) and Shannon index (evenness; 16S: r rm = 0.50, CI = 0.28 to 0.66; tuf: r rm = 20.51, CI = 20.62 to 20.38) with time after exclusion of 0-h and 8-h samples. The negative alpha diversity trend observed in the tuf data set may be explained by the reduced abundance of the taxa-rich Firmicutes phylum in nasal microbiota of growing piglets at later time points ( Fig. 1F and Fig. S3B).
The age-based dynamic changes of the microbiome compositions were further evaluated at a lower taxonomic level. Using 16S data, we found there were 22 genera markedly changed among the top 50 abundant ASVs. The relative abundance of ASVs from the genus Rothia (and Rothia nasimurium at species level) increased from 16 h to 7 days but subsequently decreased after 14 days. Decreases in abundance of Rothia was accompanied by increases of Moraxella and Streptococcus genera ( Fig. 2A). In the tuf data, Staphylococcus accounted for more than 25% of the bacterial sequences until the age of 1 day but decreased dramatically from day 2 to day 14, which agrees with the 16S data. Of the 22 identified Staphylococcus species, S. microti (6.4%), S. haemolyticus (3.2%) and S. hyicus (3.2%) were the most abundant, while S. hominis, S. simulans, S. cohnii, S. arlettae, S. epidermidis, and S. aureus each accounted for approximately 0.1% of total bacterial abundance (Fig. 2B, bottom annotation;).
The 16S rRNA and tuf gene sequencing data provided species-level resolution for some, but not all, of the ASVs. In total, 28 different species were significantly correlated with colonization of MRSA and S. aureus ( Fig. 2A; Table S1). Consistent with what we obtained at the genus level, rmcorr analysis demonstrated species such as Streptococcus agalactiae, Acinetobacter schindleri, Mannheimia varigena, Helcococcus ovis, Corynebacterium stationis, and Rothia nasimurium (r rm . 20.55; adj. P value , 0.05) to be strongly anticorrelated with MRSA colonization in 16S data ( Fig. 3A and Fig. S7). Similarly, C. stationis and M. varigena were found to be anticorrelated with S. aureus colonization (r rm . 0.55; adj. P value , 0.05). Although a low level of C. stationis was also observed in MRSA-positive samples, we exclusively observed increased C. stationis abundance in MRSA-negative samples (Fig. S7). Of note, there was a weak but insignificant correlation of the genus Corynebacterium with carriage of MRSA (r rm = 20.31; adj. P value . 0.05) and S. aureus (r rm = 20.32; adj. P value . 0.05).
To confirm the detected correlation-based associations, we performed logistic regression analysis to correlate MRSA/S. aureus colonization with the genus and species level microbiota. As expected, taxa identified as most significantly associated with MRSA/S. aureus colonization using correlation-based associations were further validated with the regression-based analysis. Genera and species found to be significantly associated with MRSA/S. aureus colonization using 16S-and tuf-based data sets are provided in Table S2.
Microbial taxa anticorrelated with MRSA/S. aureus nasal colonization tend to cooccur. As the nasal cavity is a nutrient-limited environment, the composition of nasal microbiota can be modulated by interactions between different bacterial species. Intermicrobial interactions can be a major driver of microbial community composition, and understanding such interactions can unveil important insights regarding establishment and carriage of MRSA/S. aureus in the nasal environment. Thus, we further investigated if microbial taxa identified as negatively associated with MRSA/S. aureus nasal colonization display a tendency toward cooccurrence or not. Using 16S species-level data, we identified nine Coabundant groups (CAGs), each comprising bacteria significantly correlated with each other from 16 h to day 42 (Fig. 4A). Of these, CAG 6, CAG 7, CAG 8, and CAG 9 were composed of taxa which were positively correlated (MRSA/S. aureus-positive CAGs), while CAG 1, CAG 2, CAG 3, CAG 4, and CAG 5 were composed of bacterial taxa which were negatively correlated with MRSA/S. aureus colonization (MRSA/S. aureus-negative CAGs). The constituent taxa of the CAGs not only cooccurred in terms of overall abundances, but also varied consistently over time. In particular, CAG 1/CAG 3/CAG 4/CAG 5 were anticorrelated with CAG 6/CAG 7/CAG 8 (Fig. 4B). We noted potential driver-passenger dynamics in CAG 6 whereby Moraxella spp. (marked *) is the first taxa to increase in abundance over time and is then followed by the other CAG 6 species. Using this 16S amplicon, only "unclassified Staphylococcus" in CAG 3 (marked §) was classified for this genus, but since it was negatively correlated with MRSA/S. aureus levels, we believe it is of a different species than S. aureus.
A total of three CAGs in tuf species-level data were identified, where CAG 1 and CAG 2 (containing S. aureus marked *) were positively correlated. CAG 3, on the other hand, was negatively correlated with MRSA/S. aureus nasal colonization (Fig. 5A). This CAG was the largest cluster containing taxa such as S. microti, S. simulans, E. faecium, and Streptococcus spp., and we observed no obvious driver-passenger dynamics in this potential MRSA-excluding group. Interestingly, Staphylococcus cohnii was correlated with S. aureus in CAG 2, while S. epidermidis and S. hominis were part of a separate CAG, suggesting differences in their abundance dynamics over the course of time (Fig. 5B).

DISCUSSION
Species anticorrelated to S. aureus and MRSA were identified. Detection of S. aureus in marker gene analysis can be hampered by the fact that the piglet nostrils harbor relatively small amounts of S. aureus. The low abundance of S. aureus observed in this study is in concordance with previous findings (20,21). The genus-and specieslevel resolution obtained from the sequencing data, substantiated with S. aureus-specific qPCRs and MRSA-specific culturing allowed successful identification of genera and species, which anticorrelated with S. aureus and MRSA. Interestingly, anticorrelating OTUs of Helcococcus and Acinetobacter have been described before in relation with low numbers of MRSA in pig noses (20), but there was no match with phyla or genera anticorrelated with MRSA in the study from Weese et al. (21). A limitation of our study is the small number of pigs that were analyzed, and it is possible that identified phyla or genera are not completely correlated with findings in other studies studying the pig microbiome. tuf gene sequencing identified Staphylococcus microti and Staphylococcus simulans as negatively associated with S. aureus. Such an inhibiting effect was recently described by Brown and colleagues, showing that peptides of S. simulans protected against MRSA colonization and associated skin damage in a mouse model (28). These peptides were inhibiting or disrupting of the arg-based quorum sensing of S. aureus that has been associated with colonization and virulence factor activation. Production of arg quorum sensing inhibiting peptides has been detected in multiple coagulase-negative staphylococci (CoNS), including S. simulans, from porcine nasal swabs (29). It is considered an important mechanism for bacterial interactions evoking S. aureus competition. Other competition mechanisms involved in nasal colonization, apart from the production of small molecules, include competition for adhesions sites and nutrients, antibiosis, and inducing host defenses (30).
tuf gene sequencing improves Staphylococcus species resolution. To understand the composition of the nasal microbiome and its interactions, high taxonomic resolution at the species or even strain level is needed, as identifying anticorrelating genera to MRSA could lead to misinterpretations. For example, Yan et al. showed that two species of the genus Corynebacterium, the species C. accollens and C. pseudodiphtericum might act differently on S. aureus colonization in the nasal cavity (31). They identified that these species showed either inhibition or stimulation of S. aureus growth in vitro. Therefore, in-depth analysis of individual bacterial species to find S. aureus anticorrelative species is crucial. A limitation in our study is the reliance on the V4 region of the 16S rRNA gene, as this sequence region contains low sequence diversity and is unable to discriminate S. aureus from other Staphylococcus species in microbiome analysis (26,27). However, several studies applying tuf gene sequencing have shown that this gene is discriminating of all Staphylococcus species (23,26,27) but can also monitor shifts in abundance of clinically important Staphylococcus species in the nasal microbiome (26,32,33). Using tuf gene sequencing, we identified 22 Staphylococcus species. This is in contrast to previous work where 12 Staphylococcus species were identified, with S. equorum as the most abundant in the porcine nose (32). In our study, S. microti was the most abundant staphylococcal species and was predominantly present in the first week of life. Moreover, we found that it was negatively associated with S. aureus, and its abundance decreases after day 4, when stable S. aureus colonization was established. Our species-level identification highlights the added value of complementing 16S rRNA sequencing with tuf gene sequencing, or multiple 16S rRNA gene regions (34), in microbiome studies, especially when Staphylococcus species are a target. To achieve even higher resolution and also functional information, metagenomic shotgun sequencing would be required.
Trends in the developing nasal microbiome. Here, we captured the dynamic and longitudinal development of the nasal microbiota of piglets. The identification of CAGs of bacteria also demonstrated time-dependent trends, further supporting that the porcine nasal microbiota is not stable but develops throughout time with a succession of coabundant species. The finding that the Proteobacteria and Firmicutes were the most abundant phyla agrees well with previous nasal microbiota studies of pigs (21,32,33,(35)(36)(37). A large drop observed in the relative abundance of Actinobacteria after day 7 has also been described (36). Moreover, the rise in abundance of Proteobacteria after weaning relates to Moraxella becoming the most abundant genus, and this finding is in line with the increase of Moraxella and Bergeyella upon removal of perinatal antimicrobials (35). But it is important to note that neither the piglets nor the sows received any antimicrobials for our study. Additionally, R. nasimurium from the phylum Actinobacteria has been previously described as a commensal on porcine tonsils and capable of producing the antibiotic valinomycin (38). The onset of the R. nasimurium decrease was around the time that S. aureus was detected and coincided with the decrease of the taxa from tuf-CAG 3, 16S-CAG 3, 16S-CAG 4, and 16S-CAG 5, consisting of additional anticorrelating species to S. aureus and MRSA colonization. Moreover, the taxa of 16S CAGs 6 and 7 were positively correlated with MRSA and contained the genera Oscillospira, Dorea, Peptococcus, Lactobacillus, Coprococcus, and Methanobrevibacter. This hints at a microbial shift associated with a loss of a protective effect against, or a stable colonization of, S. aureus around this time point. The question remains whether this shift is universal or an effect of host or environmental stimuli. No farm-related effects could be studied here, as the piglets were obtained from a single farm. Other environmental effects that might explain the microbial shift could be fecal input, as evident by an increase in the genus Clostridium around day 14, and other gut-related genera from 16S-CAG 6, 7, and 8. A decrease of maternal immunity after the first week of life, dietary changes approaching weaning, or applied perinatal antimicrobials are other factors that can modulate microbial shifts in the microbiome (39). This indicates that phyla and genera negatively associated with MRSA identified in silico will require further investigation with regard to their interactions with MRSA and their ecological context in the microbiome of the host.
In the human gut, the importance of an initial priming effect of natural birth on the further development of the microbiome and host immunology has been well described (40,41). As we showed that the microbiome is shaped by development of the piglets, we expect that manipulations of the microbiota in early life could later in life stabilize in the microbiome. It is important that these manipulations will not result in dysbiosis and enable colonization of pathogenic bacteria. This underlines why longitudinal investigation of a priming effect and the developing and stabilizing community in the nasal microbiome is essential. Microbes from the maternal gut, birth canal, and skin, are the first to colonize the naive nose epithelium of the newborn piglets. Some of the species present at initial time points were found throughout the study, indicating that development of the microbiome started directly at birth and stabilized over time. Detection of Archaea and anaerobic bacterial species at the later time points might indicate continuous introduction of fecal species into the nostrils of piglets. This could be a result of the rooting behavior of piglets. However, recent studies have described a large archaeal diversity in the human nose (42), and Archaea might be a stable constituent of the porcine nasal microbiome. As the number of longitudinal pig microbiome studies from birth is extremely low, more research is needed to understand the drivers of the development of the porcine nasal microbiome. Our study observed a potential early-in-life protective delay of MRSA colonization. We identified CAGs of species negatively associated with MRSA. Members of these CAGs were present at all time points. This could indicate that these species remain colonized and could establish a lower or negative MRSA presence later in life. Therefore, it is important to investigate the species negatively associated with MRSA in a larger number and more diverse set of animals and to obtain data from pigs that present a long-term stable MRSA-negative status.
Conclusion. Combining 16S rRNA and tuf marker gene sequencing with culture and qPCR-based quantification led to the identification of bacterial species negatively associated with MRSA and S. aureus in the pig nasal microbiome. The nasal microbiome developed with a time-dependent succession of coabundance groups that may indicate earlyin-life protection of S. aureus or MRSA colonization. Supplementing this study with nextgeneration sequencing free of amplification bias, such as shotgun metagenomics, will potentially lead to a higher taxonomic resolution and functional insights. The higher resolution is needed to study interactions at the strain level, enabling a better understanding of the complexities of the developing nasal microbiome, which could lead to novel strategies to reduce colonization of pathogens.

MATERIALS AND METHODS
Animal management and sampling. The study was performed in accordance with the Dutch law on research animal welfare and was approved and registered under 2014.II.05.036 by the Animal Ethical Committee of Utrecht University, the Netherlands. The study was carried out on a conventional farm where two random sows from different pens were selected. Eight landrace piglets from two litters (litter A and litter B) were sampled at 13 different time points. Piglets received colostrum and had access to solid feed ad libitum. Animals received an iron injection (200 mg per animal) at the age of 1 week as a part of normal pig-farming procedure to supplement iron deficiency. Vaccinations against mycoplasma and circovirus were performed at the age of 4 weeks. All piglets were housed in two groups of intact litters until weaning at the age of 4 weeks. As part of farm management practice, piglets from litter A were separated from the sow hours before sampling at 28 days, and piglets from litter B were moved to another pen the day after sampling at weaning. After weaning, piglets from both litters were mixed with piglets from other sows and kept in larger groups. Piglets and sows enrolled in this study did not show any illness and therefore did not receive any additional treatment or antimicrobials. A nasal swab was obtained from all piglets within the minutes after birth (t = 0 days) using a cotton swab (Medical Wire & Equipment, Wiltshire, United Kingdom). Swabs were also obtained at 8 h, 16 h, and 24 h (t = 1 day) after the first sampling, after which the piglets were sampled daily (t = 2, 3, and 4 days) and, finally, weekly until the piglets were 6 weeks old (t = 7, 14, 21, 28, 35, and 42 days). Nasal swabs were suspended in 1 ml saline supplemented with 1 mM EDTA (molecular grade; Sigma-Aldrich, the Netherlands). Suspension was subsequently subsampled in 3 aliquots for (i) microbiome analysis, (ii) real-time PCR to quantify S. aureus in general (including LA-MRSA), and (iii) bacteriological culturing to enumerate MRSA.
Quantification of S. aureus by real-time PCR. Two hundred ml of the nasal swab suspension was used to quantify S. aureus using quantitative real-time PCR (qPCR). Briefly, phocine herpes virus (PhHV) was added to the sample as an internal amplification control (43). DNA was then extracted with the High Pure PCR template preparation kit (Roche, the Netherlands) according to the manufacturer's instructions, and the sample was eluted in 50 ml elution buffer. Then, 5 ml of sample DNA was used in a real-time PCR that quantified S. aureus by targeting the femA (44) and nuc (45) genes using a predefined standard curve. Quantitative results of the PCR are reported as log CFU-equivalents (CFUeq).
Enumeration of MRSA by culturing. A 10-fold serial dilution of the nasal swab sample suspension was prepared in phosphate-buffered saline (PBS) (Gibco, the Netherlands). Next, 100 ml of each dilution (10 21 to 10 24 dilution) was plated on MRSA selective medium (Brilliance MRSA 2 agar; Oxoid, the Netherlands) and incubated at 37°C for 18 to 24 h. MRSA-suspected colonies were counted, and the number of CFU of MRSA was calculated and reported as log CFU. One MRSA-suspected colony from each sample was confirmed as LA-MRSA by targeting the ST398-specific DNA fragment C01 (46), and methicillin resistance was tested by using a mecA (44) PCR. In case the C01 gene-specific PCR was negative, S. aureus-specific PCRs targeting the femA (44) and nuc (45) genes were performed.
DNA extraction and sequencing. DNA extraction was performed using a modified version of Mag-Mini bead-beating and a magnetic bead procedure (LGC Genomics, Berlin, Germany) as described by Wyllie et al. (47). Amplicon libraries targeting the V4 region of the 16S rRNA gene were prepared using 515F (GTGCCAGCMGCCGCGGTAA) and 806R (GGACTACHVGGGTWTCTAAT) universal primers. Sequencing was performed on an Illumina MiSeq platform using v2 chemistry (2 Â 250 bp) (48). Similarly, libraries amplifying the tuf gene, a discriminatory target for Staphylococcus species were prepared using the oligonucleotides (23) tuf-F (GCCAGTTGAGGACGTATTCT) and tuf-R (CCATTTCAGTACCTTCTGGTAA), and sequencing was performed on an Illumina MiSeq platform using v3 chemistry (2 Â 300 cycle). Nontemplate DNA extraction controls were also included in the amplification and sequencing protocol to monitor potential contamination.
Taxonomy was assigned to nonchimeric sequences using the naive Bayes (NB) RDP classifier natively implemented in QIIME 2 (54). For this, the classifier was trained explicitly on the region of the gene that was sequenced and used for classification with a bootstrap confidence threshold of 80%. We used the Greengenes reference database v13.8 clustered at 99% identity for classification of 16S ASVs (55). For the tuf data, we prepared a custom reference taxonomy database by retrieving full-length bacterium-originating tuf sequences from KEGG (56) (https://www.genome.jp/dbget-bin/www_bget?ko:K02358; accessed 2019) and used it for classification of tuf ASVs using the method described for 16S data. Additionally, for 16S amplicon data, we used SPINGO for species-level classification wherever possible (57).
Initial preprocessing of the ASV table was conducted using the decontam (58) and CoDaSeq (59) packages, whereby potential reagent contaminants were identified and removed using the frequencybased method implemented in the decontam package. Next, we filtered out ASVs based on prevalence and abundance criteria using the codaseq.filter function from the CoDaSeq package. Only ASVs present in .10% of samples with a relative abundance of .0.0001 were retained for downstream analysis, which resulted in 368 ASVs for the 16S and 204 ASVs for the tuf data set. Except in the case of alpha diversity, this filtered ASV count table was used for all the downstream bioinformatic analyses.
Statistical analysis of compositional data. All statistical analyses and graphical representations were performed in R using the packages vegan (60), CoDaSeq (59), zCompositions (61), rmcorr (62), Ggplot2 (63), Heatmaply (64), and ComplexHeatmap (65). Moreover, GraPhlAn was used for visualization of phylogenetic trees generated from species-level summarized 16S and tuf data sets (66). To account for the complex compositional structure of the microbiome data and to avoid the likelihood of generating spurious correlations, we first imputed the zeros in the abundance metrices using the count zero multiplicative replacement method (cmultRepl, method = "CZM") implemented in the zCompositions package and applied a centered log-ratio transformation (CLR) using the codaSeq.clr function in the CoDaSeq package. Because the ASV table was summarized at different taxonomic levels (from phylum to species level), we used CLR transformation on each taxonomic level separately. Alpha diversity was determined using Chao1 (richness) and Shannon index (diversity), and the nonlinear association of a-diversity with time point (as numeric) was accessed by fitting the loess splines using the Ggplot2 package. The statistically significant association of time points with alpha diversity was tested using the rmcorr package. Principal-coordinate analysis (PCA) was carried out using the prcomp function in R using the Aitchison distance matrix (CLR plus Euclidean distances). Permutational multivariate analysis (PERMANOVA [67]) was performed on the Aitchison distances with 9,999 permutations to evaluate the effect of different clinical variables (i.e., time point and litter) on the nasal microbiota composition.
Association of microbiota data with metadata. Since the nasal piglet microbiota during first two initial time points (0 h and 8 h) was not stable and harbors bacteria that are commonly found in animal feces, in the uterus and cervix of the sow, or in soil, we considered them relatively unstable and excluded them (n = 15 from 16S and n = 16 from tuf) from all statistical analyses. Associations between taxa and log CFU of MRSA and log CFUeq of S. aureus were obtained using repeated measure correlation analysis from the rmcorr package (62), which determines the relationship between two continuous variables while controlling for between-individual variance. Rmcorr identifies a common regression slope and thereby estimates the association shared among all the individuals. Most popular correlation techniques, such as Pearsons correlation, assume independence of error between observations and thus cannot be used where more than one data point is obtained from individuals. Rmcorr accounts for this nonindependence among observations in repeated measurement data by removing measured variance between individuals. Similar to the Pearson correlation coefficient, the rmcorr coefficient (r rm ) ranges from 21 to 11 and reports the strength of the linear association between two variables. The rmcorr method calculates the rmcorr coefficient (r rm ), P value, and a 95% confidence interval of the rmcorr coefficient by bootstrapping the samples (n = 100). So, when there is no strong heterogeneity across subjects and parallel lines provide a good fit, the rmcorr effect size (r rm ) will be large, with tight confidence intervals. Next, in order to confirm the rmcorr correlation findings, we performed logistic regression analysis using multivariate analysis by linear models (MaAsLin2 v1.1.1) considering litter and animal ID as random effects and MRSA/S. aureus colonization events as categorical data (68). MaAsLin2 performs boosted, additive general linear models between metadata and microbial abundance. Boosting of metadata and selection of a model was performed per taxon. Microbial abundances were CLR-transformed at each taxonomic level to account for the compositional nature of the data. Multiple testing correction was carried out with the Bonferroni method where appropriate for all statistical tests (69).
Coabundance analysis. Following rmcorr correlations between each pair of species, species-level summarized taxa were clustered into the coabundant groups (CAGs) based on their CLR-transformed abundances across all the samples. Correlations were considered significant below a q value cutoff 0.05 after Benjamini-Hochberg (BH) multiple testing correction. Hierarchical clustering was performed using the Spearman distance matrix and ward linkage clustering to identify CAGs cooccurring with each other across all time points. Next, the dendrogram was cut at a height of 1.0 to generate nine and three different CAGs for the 16S and tuf data sets, respectively. Taxa comprising each CAG were plotted individually to understand longitudinal dynamics of the microbiome and its association with MRSA and S. aureus colonization.
Data availability. Sequence data are available under NCBI BioProject accession no. PRJNA687981.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.