The day‐to‐day stability of the ruminal and fecal microbiota in lactating dairy cows

Abstract In this study, we examined differences between the microbiota of the ruminal fluid (DR) and feces (DF) from five lactating dairy cows over three consecutive days using 16S rRNA gene sequence‐based analysis. Results showed significant differences between the microbial communities of the DR and DF. In particular, the relative abundance of the phyla Firmicutes and Actinobacteria was significantly lower (q < 0.001) in DR compared with DF, while the relative abundance of Bacteroidetes was significantly higher in DF than that of DR (q < 0.001). A significantly higher relative abundance of the genera Bifidobacterium, 5‐7N15, Clostridium, Epulopiscium, SMB53, Turicibacter, Dorea, Roseburia, and Akkermansia was observed in the DF, while a higher relative abundance of the genera Prevotella, Butyrivibrio, CF231, RFN20, and Succiniclasticum was observed in the DR. A further analysis using the functional prediction program PICRUSt showed that sequences belonging to the 5‐7N15, Akkermansia, Bifidobacterium, Clostridium, Dorea, Epulopiscium, Roseburia, and Turicibacter were significantly and positively correlated with glycan biosynthesis and metabolism, while CF231, Prevotella, RFN20, and Succiniclasticum were significantly and positively correlated with amino acid, lipid, carbohydrate, other amino acid, cofactors, and vitamins metabolism. No significant differences were observed across the three consecutive days in either the DR or DF ecosystems, with no significant differences in the diversity or abundance at the phylum and genus levels suggested that there is a limited day‐to‐day variability in the gut microbiota.

is available on the role of rumen fermentation in ruminants, relatively less is known about the hindgut fermentation process in dairy cows.
Because feed presented to the rumen and hindgut may differ, the microbiota in the rumen and feces are likely also distinct. Mao, Zhang, Liu, and Zhu (2015) analyzed the bacterial communities along the gastrointestinal tract in six lactating dairy cows and found significant differences in microbial community composition, including species diversity and abundance. In a different study, Liu, Zhang, Zhang, Zhu, and Mao (2016) reported decreased species diversity in fecal microbiota compared with rumen microbiota. Similarly, de Oliveira et al. (2013) found overt differences in microbial communities between the rumen and feces of Brazilian Nelore steers. The microbial composition of the ruminal fluid and feces can be further affected by other factors such as diet (McCann, Wickersham, & Loor, 2014). Thus, dissimilarities among diets and sampling sites may hinder analyses of these bacterial communities.
In addition, while some studies have examined representative microbiota samples from ruminal fluid or feces collected at a single time point, many other studies have used samples collected over three consecutive days (Shabat et al., 2016;Zhou et al., 2018).
Based on 16S rRNA sequencing techniques, Skarlupka, Kamenetsky, Jewell, and Suen (2019) found a limited day-to-day variability in the rumen microbiota from three consecutive days. However, the work carried out by Skarlupka et al. (2019) was mainly focused on the rumen microbiota and did not explore the difference in the fecal microbiota on a day-to-day basis. Thus, we characterized the microbial communities of ruminal fluid and feces from five lactating dairy cows over three consecutive days to examine differences in bacterial communities among sampling days.

| Animals and sample collection
Five five-year-old lactating Holstein dairy cows were housed in a free barn at a commercial dairy farm (Beijing, China) and were cared for according to the practices outlined in the Guide for the

Care and Use of Agricultural Animals in Agricultural Research and
Teaching (FASS, 2010). Cows were fed a total mixed ration ad libitum, formulated to meet or exceed the energy requirements of lactating dairy cows (NRC, 2001) (Table A1). The cows were adapted to the diet and barn for at least 40 days prior to sample collection. The cows had free access to fresh water. During this period, none of the cows experienced disease and did not receive any antibiotic treatment.
DR was collected via the mouth prior to morning feeding on three consecutive days using a stomach tube with a rumen vacuum sampler (A1141K, Ancitech, Winnipeg, CA). DR samples were filtered through a four-layer cheesecloth. Fecal samples were obtained from the rectum prior to morning feeding on three consecutive days using sterile gloves. All samples were immediately frozen in liquid nitrogen and stored at −80°C until further analysis.

