The experiment was conducted as per the guidelines of the Canadian Council of Animal Care (CCAC, 1993) at the Dairy Research and Technology Center of the University of Alberta (Edmonton, Alberta, Canada), and the methodology was previously described in our companion study (Villot et al., 2020). The animal use protocol was approved by the University of Alberta Animal Care Committee (Animal Use Protocol # 2645). A total of 20 singlet and naturally delivered Holstein bull calves born from May to August 2018 with a birth body weight (BW) between 35 and 55 kg were enrolled in this study. Calves were removed from the maternity pen at birth to avoid contact with the dam and moved to disinfected indoor pens (150 × 125 cm) bedded with a consistent amount of fresh straw to facilitate nesting. Animals included in this study presented healthy at birth based on general appearance, rectal temperature, umbilical appearance, and nasal discharge; none were vaccinated or treated with therapeutic fluids during the trial. Calves received two meals of standardized colostrum replacer (HeadStart, Saskatoon Colostrum Co., Ltd., Saskatoon, SK, Canada). Briefly, calves received two feedings of colostrum at 2 h after birth fed at 7% of birth BW, and at 12 h after birth fed at 3% of birth BW to receive an average amount of 180, and 120 g of IgG, respectively (Villot et al., 2020). Colostrum was bottle fed, and the remainder of the meal volume was tube-fed if the calves did not consume all of the meal within 30 min. The concentration of IgG and IgA delivered to the calves was measured in a pooled colostrum batch composed of subsamples of all the colostrum fed to the calves, which confirmed that all calves received the same amount of IgG and IgA from colostrum (Villot et al., 2019). Milk replacer was bottle-fed, and it was composed of 260 g/kg of crude protein, 160 g/kg of crude fat, and 4.58 Mcal/kg of metabolizable energy on a dry matter basis (Grober Animal Nutrition, Cambridge, Ontario, Canada). The initial volume of individual meals offered to each calf was 7.5% of their birth BW during the 2 first days of the experiment and increased to 8.5% of BW from d 3 to d 7.

A generalized randomized block design was used to evaluate the impact of SCB supplementation on bacterial colonization and its impact on the development of humoral immunity of Holstein bull calves. Calves were randomly assigned to one of two treatments according to their date of birth and initial BW: Control (CON, n = 10) fed at 5 g/d of carrier with no live yeast; and Saccharomyces cerevisiae boulardii (SCB, n = 10) fed at 5 g of live SCB per day (10 × 10⁹ CFU) from birth to 1 week of life. Both the CON and SCB treatments were fed once daily starting with the first colostrum feeding and subsequently in the morning milk replacer feeding until 1 week of age when calves were euthanized.

At 1 week of life, 3 h after the morning meal, calves were euthanized to collect tissue and digesta samples from the gastrointestinal tract. Calves were euthanized with a pentobarbital sodium injection (Euthanyl, Vetoquinol, Lavaltrie, QC, Canada) at 0.125 mL/kg of BW administered through a jugular catheter. Once the calf reached a surgical plane of anesthesia, exsanguination was performed. The esophagus and rectum were first ligated to prevent loss of digesta and the entire gut digesta were then removed with caution to avoid transfer of digesta between the different segments. Segments of the proximal jejunum, ileum, and colon, along with their digesta were collected according to Villot et al. (2020). Briefly, the proximal jejunum sample was taken 100 cm distal to the pylorus sphincter, the ileum segment was defined as 30 cm proximal to the ileocecal junction, and the colon segment was defined as 30 cm distal to the ileocecal junction. Zip-ties were placed on each end of the segments so that a sample of digesta could be collected. Intestinal digesta were removed from the sample using tweezers and placed in a 50-mL Falcon tube. The tissue was then washed in sterile PBS (pH 7.4) until clean (~3–4 washes) and placed in a sterile bag. Both gut digesta and tissue samples were immediately snap-frozen in liquid nitrogen and then stored at −80°C until further processing.

A primer set that targets the 26S rRNA gene of S. cerevisiae was used to estimate the copy numbers of SCB using real-time quantitative PCR (qPCR), according to Chang et al. (2007) and Villot et al. (2020). Briefly, the qPCR was performed in the high throughput Viia 7 Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA) using SYBR green chemistry (Fast SYBR Green Master Mix, Applied Biosystems, Foster City, CA) with specific primers (Forward: 5′- AGGAGTGCGGTTCTTTG −3′, and Reverse: 5′- TACTTACCGAGGCAAGCTACA −3′). Genomic DNA of SCB was extracted from the commercial probiotic yeast used in the study (ProTernative Milk; CNCM I-1079 strain) to generate a standard curve to estimate the gene copy number of S. cerevisiae in calf digesta from proximal jejunum, ileum, and colon. All qPCR reactions were performed using the optimized conditions described previously by Villot et al. (2020), and the quantity of S. cerevisiae was calculated from the standard curve.

