A Single Intranasal Dose of Bacterial Therapeutics to Calves Confers Longitudinal Modulation of the Nasopharyngeal Microbiota: a Pilot Study

ABSTRACT To address the emergence of antimicrobial-resistant pathogens in livestock, microbiome-based strategies are increasingly being sought to reduce antimicrobial use. Here, we describe the effects of intranasal application of bacterial therapeutics (BTs) on the bovine respiratory microbiota and used structural equation modeling to investigate the causal networks after BT application. Beef cattle received (i) an intranasal cocktail of previously characterized BT strains, (ii) an injection of metaphylactic antimicrobial (tulathromycin), or (iii) intranasal saline. Despite being transient colonizers, inoculated BT strains induced longitudinal modulation of the nasopharyngeal bacterial microbiota while showing no adverse effect on animal health. The BT-mediated changes in bacteria included reduced diversity and richness and strengthened cooperative and competitive interactions. In contrast, tulathromycin increased bacterial diversity and antibiotic resistance and disrupted bacterial interactions. Overall, a single intranasal dose of BTs can modulate the bovine respiratory microbiota, highlighting that microbiome-based strategies have potential in being utilized to mitigate bovine respiratory disease in feedlot cattle. IMPORTANCE Bovine respiratory disease (BRD) remains the most significant health challenge affecting the North American beef cattle industry and results in $3 billion in economic losses yearly. Current BRD control strategies mainly rely on antibiotics, with metaphylaxis commonly employed to mitigate BRD incidence in commercial feedlots. However, the emergence of multidrug-resistant BRD pathogens threatens to reduce the efficacy of antimicrobials. Here, we investigated the potential use of novel bacterial therapeutics (BTs) to modulate the nasopharyngeal microbiota in beef calves, which are commonly administered metaphylactic antibiotics to mitigate BRD when sourced from auction markets. By direct comparison of the BTs with an antibiotic commonly used for BRD metaphylaxis in feedlots, this study conveyed the potential use of the BTs to modulate respiratory microbiome and thereby improve resistance against BRD in feedlot cattle.

feedlots for production. During feedlot placement, cattle are most susceptible to BRD, with the majority of cases occurring within the first 60 days of feedlot placement (3). Primary viral infections or stressors, which include weaning, shipping, and comingling with new pen mates, are proposed to reduce host immunity during transition to feedlots (4). Consequently, opportunistic bacterial pathogens residing in the upper respiratory tract proliferate and translocate to the lungs, causing bronchopneumonia (5). The main BRD-associated bacteria include Mannheimia haemolytica, Pasteurella multocida, Histphilus somni, and Mycoplasma bovis (6). As a result of increased BRD susceptibility, commercial feedlots rely on antimicrobial-driven approaches to prevent BRD infections in cattle (7).
Long-acting injectable antimicrobials are commonly administered to cattle entering feedlots (i.e., metaphylaxis) for BRD prevention (8). For example, the macrolide tulathromycin was used at feedlot entry by 45.3% of feedlots in the United States for BRD mitigation (9). Metaphylactic antimicrobials treat lung infections that may be prevalent in calves entering feedlots and also prevent infection during the course of their bioactivity. It is also likely that metaphylactic antimicrobials reduce prevalence and proliferation of BRD pathogens in the upper respiratory tract, a prerequisite to lung translocation (10). However, antimicrobial resistance in BRD pathogens has increased over the last 10 years (11)(12)(13), with resistance to tulathromycin recently being detected in more than 70% of M. haemolytica and P. multocida isolated from feedlot calves (14). In addition, resistance elements have been detected in mobile elements from the BRD-associated Pasteurellaceae family, conferring multidrug resistance to antimicrobials used for both prevention and treatment of BRD in feedlot cattle (15,16). Antimicrobial resistance in BRD pathogens therefore threatens the efficacy of currently used antimicrobials in beef cattle. Alternatives to antimicrobials are therefore needed for use in novel feedlot management strategies.
The respiratory microbiota contributes to host health by providing colonization resistance against pathogens and maintaining homeostasis (17). It is hypothesized that disruption of the bovine respiratory microbiota can promote the proliferation of BRD pathogens (18). Indeed, several management factors, including transportation to a feedlot (19), diet composition (20), and antimicrobial administration (21), alter the upper respiratory tract microbiota of cattle. Recently, we observed that several lactic acid-producing bacteria (LAB) were inversely correlated with Pasteurellaceae in the nasopharynx of cattle transported to an auction market and, subsequently, a feedlot (22). Genera within the LAB order Lactobacillales were shown to be reduced in cattle that develop BRD (23,24). Additionally, isolates of LAB originating from the bovine respiratory tract have been shown to directly inhibit BRD pathogens (25,26). These data support that LAB are important community members of the bovine respiratory tract and may be integral to providing colonization resistance against BRD pathogens. Consequently, we have previously developed BTs comprised of six Lactobacillus strains that were characterized for their in vitro inhibition and exclusion of M. haemolytica and adherence to and immunomodulation of bovine turbinate cells (26). The intranasal inoculation of these BT strains was able to inhibit colonization by M. haemolytica in experimentally challenged dairy calves (27). In the present study, we further evaluated the longitudinal effects of these intranasal BTs on the nasopharyngeal microbiota of beef calves after a single intranasal dose. The effects of the BTs were also compared to those of tulathromycin, a common antimicrobial used for metaphylaxis, administered subcutaneously.