| DNA extraction, PCR amplification, and Illumina MiSeq sequencing
Total DNA was extracted from the DR and DF samples using a Qiagen DNA Extraction Kit (Qiagen, Hilden, Germany). DNA samples were quantified using a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and then stored at −80°C until being used as a template for the PCR assays. Barcoded primers 343F (5ʹ-GATCCTACGGGAGGCAGCA-3ʹ) and 534R (5ʹ-GCTTACCGCGGCTGCTGGC-3ʹ) were used to amplify the V3-V4 region of the 16S rRNA gene (Ji et al., 2017). PCR reactions for each sample were carried out on a Mastercycler Gradient (Eppendorf, Germany) using 25 μl reaction volumes containing 2 μL template DNA, 5 μM of each forward and reverse primers, 12.5 μL 2 × Taq PCR MasterMix (KAPA Biosystems, Wilmington, USA), 3 μL BSA(2ng/μl), and 5.5 μL ddH 2 O. PCR assays consist of an initial denaturing step at 95°C for 5 min followed by 32 cycles of 95°C for 45 s, 55°C for 50 s, and 72°C for 45 s with a final extension at 72°C for 10 min. Samples were run on a 2% agarose gel using electrophoresis, and those samples with a band appearing between 200 and 210 bp were extracted and purified using a QIAquick Gel Extraction Kit (QIAGEN, Germany).
Following qualification and quantification, purified amplicons were pooled at equal molarity and sequenced on an Illumina MiSeq platform using a MiSeq Reagent Kit V3 (600-cycle, Illumina, San Diego, CA, USA) according to standard protocols (Caporaso et al., 2012).

| Taxonomic classification and diversity analysis
Shorter reads (lower than 200 bp), low-quality sequences (scores lower than 20), ambiguous bases, and chimeras were removed by USEARCH in the QIIME1 pipeline (version 1.5.0) (Caporaso et al., 2011). Clean, paired-end sequence reads were merged using FLASH (version 1.2.7) (Magoc & Salzberg, 2011), and 16S rRNA sequences were classified using UCLUST (version 1.2.22) against the SILVA bacterial database (SILVA version 119, released in July, 2014) (Pruesse et al., 2007). Operational taxonomic units (OTUs) determined at 97% similarity was carried out using UCLUST (Edgar, 2010). All singletons and doubletons were removed using UCLUST to generate a representative OTU table. To ensure even sequencing depth across samples, the number of tags per sample was randomly subsampled to 15,900 for bacterial community analysis.

| Statistical analysis
Alpha diversity analysis was performed using QIIME including Chao 1 index calculation, determination of the number of OTUs, phylogenetic diversity whole tree analysis, and Shannon's diversity index were calculated from the OTU table. Beta diversity indices were measured based on unweighted UniFrac distances and displayed using principal coordinate analysis (PCoA) (Lozupone & Knight, 2005). Analysis of similarities (ANOSIM) was performed to determine the overall difference in microbial communities by sampling sites using the Vegan package (Oksanen et al., 2015) in R (version 3.1.2). Functional classification was predicted using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt, https ://hutte nhower.sph.harva rd.edu/ galaxy) from 16S rRNA sequencing, and Kyoto Encyclopedia of Gene and Genomes (KEGG) functional composition profiles were generated (Langille et al., 2013). Spearman's rank correlation analysis was used to analyze the correlation between predicted functional profiles and the genus (relative abundance of more than 1%).
Specific bacterial taxa were tested using the Kruskal-Wallis analysis to assess differences in relative abundance on different sampling days within the same community, while the Wilcoxon rank-sum test was used to determine differences in relative abundance at the phylum and genus levels between communities. Statistical analysis was performed in R (version 3.1.2), with statistical significance accepted at p < .05. All P-values from the multiple comparison analyses were adjusted based on the false discovery rate (Benjamini & Hochberg, 1995), with statistical significance accepted at adjusted q-values < 0.05.