Data on sIgA concentration in intestinal digesta, as well as the analysis procedures, are presented in Villot et al. (2019). Briefly, sIgA was measured in intestinal digesta from jejunum, ileum, and colon on d 7. The concentration of sIgA was analyzed using a sandwich ELISA technique with a commercial kit (Bethyl Laboratories Inc., TX), as per the manufacturer’s recommendations and as described by Villot et al. (2020). All samples were tested in duplicate, and samples with an intra and inter-assay CV greater than 10% were repeated.

The concentrations of short-chain fatty acids (SCFAs; acetic, propionic, butyric, and caproic) were determined via gas chromatography, as described by Schlau et al. (2012). Briefly, digesta samples were thawed and centrifuged at 13,000 rpm for 5 min at 4°C and the supernatant mixed with 25% meta-phosphoric acid (4:1; vol/vol). The samples were then centrifuged again at 13,000 rpm for 5 min at 4°C, and the supernatant was mixed with internal standard solution (5:1; vol/vol) before incubation at −20°C overnight. After incubation, samples were centrifuged at 19,000 × g at 4°C for 5 min and the supernatant transferred to vials. Samples were analyzed by gas chromatography (Varian, Palo Alto, CA) using a flame ionization detector and a capillary column (CP-WAX 58 FFAP 25 m 0.53 mm, Varian CP7767, Varian Analytical Instruments, Walnut Creek, CA).

Extraction of microbial genomic DNA was performed similarly to that previously reported by Yu and Morrison (2004). Briefly, the digesta and tissue samples (~0.3 g ± 0.1 g) were washed twice with Tris-EDTA buffer. After the addition of cell-lysis buffer containing 4% SDS, samples were subjected to physical disruption at 5,000 rpm for 3 min using Biospec Mini Beads Beater 8 (BioSpec, Bartlesville, OK), followed by incubation at 70°C for 15 min and centrifugation for 5 min at 16,000 × g at 10°C. The bead beating, incubation, and centrifugation steps were repeated, and impurities were removed from the supernatant using 10 M ammonium acetate, followed by DNA precipitation using isopropanol. After precipitation, DNA was further purified using QIAmp fast DNA stool mini kit (Qiagen Inc., Germantown, MD). The quantity and quality of DNA were evaluated using the NanoDrop 1000 spectrophotometer, and DNA was stored at −20°C until further use.

For library preparation and sequencing, a fragment of the 16S rRNA gene spanning the V1–V3 hypervariable region was amplified by PCR using dual indexed {universal bacterial primers [9F (5′-GAGTTTGATCMTGGCTCAG-3′); 515R (5′- CCGCGGCKGCTGGCAC- 3′); Kroes et al., 1999]}. The amplicons with targeted size (~480 bp) were purified from 1% agarose gel using a QIAEX II gel extraction kit (Qiagen Inc.). The quality and quantity of purified PCR products were checked using a NanoDrop 1000 (NanoDrop Technologies, Wilmington, DE) to ensure that the concentration of DNA from all samples was higher than 25 ng/μL. Sequencing libraries were prepared using the Nextera XT Index (Illumina) and sequenced on an Illumina MiSeq platform (2 × 300 bp) at Genome Quebec (McGill University, Montreal, QC, Canada).