RESULTS
Calf health and weight gain. Calves were monitored daily for clinical signs of BRD throughout the study. Elevated rectal temperature ($39.7°C) was detected in five calves (three from the control group [CTRL] and two from the bacterial therapeutics group [BT]) (see Table S1 in the supplemental material) that recovered in response to a single injection of the antibiotic Micotil (tilmicosin). The remainder of the experimental calves were healthy during the course of study. The calves were weighed first at arrival (day 21) which was 24 h prior to treatment administration and subsequently on a biweekly basis for the 42 days of study. The average daily gain was not different among treatments (P = 0.506) (Fig. S1).
Prevalence of BRD-associated pathogens determined by NP swab culturing. The presence of M. haemolytica, P. multocida, and H. somni was evaluated by culturing the nasopharyngeal (NP) swabs collected during the first 28 days of the study (Fig. 1). Overall, M. haemolytica had the highest prevalence of the three pathogens on day 21 (24 to 40% across treatments). No significant difference was detected between treatment groups (P . 0.05) at any sampling time point for M. haemolytica; however, there was a tendency (P = 0.06) for reduced prevalence in metaphylaxis (MP)-treated calves on day 7. For P. multocida, prevalence in CTRL (10 to 40%) and BT (5 to 36%) calves was similar (P . 0.05). In MP calves, prevalence of P. multocida was lower (range, 0 to 5%) on days 7 and 14 than in CTRL and BT calves (P , 0.05). Overall, the prevalence of H. somni remained low throughout the study for all treatment groups and was not different between treatment groups at any sampling point (P . 0.05). Only the MP group had colonization rates of 0% for P. multocida (days 1, 4, and 14) and M. haemolytica (day 1), which occurred after metaphylactic treatment.
Total bacteria and Lactobacillus in NP swabs determined by quantitative PCR. The abundances of total bacteria and Lactobacillus in nasal swabs was estimated by quantifying the gene copy numbers of general and Lactobacillus-specific 16S rRNA (Fig. 2). An interaction between treatment and time affected the total bacterial number (P = 0.02) ( Fig. 2A). Compared to BT and CTRL calves, tulathromycin injection reduced bacteria in NP swabs on days 4, 7, and 42 (P , 0.05). Total bacteria in NP swabs of BT and CTRL calves were not different over the course of the study (P . 0.05).
The estimated number of Lactobacillus in NP swabs was affected by the interaction of treatment and time (P = 0.01). The mean total Lactobacillus 16S rRNA copy number per swab from the BT group increased (P , 0.01) from day 21 to day 1 (24 h post-BT inoculation) (Fig. 2B)  (ii) Community structure of the NP microbiota. Permutational multivariate analysis of variance (PERMANOVA) revealed that an interaction of treatment and time had an effect (R 2 = 0.04, P = 0.001) on the microbial structure of the NP microbiota (Fig. 3). However, time had a larger effect on microbial structure (R 2 = 0.141, P = 0.001) than treatment (R 2 = 0.034, P = 0.001). As indicated by the detrended correspondence analysis (DCA) plots (Fig. 3), the microbiota tended to cluster by early (days 21 to 2), mid-(days 4 to 7), and late (days 14 to 28) time points. Clustering according to treatment was most evident on days 28 and 42.
Alpha diversity, as assessed by microbial richness and Shannon diversity index, revealed that both indices were affected by a treatment and time interaction (P , 0.05) (Fig. 4). After BT inoculation, the NP microbiota of BT calves was reduced in richness throughout the study compared to CTRL and MP calves (P # 0.014), except for on day 7 when richness was similar in BT and CTRL groups (P = 0.213). In contrast, MP calves had increased richness on days 7 to 42 compared to BT calves and increased richness on days 7 and 14 compared to CTRL calves (P # 0.001). Similar to richness, the Shannon diversity index was reduced in BT calves from days 7 to 42 compared to CTRL and MP treatments (P # 0.0001). The Shannon diversity was greater in MP calves than to BT and CTRL calves, but only on days 1 and 7 (P # 0.016).
(iii) Composition of the NP microbiota. Across time and treatment groups, a total of 14 different bacterial phyla were identified, among which Proteobacteria (36.4%), Tenericutes (22.2%), Firmicutes (17.4%), Actinobacteria (12.9%), and Bacteroidetes (9.9%) were the most relatively abundant and together constituted 98.8% of the sequences. The diversity of genera within each phylum varied, with the relative abundance of a single genus ranging from ,1% to 100% of a phylum. Overall, the 10 most relatively abundant genera across treatments and time included Mycoplasma (22.2%), Moraxella (18.8%),  Changes in microbial composition following BT and tulathromycin treatment. (i) Changes in the five most relatively abundant phyla. Noticeable changes in NP microbial composition at the phylum level were observed in response to treatment and time effects (Fig. 5). The relative abundance of Proteobacteria in CTRL calves varied over the 42 days of study, with a gradual increase in the first 7 days followed by a decline in the remaining 5 weeks of the study. In BT calves, the relative abundance of Proteobacteria was similar to CTRL calves except for on day 28, when it was increased (P = 0.02). In MP calves, however, Proteobacteria had a lower relative abundance than CTRL calves on days 4, 7, and 14 (P # 0.008). The relative abundance of Firmicutes did not differ between BT and CTRL groups at any sampling time (P . 0.05). However, Firmicutes was increased in MP calves within the first 14 days of antibiotic injection compared to CTRL calves (P , 0.05). Similarly, Bacteroidetes became significantly enriched in MP calves on days 7 and 14 relative to the CTRL group (P # 0.022). Compared to CTRL calves, the abundance of Actinobacteria was reduced in BT calves on day 42 (P = 0.011), while it was increased in MP calves on day 1 (P = 0.038). The relative abundance of Tenericutes was not affected by treatment (P . 0.05).
(ii) Changes in lactic acid-producing bacteria at the family level. The LAB taxonomically belong to the Lactobacillales order (28). Five different LAB families, including Aerococcaceae, Carnobacteriaceae, Enterococcacease, Lactobacillaceae, and Streptococcaceae, were detected in the present study (Fig. S2). Overall, no treatment effects on the abundance of these families were observed, although the relative abundance of Lactobacillaceae was .10% in four calves from the BT group.
(iii) Changes in microbial composition at the genus level. Beta-binomial regression analysis allowed for identification of compositional differences of NP genera between treatment groups across time. In total, we identified 28 genera within the BT and MP groups whose change in relative abundance from day 21 (baseline) was Intranasal bacterial Therapeutics and Microbiota mSystems significantly different from changes that occurred in the CTRL group (P , 0.05). As shown in Fig. 6, 4 of the 10 most abundant genera (Ruminococcaceae_UCG-005, Psychrobacter, Jeotgalicoccus, and Planococcus) were included among these taxa that differed between treatment groups. In general, most of the 28 taxa became less abundant in the BT group following BT inoculation. Among these, Ruminococcaceae_NK4A214_group, Paeniclostridium, Lachnospiraceae_UCG-005_group, and Cellvibrio experienced a consistent decline during the entire post-BT inoculation period. For some taxa, the magnitude of change in relative abundance from the baseline varied among sampling times. For BT calves, the most significant reduction from baseline (day 21) was observed for Ruminococcaceae_NK4A214_group (day 2), Lachnospiraceae_NK4A136_group (day 1), Facklamia (day 28), Celllvibrio (day 2), and Acetitomaculum (day 4). Lactobacillus was the only taxon in the BT group that experienced an increase in relative abundance that was .1, occurring on day 1. In contrast to the BT group, most of the 28 genera became enriched in the MP calves, with several taxa having increases in relative abundances which were .2. The most immediate and consistent enrichment following tulathromycin injection was observed for Jeotgalibaca.
(iv) Changes in relative abundance of Lactobacillus spp. To identify whether the significant enrichment of genus Lactobacillus observed on day 1 (Fig. 2B) was due to the inoculated BT Lactobacillus strains, species-level taxonomic identification was performed on sequencing data from the Lactobacillus genus (Fig. S3). Within Lactobacillus, (v) Changes in relative abundance of BRD-associated genera. The relative abundances of Mannheimia, Pasteurella, Histophilus, and Mycoplasma, which encompass BRDassociated pathogens, were not affected by treatment (P . 0.05) (Fig. S4).
Microbial interactions and dynamics of the NP microbiota. (i) Interaction network structure among all observed genera. To evaluate overall dynamics of microbial communities, ecological modeling was used to analyze the interaction of all genera. As shown in network plots (Fig. 7), distinct microbial interaction network structures were observed between CTRL, BT, and MP groups. Compared to the CTRL group, the  interaction network of microbiota from BT calves was more complex with a greater number of genera-genera interactions. In contrast, there was a large decrease in genera-genera interactions among the microbial community of calves in the MP group, with only 22 genera identified in the network model. Even among these 22 genera, the interactions were connected by two separate hubs, indicating that the tulathromycin injection diminished the interaction network among the NP microbial community.
(ii) Structure of causal relationships among 16 selected genera. Based on the observed distinct changes in species interaction network among NP microbiota in response to both BT and MP treatments, we further evaluated the relationship of 16 targeted genera using causal structure-based path modeling. Path analysis is a member of the structural equation modeling tools that enables the identification of causal relationships between measured variables (29). This path model analysis provides greater detail on the interaction between selected bacterial genera, as it predicts causality and direction of the causality, and accounts for the time effect and unmeasured effects. Three path models (CTRL, BT, and MP) are described in Table 1 and are depicted as diagrams of causal relationships of the relative abundances in Fig. S5.
Details on model fit statistics and the number of iterations performed by the CALIS procedure are shown in Table 1. The P values for the chi-square statistic for the modified path models were .0.05, indicating good model fitting. The root mean square error of approximation (RMSEA) was ,0.05, and the 90% confidence limits of the RMSEA were also ,0.05. Therefore, the null hypothesis that the modified path models closely fit the data was retained. The values of the Tucker-Lewis index were also .0.9, where the value of a true model would be expected to be 1.
The variables that constituted the subsets of exogenous (independent) and endogenous (dependent) variables varied in the three path models ( Table 1). The day was initially hypothesized to be exogenous, and no path or covariance was added to change the position in the model of day from being an independent explanatory variable. The process of model modification positioned the abundance of Histophilus as an exogenous variable in the modified path models, and it was a predictor variable with a statistically significant standardized total effect on Alloprevotella (b model 1 Hist;Allo = 0.18 6 0.07, P = 0.005) in model 1 (of the BT group) (Table S2), while it was not a predictor variable in model 2 (CTRL) or model 3 (MP) (Fig. S5).
(iii) Robust relationships. Path relationships that were statistically significant in the three modified models were considered to be robust. As such, the negative and statistically significant (all P , 0.0001) values of the standardized direct effect from day to Atopostipes in the modified path models represent a robust causal relationship. Therefore, under the conditions of this experiment, the relative abundance of Atopostipes is expected to decrease monotonically over time. In contrast, the standardized direct effect from day to Psychrobacter was positive in the modified path models, and therefore, the relative abundance of Psychrobacter is expected to increase monotonically over time under conditions similar to the experiment. Robust causal relationships were also identified between observed genera. For example, the standardized positive and statistically Rumi;Jeot = 0.50 6 0.06). These positive path relationships can be interpreted as commensal biological relationships where the bacteria that are a source of variance experience no benefit or harm, but the bacteria that extract variance can benefit by increasing monotonically in relative abundance.
(iv) Co-occurrence of path relationships for treatment groups. Path relationships among the observed genera were identified in model 1 (BT group) and model 2 (CTRL group). For the BT (model 1) and CTRL (model 2) groups, there were statistically significant (P # 0.0001) and positive standardized total effects from the Christensenellaceae_R7_group  Rumi;Rike = 0.78 6 0.06).
In the MP group (model 3) (Fig. S5C), there were fewer path relationships that cooccurred in the BT (model 1) and CTRL (model 2) groups. The causal relationships between genera that were unidirectional for the CTRL and BT groups were frequently bidirectional in the MP group. The relative abundances of seven genera become exogenous variables in model 3 (Table 1). One genus in the MP group (Atopostipes) exhibited standardized total effects that co-occurred in the BT group. The standardized total effects from Atopostipes were positive and statistically significant (P # 0.02) compared to Acetitomaculum (b Atop;Jeot = 0.10 6 0.04). The standardized total effects from the relative abundance of Lactobacillus to the relative abundances of other bacteria are shown in Table S3. In this case, the effects from Lactobacillus to Acetitomaculum, Acinetobacter, Alloprevotella, and Jeotgalibaca were positive in both the BT and CTRL groups. In model 3 (MP), the relative abundance of Lactobacillus did not affect the relative abundances of the selected genera. Based on CTRL data (model 2), Lactobacillus was modeled to increase monotonically the relative abundance of Acetitomaculum, which inhibits Mannheimia (Fig. S5). The effect of Lactobacillus on Mannheimia was therefore indirect in model 2. In the BT group (model 1), Lactobacillus had a direct positive effect and an indirect negative effect on Mannheimia. The indirect inhibition of Mannheimia was mediated by Phascolarctobacterium, which inhibited Mannheimia. The standardized total effects from Lactobacillus to Mycoplasma and Pasteurella were negative in the CTRL group, but these effects were not statistically significant in the BT group.
Antimicrobial resistance determinants in the NP microbiota. The macrolide resistance gene msr(E) increased in the NP microbiome of MP calves during the first 28 days of the study (P , 0.01) (Fig. 8A). The abundance of msr(E) was greater in MP calves than CTRL and BT calves during the last 2 weeks of study (P , 0.05). The Intranasal bacterial Therapeutics and Microbiota mSystems abundance of the tetracycline-resistant gene, tet(H), was not affected by treatment or time (P . 0.05) (Fig. 8B).

DISCUSSION
Studies have suggested that mutualistic and antagonistic interactions take place within the microbial community of the bovine respiratory tract (24,25). These interactions may contribute positively or negatively to microbiota-mediated colonization resistance against respiratory pathogens (18,30). Specifically, in the nasopharynx, LAB have been negatively associated with Pasteurellaceae, and certain LAB strains were capable of directly inhibiting BRD bacterial pathogens (25,26). These data suggested that LAB have potential as BTs for mitigating BRD-associated pathogens. In the present study, we evaluated the effect of inoculating 6 previously characterized Lactobacillus strains (26) directly into the upper respiratory tract of cattle on the microbiota and compared those effects to a common metaphylactic antimicrobial, tulathromycin. Throughout the study, all calves remained healthy, with the exception of 5 that were treated for BRD. However, the rates of BRD were not affected by treatment, and weight gain was similar across all treatments, indicating that the BTs did not adversely affect inoculated calves.
Colonization by BTs. A single application of BTs was selected in order to fit within modern management systems of beef cattle, which employ a single processing event upon arrival at feedlots. Despite the BT strains originating from the nasopharynx of feedlot cattle and displaying strong in vitro adhesion to bovine turbinate cells (26); real-time PCR and 16S rRNA SV analyses indicated that colonization by BT strains was transient, lasting up to 48 h. Similar findings from gastrointestinal studies support the difficulty that exogenous strains have in colonizing established microbial communities (31,32). The presence of similar indigenous species can affect colonization by exogenous bacteria (33), potentially by limiting resource availability (34,35). For example, colonization by Bifidobacterium longum (AH1206) was impeded in the gastrointestinal tract when endogenous B. longum was already present (35). Indeed, Lactobacillus species were detected prior to inoculation of BTs, though these were mainly assigned to the operational taxonomic unit (OTU) L. fermentum/L. mucosae. Thus, inoculating the BTs earlier in the life of cattle may increase colonization potential if done prior to establishment of similar Lactobacillus spp.
Longitudinal effects of BTs and tulathromycin on the respiratory microbiota. Changes in the composition of the NP microbiota were observed after BT administration. Despite BT treatment resulting in both increases and decreases of 28 taxa that significantly changed from baseline levels prior to inoculation (Fig. 6), the strongest BT effects on these taxa were inhibitory. The BT strains have previously been shown to produce lactate or hydrogen peroxide or encode bacteriocins (26), which may have led to inhibition of resident bacteria. These direct impacts, however, would have only occurred within 48 h of the BTs colonizing the nasopharynx. It is interesting to note despite being transient, a single administration of the BTs had a long-term impact on structure and diversity of the NP microbiota. Both a reduction in richness and diversity was observed up to 42 days after BT inoculation. Information regarding the effect of probiotics or BTs on respiratory microbiota in animal models is limited; however, similar effects have been reported for probiotics targeting the gastrointestinal microbiota. For example, Zhang and colleagues observed that the gut microbial community was impacted up to 2 weeks after administration of a Lactobacillus casei probiotic strain (36). Other studies have reported prolonged effects on the gut microbiota for up to 1 (33) and 5 (37) months following cessation of a probiotic cocktail containing 11 strains fed to humans. The underlying mechanisms by which the BT Lactobacillus strains exerted longterm modulation of the microbiota are difficult explain but are likely related to their initial effects on community members.
Bacterial communities in most niches form complex ecological interaction webs, and such interactions are important in maintaining microbiome homeostasis and a symbiotic relationship between microbe and host (38). We therefore evaluated the microbiome-wide community networks to gain insights into the changes in microbial community interactions in NP microbiota in response to BTs and antibiotic administration. The interaction network among all observed genera was predicted based on interaction networks from sequencing data. After observing the distinct interaction network structure among treatment groups, we decided to further investigate the causal networks among targeted taxa using structure equation modeling. Mainali et al. were the first to apply causal models to detect interaction networks in the human microbiome using conditional Granger causality (39). These authors argued that the causal models may provide more accurate prediction of interaction networks among the microbial community than standard correlation/network analysis, as correlation is neither necessary nor sufficient to establish causation, and environmental filtering can lead to a correlation between noninteracting taxa.
The path model of the BT group revealed that a moderate degree of alterations had taken place in the causal relationship structure among observed genera, as seen by the changes in the magnitude of the direct and indirect effects from one genus to another. It is possible that the indirect effects of the BTs are what caused the prolonged effects on community diversity. Given the computational complexity of the path models, they were limited to 16 genera. However, it was interesting that ecological network analysis of all observed genera showed a more complex network for the BT calves than both CTRL and MP calves. Yang and colleagues suggested that probiotics (Paracccus marcusii DB11 and Bacillus cereus G19) promote intestinal microbiota homeostasis by enhancing species-species interactions and increasing the number of connecters and/or module hubs within the network (40). Thus, BT administration may have promoted NP microbiota homeostasis by strengthening and promoting speciesspecies interactions. This, in turn, may have reduced colonization potential by new bacteria, which is supported by the decreased richness observed for BT calves. While we did not sample calves prior to arrival, several studies have shown that diversity of the NP microbiota increases within days after feedlot placement, and it has been suggested to be linked to susceptibility to BRD (19). It would therefore be interesting to measure the effect of BT inoculation in calves prior to transport in future studies to evaluate whether they promote stabilization of the respiratory microbiota during transport.
Tulathromycin also altered the NP microbiota structure, diversity, and composition compared to CTRL calves. Measured by real-time PCR, the total bacteria per NP swab was reduced in MP calves on days 4 and 7. This was likely due to inhibition of members of the phylum Proteobacteria, which decreased in MP calves. Interestingly, despite this reduction, Shannon diversity and richness were increased in MP calves on day 7 compared to the other treatments. While Holman et al. (21) did not see an increase in diversity of the NP microbiota of calves administered tulathromycin, they did observe an increase in diversity in calves administered oxytetracycline. In addition, like our study, Holman and colleagues observed an increase in Ruminococcaceae_UCG-005, Rikenellaceae_RC9_gut_group, Phascolarctobacterium, Facklamia, Jeotgalibaca, and Acinetobacter following tulathromycin injection (21), supporting that those injectable antimicrobials may lead to an increase in microbial richness and diversity of the respiratory microbiota.
Most often, the diversity of the microbial community is claimed to be positively associated with the stability of microbiota (41) and health (42). Typically, however, studies show positive associations between bacterial diversity and health related to the gastrointestinal tract (43,44). However, it has also been argued that diversity in host-associated microbial communities may not always be associated with microbial community stability and health. For example, higher bacterial diversity and richness were observed in the upper respiratory tract of children with invasive pneumococcal disease than in healthy children (45). It was proposed by the authors that the higher diversity and richness of the NP microbiota was associated with impaired immune response. It was interesting in our study that antimicrobial treatment increased diversity of the respiratory tract, whereas antimicrobial administration can reduce diversity in the gastrointestinal tract. Perhaps this is a reflection of increased exposure to exogenous bacteria that the respiratory tract faces compared to the digestive tract.
Functional properties of the gut microbiota, including colonization resistance at the community level, are believed to be maintained by the active interactions among genetically distinct and diverse microbial species, which allows the microbial community to perform complex metabolic activities. A single monospecies population or multispecies but with no interconnectivity could not provide such collective functions of microbiota (46). In addition, the functional activities and stability of a microbiota are influenced by the positive and negative feedback loops generated as a result of the cooperation (47) and competition (48) among the different microbial species (49). Also, evolution of species-species interactions has been reported to determine microbial community productivity in new environments (50). Although proportionally equal positive and negative interactions between the bacterial species within NP microbiota were observed in both CTRL and BT groups, such cooperative and competitive interactions were more intensive in the BT group than in the CTRL group. In contrast, there were only positive interactions observed between these fewer species that remained in the network model as interconnected species in the MP group.
The stability of a mammalian microbiota depends on how the species interacts with one another (51). Weak and competitive interactions are stabilizing, and they limit positive feedback loops and the possibility that if one species decreases, it will result in the decrease of others. Cooperative interactions determine the productivity of microbiome, which is the efficiency of converting resources into energy (51). In the MP group, there was only a cooperative interaction observed in NP microbiota, and the competitive interaction was missing, suggesting that NP microbiota in MP cattle would most likely experience dysbiosis, whereas in BT calves, both cooperative and competitive interactions were present with higher magnitude. This indicated that the NP microbiota in BT calves was most likely stable with normal community functions. It has been argued that the asymmetry of an unhealthy microbiome can relate to nonneutral states created by strong stressors, reducing host ability to contain certain bacteria and resulting in overgrowth of abundance and multiplication of species (52). We observed that the abundance of most of the observed significant taxa (n = 28) was enriched in calves that received antibiotic. The overgrowth of these taxa therefore might be due to the nonneutral state of microbiota induced by antibiotic.
While the increase in diversity following tulathromycin was unexpected, we hypothesize that it was related to the effects of the antimicrobial on the community network. For MP calves, 7 genera were determined to be exogenous to the developed path model, which had fewer interactions than those developed for the BT and CTRL groups of calves. In addition, the ecological network was also far less complex, and together, these data indicate a substantial perturbation of the ecological network in MP calves. Likewise, Yang et al. (40) reported that the use of the antimicrobial florfenicol resulted in deterioration of the ecological network among intestinal microbiota of sea cucumbers, leading to the homeostatic collapse of microbiota. Thus, the deterioration of the species interactions by tulathromycin may have made the NP microbiota in MP calves more permissive to exogenous bacteria colonization, increasing diversity.
While the impact this may have on cattle health was not possible to evaluate in the present study, it is an area that warrants further investigation. Despite the use of metaphylactic antimicrobials, the overall rate of BRD in feedlots has not declined over the last 40 years (4). In fact, some feedlots are observing a change in BRD incidence, with cases occurring later in the feeding period (E. Janzen, University of Calgary, personal communication). While unsupported, it is tempting to speculate that metaphylaxis treatment may prevent short-term BRD incidence through direct inhibition of pathogens, but the causal deterioration of keystone taxa (53) may increase the risk of longer-term respiratory complications. In contrast, we observed that BTs promoted microbiome homeostasis by stabilizing the NP bacterial community structure and genera interactions.

Longitudinal effects of BTs and tulathromycin on BRD-associated pathogens.
The MP treatment reduced the prevalence of culturable M. haemolytica. In support of this, tulathromycin injection has previously been shown to reduce NP colonization of M. haemolytica in feedlot cattle (10). In contrast, despite showing strong inhibitory properties against M. haemolytica in vitro, the BTs did not reduce prevalence of M. haemolytica. However, the BT strains were tested for competitive exclusion in vitro (26). The fact that 24% of calves in the BT group were M. haemolytica positive prior to administration suggests that the BT strains may not be effective for displacement of M. haemolytica. Alternatively, the BT strains may have a more pronounced effect at inhibiting M. haemolytica when this pathogen in present in greater concentrations.
Despite MP calves being the only group to have several time points where no M. haemolytica or P. multocida could be cultivated from NP swabs, sequence analysis of Mannheimia, Pasteurella, Histophilus, and Mycoplasma revealed no differences in relative abundance of these BRD-associated genera. For Mannheimia, Pasteurella, and Histophilus, this was likely a result of interanimal variation and overall limited relative abundance of these genera. While a previous study did show that tulathromycin injection reduced NP Pasteurella, the calves in that study were ranch derived and had higher levels of Pasteurella at feedlot entry (21). The limited number of BRD cases and reduced relative abundance of BRD-associated genera in our study was likely attributed to the calves being at lower risk for BRD, having come from a single source and being placed into individual stalls instead of being comingled.
It was interesting that Histophilus was exogenous to microbial path models for all three treatments. H. somni has previously been shown to have a 32.5 times greater chance of being isolated from feedlot calves after 40 days on feed (14). This suggests that genera that were not included in the models, or external factors associated with feedlots, may be related to its colonization of cattle. Lactobacillus had a direct positive effect on Mannheimia for BT calves, though indirect negative effects were also observed. This finding is difficult to explain but may be related to limiting model analysis to 16 genera. Although the relative abundance of Mannheimia did not increase as a result of BT treatment, future studies with Lactobacillus-based BTs should be tested for their effects on Mannheimia. In contrast, for CTRL calves, Lactobacillus had overall negative effects on Mannheimia and Pasteurella, though they were not direct. Thus, this supports previous data showing the importance of indigenous Lactobacillus in bovine respiratory health (22,24).
Longitudinal effects of BTs and tulathromycin on antimicrobial resistance determinants. The macrolide resistance gene msr(E) encodes a multidrug efflux pump conferring resistance to macrolides, including tulathromycin (54). This gene has been detected in BRD pathogens and has been detected within integrative conjugate elements (15,16). Given the increase in macrolide resistance in feedlot BRD pathogens over the last 10 years (11,14), evaluation of resistance genes in respiratory bacteria is important to maintaining proper selection of antimicrobials. Despite altering the microbiota, BT-treated calves did not select for bacteria that carried msr(E) or tet(H). In contrast, msr(E) increased in MP calves, showing that tulathromycin administration selected for bacteria carrying this resistance gene. Similarly, tulathromycin use has previously been associated with increased antimicrobial resistance determinants in NP microbiomes of feedlot cattle (21,55).
Conclusion. A single dose of intranasal BTs, which were developed from bovine respiratory commensal Lactobacillus spp., induced longitudinal modulation of the NP microbiota in auction market beef calves, with no adverse effects on animal health and growth performance. While no differences in the relative abundances of BRD-associated genera were observed between treatments, the ecological networks of NP bacteria from BT-treated calves became more integrated. It was proposed that this resulted in a more stable microbiome with increased resilience against exogenous microorganisms. In contrast, disruption of the microbiome after tulathromycin treatment reduced resilience, leading to increased diversity in MP calves. Overall, this study showed that the bovine respiratory microbiota can be altered by administration of BTs and may therefore provide new opportunities to enhance microbiome-mediated respiratory resistance against BRD pathogens. We also showed the usefulness of employing network and structure equation modeling to evaluate the impact of BTs on host microbiota. Future studies should consider the optimal administration of BTs, as improved resilience against BRD pathogens may result from administration prior to calves being shipped to feedlots.

MATERIALS AND METHODS
Animals and experimental design. Animals used in this study were cared for in agreement with the Canadian Council for Animal Care guidelines (56). All the procedures and protocols with respect to animal handling and sampling were reviewed and approved by the Animal Care Committee at the Lethbridge Research and Development Centre, Agriculture and Agri-Food Canada (Lethbridge, AB, Canada).
Sixty cross-bred beef heifers, approximately 6 months old (initial body weight [BW] = 266 6 13 kg) and originating from a single cow-calf ranch, were purchased from a local auction market and transported to the Lethbridge Research and Development Centre feedlot (,10 km distance). Upon arrival, the calves were weighed, and NP swabs were collected (day 21). Calves were then blocked by weight and randomly assigned to three treatment groups (n = 20 per treatment) as follows: (i) the BT group received an intranasal cocktail of six Lactobacillus strains suspended in phosphate-buffered saline (PBS) in equal concentrations (3 Â 10 9 CFU per nostril), (ii) the metaphylaxis (MP) group received a subcutaneous injection of tulathromycin (2.5 mg/kg BW), and (iii) the control (CTRL) group received intranasal PBS without bacteria. Each of the 60 calves was housed in an individual pen throughout the study and was fed once daily a diet containing 75% barley silage, 22.5% dry roll barley, and 2.5% standard feedlot supplement. Calves had free access to drinking water.
Preparation of BT inoculum. The BT cocktail was a mixture of six Lactobacillus strains (1 Â 10 9 CFU mL 21 ), L. amylovorus (isolate 72B), L. buchneri (63A and 86D), L. curvatus (103C), and L. paracasei (3E and 57A). These isolates were inoculated on Lactobacillus De Man, Rogosa, and Sharpe (MRS) agar (Dalynn Biologicals, Calgary, AB, Canada) and incubated for 48 h at 37°C in 10% CO 2 . One day prior to nasal inoculation, a single colony of each strain was inoculated into 5 mL MRS broth and incubated at 37°C with agitation at 200 rpm. After 18 h of incubation, each bacterial culture was centrifuged at 7,600 Â g for 10 min, the supernatant was discarded, and the pellet was resuspended with prewarmed (37°C) PBS to achieve a target concentration of 1 Â 10 9 CFU per mL 21 using preestablished optical density at 600 nm (OD 600 ) values. The BT cocktail was prepared by mixing the six Lactobacillus isolates at equal ratios in PBS. One hour prior to inoculation, 3 mL of the BT cocktail was loaded in a sterile 10-mL syringe, and the syringe tip was covered with sterile needle to prevent any leakage. For the control group, 3 mL PBS was loaded similarly into a syringe.
Administration of BT cocktail and tulathromycin. On day 0, calves were restrained in a squeeze chute and administered treatments. Sterile laryngotracheal mucosal atomization devices (LMA MADgic laryngotracheal mucosal atomization device without syringe; catalog no. MADAY 700; Teleflex, Morrisville, NC) were fitted to loaded syringes, and the atomization device was inserted into each nostril of calves (approximately 15 cm) and sprayed until the syringe was empty. One atomization device was used for the two nostrils of each calf. A total of 6 mL of BT inoculum (3 mL per nasal cavity) was administered to calves in the BT treatment group. For control group animals, PBS was sprayed into the nostrils, similar to the BT inoculum (3 mL per nostril). The metaphylaxis group received a single subcutaneous injection of long-acting tulathromycin (2.5 mg/kg body weight).
Nasopharyngeal swab sampling and processing. In addition to sampling at feedlot arrival (day 21, 24 h prior to treatment administration), NP samples were collected from the right nostril of each calf in the study on days 1, 2, 4, 7, 14, 28, and 42, for a total of 8 sampling days. The NP sampling procedures were described previously (19). Prior to sampling, the right nostril was wiped clean with 70% ethanol. Extended guarded swabs (27 cm) with a rayon bud (catalog no. MW124; Medical Wire & Equipment, Corsham, England) were used for sampling. Swabs were taken while the animals were restrained in a squeeze chute. Swab tips were then cut and placed in a sterile 1.5-mL tube on ice. Samples were transported to the lab and processed within 1 h of collection. At the lab, the swab tip was transferred into a cryovial containing 1 mL brain heart infusion (BHI) with 20% glycerol and vortexed.
(i) Isolation and detection of bovine respiratory pathogens. Aliquots of swab suspension were plated for isolation and detection of BRD-associated pathogens, including M. haemolytica, P. multocida, and H. somni. Culturing, isolation, and PCR identification of the pathogenic isolates were described previously (19,22). The remaining swab suspensions in BHI glycerol stock were stored at 280°C for DNA extraction.
After quality check with FastQC 0.11.5 and MultiQC 1.0 (57), primers and low-quality sequences were trimmed off the raw sequence reads using cutadapt 1.14 (58). The trimmed reads were used to construct amplicon sequence variants (ASVs) using dada2 1.10.0 (59) in R 3.5.1 (60). Unless otherwise stated, all dada2 functions were used with default parameters. Reads were first filtered with dada2::filterAndTrim with a maximum expected error of 1. Error rates were learned for the forward and reverse reads separately, and these error rates were used to infer exact sequences (error correct) for each sample from dereplicated and trimmed reads using pooled=TRUE for the dada2::dada. Following this, the forward and reverse reads were merged using dada2::mergePairs. Chimeras were removed with dada2::removeBimeraDenovo, and taxonomy was assigned using the naive Bayesian classifier (61) as implemented in dada2::assignTaxonomy trained with the Silva training set version 132 (https://zenodo.org/record/1172783#.Y9kyanbMK3A). Specieslevel assignment was done with dada2::addSpecies that uses exact matching to assign species where possible. ASVs were aligned with ssu-align 0.1.1 (62).
Quantification of bacterial, Lactobacillus spp., and antibiotic resistance determinants using quantitative PCR. Real-time PCR was performed to quantify copies of 16S rRNA genes in DNA from nasopharyngeal swabs that were specific to all bacteria, or Lactobacillus, to estimate the abundances of total bacteria and Lactobacillus genera, respectively. Total 16S rRNA gene copies were amplified using primers 515F and 806R described above. Lactobacillus-specific 16S rRNA gene copies were amplified using a genus-specific primer described previously (63). In addition to 16S rRNA genes, the macrolide resistance gene msr(E) and tetracycline resistance gene tet(H) were quantified from the DNA extracted from nasal swabs. The primers for msr(E) and tet(H) were reported previously by Klima et al. (15) and Zhu et al. (64), respectively.
To generate standards for PCR, amplicons of each target gene were cloned into competent E. coli cells using the TOPO cloning reaction kit (Invitrogen) according to the instructions of the manufacturer. Plasmids containing amplicon inserts were purified by the QIAprep Spin miniprep kit (Qiagen, Hilden, Germany) and then serially diluted. Each real-time PCR mixture (25 mL) contained 1 Â iQ SYBR green supermix (Bio-Rad Laboratories Inc.), 0.4 mM each primer, 0.1 mg/mL bovine serum albumin (BSA; New England Biolabs, Pickering, ON, Canada), and 25 ng of DNA. For each PCR, the amount of DNA extracted from the NP swabs was normalized to 10 ng/mL. The quantification of target genes was performed on a CFX96 Touch real-time PCR detection system (Bio-Rad Laboratories Inc.) with the following conditions: an initial denaturation at 95°C for 3 min, followed by 40 cycles at 95°C for 25 s, 50°C for 30 s, and then 72°C for 45 s. For quantification of total and Lactobacillus-specific 16S rRNA gene copies and the resistant gene copies, standards were prepared for each gene using the respective pDrive plasmid containing inserted amplicons and concentrations of 10 6 , 10 5 , 10 4 , 10 3 , and 10 2 copies per reaction (in duplicate). Melt curve analysis was performed on all PCRs to ensure specific amplification. The temperature range was 60°C to 95°C, and fluorescence was measured at 0.5°C intervals.
Statistical analysis. Statistical analysis of sequence data was done with R 3.6.0 with phyloseq 1.28.0 (57) and vegan 2.5.5 (58). Plots were created with ggplot2 3.1.1. Samples with ,1,000 sequences (n = 1) were excluded from further analyses. Sequences matching mitochondria or chloroplast were removed along with any sequences that were not assigned to Bacteria. A filtered copy of the ASV sequence table was created that retained ASVs present (count $ 2) in at least 1% of the samples. Samples with 0 reads remaining after the filter was applied were removed from the ASV table. This served to reduce noise for downstream analysis. The unfiltered version of the sequence table was used for alpha diversity, which was assessed with the Shannon diversity index using scikit-bio 0.3.1 (65). Richness was estimated with a Poisson model using breakaway 4.6.8 (66). To calculate beta diversity, the ASV counts (filtered table) were normalized with a variance-stabilizing transform (using DESeq2 1.24.0) (67) with size factors calculated using GMPR (68), and then sample-sample distances were determined with the Bray-Curtis metric and visualized with detrended correspondence analysis (DCA). Permutational multivariate analysis of variance (PERMANOVA) with 10,000 permutations was used to determine the effect of treatment and sampling time on the microbial community structure. For statistical testing, the model included an interaction between day and treatment to detect differences between treatment groups over time. The beta function from breakaway 4.6.8 (69), which models both observed and unobserved diversity, was used to test alpha diversity. Corncob 0.1.0 (70) was used to test for differentially abundant taxa by fitting a model with and without the interaction term. This was done to identify taxa that showed a change from baseline in the BT and MP groups that differed significantly from the CTRL group. Corncob uses a betabinomial-based model that controls for correlated observations (taxa) and overdispersion.
Ecological network modeling was performed to evaluate the directed microbial interactions among the nasopharyngeal microbial communities using BEEM-static in R. Briefly, based on the abundance profile of all observed genera derived from 8 sampling time points, generalized Lotka-Volterra models (gLVMs), coupling biomass estimation and model inference in an expectation maximization-like algorithm (BEEM), were used to construct three (CTRL, BT, MP groups) ecological network models as described by Li et al. (71). The interaction network inferred by BEEM-statistic was visualized by the plots generated using the graph package of R.
Generalized linear models (SAS Proc GLIMMIX) were built on the prevalence data of M. haemolytica, P. multocida, and H. somni determined by culturing. The Bernoulli/binary distributions were selected for the models. Separate F-tests of the treatment effect were produced, using the SLICE option in the LSMEANS statement, for the respective days of observations. The relative abundance of phyla and quantitative PCR (qPCR) results data were analyzed using the GLIMMIX procedure in SAS (SAS 9.4, SAS Institute Inc., Cary, NC). The individual calf, treatment, and time were included in the CLASS statement. The models were "generalized" due to the specification of response distributions that were not Gaussian normal. Models were "mixed" due to the inclusion of fixed effects (treatment-nested-in-time and time) and random effects (individual). Variance heterogeneity was modeled using a "RANDOM _RESIDUAL_/ GROUP = TreatmentÂTime" statement. Response distributions and structures of the variance-covariance matrix were selected for each genus based on the model fit statistics, i.e., the Bayesian information criterion (BIC). Preliminary models that specified the beta-binomial distribution did not converge. Therefore, alternative distributions were tested, gamma, inverse Gaussian, lognormal, shifted t, Gaussian normal, exponential, and geometric. Statistical significance was declared at P values of ,0.05, and a trend toward significance was considered at values of 0.05 , P , 0.10.
Path analysis was performed to model the interrelationships of selected bacterial species within the NP microbial communities. The construction of an initial path model was based on theories regarding the causal relationships. These path diagrams were used to illustrate the strength and direction of the causal relationships between variables. To identify biological interactions within the NP microbiota and the changes in biological interactions in response to the BT and MP treatment, 16 genera were selected, 12 that exhibited the greatest change in relative abundance in BT and MP groups relative to the CTRL group over the course of study and 4 BRD-associated genera (Mannheimia, Pasteurella, Histophilus, and Mycoplasma). The relative abundance data for these 16 genera were sorted into three subsets based on the treatment groups (BT, CTRL, and MP). The methods of path modeling followed those of Schwinghamer et al. (29). The CORR procedure in SAS (SAS 9.4; SAS Institute Inc.) was used to calculate the matrices of Spearman rank-based correlation coefficients that were used as input for path modeling with SAS Proc CALIS. The initial hypothetical model was: Three modified path models were developed based on the addition and subtraction of paths and covariance terms, based on Wald statistics and LaGrange multiplier statistics that were calculated using the MODIFICATION option in the PROC CALIS statement. Modifications were selected based on lower (better) values of Schwarz's BIC and a value equal to zero for the stability criterion of reciprocal causation. A P value of ,0.05 was considered significant.

Ruminococcaceae UCG005
Data availability. Raw sequence data are available from the NCBI Sequence Read Archive under BioProject accession number PRJNA830920.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.