| Differences in bacterial community diversity and composition between sampling sites and over different days
Our 16S rRNA gene sequencing produced an average of 22,332 ± 4,411 high-quality sequences for each of the 30 samples in our study. Based on a sequence similarity cutoff of 97%, a total of 1,679 OTUs were detected across all samples. Rarefaction analysis was conducted to assess OTU coverage, producing a Good's coverage value > 0.97 ( Figure A1) for each sample, implying that the sequence coverage was sufficient. Chao 1, number of OTUs, phylogenetic diversity whole tree analysis, and Shannon diversity index values for DR were significantly higher than those obtained for DF (q < 0.001).
Importantly, no differences in these values were observed for either community across different sampling days (q > 0.05) ( Table 1).
We then used a PCoA to examine difference in microbial community structure. While there was an obvious separation of the DR and DF bacterial communities, samples from the same community col-

| Bacterial community composition
Only 497 of the 1,679 identified OTUs were shared between the DF and DR bacterial communities ( Figure A2). Overall, 883 OTUs were specific to the DR samples and 386 OTUs were specific to the DF samples ( Figure A2).
A total of 17 phyla were identified across all samples. The bacterial communities from all samples were dominated by the genera Firmicutes and Bacteroidetes (Figure 3a). In the DF bacterial community, Firmicutes was the predominant phylum, with a relative abundance of up to 80.79%, followed by Actinobacteria (11.97%) and Bacteroidetes (5.51%) (Figure 3a). The predominant phylum in the DR microbiota was Bacteroidetes (58.82%), followed by Firmicutes

| Differences in community composition between the DF and DR samples
The relative abundance of Firmicutes and Actinobacteria was significantly increased in the DF (q < 0.05) ( Figure 3b and Table A2), while the DR microbiota contained a higher relative abundance of Bacteroidetes (q < 0.001). At the phylum level, no differences were observed in either the DR or DF communities over the three consecutive sampling days except for in the relative abundance of Bacteroidetes, Firmicutes, and Actinobacteria (q < 0.001). At the phylum level, no differences were observed over three consecutive sampling days in DR bacterial communities and in DF bacterial communities (q > 0.05).
At the genus level, the relative abundance of Bifidobacterium, 5-7N15, Clostridium, Epulopiscium, SMB53, Turicibacter, Dorea, Roseburia, and Akkermansia was significantly higher (q < 0.01) in the DF community compared with the DR community, while the relative abundance of Prevotella, CF231, Butyrivibrio, RFN20, and Succiniclasticum was significantly lower (q < 0.05) in DF community than DR community (Figure 4b and Table A3). There were no changes in the relative abundance of genera among different sampling days for the DR and DF microbial communities.

| Divergence of predicted microbial functions in the DR and DF Groups
We then performed a PICRUSt analysis and found 22 KEGG pathways (level 2), 18 of which were significantly different between the DR and DF samples. Lipid metabolism (10.63%), biosynthesis of other secondary metabolites (7.40%), xenobiotics biodegradation and metabolism (4.84%), amino acid metabolism (3.74%), and signaling molecules and interaction (3.59%) were identified as the top 5 predicted functions for the DR microbiota, which were also the top predicted function for DF microbiota. Cell communication, cell growth and death, and glycan biosynthesis and metabolism were significantly higher in DF compared to the DR group (q < 0.05), whereas lipid metabolism, amino acid metabolism, carbohydrate metabolism, enzyme families, metabolism of cofactors and vitamins, and metabolism of other amino acids were significantly higher in DR compared to the DR group (q < 0.05, Figure 5a).

| Relationship between bacterial community and predicted function
To examine the relationship between DR and DF communities, Spearman's rank correlation was used to identify linkages between predicted function and taxonomy (Figure 5b). We found that the relative abundance of Prevotella, Succiniclasticum, CF231, and RFN20 spp. was positively correlated with amino acid metabolism, biosynthesis of other secondary metabolites, carbohydrate metabolism, digestive system, enzyme families, folding sorting and degradation, lipid metabolism, membrane transport, metabolism of cofactors and vitamins, metabolism of other amino acids, metabolism of terpenoids and polyketides, nucleotide metabolism, replication and repair, signaling molecules, and interaction and transcription. We also found that the relative abundance of 5-7N15, Akkermansia, Bifidobacterium, Clostridium, Dorea, Epulopiscum, Roseburia, and Turicibacter was negatively correlated with those predicted function (q < 0.05, Figure 5b).