The adapter sequences were trimmed from the raw fastq files, and the trimmed reads were demultiplexed according to the samples using the bcl2fastq2 conversion software version 2.20.0. (Illumina). For bioinformatic analysis, the sorted reads were imported and processed using the Quantitative Insight into Microbial Ecology (QIIME2) package version 2021.2 (Bolyen et al., 2019). First, low-quality reads (Phred score < 20) and short (<100 bp) reads were filtered out. This was followed by denoising and merging using the plugin DADA2 to generate an amplicon sequence variant (ASV) feature table. Sequences shorter than 400 bp as well as chimeric sequences and singleton ASVs were excluded from further analyses. Taxonomic classification was performed in QIIME2 using a scikit-learn naïve Bayes classifier trained for the V1-V3 region. Next, taxonomic assignment was performed using the trained classifier using the SILVA database (SILVA release 132) with 99% identity for representative bacterial sequences. Alpha and beta diversity were performed using QIIME2 with a sampling depth of 15,000 for the digesta samples, and 7,000 for the tissue samples. Alpha-diversity analysis was conducted with standard diversity metrics accessed via QIIME2, including Chao1, Shannon index, and Faith index. A nonparametric analysis of variance (Kruskal-Wallis test) was used to test differences in α-diversity among treatment groups and to calculate p-values within QIIME2. Beta-diversity was measured by calculating the weighted and unweighted UniFrac distances using QIIME2. Non-metric multidimensional scaling (NMDS) was applied on the resulting distance matrices, and analysis of similarity (ANOSIM) were used to calculate p-values and to test differences in β-diversity among treatment groups for significance. Furthermore, the linear discriminant analysis of effect size (LEfSe; Segata et al., 2011) was used to evaluate ASVs differences between treatments. The data were displayed using linear discriminant analysis score plots for significant ASVs after Kruskal–Wallis and Wilcoxon tests. A size-effect threshold of 3.5 on the logarithmic LDA score was used to identify discriminating taxa. Finally, the relative abundance of discriminant taxa from the LEfSe analysis were compared between groups using among the groups, and Mann–Whitney U test was used to identify taxa that were significantly different.

Microbiota network analysis was performed at the phylum and at the genus level and separately for each treatment group. Genera with less than 0.1% mean relative abundance and less than 25% prevalence among the different samples were excluded. The ASV relative abundance table obtained from QIIME was imported into R package Phyloseq and converted into an adjacency matrix through the use of R package SPIEC-EASI (Kurtz et al., 2015). Next, the Meinshausen-Bühlman neighborhood selection method was used to calculate the conditional dependence between the different ASV pairs. This approach was developed to perform network analysis for compositional and sparse data, such as relative abundance of ASV data, which prevents indirect and spurious associations from being generated (Kurtz et al., 2015). The obtained matrices were visualized as networks at both Phylum and Genus levels. The networks model the relationships between different ASVs, where nodes represent ASVs, and edges represent their co-occurrence across different samples. In order to evaluate ASV centrality in the networks, their hub scores were calculated and extracted from the networks using the SPIEC-EASI hub_score() function, and ranked. The 10 species with the highest hub scores were defined as the keystone species in each group.

Lastly, non-parametric Spearman rank correlation coefficient analysis was used to test the relationship between secretory IgA in ileum digesta (sIgA), short chain fatty acids profile (acetate, propionate, caproate, and total SCFAs), alpha diversity measures (species richness, evenness, and phylogenetic diversity), with the bacterial communities differentially abundant between treatments and Saccharomyces cerevisiae gene copy numbers present in ileum digesta. Physiological data were published previously by Villot et al. (2020) in a companion paper. The resulting correlation matrix was visualized in a heatmap format generated by the corrplot package of R [Corrplot: visualization of a correlation matrix; R package version 4.1.2. 2021].

Supplementation of newborn calves with SCB increased species richness at 1 week of age in the ileum digesta (p = 0.01; Figure 1A) and tended to increase species richness in colon digesta (p = 0.08; Figure 1A), as denoted by the Chao1 index. No differences in evenness (Shannon index) of the bacterial microbiota of calves were observed between treatments in any of the locations in the digesta samples (p > 0.23; Figure 1A). Additionally, SCB supplementation increased the phylogenetic diversity (PD) in the ileum digesta (p < 0.01; Figure 1A), as denoted by the PD faith index. In the tissue samples, no difference between groups were observed in species richness (p > 0.34; Figure 1B), evenness (p > 0.21; Figure 1B), or PD (p > 0.17; Figure 1B). To compare bacterial communities between groups, we used unweighted UniFrac distance matrixes in both digesta and tissue samples of the different segments of the intestine and were visualized by NMDS (Figure 2). These results revealed that the bacterial profiles generated from ileum digesta-associated communities tended to separate into distinct clusters according to treatment (ANOSIM R = 0.20, p = 0.02; Figure 2B), with no differences observed in the other regions of the intestine both for mucosa and digesta.

Copy numbers of 26S rRNA S. cerevisiae gene were detected in the digesta of all intestinal compartments measured. Supplementation of SCB in milk significantly increased the density of S. cerevisiae per g of digesta in jejunum, ileum, and colon compared with non-supplemented calves (p < 0.001; Table 1).