| D ISCUSS I ON
The gastrointestinal tract of ruminant animals harbors a large number of symbiotic microbes that allow the host to acquire nutrients    reported to be more highly abundant in the feces of cattle (Callaway et al., 2010;Mao, Huo, & Zhu, 2013) and is thought to be associated with dermatitis in cattle and contagious dermatitis in sheep (Evans, Brown, Hartley, Smith, & Carter, 2012;Sadet, Martin, Meunier, & Morgavi, 2007).
Thirteen genera (relative abundance of each sample > 1%) were differentially abundant between the DR and DF (Table A3) The Spearman correlation between the genus level relative abundance more than 1% and the selected predicted KEGG pathways. *0.01 < q < 0.05, **0.01 < q < 0.001, ***0.001 < q < 0.0001, ****0.0001 < q < 0.00001. The reads color genus name means the higher relative abundance in DR than those of DF  (Flint & Bayer, 2008;Dodd, Mackie, & Cann, 2011). Butyrate and propionate are the most important SCFAs in dairy cattle. Members of the genus Butyrivibrio are butyrate producers (Mohammed et al., 2014), and Succiniclasticum spp. can produce propionate from succinate (van Gylswyk, 1995). This may partly explain the significantly higher abundance of Butyrivibrio and Succiniclasticum in the DR community (Figure 4b and Table A3). Two unclassified bacteria were also identified in the present and the other studies (Jewell, McCormick, Odt, Weimer, & Suen, 2015;Zhu et al., 2018), including CF231 and RFN20, which were significantly abundant in the DR community.
In addition, in the present study, the predicted function of the lipid, carbohydrate, and amino acid metabolism by PICRUSt were positive and significantly correlated with Prevotella, CF231, RFN20, and Succiniclasticum (Figure 5a,b). Butyrivibrio has strongly correlated with predicted lipid, carbohydrate, cofactors and vitamins, and other amino acid and metabolism (Figure 5b).
The SCFAs producer bacteria and metabolism function were both enriched in the DR community suggested rumen bacteria was related to the diet and rumen function of plant digestion and VFAs production.
Another important finding of this study was the lack of significant variation among samples collected over three consecutive days from each of the gut communities, indicating that dairy cows have two relatively independent and stable microbial communities in the gut when clinically healthy. This is consistent with the results of Skarlupka et al. (2019); they compared the rumen liquid and solid community from three consecutive days and found that the rumen community is limited day-to-day variability in the rumen microbiota.
But in our study, we found that feces bacterial community were also like stable on a day-to-day basis, indicating that the gut microbial community is likely stable on a day-to-day basis.

| CON CLUS IONS
Our results demonstrate striking differences in the composition of bacterial communities from the DR and DF of lactating dairy cattle, indicating that two relatively independent and stable microbial communities were existing in the gut of dairy cows. Furthermore, no significant differences in the microbiota among samples collected over three consecutive days from either the DR or DF indicated that the gut microbial community is likely stable on a day-to-day basis. Part of this work was assisted by Zhongdi Dairy Farm (Beijing, China).

CO N FLI C T O F I NTE R E S T S
None declared.

E TH I C S S TATEM ENT
The experimental protocol was reviewed and approved by the Ethical Committee of the College of Animal Science and Technology of China Agricultural University. Abbreviations: DF, feces samples; DR, rumen fluid samples.

DATA AVA I L A B I L I T Y S TAT E M E N T
*q value < 0.05, **q value < 0.01, ***q value < 0.001, ns means q > 0.05.

TA B L E A 2
The relative abundance of major bacterial phyla (>1%) between rumen fluid and feces Abbreviations: DF, feces samples; DR, rumen fluid samples.