We performed a linear discriminant analysis coupled with effect size measurements (LEfSe) in ileal digesta samples given that it was the only location and sample type where differences in species richness, phylogenetic diversity, and clustering of the bacterial community according to treatment were observed (Figure 1A). Results from LEfSe analysis were used to generate a circular cladogram for the ileal digesta microbiota and to identify the differentially abundant taxa between treatments (Figure 3). The LEfSe analysis in ileum digesta samples indicated that the phyla Actinobacteriota, including the families Corynebacteriaceae, Eggerthellaceae, as well as the families Bacillaceae, Erysipelotrichaceae, Aerococcaceae, Eubacteriaceae, and Ruminococcaceae were among the discriminatory taxa between CON and SCB groups. Next, the relative abundances of the discriminatory taxa between treatments were compared. The relative abundance of the families Ruminococcaceae, Bacillaceae, Eubacteriaceae, Corynebacteriaceae, Eggerthellaceae, Erysipelotrichaceae, Aerococcaceae, and Eubacteriaceae (p < 0.04; Figure 4A), as well as members of the genus Eggerthella, Geobacillus, and Anoxybacillus were significantly increased in the SCB group (p < 0.05; Figure 4B) compared with CON.

To further assess community ecological parameters, we performed a co-occurrence network analysis using SPIEC-EASI to infer microbial ecological networks differences between the CON and SCB groups in ileal digesta. The resulting network revealed that probiotic supplementation increased the microbial interconnectivity between nodes in the network as denoted by the increase in total edges within the SCB network (417 edges) compared to CON (176 edges), and the number of total nodes (SCB = 198 vs. CON = 135) present in the network (Figure 5). In addition, the keystone species (composed of the 10 taxa most central to the microbiota networks according to the hub scores) were different between groups and show that SCB increased the number of different taxa at the genus level that are central to the microbiota network (Figure 6). For instance, the keystone species in CON group consisted primarily of bacteria from the genus Lactobacillus, Escherichia Shigella, and Veillonella (Figure 6A, Hub Score > 0.6), whereas the SCB group consisted of Streptococcus, Lactococcus, Lactobacillus, Pelomonas, Actinomyces, Eubacterium, Lachnoclostridium (Figure 6B, Hub Score > 0.6). Furthermore, the average hub scores were higher for those keystone species in SCB compared with CON groups (0.75 vs. 0.60), suggesting that these keystone species were not only different between treatment groups, but were also more interconnected within their ecological networks.

The concentrations of sIgA was measured in the digesta of distal jejunum, Ileum, Colon, and on feces. Supplementation of SCB significantly increased sIgA concentrations in the ileum, and colon digesta (p < 0.001; Table 2). Concentrations of total SCFAs, acetate, propionate, and caproic acid in proximal jejunum, ileum, and colon were not affected by treatment (p ≥ 0.25; Figure 7; Supplementary Table S3), nor the acetate to propionate ratio (p = 0.80; Supplementary Table 3). However, the acetate concentration was higher in the colon (p < 0.01; Figure 7A; Supplementary Table S3) compared to the ileum and proximal jejunum. Furthermore, when doing pairwise comparison of treatment effects by gut location and adjusted by Tukey for multiple comparisons, the concentration of caproic acid was higher in the ileum for SCB compared to CON (p < 0.03; Figure 7C; Supplementary Table S3).

The relative abundance of the discriminatory taxa obtained in LEfSe in addition to the Saccharomyces cerevisiae gene copy numbers in ileum digesta was then correlated to sIgA, as well as short chain fatty acids profile, and alpha diversity indices in the same segment and sample type using spearman rank correlations (Figure 8). The gene copy numbers of SCB were positively correlated with the concentration of sIgA in ileum digesta (p < 0.01; Spearman’ ρ = 0.58). In addition, the relative abundance of the family Corynebacteriaceae (p < 0.01; Spearman’ ρ = 0.73), and Eubacteriaceae (p < 0.01; Spearman’ ρ = 0.65) were positively correlated with the concentrations of sIgA in the ileum digesta. Additionally, the relative abundance of bacterial families Ruminococcaceae (p < 0.01; Spearman’ ρ = 0.62), Aerococcaceae (p < 0.01; Spearman’ ρ = 0.59), and Corynebacteriaceae (p = 0.03; Spearman’ ρ = 0.48), and genus Eggerthella (p = 0.01; Spearman’ ρ = 0.55) were positively correlated with species richness (Chao1 index). Total SCFAs, acetate, propionate, and their ratio were not correlated with SCB gene copy numbers nor with any of the differentially abundant taxa between groups (Figure 8).